DGtal  1.5.beta
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  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 // ------------------------------------------------------------------------
187 template < typename TPredicate >
188 inline
189 int 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 // ------------------------------------------------------------------------
214 template < typename TPredicate >
215 inline
216 std::pair<typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay,
217  typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay>
218 DGtal::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 // ------------------------------------------------------------------------
232 template < typename TPredicate >
233 inline
234 std::vector<typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay>
235 DGtal::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 // ------------------------------------------------------------------------
298 template < typename TPredicate >
299 inline
300 bool
301 DGtal::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 // ------------------------------------------------------------------------
336 template < typename TPredicate >
337 inline
338 typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay
339 DGtal::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 // ------------------------------------------------------------------------
376 template < typename TPredicate >
377 inline
378 typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay
379 DGtal::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 // ------------------------------------------------------------------------
395 template < typename TPredicate >
396 inline
397 typename DGtal::PlaneProbingR1Neighborhood<TPredicate>::PointOnProbingRay
398 DGtal::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  */
433 template <typename TPredicate>
434 inline
435 void
436 DGtal::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  */
445 template <typename TPredicate>
446 inline
447 bool
448 DGtal::PlaneProbingR1Neighborhood<TPredicate>::isValid() const
449 {
450  return true;
451 }
452 
453 
454 
455 ///////////////////////////////////////////////////////////////////////////////
456 // Implementation of inline functions //
457 
458 template <typename TPredicate>
459 inline
460 std::ostream&
461 DGtal::operator<< ( std::ostream & out,
462  const PlaneProbingR1Neighborhood<TPredicate> & object )
463 {
464  object.selfDisplay( out );
465  return out;
466 }
467 
468 // //
469 ///////////////////////////////////////////////////////////////////////////////
470 
471