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.
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.
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/>.
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
24 * Implementation of inline methods defined in ChamferNorm2D.h
26 * This file is part of the DGtal library.
30//////////////////////////////////////////////////////////////////////////////
32#include "DGtal/arithmetic/IntegerComputer.h"
33//////////////////////////////////////////////////////////////////////////////
35///////////////////////////////////////////////////////////////////////////////
36// IMPLEMENTATION of inline methods.
37///////////////////////////////////////////////////////////////////////////////
39///////////////////////////////////////////////////////////////////////////////
40template <typename TSpace>
42DGtal::experimental::ChamferNorm2D<TSpace>::~ChamferNorm2D()
45///////////////////////////////////////////////////////////////////////////////
46template <typename TSpace>
48DGtal::experimental::ChamferNorm2D<TSpace>::ChamferNorm2D(const unsigned int N)
51 myDirections.reserve(N);
55 IntegerComputer<typename Vector::Component> IC;
57 myDirections.push_back( Vector(0,-1));
59 for(int i = 1; i <= NN; i++)
60 for(int j = -NN; j <= NN; j++)
62 typename Vector::Component gc = IC.gcd( i,j );
64 myDirections.push_back( Vector( i,j) ) ;
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));
79 std::sort(myDirections.begin(), myDirections.end(), myLessThanAngular);
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"));
90///////////////////////////////////////////////////////////////////////////////
91template <typename TSpace>
93DGtal::experimental::ChamferNorm2D<TSpace>::ChamferNorm2D(const Directions &aDirSet,
94 const Directions &aNormalDirSet,
96 myDirections(aDirSet), myNormals(aNormalDirSet), myNorm(norm)
99///////////////////////////////////////////////////////////////////////////////
100template <typename TSpace>
103DGtal::experimental::ChamferNorm2D<TSpace>::selfDisplay ( std::ostream & out ) const
105 out << "[ChamferNorm2D] mask size= "<<myDirections.size();
107///////////////////////////////////////////////////////////////////////////////
108template <typename TSpace>
111DGtal::experimental::ChamferNorm2D<TSpace>::isValid() const
115///////////////////////////////////////////////////////////////////////////////
116template <typename TSpace>
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
123 ASSERT(Lmin[aDimension] < Lmax[aDimension]);
124 ASSERT(Lmin[(aDimension+1) %2] == Lmax[(aDimension+1)%2]);
125 boost::ignore_unused_variable_warning(Lmax);
130 if ((aQ[0] - aP[0]) == NumberTraits<Abscissa>::ZERO)
137 slope = (aQ[1]-aP[1])/static_cast<double>((aQ[0] - aP[0]));
140 y = slope * Lmin[0] + aP[1] - slope*aP[0];
141 return static_cast<Abscissa>( floor(y));
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;
149 if ((aQ[1] - aP[1]) == NumberTraits<Abscissa>::ZERO)
157 return static_cast<Abscissa>( floor(x));
160///////////////////////////////////////////////////////////////////////////////
161template <typename TSpace>
163typename DGtal::experimental::ChamferNorm2D<TSpace>::Vector
164DGtal::experimental::ChamferNorm2D<TSpace>::getNormalFromCone(ConstIterator aCone) const
166 size_t pos = aCone - myDirections.begin();
167 return myNormals[pos];
169///////////////////////////////////////////////////////////////////////////////
170template <typename TSpace>
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
177 ASSERT(Lmin[aDimension] < Lmax[aDimension]);
178 ASSERT(Lmin[(aDimension+1) %2] == Lmax[(aDimension+1)%2]);
179 boost::ignore_unused_variable_warning(Lmax);
184 if ((aQ[0] - aP[0]) == NumberTraits<Abscissa>::ZERO)
191 slope = (aQ[1]-aP[1])/static_cast<double>((aQ[0] - aP[0]));
194 y = slope * Lmin[0] + aP[1] - slope*aP[0];
195 return static_cast<Abscissa>( ceil(y));
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;
203 if ((aQ[1] - aP[1]) == NumberTraits<Abscissa>::ZERO)
211 return static_cast<Abscissa>( ceil(x));
214///////////////////////////////////////////////////////////////////////////////
215template <typename TSpace>
217typename DGtal::experimental::ChamferNorm2D<TSpace>::ConstIterator
218DGtal::experimental::ChamferNorm2D<TSpace>::getCone(const Vector &aDirection,
219 ConstIterator aBegin,
220 ConstIterator aEnd) const
222 return std::upper_bound(aBegin, aEnd , canonicalRay(aDirection), myLessThanAngular) -1;
224///////////////////////////////////////////////////////////////////////////////
225template <typename TSpace>
227typename DGtal::experimental::ChamferNorm2D<TSpace>::Vector
228DGtal::experimental::ChamferNorm2D<TSpace>::canonicalRay(const Vector &aRay) const
230 return Vector( (aRay[0]<0)? -aRay[0]:aRay[0], aRay[1] );
232///////////////////////////////////////////////////////////////////////////////
233template <typename TSpace>
235typename DGtal::experimental::ChamferNorm2D<TSpace>::RawValue
236DGtal::experimental::ChamferNorm2D<TSpace>::rawDistance(const Point &P,
237 const Point &Q) const
239 Vector ray = canonicalRay(Q - P);
242 ConstIterator it = getCone(ray, myDirections.begin(), myDirections.end());
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];
249///////////////////////////////////////////////////////////////////////////////
250template <typename TSpace>
252typename DGtal::experimental::ChamferNorm2D<TSpace>::Value
253DGtal::experimental::ChamferNorm2D<TSpace>::operator()(const Point &P,
254 const Point &Q) const
256 return this->rawDistance(P,Q)/myNorm;
258///////////////////////////////////////////////////////////////////////////////
259template <typename TSpace>
261typename DGtal::experimental::ChamferNorm2D<TSpace>::ConstIterator
262DGtal::experimental::ChamferNorm2D<TSpace>::shrinkPSubMask(ConstIterator aBegin,
264 const Point &aP, const Point &aQ,
265 const Point &Lmin, const Point &Lmax,
266 const Dimension aDimension,
268 Point &nextMidPoint) const
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]);
277 trace.info()<<"Checking "<<*aBegin<<" -- "<< *aEnd << " adimension="<<aDimension << " P="<<aP<<" Q="<<aQ<<std::endl;
280 if ((aEnd - aBegin) > 1)
283 ConstIterator mid = aBegin + (aEnd - aBegin) / 2;
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);
293 midPoint[aDimension] = upperMid;
295 nextMidPoint[aDimension] = lowerNext;
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);
302 //Distance w.r.t. P, O(1) per distance
303 Vector normalDir = getNormalFromCone(mid);
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];
309 bool PcloserMid = (dmidP < dmidQ);
310 bool PcloserNextMid = (dnextmidP < dnextmidQ);
312 //We have localized the cone
313 if (PcloserMid != PcloserNextMid)
318 //We decide to which direction we need to recurse
319 if (PcloserMid && PcloserNextMid)
321 if (aP[aDimension] < aQ[aDimension])
324 ConstIterator c = shrinkPSubMask(mid,aEnd,aP, aQ, Lmin, Lmax, aDimension, midPoint, nextMidPoint);
330 ConstIterator c = shrinkPSubMask(aBegin, mid, aP, aQ, Lmin, Lmax, aDimension, midPoint, nextMidPoint);
337 if (aP[aDimension] < aQ[aDimension])
340 ConstIterator c = shrinkPSubMask(aBegin, mid, aP, aQ, Lmin, Lmax, aDimension, midPoint, nextMidPoint);
346 ConstIterator c = shrinkPSubMask(mid,aEnd,aP, aQ, Lmin, Lmax, aDimension, midPoint, nextMidPoint);
353 ///////////////////////////////////////////////////////////////////////////////
354template <typename TSpace>
356typename DGtal::experimental::ChamferNorm2D<TSpace>::ConstIterator
357DGtal::experimental::ChamferNorm2D<TSpace>::shrinkP(ConstIterator aBegin,
359 const Point &aP, const Point &aQ,
360 const Point &Lmin, const Point &Lmax,
361 const Dimension aDimension,
363 Point &nextMidPoint) const
365 ASSERT(Lmin[aDimension] < Lmax[aDimension]);
366 ASSERT(Lmin[(aDimension+1) %2] == Lmax[(aDimension+1)%2]);
368 //Special case P in L
369 if (aP[(aDimension +1)%2] == Lmin[(aDimension+1)%2])
371 //Special Case P on L
374 if (aP[aDimension] > aQ[aDimension])
378 //We retrun the last cone (t,t+1) such that t+1 == aEnd
379 ConstIterator preEnd = aEnd - 1;
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]);
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];
409 P[0] = Lm[0] - (P[0] - Lm[0]);
410 Q[0] = Lm[0] - (Q[0] - Lm[0]);
413 c = shrinkPSubMask(aBegin, aEnd, P, Q, Lm, LM, 1,midPoint, nextMidPoint);
415 //Setting back the midpoint
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];
425 // 1str half space (x>0, verticaL)
426 c = shrinkPSubMask(aBegin, aEnd, P, Q, Lm, LM, 1,midPoint, nextMidPoint);
429///////////////////////////////////////////////////////////////////////////////
430template <typename TSpace>
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
437 Point midPointP,nextMidPointP;
438 Point midPointQ,nextMidPointQ;
440 ConstIterator itBeg = myDirections.begin();
441 ConstIterator itEnd = myDirections.end();
442 ConstIterator coneP = shrinkP(itBeg, itEnd, P, Q, startingPoint, endPoint, aDimension, midPointP, nextMidPointP);
444 ConstIterator coneQ = shrinkP(itBeg, itEnd, Q, P, startingPoint,endPoint,aDimension, midPointQ, nextMidPointQ);
447 Abscissa voroAbscissa;
448 Vector normalP = getNormalFromCone(coneP);
449 Vector normalQ = getNormalFromCone(coneQ);
453 //Symmetry w.r.t. L is necessary
454 if (P[0] > startingPoint[0])
456 if (Q[0] > startingPoint[0])
459 if (normalP[1] == normalQ[1])
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;
466 if (coneP == itBeg) // and thus coneQ=itBeg
468 //Return lowest valid cone abscissa
469 if ((nextMidPointP != P) && (nextMidPointQ !=Q)) //P on
470 return (nextMidPointQ[1]<nextMidPointP[1])? nextMidPointQ[1]:nextMidPointP[1];
472 if (nextMidPointP == P)
473 return nextMidPointQ[1];
475 return nextMidPointP[1];
478 if ((coneP +1) == itEnd) // and thus (coneQ+1)=itEnd
480 //Return highest valid cone abscissa
481 if ((midPointP != P) && (midPointQ !=Q)) //P on
482 return (midPointQ[1]> midPointP[1])? midPointQ[1] : midPointP[1];
486 //midPointP===P trick
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
499 return nextMidPointP[1];
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] ) ));
507 //Symmetry w.r.t. L if necessary
510 Abscissa tmp = normalP[0];
511 normalP[0] = normalP[1];
515 normalQ[0] = normalQ[1];
518 if (P[1] > startingPoint[1])
520 if (Q[1] > startingPoint[1])
523 if (normalP[0] == normalQ[0])
525 //same cone, Voro edge should be on cone extermities
526 RawValue dpm = this->rawDistance(P, midPointP);
527 RawValue dqm = this->rawDistance(Q, midPointQ);
530 else return nextMidPointP[0];
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]) ));
537///////////////////////////////////////////////////////////////////////////////
538template <typename TSpace>
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
545 ASSERT( u[aDimension] < v[aDimension]);
546 ASSERT( u[aDimension] < w[aDimension]);
547 ASSERT( v[aDimension] < w[aDimension]);
549 Abscissa x_uv = getLowerVoronoiEdgeAbscissa(u, v, startingPoint, endPoint, aDimension);
551 if ((x_uv > endPoint[aDimension]))
557 Abscissa x_vw = getLowerVoronoiEdgeAbscissa(v, w, startingPoint, endPoint, aDimension);
558 if ((x_vw < startingPoint[aDimension]))
564 return (x_uv > x_vw);
566///////////////////////////////////////////////////////////////////////////////
567// Implementation of inline functions //
568template <typename TSpace>
571DGtal::operator<< ( std::ostream & out,
572 const DGtal::experimental::ChamferNorm2D<TSpace> & object )
574 object.selfDisplay( out );
578///////////////////////////////////////////////////////////////////////////////