DGtal  1.2.0
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 // ------------------------------------------------------------------------
42 template < typename TPredicate >
43 inline
44 DGtal::PlaneProbingR1Neighborhood<TPredicate>::
45 PlaneProbingR1Neighborhood(Predicate const& aPredicate, Point const& aQ, Triangle const& aM)
46  : DGtal::PlaneProbingRNeighborhood<TPredicate>(aPredicate, aQ, aM)
47 {}
48 
49 // ------------------------------------------------------------------------
50 template < typename TPredicate >
51 inline
52 DGtal::PlaneProbingR1Neighborhood<TPredicate>::
53 ~PlaneProbingR1Neighborhood()
54 {}
55 
56 ///////////////////////////////////////////////////////////////////////////////
57 // ----------------------- Plane Probing services ------------------------------
58 
59 // ------------------------------------------------------------------------
60 template < typename TPredicate >
61 inline
62 typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::HexagonState
63 DGtal::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 // ------------------------------------------------------------------------
180 template < typename TPredicate >
181 inline
182 int 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 // ------------------------------------------------------------------------
207 template < typename TPredicate >
208 inline
209 std::pair<typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay,
210  typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay>
211 DGtal::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 // ------------------------------------------------------------------------
225 template < typename TPredicate >
226 inline
227 std::vector<typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay>
228 DGtal::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 // ------------------------------------------------------------------------
291 template < typename TPredicate >
292 inline
293 bool
294 DGtal::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 // ------------------------------------------------------------------------
329 template < typename TPredicate >
330 inline
331 typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay
332 DGtal::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 // ------------------------------------------------------------------------
369 template < typename TPredicate >
370 inline
371 typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay
372 DGtal::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 // ------------------------------------------------------------------------
388 template < typename TPredicate >
389 inline
390 typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay
391 DGtal::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  */
426 template <typename TPredicate>
427 inline
428 void
429 DGtal::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  */
438 template <typename TPredicate>
439 inline
440 bool
441 DGtal::PlaneProbingR1Neighborhood<TPredicate>::isValid() const
442 {
443  return true;
444 }
445 
446 
447 
448 ///////////////////////////////////////////////////////////////////////////////
449 // Implementation of inline functions //
450 
451 template <typename TPredicate>
452 inline
453 std::ostream&
454 DGtal::operator<< ( std::ostream & out,
455  const PlaneProbingR1Neighborhood<TPredicate> & object )
456 {
457  object.selfDisplay( out );
458  return out;
459 }
460 
461 // //
462 ///////////////////////////////////////////////////////////////////////////////
463 
464