DGtal 1.4.2
Loading...
Searching...
No Matches
PlaneProbingR1Neighborhood.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
19 * @author Jocelyn Meyron (\c jocelyn.meyron@liris.cnrs.fr )
20 * Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
21 *
22 * @date 2020/12/04
23 *
24 * Implementation of inline methods defined in PlaneProbingR1Neighborhood.h
25 *
26 * This file is part of the DGtal library.
27 */
28
29
30//////////////////////////////////////////////////////////////////////////////
31#include <cstdlib>
32//////////////////////////////////////////////////////////////////////////////
33
34///////////////////////////////////////////////////////////////////////////////
35// IMPLEMENTATION of inline methods.
36///////////////////////////////////////////////////////////////////////////////
37
38///////////////////////////////////////////////////////////////////////////////
39// ----------------------- Standard services ------------------------------
40
41// ------------------------------------------------------------------------
42template < typename TPredicate >
43inline
44DGtal::PlaneProbingR1Neighborhood<TPredicate>::
45PlaneProbingR1Neighborhood(Predicate const& aPredicate, Point const& aQ, Triangle const& aM)
46 : DGtal::PlaneProbingRNeighborhood<TPredicate>(aPredicate, aQ, aM)
47{}
48
49// ------------------------------------------------------------------------
50template < typename TPredicate >
51inline
52DGtal::PlaneProbingR1Neighborhood<TPredicate>::
53~PlaneProbingR1Neighborhood()
54{}
55
56///////////////////////////////////////////////////////////////////////////////
57// ----------------------- Plane Probing services ------------------------------
58
59// ------------------------------------------------------------------------
60template < typename TPredicate >
61inline
62typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::HexagonState
63DGtal::PlaneProbingR1Neighborhood<TPredicate>::hexagonState ()
64{
65 this->myCandidates.clear();
66
67 int code = getNeighborhoodCode();
68
69 switch (code) {
70 // One vertex of the H-neighborhood is in P
71 case 1:
72 this->myCandidates.push_back(PointOnProbingRay({ 0, 1, 2 }, 0));
73 break;
74
75 case 2:
76 this->myCandidates.push_back(PointOnProbingRay({ 0, 2, 1 }, 0));
77 break;
78
79 case 4:
80 this->myCandidates.push_back(PointOnProbingRay({ 1, 2, 0 }, 0));
81 break;
82
83 case 8:
84 this->myCandidates.push_back(PointOnProbingRay({ 1, 0, 2 }, 0));
85 break;
86
87 case 16:
88 this->myCandidates.push_back(PointOnProbingRay({ 2, 0, 1 }, 0));
89 break;
90
91 case 32:
92 this->myCandidates.push_back(PointOnProbingRay({ 2, 1, 0 }, 0));
93 break;
94
95 // Two vertices of the H-neighborhood are in P
96 case 6:
97 this->myCandidates.push_back(PointOnProbingRay({ 0, 2, 1 }, 0));
98 this->myCandidates.push_back(PointOnProbingRay({ 1, 2, 0 }, 0));
99 break;
100
101 case 24:
102 this->myCandidates.push_back(PointOnProbingRay({ 1, 0, 2 }, 0));
103 this->myCandidates.push_back(PointOnProbingRay({ 2, 0, 1 }, 0));
104 break;
105
106 case 33:
107 this->myCandidates.push_back(PointOnProbingRay({ 2, 1, 0 }, 0));
108 this->myCandidates.push_back(PointOnProbingRay({ 0, 1, 2 }, 0));
109 break;
110
111 case 3:
112 this->myCandidates.push_back(closestRayPoint(candidateRay(0)));
113 break;
114
115 case 12:
116 this->myCandidates.push_back(closestRayPoint(candidateRay(1)));
117 break;
118
119 case 48:
120 this->myCandidates.push_back(closestRayPoint(candidateRay(2)));
121 break;
122
123 // Three vertices of the H-neighborhood are in P
124 case 35:
125 {
126 std::pair<PointOnProbingRay, PointOnProbingRay> candidate = candidateRay(0);
127 std::vector<PointOnProbingRay> shortList = { candidate.second, PointOnProbingRay({ 2, 1, 0 }, 0) };
128 PointOnProbingRay closest = this->closestPointInList(shortList);
129 //PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 2, 1, 0 }, 0) });
130 this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest)));
131 }
132 break;
133
134 case 7:
135 {
136 std::pair<PointOnProbingRay, PointOnProbingRay> candidate = candidateRay(0);
137 std::vector<PointOnProbingRay> shortList = { candidate.second, PointOnProbingRay({ 1, 2, 0 }, 0) };
138 PointOnProbingRay closest = this->closestPointInList(shortList);
139 this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest)));
140 }
141 break;
142
143 case 14:
144 {
145 std::pair<PointOnProbingRay, PointOnProbingRay> candidate = candidateRay(1);
146 std::vector<PointOnProbingRay> shortList = { candidate.second, PointOnProbingRay({ 0, 2, 1 }, 0) };
147 PointOnProbingRay closest = this->closestPointInList(shortList);
148 this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest)));
149 }
150 break;
151
152 case 28:
153 {
154 std::pair<PointOnProbingRay, PointOnProbingRay> candidate = candidateRay(1);
155 std::vector<PointOnProbingRay> shortList = { candidate.second, PointOnProbingRay({ 2, 0, 1 }, 0) };
156 PointOnProbingRay closest = this->closestPointInList(shortList);
157 this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest)));
158 }
159 break;
160
161 case 56:
162 {
163 std::pair<PointOnProbingRay, PointOnProbingRay> candidate = candidateRay(2);
164 std::vector<PointOnProbingRay> shortList = { candidate.second, PointOnProbingRay({ 1, 0, 2 }, 0) };
165 PointOnProbingRay closest = this->closestPointInList(shortList);
166 this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest)));
167 }
168 break;
169
170 case 49:
171 {
172 std::pair<PointOnProbingRay, PointOnProbingRay> candidate = candidateRay(2);
173 std::vector<PointOnProbingRay> shortList = { candidate.second, PointOnProbingRay({ 0, 1, 2 }, 0) };
174 PointOnProbingRay closest = this->closestPointInList(shortList);
175 this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest)));
176 }
177 break;
178
179 default:
180 break;
181 }
182
183 return this->classify(myState);
184}
185
186// ------------------------------------------------------------------------
187template < typename TPredicate >
188inline
189int DGtal::PlaneProbingR1Neighborhood<TPredicate>::getNeighborhoodCode () const
190{
191 int code = 0;
192
193 myState = { false, false, false, false, false, false };
194 for (int i = 0; i < 6; ++i) {
195 PointOnProbingRay const& r = this->myNeighborhood[i];
196
197 // We skip the ray if it is not part of the rays that should be considered at this step
198 if (! this->isNeighbor(r))
199 {
200 continue;
201 }
202
203 myState[i] = this->myPredicate(this->absolutePoint(r));
204 if (myState[i])
205 {
206 code += (1 << i);
207 }
208 }
209
210 return code;
211}
212
213// ------------------------------------------------------------------------
214template < typename TPredicate >
215inline
216std::pair<typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay,
217 typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay>
218DGtal::PlaneProbingR1Neighborhood<TPredicate>::candidateRay (Index const& index) const
219{
220 Index indexM1 = (index + 2) % 3, indexM2 = (index + 1) % 3;
221 PointOnProbingRay r1({ index, indexM1, indexM2 }, 0),
222 r2({ index, indexM2, indexM1 }, 0);
223
224 if (detail::squaredNorm(this->myM[indexM1]) >= detail::squaredNorm(this->myM[indexM2])) {
225 return std::make_pair(r1, r2.getBase());
226 } else {
227 return std::make_pair(r2, r1.getBase());
228 }
229}
230
231// ------------------------------------------------------------------------
232template < typename TPredicate >
233inline
234std::vector<typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay>
235DGtal::PlaneProbingR1Neighborhood<TPredicate>::intersectSphereRay (PointOnProbingRay const& aPoint, PointOnProbingRay const& aRay) const
236{
237 using NumberTraits = DGtal::NumberTraits<Integer>;
238
239 PointOnProbingRay startingPoint = aRay.getBase();
240 assert(this->myPredicate(this->absolutePoint(startingPoint)));
241
242 Point origin = -this->myM[aRay.sigma(0)],
243 y = this->myM[aRay.sigma(1)],
244 x = this->myM[aRay.sigma(2)];
245
246 Integer a = detail::squaredNorm(x),
247 b = detail::distToSphere<Point>({ -this->myM[0], -this->myM[1], -this->myM[2],
248 this->relativePoint(aPoint), origin + x }) - a + 2*x.dot(y),
249 c = detail::distToSphere<Point>({ -this->myM[0], -this->myM[1], -this->myM[2],
250 this->relativePoint(aPoint), origin + y });
251 Integer delta = b*b - 4*a*c;
252
253 std::vector<PointOnProbingRay> res;
254
255 if (delta < 0) {
256 return res;
257 }
258
259 Integer i1 = std::ceil(double(-b - sqrt(delta)) / (2*a) ),
260 i2 = std::floor(double(-b + sqrt(delta)) / (2*a) );
261 Integer zero = NumberTraits::ZERO;
262 if (i1 <= i2 && i2 >= zero) {
263 i1 = std::max(zero, i1);
264
265 PointOnProbingRay p1(aRay.sigma(), i1), p2(aRay.sigma(), i2);
266
267 // Case of cospherical points
268 if (detail::distToSphere<Point>({ -this->myM[0], -this->myM[1], -this->myM[2],
269 this->relativePoint(aPoint) , this->relativePoint(p1) }) == 0) {
270 if (this->absolutePoint(p1) >= this->absolutePoint(aPoint)) {
271 i1++;
272 }
273 }
274
275 if (detail::distToSphere<Point>({ -this->myM[0], -this->myM[1], -this->myM[2],
276 this->relativePoint(aPoint) , this->relativePoint(p2) }) == 0) {
277 if (this->absolutePoint(p2) >= this->absolutePoint(aPoint)) {
278 i2--;
279 }
280 }
281
282 if (i1 > i2) {
283 return res;
284 }
285
286 p1 = PointOnProbingRay(aRay.sigma(), i1);
287 p2 = PointOnProbingRay(aRay.sigma(), i2);
288
289 // Add extremal points to the list
290 res.push_back(p1);
291 res.push_back(p2);
292 }
293
294 return res;
295}
296
297// ------------------------------------------------------------------------
298template < typename TPredicate >
299inline
300bool
301DGtal::PlaneProbingR1Neighborhood<TPredicate>::isValidIntersectSphereRay (PointOnProbingRay const& aPoint, PointOnProbingRay const& aRay,
302 std::vector<PointOnProbingRay> const& aLst) const
303{
304 PointOnProbingRay startingPoint = aRay.getBase();
305 assert(this->myPredicate(this->absolutePoint(startingPoint)));
306
307 Point refPoint = this->relativePoint(aPoint);
308 bool res = true;
309
310 if (aLst.size() > 0) {
311 PointOnProbingRay currentX = startingPoint;
312 while (! this->isSmallest(refPoint, this->relativePoint(currentX))) {
313 currentX = currentX.next(1);
314 }
315
316 PointOnProbingRay first = currentX;
317 while (this->isSmallest(refPoint, this->relativePoint(currentX))) {
318 currentX = currentX.next(1);
319 }
320
321 PointOnProbingRay last = currentX.previous(1);
322 return (aLst.size() == 2) && (first == aLst[0]) && (last == aLst[1]);
323 } else {
324 int n = 15;
325 PointOnProbingRay currentX = startingPoint;
326 while (res && currentX.position() < n) {
327 res = !this->isSmallest(refPoint, this->relativePoint(currentX));
328 currentX = currentX.next(1);
329 }
330 }
331
332 return res;
333}
334
335// ------------------------------------------------------------------------
336template < typename TPredicate >
337inline
338typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay
339DGtal::PlaneProbingR1Neighborhood<TPredicate>::closestPointOnRayConstant (PointOnProbingRay const& aRay) const
340{
341 using NumberTraits = DGtal::NumberTraits<Integer>;
342
343 PointOnProbingRay startingPoint = aRay.getBase();
344 assert(this->myPredicate(this->absolutePoint(startingPoint)));
345
346 Point origin = -this->myM[aRay.sigma(0)],
347 u = this->myM[aRay.sigma(1)],
348 v = this->myM[aRay.sigma(2)];
349
350 Integer a = detail::squaredNorm(v), b = 3*a,
351 c = 2 * u.dot(v) + detail::distToSphere<Point>({ -this->myM[0], -this->myM[1], -this->myM[2],
352 origin + u , origin + v });
353 Integer delta = b*b - 4*a*c;
354
355 if (delta >= 0) {
356 Integer index = std::ceil(double(-b + sqrt(delta)) / (2*a));
357 index = std::max(NumberTraits::ZERO, index);
358
359 // Case of cospherical points
360 PointOnProbingRay p1(aRay.sigma(), index), p2(aRay.sigma(), index+1);
361
362 if (detail::distToSphere<Point>({ -this->myM[0], -this->myM[1], -this->myM[2],
363 this->relativePoint(p1), this->relativePoint(p2) }) == 0) {
364 if (this->absolutePoint(p1) >= this->absolutePoint(p2)) {
365 index++;
366 }
367 }
368
369 return PointOnProbingRay(aRay.sigma(), index);
370 } else {
371 return aRay.getBase();
372 }
373}
374
375// ------------------------------------------------------------------------
376template < typename TPredicate >
377inline
378typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay
379DGtal::PlaneProbingR1Neighborhood<TPredicate>::closestPointOnRayLinear (PointOnProbingRay const& aRay) const
380{
381 PointOnProbingRay startingPoint = aRay.getBase();
382 assert(this->myPredicate(this->absolutePoint(startingPoint)));
383
384 PointOnProbingRay previousX = startingPoint;
385 PointOnProbingRay currentX = previousX.next(1);
386 while (this->isSmallest(this->relativePoint(previousX), this->relativePoint(currentX))) {
387 previousX = currentX;
388 currentX = previousX.next(1);
389 }
390
391 return previousX;
392}
393
394// ------------------------------------------------------------------------
395template < typename TPredicate >
396inline
397typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay
398DGtal::PlaneProbingR1Neighborhood<TPredicate>::closestRayPoint (std::pair<PointOnProbingRay, PointOnProbingRay> const& aRayPoint) const
399{
400 PointOnProbingRay const& ray = aRayPoint.first;
401 PointOnProbingRay const& point = aRayPoint.second;
402
403 PointOnProbingRay res = point;
404 std::vector<PointOnProbingRay> intersections = intersectSphereRay(point, ray);
405 assert(isValidIntersectSphereRay(point, ray, intersections));
406
407 if (intersections.size() > 0) {
408 PointOnProbingRay y1 = intersections[0], y2 = intersections[1];
409 if (this->myPredicate(this->absolutePoint(y1))) {
410 PointOnProbingRay y = closestPointOnRayConstant(ray);
411 assert(y == closestPointOnRayLinear(ray));
412 if (y1 <= y && y <= y2) {
413 if (this->myPredicate(this->absolutePoint(y))) {
414 res = y;
415 } else {
416 res = this->closestPointOnRayLogWithPredicate(y1);
417 assert(res == this->closestPointOnRayLinearWithPredicate(ray.getBase()));
418 }
419 }
420 }
421 }
422
423 return res;
424}
425
426///////////////////////////////////////////////////////////////////////////////
427// Interface - public :
428
429/**
430 * Writes/Displays the object on an output stream.
431 * @param out the output stream where the object is written.
432 */
433template <typename TPredicate>
434inline
435void
436DGtal::PlaneProbingR1Neighborhood<TPredicate>::selfDisplay ( std::ostream & out ) const
437{
438 out << "[PlaneProbingR1Neighborhood]";
439}
440
441/**
442 * Checks the validity/consistency of the object.
443 * @return 'true' if the object is valid, 'false' otherwise.
444 */
445template <typename TPredicate>
446inline
447bool
448DGtal::PlaneProbingR1Neighborhood<TPredicate>::isValid() const
449{
450 return true;
451}
452
453
454
455///////////////////////////////////////////////////////////////////////////////
456// Implementation of inline functions //
457
458template <typename TPredicate>
459inline
460std::ostream&
461DGtal::operator<< ( std::ostream & out,
462 const PlaneProbingR1Neighborhood<TPredicate> & object )
463{
464 object.selfDisplay( out );
465 return out;
466}
467
468// //
469///////////////////////////////////////////////////////////////////////////////
470
471