DGtal 1.4.0
Loading...
Searching...
No Matches
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///////////////////////////////////////////////////////////////////////////////
40template <typename TSpace>
41inline
42DGtal::experimental::ChamferNorm2D<TSpace>::~ChamferNorm2D()
43{
44}
45///////////////////////////////////////////////////////////////////////////////
46template <typename TSpace>
47inline
48DGtal::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///////////////////////////////////////////////////////////////////////////////
91template <typename TSpace>
92inline
93DGtal::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///////////////////////////////////////////////////////////////////////////////
100template <typename TSpace>
101inline
102void
103DGtal::experimental::ChamferNorm2D<TSpace>::selfDisplay ( std::ostream & out ) const
104{
105 out << "[ChamferNorm2D] mask size= "<<myDirections.size();
106}
107///////////////////////////////////////////////////////////////////////////////
108template <typename TSpace>
109inline
110bool
111DGtal::experimental::ChamferNorm2D<TSpace>::isValid() const
112{
113 return true;
114}
115///////////////////////////////////////////////////////////////////////////////
116template <typename TSpace>
117inline
118typename DGtal::experimental::ChamferNorm2D<TSpace>::Abscissa
119DGtal::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///////////////////////////////////////////////////////////////////////////////
161template <typename TSpace>
162inline
163typename DGtal::experimental::ChamferNorm2D<TSpace>::Vector
164DGtal::experimental::ChamferNorm2D<TSpace>::getNormalFromCone(ConstIterator aCone) const
165{
166 size_t pos = aCone - myDirections.begin();
167 return myNormals[pos];
168}
169///////////////////////////////////////////////////////////////////////////////
170template <typename TSpace>
171inline
172typename DGtal::experimental::ChamferNorm2D<TSpace>::Abscissa
173DGtal::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///////////////////////////////////////////////////////////////////////////////
215template <typename TSpace>
216inline
217typename DGtal::experimental::ChamferNorm2D<TSpace>::ConstIterator
218DGtal::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///////////////////////////////////////////////////////////////////////////////
225template <typename TSpace>
226inline
227typename DGtal::experimental::ChamferNorm2D<TSpace>::Vector
228DGtal::experimental::ChamferNorm2D<TSpace>::canonicalRay(const Vector &aRay) const
229{
230 return Vector( (aRay[0]<0)? -aRay[0]:aRay[0], aRay[1] );
231}
232///////////////////////////////////////////////////////////////////////////////
233template <typename TSpace>
234inline
235typename DGtal::experimental::ChamferNorm2D<TSpace>::RawValue
236DGtal::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///////////////////////////////////////////////////////////////////////////////
250template <typename TSpace>
251inline
252typename DGtal::experimental::ChamferNorm2D<TSpace>::Value
253DGtal::experimental::ChamferNorm2D<TSpace>::operator()(const Point &P,
254 const Point &Q) const
255{
256 return this->rawDistance(P,Q)/myNorm;
257}
258///////////////////////////////////////////////////////////////////////////////
259template <typename TSpace>
260inline
261typename DGtal::experimental::ChamferNorm2D<TSpace>::ConstIterator
262DGtal::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 ///////////////////////////////////////////////////////////////////////////////
354template <typename TSpace>
355inline
356typename DGtal::experimental::ChamferNorm2D<TSpace>::ConstIterator
357DGtal::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///////////////////////////////////////////////////////////////////////////////
430template <typename TSpace>
431inline
432typename DGtal::experimental::ChamferNorm2D<TSpace>::Abscissa
433DGtal::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///////////////////////////////////////////////////////////////////////////////
538template <typename TSpace>
539inline
540bool
541DGtal::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 //
568template <typename TSpace>
569inline
570std::ostream&
571DGtal::operator<< ( std::ostream & out,
572 const DGtal::experimental::ChamferNorm2D<TSpace> & object )
573{
574 object.selfDisplay( out );
575 return out;
576}
577// //
578///////////////////////////////////////////////////////////////////////////////