DGtal  1.5.beta
ChamferNorm2D.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file ChamferNorm2D.ih
19  * @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr )
20  * Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
21  *
22  * @date 2013/12/18
23  *
24  * Implementation of inline methods defined in ChamferNorm2D.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 #include "DGtal/arithmetic/IntegerComputer.h"
33 //////////////////////////////////////////////////////////////////////////////
34 
35 ///////////////////////////////////////////////////////////////////////////////
36 // IMPLEMENTATION of inline methods.
37 ///////////////////////////////////////////////////////////////////////////////
38 
39 ///////////////////////////////////////////////////////////////////////////////
40 template <typename TSpace>
41 inline
42 DGtal::experimental::ChamferNorm2D<TSpace>::~ChamferNorm2D()
43 {
44 }
45 ///////////////////////////////////////////////////////////////////////////////
46 template <typename TSpace>
47 inline
48 DGtal::experimental::ChamferNorm2D<TSpace>::ChamferNorm2D(const unsigned int N)
49 {
50  myNorm = 1.0;
51  myDirections.reserve(N);
52  myNormals.reserve(N);
53  int NN = N;
54 
55  IntegerComputer<typename Vector::Component> IC;
56 
57  myDirections.push_back( Vector(0,-1));
58 
59  for(int i = 1; i <= NN; i++)
60  for(int j = -NN; j <= NN; j++)
61  {
62  typename Vector::Component gc = IC.gcd( i,j );
63  if (gc == 1)
64  myDirections.push_back( Vector( i,j) ) ;
65  }
66 
67  //3-4 mask embedding
68  unsigned int nbperSlots = static_cast<unsigned int>(myDirections.size() / 4);
69  for(unsigned int k=0; k < nbperSlots; k++)
70  myNormals.push_back( Vector(1,-3));
71  for(unsigned int k=0; k < nbperSlots; k++)
72  myNormals.push_back( Vector(3,-1));
73  for(unsigned int k=0; k < nbperSlots; k++)
74  myNormals.push_back( Vector(3,1));
75  for(unsigned int k=0; k < nbperSlots; k++)
76  myNormals.push_back( Vector(1,3));
77 
78  //Sorting
79  std::sort(myDirections.begin(), myDirections.end(), myLessThanAngular);
80 
81 #ifdef DEBUG_VERBOSE
82  //Display
83  trace.info()<< "Constructing mask: "<<std::endl;
84  std::copy( myDirections.begin(), myDirections.end(), std::ostream_iterator<Vector>(std::cout, "\n"));
85  trace.info()<<std::endl;
86  std::copy( myNormals.begin(), myNormals.end(), std::ostream_iterator<Vector>(std::cout, "\n"));
87 #endif
88 
89 }
90 ///////////////////////////////////////////////////////////////////////////////
91 template <typename TSpace>
92 inline
93 DGtal::experimental::ChamferNorm2D<TSpace>::ChamferNorm2D(const Directions &aDirSet,
94  const Directions &aNormalDirSet,
95  const Value norm):
96  myDirections(aDirSet), myNormals(aNormalDirSet), myNorm(norm)
97 {
98 }
99 ///////////////////////////////////////////////////////////////////////////////
100 template <typename TSpace>
101 inline
102 void
103 DGtal::experimental::ChamferNorm2D<TSpace>::selfDisplay ( std::ostream & out ) const
104 {
105  out << "[ChamferNorm2D] mask size= "<<myDirections.size();
106 }
107 ///////////////////////////////////////////////////////////////////////////////
108 template <typename TSpace>
109 inline
110 bool
111 DGtal::experimental::ChamferNorm2D<TSpace>::isValid() const
112 {
113  return true;
114 }
115 ///////////////////////////////////////////////////////////////////////////////
116 template <typename TSpace>
117 inline
118 typename DGtal::experimental::ChamferNorm2D<TSpace>::Abscissa
119 DGtal::experimental::ChamferNorm2D<TSpace>::getLowerRayIntersection(const Vector &aP, const Vector &aQ,
120  const Point &Lmin, const Point &Lmax,
121  const Dimension aDimension) const
122 {
123  ASSERT(Lmin[aDimension] < Lmax[aDimension]);
124  ASSERT(Lmin[(aDimension+1) %2] == Lmax[(aDimension+1)%2]);
125  boost::ignore_unused_variable_warning(Lmax);
126 
127  double slope, y;
128  if (aDimension == 1)
129  {
130  if ((aQ[0] - aP[0]) == NumberTraits<Abscissa>::ZERO)
131  {
132  if (aQ[1] > aP[1])
133  return myInfinity;
134  else
135  return -myInfinity;
136  }
137  slope = (aQ[1]-aP[1])/static_cast<double>((aQ[0] - aP[0]));
138  y=0.0;
139 
140  y = slope * Lmin[0] + aP[1] - slope*aP[0];
141  return static_cast<Abscissa>( floor(y));
142  }
143  else
144  {
145  ASSERT( (aQ[0] - aP[0]) != NumberTraits<Abscissa>::ZERO);
146  slope = (aQ[1]-aP[1])/static_cast<double>((aQ[0] - aP[0]));
147  double x = (Lmin[1] - aP[1] + slope*aP[0]) / slope;
148 
149  if ((aQ[1] - aP[1]) == NumberTraits<Abscissa>::ZERO)
150  {
151  if (aQ[0] > aP[0])
152  return myInfinity;
153  else
154  return -myInfinity;
155  }
156 
157  return static_cast<Abscissa>( floor(x));
158  }
159 }
160 ///////////////////////////////////////////////////////////////////////////////
161 template <typename TSpace>
162 inline
163 typename DGtal::experimental::ChamferNorm2D<TSpace>::Vector
164 DGtal::experimental::ChamferNorm2D<TSpace>::getNormalFromCone(ConstIterator aCone) const
165 {
166  size_t pos = aCone - myDirections.begin();
167  return myNormals[pos];
168 }
169 ///////////////////////////////////////////////////////////////////////////////
170 template <typename TSpace>
171 inline
172 typename DGtal::experimental::ChamferNorm2D<TSpace>::Abscissa
173 DGtal::experimental::ChamferNorm2D<TSpace>::getUpperRayIntersection(const Vector &aP, const Vector &aQ,
174  const Point &Lmin, const Point &Lmax,
175  const Dimension aDimension) const
176 {
177  ASSERT(Lmin[aDimension] < Lmax[aDimension]);
178  ASSERT(Lmin[(aDimension+1) %2] == Lmax[(aDimension+1)%2]);
179  boost::ignore_unused_variable_warning(Lmax);
180 
181  double slope, y;
182  if (aDimension == 1)
183  {
184  if ((aQ[0] - aP[0]) == NumberTraits<Abscissa>::ZERO)
185  {
186  if (aQ[1] > aP[1])
187  return myInfinity;
188  else
189  return -myInfinity;
190  }
191  slope = (aQ[1]-aP[1])/static_cast<double>((aQ[0] - aP[0]));
192  y=0.0;
193 
194  y = slope * Lmin[0] + aP[1] - slope*aP[0];
195  return static_cast<Abscissa>( ceil(y));
196  }
197  else
198  {
199  ASSERT( (aQ[0] - aP[0]) != NumberTraits<Abscissa>::ZERO);
200  slope = (aQ[1]-aP[1])/static_cast<double>((aQ[0] - aP[0]));
201  double x = (Lmin[1] - aP[1] + slope*aP[0]) / slope;
202 
203  if ((aQ[1] - aP[1]) == NumberTraits<Abscissa>::ZERO)
204  {
205  if (aQ[0] > aP[0])
206  return myInfinity;
207  else
208  return -myInfinity;
209  }
210 
211  return static_cast<Abscissa>( ceil(x));
212  }
213 }
214 ///////////////////////////////////////////////////////////////////////////////
215 template <typename TSpace>
216 inline
217 typename DGtal::experimental::ChamferNorm2D<TSpace>::ConstIterator
218 DGtal::experimental::ChamferNorm2D<TSpace>::getCone(const Vector &aDirection,
219  ConstIterator aBegin,
220  ConstIterator aEnd) const
221 {
222  return std::upper_bound(aBegin, aEnd , canonicalRay(aDirection), myLessThanAngular) -1;
223 }
224 ///////////////////////////////////////////////////////////////////////////////
225 template <typename TSpace>
226 inline
227 typename DGtal::experimental::ChamferNorm2D<TSpace>::Vector
228 DGtal::experimental::ChamferNorm2D<TSpace>::canonicalRay(const Vector &aRay) const
229 {
230  return Vector( (aRay[0]<0)? -aRay[0]:aRay[0], aRay[1] );
231 }
232 ///////////////////////////////////////////////////////////////////////////////
233 template <typename TSpace>
234 inline
235 typename DGtal::experimental::ChamferNorm2D<TSpace>::RawValue
236 DGtal::experimental::ChamferNorm2D<TSpace>::rawDistance(const Point &P,
237  const Point &Q) const
238 {
239  Vector ray = canonicalRay(Q - P);
240 
241  //The cone
242  ConstIterator it = getCone(ray, myDirections.begin(), myDirections.end());
243 
244  // v = dir_i/w_i
245  Vector normalDir = getNormalFromCone(it);
246  //distance as the scalar product with the associated normal
247  return ray[0]*normalDir[0] + ray[1]*normalDir[1];
248 }
249 ///////////////////////////////////////////////////////////////////////////////
250 template <typename TSpace>
251 inline
252 typename DGtal::experimental::ChamferNorm2D<TSpace>::Value
253 DGtal::experimental::ChamferNorm2D<TSpace>::operator()(const Point &P,
254  const Point &Q) const
255 {
256  return this->rawDistance(P,Q)/myNorm;
257 }
258 ///////////////////////////////////////////////////////////////////////////////
259 template <typename TSpace>
260 inline
261 typename DGtal::experimental::ChamferNorm2D<TSpace>::ConstIterator
262 DGtal::experimental::ChamferNorm2D<TSpace>::shrinkPSubMask(ConstIterator aBegin,
263  ConstIterator aEnd,
264  const Point &aP, const Point &aQ,
265  const Point &Lmin, const Point &Lmax,
266  const Dimension aDimension,
267  Point &midPoint,
268  Point &nextMidPoint) const
269 {
270  ASSERT(Lmin[aDimension] < Lmax[aDimension]);
271  ASSERT(Lmin[(aDimension+1) %2] == Lmax[(aDimension+1)%2]);
272  ASSERT(aP[(aDimension+1) %2] <=Lmin[(aDimension+1) %2]);
273  ASSERT(aDimension == 1);
274  ASSERT(aP[aDimension] != aQ[aDimension]);
275 
276 #ifdef DEBUG_VERBOSE
277  trace.info()<<"Checking "<<*aBegin<<" -- "<< *aEnd << " adimension="<<aDimension << " P="<<aP<<" Q="<<aQ<<std::endl;
278 #endif
279 
280  if ((aEnd - aBegin) > 1)
281  {
282  //Mid
283  ConstIterator mid = aBegin + (aEnd - aBegin) / 2;
284 
285  //beginPoint, endPoint
286  Abscissa upperMid,lowerNext;
287  midPoint = aP + (*mid);
288  nextMidPoint = aP + *(mid+1);
289  upperMid = getUpperRayIntersection(aP, midPoint, Lmin, Lmax, aDimension);
290  lowerNext = getLowerRayIntersection(aP, nextMidPoint, Lmin, Lmax, aDimension);
291 
292  midPoint = Lmin;
293  midPoint[aDimension] = upperMid;
294  nextMidPoint = Lmin;
295  nextMidPoint[aDimension] = lowerNext;
296 
297  //Distances w.r.t. Q, O(log(n)) per distance
298  RawValue dmidQ = this->rawDistance(aQ, midPoint);
299  RawValue dnextmidQ = this->rawDistance(aQ, nextMidPoint);
300 
301 
302  //Distance w.r.t. P, O(1) per distance
303  Vector normalDir = getNormalFromCone(mid);
304 
305  RawValue dmidP = (midPoint[0]-aP[0])*normalDir[0] + (midPoint[1]-aP[1])*normalDir[1];
306  normalDir = getNormalFromCone(mid);
307  RawValue dnextmidP = (nextMidPoint[0]-aP[0])*normalDir[0] + (nextMidPoint[1]-aP[1])*normalDir[1];
308 
309  bool PcloserMid = (dmidP < dmidQ);
310  bool PcloserNextMid = (dnextmidP < dnextmidQ);
311 
312  //We have localized the cone
313  if (PcloserMid != PcloserNextMid)
314  {
315  return mid;
316  }
317 
318  //We decide to which direction we need to recurse
319  if (PcloserMid && PcloserNextMid)
320  {
321  if (aP[aDimension] < aQ[aDimension])
322  {
323  //Recurse up
324  ConstIterator c = shrinkPSubMask(mid,aEnd,aP, aQ, Lmin, Lmax, aDimension, midPoint, nextMidPoint);
325  return c;
326  }
327  else
328  {
329  //Recurse down
330  ConstIterator c = shrinkPSubMask(aBegin, mid, aP, aQ, Lmin, Lmax, aDimension, midPoint, nextMidPoint);
331  return c;
332  }
333  }
334  else
335  {
336  // mid and next in Q
337  if (aP[aDimension] < aQ[aDimension])
338  {
339  //Recurse down
340  ConstIterator c = shrinkPSubMask(aBegin, mid, aP, aQ, Lmin, Lmax, aDimension, midPoint, nextMidPoint);
341  return c;
342  }
343  else
344  {
345  //Recurse up
346  ConstIterator c = shrinkPSubMask(mid,aEnd,aP, aQ, Lmin, Lmax, aDimension, midPoint, nextMidPoint);
347  return c;
348  }
349  }
350  }
351  return aBegin;
352 }
353  ///////////////////////////////////////////////////////////////////////////////
354 template <typename TSpace>
355 inline
356 typename DGtal::experimental::ChamferNorm2D<TSpace>::ConstIterator
357 DGtal::experimental::ChamferNorm2D<TSpace>::shrinkP(ConstIterator aBegin,
358  ConstIterator aEnd,
359  const Point &aP, const Point &aQ,
360  const Point &Lmin, const Point &Lmax,
361  const Dimension aDimension,
362  Point &midPoint,
363  Point &nextMidPoint) const
364 {
365  ASSERT(Lmin[aDimension] < Lmax[aDimension]);
366  ASSERT(Lmin[(aDimension+1) %2] == Lmax[(aDimension+1)%2]);
367 
368  //Special case P in L
369  if (aP[(aDimension +1)%2] == Lmin[(aDimension+1)%2])
370  {
371  //Special Case P on L
372  midPoint = aP;
373  nextMidPoint = aP;
374  if (aP[aDimension] > aQ[aDimension])
375  return aBegin;
376  else
377  {
378  //We retrun the last cone (t,t+1) such that t+1 == aEnd
379  ConstIterator preEnd = aEnd - 1;
380  return preEnd;
381  }
382  }
383 
384  Point P,Q, Lm, LM;
385  ConstIterator c;
386  P = aP;
387  Q = aQ;
388  Lm = Lmin;
389  LM = Lmax;
390 
391  if (aDimension == 1)
392  if (aP[0] > Lmin[0])
393  {
394  //Putting P to the left of L
395  P[0] = Lmin[0] - (aP[0] - Lmin[0]);
396  Q[0] = Lmin[0] - (aQ[0] - Lmin[0]);
397  }
398 
399  if (aDimension == 0)
400  {
401  //Symmetry x<->y
402  Q[0] = aQ[1] ; Q[1] = aQ[0];
403  P[0] = aP[1] ; P[1] = aP[0];
404  Lm[0] = Lmin[1] ; Lm[1] = Lmin[0];
405  LM[0] = Lmax[1] ; LM[1] = Lmax[0];
406  //Swapping x<->y
407  if (P[0] > Lm[0])
408  {
409  P[0] = Lm[0] - (P[0] - Lm[0]);
410  Q[0] = Lm[0] - (Q[0] - Lm[0]);
411  }
412 
413  c = shrinkPSubMask(aBegin, aEnd, P, Q, Lm, LM, 1,midPoint, nextMidPoint);
414 
415  //Setting back the midpoint
416  //Transform back
417  Point mmid = midPoint;
418  midPoint[0] = mmid[1] ; midPoint[1] = mmid[0];
419  Point nextmid = nextMidPoint;
420  nextMidPoint[0] = nextmid[1] ; nextMidPoint[1] = nextmid[0];
421 
422  return c;
423  }
424 
425  // 1str half space (x>0, verticaL)
426  c = shrinkPSubMask(aBegin, aEnd, P, Q, Lm, LM, 1,midPoint, nextMidPoint);
427  return c;
428 }
429 ///////////////////////////////////////////////////////////////////////////////
430 template <typename TSpace>
431 inline
432 typename DGtal::experimental::ChamferNorm2D<TSpace>::Abscissa
433 DGtal::experimental::ChamferNorm2D<TSpace>::getLowerVoronoiEdgeAbscissa(const Point &P, const Point &Q,
434  const Point &startingPoint, const Point &endPoint,
435  const Dimension aDimension) const
436 {
437  Point midPointP,nextMidPointP;
438  Point midPointQ,nextMidPointQ;
439 
440  ConstIterator itBeg = myDirections.begin();
441  ConstIterator itEnd = myDirections.end();
442  ConstIterator coneP = shrinkP(itBeg, itEnd, P, Q, startingPoint, endPoint, aDimension, midPointP, nextMidPointP);
443 
444  ConstIterator coneQ = shrinkP(itBeg, itEnd, Q, P, startingPoint,endPoint,aDimension, midPointQ, nextMidPointQ);
445 
446 
447  Abscissa voroAbscissa;
448  Vector normalP = getNormalFromCone(coneP);
449  Vector normalQ = getNormalFromCone(coneQ);
450 
451  if (aDimension == 1)
452  {
453  //Symmetry w.r.t. L is necessary
454  if (P[0] > startingPoint[0])
455  normalP[0] *= -1;
456  if (Q[0] > startingPoint[0])
457  normalQ[0] *= -1;
458 
459  if (normalP[1] == normalQ[1])
460  {
461  //same cone, Voro edge should be on cone extermities
462  //We first recompute the cone edges
463  Point midPointPP, nextMidPointPP;
464  Point midPointQQ, nextMidPointQQ;
465 
466  if (coneP == itBeg) // and thus coneQ=itBeg
467  {
468  //Return lowest valid cone abscissa
469  if ((nextMidPointP != P) && (nextMidPointQ !=Q)) //P on
470  return (nextMidPointQ[1]<nextMidPointP[1])? nextMidPointQ[1]:nextMidPointP[1];
471 
472  if (nextMidPointP == P)
473  return nextMidPointQ[1];
474 
475  return nextMidPointP[1];
476  }
477 
478  if ((coneP +1) == itEnd) // and thus (coneQ+1)=itEnd
479  {
480  //Return highest valid cone abscissa
481  if ((midPointP != P) && (midPointQ !=Q)) //P on
482  return (midPointQ[1]> midPointP[1])? midPointQ[1] : midPointP[1];
483 
484  if (midPointP == P)
485  {
486  //midPointP===P trick
487  return midPointQ[1];
488  }
489  return midPointP[1];
490  }
491 
492  RawValue dpm = this->rawDistance(P, midPointP);
493  RawValue dqm = this->rawDistance(Q, midPointQ);//
494  //same cone, Voro edge should be on cone extermities
495  //SPECIAL NORMAL CONE CASE
496  if (dpm == dqm)
497  return midPointP[1];
498  else
499  return nextMidPointP[1];
500 
501  }
502  else
503  voroAbscissa = static_cast<Abscissa>(floor((double) (P[1]*normalP[1] - Q[1]*normalQ[1] - (startingPoint[0] - P[0])*normalP[0] + (startingPoint[0] - Q[0])*normalQ[0]) /(double)(normalP[1] - normalQ[1] ) ));
504  }
505  else
506  {
507  //Symmetry w.r.t. L if necessary
508  //Symmetry x<->y
509 
510  Abscissa tmp = normalP[0];
511  normalP[0] = normalP[1];
512  normalP[1] = tmp;
513 
514  tmp = normalQ[0];
515  normalQ[0] = normalQ[1];
516  normalQ[1] = tmp;
517 
518  if (P[1] > startingPoint[1])
519  normalP[1] *= -1;
520  if (Q[1] > startingPoint[1])
521  normalQ[1] *= -1;
522 
523  if (normalP[0] == normalQ[0])
524  {
525  //same cone, Voro edge should be on cone extermities
526  RawValue dpm = this->rawDistance(P, midPointP);
527  RawValue dqm = this->rawDistance(Q, midPointQ);
528  if (dpm == dqm)
529  return midPointP[0];
530  else return nextMidPointP[0];
531 
532  }
533  else voroAbscissa = static_cast<Abscissa>(floor((double) (P[0]*normalP[0] - Q[0]*normalQ[0] - (startingPoint[1] - P[1])*normalP[1] + (startingPoint[1] - Q[1])*normalQ[1]) /(double)(normalP[0] - normalQ[0]) ));
534  }
535  return voroAbscissa;
536 }
537 ///////////////////////////////////////////////////////////////////////////////
538 template <typename TSpace>
539 inline
540 bool
541 DGtal::experimental::ChamferNorm2D<TSpace>::hiddenBy(const Point &u, const Point &v, const Point &w,
542  const Point &startingPoint, const Point &endPoint,
543  const Dimension aDimension) const
544 {
545  ASSERT( u[aDimension] < v[aDimension]);
546  ASSERT( u[aDimension] < w[aDimension]);
547  ASSERT( v[aDimension] < w[aDimension]);
548 
549  Abscissa x_uv = getLowerVoronoiEdgeAbscissa(u, v, startingPoint, endPoint, aDimension);
550 
551  if ((x_uv > endPoint[aDimension]))
552  {
553  // outside1 ==> true
554  return true;
555  }
556 
557  Abscissa x_vw = getLowerVoronoiEdgeAbscissa(v, w, startingPoint, endPoint, aDimension);
558  if ((x_vw < startingPoint[aDimension]))
559  {
560  //outside ==> true
561  return true;
562  }
563 
564  return (x_uv > x_vw);
565 }
566 ///////////////////////////////////////////////////////////////////////////////
567 // Implementation of inline functions //
568 template <typename TSpace>
569 inline
570 std::ostream&
571 DGtal::operator<< ( std::ostream & out,
572  const DGtal::experimental::ChamferNorm2D<TSpace> & object )
573 {
574  object.selfDisplay( out );
575  return out;
576 }
577 // //
578 ///////////////////////////////////////////////////////////////////////////////