DGtal 1.3.0
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 PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 2, 1, 0 }, 0) });
128 this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest)));
129 }
130 break;
131
132 case 7:
133 {
134 std::pair<PointOnProbingRay, PointOnProbingRay> candidate = candidateRay(0);
135 PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 1, 2, 0 }, 0) });
136 this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest)));
137 }
138 break;
139
140 case 14:
141 {
142 std::pair<PointOnProbingRay, PointOnProbingRay> candidate = candidateRay(1);
143 PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 0, 2, 1 }, 0) });
144 this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest)));
145 }
146 break;
147
148 case 28:
149 {
150 std::pair<PointOnProbingRay, PointOnProbingRay> candidate = candidateRay(1);
151 PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 2, 0, 1 }, 0) });
152 this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest)));
153 }
154 break;
155
156 case 56:
157 {
158 std::pair<PointOnProbingRay, PointOnProbingRay> candidate = candidateRay(2);
159 PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 1, 0, 2 }, 0) });
160 this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest)));
161 }
162 break;
163
164 case 49:
165 {
166 std::pair<PointOnProbingRay, PointOnProbingRay> candidate = candidateRay(2);
167 PointOnProbingRay closest = this->closestPointInList({ candidate.second, PointOnProbingRay({ 0, 1, 2 }, 0) });
168 this->myCandidates.push_back(closestRayPoint(std::make_pair(candidate.first, closest)));
169 }
170 break;
171
172 default:
173 break;
174 }
175
176 return this->classify(myState);
177}
178
179// ------------------------------------------------------------------------
180template < typename TPredicate >
181inline
182int DGtal::PlaneProbingR1Neighborhood<TPredicate>::getNeighborhoodCode () const
183{
184 int code = 0;
185
186 myState = { false, false, false, false, false, false };
187 for (int i = 0; i < 6; ++i) {
188 PointOnProbingRay const& r = this->myNeighborhood[i];
189
190 // We skip the ray if it is not part of the rays that should be considered at this step
191 if (! this->isNeighbor(r))
192 {
193 continue;
194 }
195
196 myState[i] = this->myPredicate(this->absolutePoint(r));
197 if (myState[i])
198 {
199 code += (1 << i);
200 }
201 }
202
203 return code;
204}
205
206// ------------------------------------------------------------------------
207template < typename TPredicate >
208inline
209std::pair<typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay,
210 typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay>
211DGtal::PlaneProbingR1Neighborhood<TPredicate>::candidateRay (int index) const
212{
213 int indexM1 = (index - 1 + 3) % 3, indexM2 = (index - 2 + 3) % 3;
214 PointOnProbingRay r1({ index, indexM1, indexM2 }, 0),
215 r2({ index, indexM2, indexM1 }, 0);
216
217 if (detail::squaredNorm(this->myM[indexM1]) >= detail::squaredNorm(this->myM[indexM2])) {
218 return std::make_pair(r1, r2.getBase());
219 } else {
220 return std::make_pair(r2, r1.getBase());
221 }
222}
223
224// ------------------------------------------------------------------------
225template < typename TPredicate >
226inline
227std::vector<typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay>
228DGtal::PlaneProbingR1Neighborhood<TPredicate>::intersectSphereRay (PointOnProbingRay const& aPoint, PointOnProbingRay const& aRay) const
229{
230 using NumberTraits = DGtal::NumberTraits<Integer>;
231
232 PointOnProbingRay startingPoint = aRay.getBase();
233 assert(this->myPredicate(this->absolutePoint(startingPoint)));
234
235 Point origin = -this->myM[aRay.sigma(0)],
236 y = this->myM[aRay.sigma(1)],
237 x = this->myM[aRay.sigma(2)];
238
239 Integer a = detail::squaredNorm(x),
240 b = detail::distToSphere<Point>({ -this->myM[0], -this->myM[1], -this->myM[2],
241 this->relativePoint(aPoint), origin + x }) - a + 2*x.dot(y),
242 c = detail::distToSphere<Point>({ -this->myM[0], -this->myM[1], -this->myM[2],
243 this->relativePoint(aPoint), origin + y });
244 Integer delta = b*b - 4*a*c;
245
246 std::vector<PointOnProbingRay> res;
247
248 if (delta < 0) {
249 return res;
250 }
251
252 Integer i1 = std::ceil(double(-b - sqrt(delta)) / (2*a) ),
253 i2 = std::floor(double(-b + sqrt(delta)) / (2*a) );
254 Integer zero = NumberTraits::ZERO;
255 if (i1 <= i2 && i2 >= zero) {
256 i1 = std::max(zero, i1);
257
258 PointOnProbingRay p1(aRay.sigma(), i1), p2(aRay.sigma(), i2);
259
260 // Case of cospherical points
261 if (detail::distToSphere<Point>({ -this->myM[0], -this->myM[1], -this->myM[2],
262 this->relativePoint(aPoint) , this->relativePoint(p1) }) == 0) {
263 if (this->absolutePoint(p1) >= this->absolutePoint(aPoint)) {
264 i1++;
265 }
266 }
267
268 if (detail::distToSphere<Point>({ -this->myM[0], -this->myM[1], -this->myM[2],
269 this->relativePoint(aPoint) , this->relativePoint(p2) }) == 0) {
270 if (this->absolutePoint(p2) >= this->absolutePoint(aPoint)) {
271 i2--;
272 }
273 }
274
275 if (i1 > i2) {
276 return res;
277 }
278
279 p1 = PointOnProbingRay(aRay.sigma(), i1);
280 p2 = PointOnProbingRay(aRay.sigma(), i2);
281
282 // Add extremal points to the list
283 res.push_back(p1);
284 res.push_back(p2);
285 }
286
287 return res;
288}
289
290// ------------------------------------------------------------------------
291template < typename TPredicate >
292inline
293bool
294DGtal::PlaneProbingR1Neighborhood<TPredicate>::isValidIntersectSphereRay (PointOnProbingRay const& aPoint, PointOnProbingRay const& aRay,
295 std::vector<PointOnProbingRay> const& aLst) const
296{
297 PointOnProbingRay startingPoint = aRay.getBase();
298 assert(this->myPredicate(this->absolutePoint(startingPoint)));
299
300 Point refPoint = this->relativePoint(aPoint);
301 bool res = true;
302
303 if (aLst.size() > 0) {
304 PointOnProbingRay currentX = startingPoint;
305 while (! this->isSmallest(refPoint, this->relativePoint(currentX))) {
306 currentX = currentX.next(1);
307 }
308
309 PointOnProbingRay first = currentX;
310 while (this->isSmallest(refPoint, this->relativePoint(currentX))) {
311 currentX = currentX.next(1);
312 }
313
314 PointOnProbingRay last = currentX.previous(1);
315 return (aLst.size() == 2) && (first == aLst[0]) && (last == aLst[1]);
316 } else {
317 int n = 15;
318 PointOnProbingRay currentX = startingPoint;
319 while (res && currentX.index() < n) {
320 res = !this->isSmallest(refPoint, this->relativePoint(currentX));
321 currentX = currentX.next(1);
322 }
323 }
324
325 return res;
326}
327
328// ------------------------------------------------------------------------
329template < typename TPredicate >
330inline
331typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay
332DGtal::PlaneProbingR1Neighborhood<TPredicate>::closestPointOnRayConstant (PointOnProbingRay const& aRay) const
333{
334 using NumberTraits = DGtal::NumberTraits<Integer>;
335
336 PointOnProbingRay startingPoint = aRay.getBase();
337 assert(this->myPredicate(this->absolutePoint(startingPoint)));
338
339 Point origin = -this->myM[aRay.sigma(0)],
340 u = this->myM[aRay.sigma(1)],
341 v = this->myM[aRay.sigma(2)];
342
343 Integer a = detail::squaredNorm(v), b = 3*a,
344 c = 2 * u.dot(v) + detail::distToSphere<Point>({ -this->myM[0], -this->myM[1], -this->myM[2],
345 origin + u , origin + v });
346 Integer delta = b*b - 4*a*c;
347
348 if (delta >= 0) {
349 Integer index = std::ceil(double(-b + sqrt(delta)) / (2*a));
350 index = std::max(NumberTraits::ZERO, index);
351
352 // Case of cospherical points
353 PointOnProbingRay p1(aRay.sigma(), index), p2(aRay.sigma(), index+1);
354
355 if (detail::distToSphere<Point>({ -this->myM[0], -this->myM[1], -this->myM[2],
356 this->relativePoint(p1), this->relativePoint(p2) }) == 0) {
357 if (this->absolutePoint(p1) >= this->absolutePoint(p2)) {
358 index++;
359 }
360 }
361
362 return PointOnProbingRay(aRay.sigma(), index);
363 } else {
364 return aRay.getBase();
365 }
366}
367
368// ------------------------------------------------------------------------
369template < typename TPredicate >
370inline
371typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay
372DGtal::PlaneProbingR1Neighborhood<TPredicate>::closestPointOnRayLinear (PointOnProbingRay const& aRay) const
373{
374 PointOnProbingRay startingPoint = aRay.getBase();
375 assert(this->myPredicate(this->absolutePoint(startingPoint)));
376
377 PointOnProbingRay previousX = startingPoint;
378 PointOnProbingRay currentX = previousX.next(1);
379 while (this->isSmallest(this->relativePoint(previousX), this->relativePoint(currentX))) {
380 previousX = currentX;
381 currentX = previousX.next(1);
382 }
383
384 return previousX;
385}
386
387// ------------------------------------------------------------------------
388template < typename TPredicate >
389inline
390typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay
391DGtal::PlaneProbingR1Neighborhood<TPredicate>::closestRayPoint (std::pair<PointOnProbingRay, PointOnProbingRay> const& aRayPoint) const
392{
393 PointOnProbingRay const& ray = aRayPoint.first;
394 PointOnProbingRay const& point = aRayPoint.second;
395
396 PointOnProbingRay res = point;
397 std::vector<PointOnProbingRay> intersections = intersectSphereRay(point, ray);
398 assert(isValidIntersectSphereRay(point, ray, intersections));
399
400 if (intersections.size() > 0) {
401 PointOnProbingRay y1 = intersections[0], y2 = intersections[1];
402 if (this->myPredicate(this->absolutePoint(y1))) {
403 PointOnProbingRay y = closestPointOnRayConstant(ray);
404 assert(y == closestPointOnRayLinear(ray));
405 if (y1 <= y && y <= y2) {
406 if (this->myPredicate(this->absolutePoint(y))) {
407 res = y;
408 } else {
409 res = this->closestPointOnRayLogWithPredicate(y1);
410 assert(res == this->closestPointOnRayLinearWithPredicate(ray.getBase()));
411 }
412 }
413 }
414 }
415
416 return res;
417}
418
419///////////////////////////////////////////////////////////////////////////////
420// Interface - public :
421
422/**
423 * Writes/Displays the object on an output stream.
424 * @param out the output stream where the object is written.
425 */
426template <typename TPredicate>
427inline
428void
429DGtal::PlaneProbingR1Neighborhood<TPredicate>::selfDisplay ( std::ostream & out ) const
430{
431 out << "[PlaneProbingR1Neighborhood]";
432}
433
434/**
435 * Checks the validity/consistency of the object.
436 * @return 'true' if the object is valid, 'false' otherwise.
437 */
438template <typename TPredicate>
439inline
440bool
441DGtal::PlaneProbingR1Neighborhood<TPredicate>::isValid() const
442{
443 return true;
444}
445
446
447
448///////////////////////////////////////////////////////////////////////////////
449// Implementation of inline functions //
450
451template <typename TPredicate>
452inline
453std::ostream&
454DGtal::operator<< ( std::ostream & out,
455 const PlaneProbingR1Neighborhood<TPredicate> & object )
456{
457 object.selfDisplay( out );
458 return out;
459}
460
461// //
462///////////////////////////////////////////////////////////////////////////////
463
464
Generated on Sun Nov 27 2022 15:15:39 for DGtal by doxygen 1.9.5