DGtal  1.5.beta
PlaneProbingTetrahedronEstimator.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 PlaneProbingTetrahedronEstimator.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cassert>
32 #include <cstdlib>
33 #include <array>
34 #include <vector>
35 #include "DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h"
36 #include "DGtal/geometry/surfaces/estimation/PlaneProbingHNeighborhood.h"
37 #include "DGtal/geometry/surfaces/estimation/PlaneProbingRNeighborhood.h"
38 #include "DGtal/geometry/surfaces/estimation/PlaneProbingR1Neighborhood.h"
39 #include "DGtal/geometry/surfaces/estimation/PlaneProbingLNeighborhood.h"
40 //////////////////////////////////////////////////////////////////////////////
41 
42 ///////////////////////////////////////////////////////////////////////////////
43 // IMPLEMENTATION of inline methods.
44 ///////////////////////////////////////////////////////////////////////////////
45 
46 namespace DGtal
47 {
48  namespace detail
49  {
50  // Helper class to choose the PlaneProbingNeighborhood class at compile-time
51  // (used in the constructor of PlaneProbingTetrahedronEstimator)
52  template < typename TPredicate, DGtal::ProbingMode mode >
53  struct PlaneProbingNeighborhoodSelector
54  {
55  static DGtal::PlaneProbingNeighborhood<TPredicate>*
56  select(TPredicate const&,
57  typename DGtal::PlaneProbingNeighborhood<TPredicate>::Point const&,
58  typename DGtal::PlaneProbingNeighborhood<TPredicate>::Triangle const&);
59  };
60 
61  template < typename TPredicate >
62  struct PlaneProbingNeighborhoodSelector<TPredicate, DGtal::ProbingMode::H>
63  {
64  static DGtal::PlaneProbingNeighborhood<TPredicate>*
65  select(TPredicate const& aPredicate,
66  typename DGtal::PlaneProbingNeighborhood<TPredicate>::Point const& aQ,
67  typename DGtal::PlaneProbingNeighborhood<TPredicate>::Triangle const& aM)
68  {
69  return new DGtal::PlaneProbingHNeighborhood<TPredicate>(aPredicate, aQ, aM);
70  }
71  };
72 
73  template < typename TPredicate >
74  struct PlaneProbingNeighborhoodSelector<TPredicate, DGtal::ProbingMode::R>
75  {
76  static DGtal::PlaneProbingNeighborhood<TPredicate>*
77  select(TPredicate const& aPredicate,
78  typename DGtal::PlaneProbingNeighborhood<TPredicate>::Point const& aQ,
79  typename DGtal::PlaneProbingNeighborhood<TPredicate>::Triangle const& aM)
80  {
81  return new DGtal::PlaneProbingRNeighborhood<TPredicate>(aPredicate, aQ, aM);
82  }
83  };
84 
85  template < typename TPredicate >
86  struct PlaneProbingNeighborhoodSelector<TPredicate, DGtal::ProbingMode::R1>
87  {
88  static DGtal::PlaneProbingNeighborhood<TPredicate>*
89  select(TPredicate const& aPredicate,
90  typename DGtal::PlaneProbingNeighborhood<TPredicate>::Point const& aQ,
91  typename DGtal::PlaneProbingNeighborhood<TPredicate>::Triangle const& aM)
92  {
93  return new DGtal::PlaneProbingR1Neighborhood<TPredicate>(aPredicate, aQ, aM);
94  }
95  };
96 
97  template < typename TPredicate >
98  struct PlaneProbingNeighborhoodSelector<TPredicate, DGtal::ProbingMode::L>
99  {
100  static DGtal::PlaneProbingNeighborhood<TPredicate>*
101  select(TPredicate const& aPredicate,
102  typename DGtal::PlaneProbingNeighborhood<TPredicate>::Point const& aQ,
103  typename DGtal::PlaneProbingNeighborhood<TPredicate>::Triangle const& aM)
104  {
105  return new DGtal::PlaneProbingLNeighborhood<TPredicate>(aPredicate, aQ, aM);
106  }
107  };
108 
109  } // namespace detail
110 } // namespace DGtal
111 
112 ///////////////////////////////////////////////////////////////////////////////
113 // ----------------------- Standard services ------------------------------
114 
115 // ------------------------------------------------------------------------
116 template < typename TPredicate, DGtal::ProbingMode mode >
117 inline
118 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::
119 PlaneProbingTetrahedronEstimator (Point const& aPoint, Triangle const& aM, Predicate const& aPredicate)
120  : myM(aM), myPredicate(aPredicate), myS(aM[0] + aM[1] + aM[2]), myQ(aPoint + myS)
121 {
122  myNeighborhood = DGtal::detail::PlaneProbingNeighborhoodSelector<TPredicate, mode>::select(myPredicate, myQ, myM);
123 }
124 
125 // ------------------------------------------------------------------------
126 template < typename TPredicate, DGtal::ProbingMode mode >
127 inline
128 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::
129 ~PlaneProbingTetrahedronEstimator ()
130 {
131  delete myNeighborhood;
132  myNeighborhood = nullptr;
133 }
134 
135 ///////////////////////////////////////////////////////////////////////////////
136 // ----------------------- Plane probing services ------------------------------
137 
138 // ------------------------------------------------------------------------
139 template < typename TPredicate, DGtal::ProbingMode mode >
140 inline
141 typename DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::Vector const&
142 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::m (int aIndex) const
143 {
144  assert(aIndex == 0 || aIndex == 1 || aIndex == 2);
145  return myM[aIndex];
146 }
147 
148 // ------------------------------------------------------------------------
149 template < typename TPredicate, DGtal::ProbingMode mode >
150 inline
151 typename DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::Point const&
152 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::q () const
153 {
154  return myQ;
155 }
156 
157 // ------------------------------------------------------------------------
158 template < typename TPredicate, DGtal::ProbingMode mode >
159 inline
160 typename DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::Point
161 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::getOrigin () const
162 {
163  return myQ - (myM[0] + myM[1] + myM[2]);
164 }
165 
166 // ------------------------------------------------------------------------
167 template < typename TPredicate, DGtal::ProbingMode mode >
168 inline
169 typename DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::Triangle
170 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::vertices () const
171 {
172  return { myQ - myM[0], myQ - myM[1], myQ - myM[2] };
173 }
174 
175 // ------------------------------------------------------------------------
176 template < typename TPredicate, DGtal::ProbingMode mode >
177 inline
178 std::pair<typename DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::Vector,
179  typename DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::Vector>
180 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::getBasis () const
181 {
182  using DGtal::detail::squaredNorm;
183 
184  Vector u = myM[1] - myM[0],
185  v = myM[2] - myM[1],
186  w = myM[0] - myM[2];
187  ASSERT(w == -u - v);
188 
189  if (squaredNorm(u) < squaredNorm(v)) {
190  if (squaredNorm(u) < squaredNorm(w)) {
191  if (squaredNorm(-w) < squaredNorm(v)) {
192  return std::make_pair(u, -w);
193  } else {
194  return std::make_pair(u, v);
195  }
196  } else {
197  if (squaredNorm(-v) < squaredNorm(u)) {
198  return std::make_pair(w, -v);
199  } else {
200  return std::make_pair(w, u);
201  }
202  }
203  } else {
204  if (squaredNorm(v) < squaredNorm(w)) {
205  if (squaredNorm(-u) < squaredNorm(w)) {
206  return std::make_pair(v, -u);
207  } else {
208  return std::make_pair(v, w);
209  }
210  } else {
211  if (squaredNorm(-v) < squaredNorm(u)) {
212  return std::make_pair(w, -v);
213  } else {
214  return std::make_pair(w, u);
215  }
216  }
217  }
218 }
219 
220 // ------------------------------------------------------------------------
221 template < typename TPredicate, DGtal::ProbingMode mode >
222 inline
223 bool
224 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::isReduced () const
225 {
226  auto basis = getBasis();
227  return DGtal::detail::isBasisReduced(basis.first, basis.second);
228 }
229 
230 // ------------------------------------------------------------------------
231 template < typename TPredicate, DGtal::ProbingMode mode >
232 inline
233 typename DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::Vector
234 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::getNormal () const
235 {
236  auto basis = getBasis();
237  return basis.first.crossProduct(basis.second);
238 }
239 
240 // ------------------------------------------------------------------------
241 template < typename TPredicate, DGtal::ProbingMode mode >
242 inline
243 std::pair<bool, typename DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::UpdateOperation>
244 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::advance (std::vector<PointOnProbingRay> const& aNeighbors)
245 {
246  myNeighborhood->setNeighbors(aNeighbors);
247  HexagonState state = hexagonState();
248 
249  if (state == HexagonState::Planar)
250  {
251  UpdateOperation op = myNeighborhood->closestCandidate();
252  update(op);
253  ASSERT(isValid());
254  return { true, op };
255  }
256 
257  return { false, {} };
258 }
259 
260 // ------------------------------------------------------------------------
261 template < typename TPredicate, DGtal::ProbingMode mode >
262 inline
263 std::pair<bool, typename DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::UpdateOperation>
264 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::advance ()
265 {
266  return advance({});
267 }
268 
269 // ------------------------------------------------------------------------
270 template < typename TPredicate, DGtal::ProbingMode mode >
271 inline
272 typename DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::Quantity
273 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::compute (std::vector<PointOnProbingRay> const& aNeighbors)
274 {
275  while (advance(aNeighbors).first) {}
276 
277  return getNormal();
278 }
279 
280 // ------------------------------------------------------------------------
281 template < typename TPredicate, DGtal::ProbingMode mode >
282 inline
283 typename DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::Quantity
284 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::compute ()
285 {
286  return compute({});
287 }
288 
289 // ------------------------------------------------------------------------
290 template < typename TPredicate, DGtal::ProbingMode mode >
291 inline
292 typename DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::Neighborhood::HexagonState
293 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::hexagonState () const
294 {
295  return myNeighborhood->hexagonState();
296 }
297 
298 // ------------------------------------------------------------------------
299 template < typename TPredicate, DGtal::ProbingMode mode >
300 inline
301 void
302 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::translateQ (Vector const& aTranslation)
303 {
304  myQ += aTranslation;
305 }
306 
307 // ------------------------------------------------------------------------
308 template < typename TPredicate, DGtal::ProbingMode mode >
309 inline
310 void
311 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::translateQ ()
312 {
313  UpdateOperation const& lastOp = myOperations[myOperations.size() - 1];
314  Point translation = lastOp.coeffs[1] * myM[lastOp.sigma[1]] + lastOp.coeffs[2] * myM[lastOp.sigma[2]];
315 
316  translateQ(translation);
317 }
318 
319 ///////////////////////////////////////////////////////////////////////////////
320 // ------------------------- Internals ------------------------------------
321 
322 // ------------------------------------------------------------------------
323 template < typename TPredicate, DGtal::ProbingMode mode >
324 inline
325 void
326 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::applyOperation (UpdateOperation const& aOp)
327 {
328  myM[aOp.sigma[0]] = aOp.coeffs[0] * myM[aOp.sigma[0]] + aOp.coeffs[1] * myM[aOp.sigma[1]] + aOp.coeffs[2] * myM[aOp.sigma[2]];
329 }
330 
331 // ------------------------------------------------------------------------
332 template < typename TPredicate, DGtal::ProbingMode mode >
333 inline
334 void
335 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::update (UpdateOperation const& aOp)
336 {
337  applyOperation(aOp);
338  myOperations.push_back(aOp);
339 }
340 
341 ///////////////////////////////////////////////////////////////////////////////
342 // Interface - public :
343 
344 // ------------------------------------------------------------------------
345 template <typename TPredicate, DGtal::ProbingMode mode>
346 inline
347 void
348 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::selfDisplay ( std::ostream & out ) const
349 {
350  out << "[PlaneProbingTetrahedronEstimator]";
351 }
352 
353 // ------------------------------------------------------------------------
354 template <typename TPredicate, DGtal::ProbingMode mode>
355 inline
356 bool
357 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::isValid() const
358 {
359  return ( isUnimodular() && isProjectedInside() );
360 }
361 
362 // ------------------------------------------------------------------------
363 template <typename TPredicate, DGtal::ProbingMode mode>
364 inline
365 bool
366 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::isInside() const
367 {
368  Triangle v = vertices();
369  for (int i = 0; i < 3; ++i) {
370  if (! myPredicate(v[i])) {
371  return false;
372  }
373  }
374  return true;
375 }
376 
377 // ------------------------------------------------------------------------
378 template <typename TPredicate, DGtal::ProbingMode mode>
379 inline
380 bool
381 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::isUnimodular() const
382 {
383  Point nest = getNormal();
384  for (int i = 0; i < 3; ++i) {
385  if (myM[i].dot(nest) != 1) {
386  return false;
387  }
388  }
389  return true;
390 }
391 
392 // ------------------------------------------------------------------------
393 template < typename TPredicate, DGtal::ProbingMode mode >
394 inline
395 bool
396 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::isProjectedInside (Triangle const& aTriangle) const
397 {
398  Triangle vec = { myQ - aTriangle[0], myQ - aTriangle[1], myQ - aTriangle[2] };
399 
400  Point s = myM[0] + myM[1] + myM[2];
401  std::array<bool, 3> res, res_not;
402  for (int i = 0; i < 3; ++i)
403  {
404  int im1 = (i - 1 + 3) % 3;
405  Point nk = vec[im1].crossProduct(vec[i]);
406 
407  res[i] = (nk.dot(-s) <= 0);
408  res_not[i] = !res[i];
409  }
410 
411  return (res[0] && res[1] && res[2]) || (res_not[0] && res_not[1] && res_not[2]);
412 }
413 
414 // ------------------------------------------------------------------------
415 template < typename TPredicate, DGtal::ProbingMode mode >
416 inline
417 bool
418 DGtal::PlaneProbingTetrahedronEstimator<TPredicate, mode>::isProjectedInside () const
419 {
420  return isProjectedInside(vertices());
421 }
422 
423 
424 ///////////////////////////////////////////////////////////////////////////////
425 // Implementation of inline functions //
426 
427 template <typename TPredicate, DGtal::ProbingMode mode>
428 inline
429 std::ostream&
430 DGtal::operator<< ( std::ostream & out,
431  const PlaneProbingTetrahedronEstimator<TPredicate, mode> & object )
432 {
433  object.selfDisplay( out );
434  return out;
435 }
436 
437 // //
438 ///////////////////////////////////////////////////////////////////////////////
439 
440