DGtal  1.5.beta
MaximalSegmentSliceEstimation.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/15
23  *
24  * Implementation of inline methods defined in MaximalSegmentSliceEstimation.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 TSurface >
43 inline
44 DGtal::MaximalSegmentSliceEstimation<TSurface>::
45 MaximalSegmentSliceEstimation ()
46 {}
47 
48 // ------------------------------------------------------------------------
49 template < typename TSurface >
50 inline
51 DGtal::MaximalSegmentSliceEstimation<TSurface>::
52 MaximalSegmentSliceEstimation (ConstAlias<Surface> aSurface)
53  : mySurface(aSurface)
54 {
55 }
56 
57 // ------------------------------------------------------------------------
58 template < typename TSurface >
59 inline
60 DGtal::MaximalSegmentSliceEstimation<TSurface>::
61 MaximalSegmentSliceEstimation ( const MaximalSegmentSliceEstimation & other )
62  : mySurface(other.mySurface), myH(other.myH)
63 {
64 }
65 
66 // ------------------------------------------------------------------------
67 template < typename TSurface >
68 inline
69 DGtal::MaximalSegmentSliceEstimation<TSurface>&
70 DGtal::MaximalSegmentSliceEstimation<TSurface>::operator= ( const MaximalSegmentSliceEstimation & other )
71 {
72  if (this != &other)
73  {
74  mySurface = other.mySurface;
75  myH = other.myH;
76  }
77 
78  return *this;
79 }
80 
81 // ------------------------------------------------------------------------
82 template < typename TSurface >
83 inline
84 DGtal::MaximalSegmentSliceEstimation<TSurface>::
85 ~MaximalSegmentSliceEstimation ()
86 {}
87 
88 // ----------------- model of CSurfelLocalEstimator -----------------------
89 
90 // ------------------------------------------------------------------------
91 template < typename TSurface >
92 template < typename SurfelConstIterator >
93 inline
94 void DGtal::MaximalSegmentSliceEstimation<TSurface>::
95 init (Scalar const& h, SurfelConstIterator /* itb */, SurfelConstIterator /* ite */)
96 {
97  myH = h;
98 }
99 
100 // ------------------------------------------------------------------------
101 template < typename TSurface >
102 template < typename SurfelConstIterator >
103 inline
104 typename DGtal::MaximalSegmentSliceEstimation<TSurface>::Quantity
105 DGtal::MaximalSegmentSliceEstimation<TSurface>::
106 eval (SurfelConstIterator it) const
107 {
108  ASSERT(mySurface != nullptr);
109 
110  KSpace const& K = space();
111  Surfel s = *it;
112 
113  // q = orthogonal direction
114  // q1, q2 = directions that span the plane of the surfel
115  Dimension q = K.sOrthDir(s);
116  typename KSpace::DirIterator q1 = K.sDirs(s), q2 = q1;
117  ++q2;
118 
119  // Plane (q1, q)
120  Tracker *tracker1 = mySurface->container().newTracker(s);
121  SurfaceSlice slice1(tracker1, *q1);
122  delete tracker1;
123 
124  Point2 n_left1, n_right1;
125  std::tie(n_left1, n_right1) = getExtremalPoints(slice1, *q1, q);
126 
127  // Plane (q2, q)
128  Tracker *tracker2 = mySurface->container().newTracker(s);
129  SurfaceSlice slice2(tracker2, *q2);
130  delete tracker2;
131 
132  Point2 n_left2, n_right2;
133  std::tie(n_left2, n_right2) = getExtremalPoints(slice2, *q2, q);
134 
135  // Find the components whose normals are in different quadrants => flat direction
136  auto signComponent = [](double x) -> int {
137  return (x >= 0) ? 1 : -1;
138  };
139 
140  auto differentQuadrant = [&signComponent](Point2 const& u, Point2 const& v) -> bool {
141  int sxu = signComponent(u[0]), syu = signComponent(u[1]);
142  int sxv = signComponent(v[0]), syv = signComponent(v[1]);
143  return (sxu != sxv) || (syu != syv);
144  };
145 
146  std::vector<int> flatDirections;
147  if (differentQuadrant(n_right1, n_left1))
148  {
149  flatDirections.push_back(*q1);
150  }
151 
152  if (differentQuadrant(n_right2, n_left2))
153  {
154  flatDirections.push_back(*q2);
155  }
156 
157  // Project the normal in the correct plane
158  RealPoint2 n1 = 0.5 * (n_right1 + n_left1);
159  RealPoint2 n2 = 0.5 * (n_right2 + n_left2);
160 
161  RealPoint n1_3d = RealPoint::zero, n2_3d = RealPoint::zero;
162 
163  n1_3d[*q1] = n1[0];
164  n1_3d[q] = n1[1];
165 
166  n2_3d[*q2] = n2[0];
167  n2_3d[q] = n2[1];
168 
169  // Orient it according to the trivial normal of the surfel
170  RealVector n = n1_3d.crossProduct(n2_3d);
171  Point trivial = trivialNormal(s);
172  if (trivial.dot(n) < 0)
173  {
174  n = -n;
175  }
176 
177  // Set to zero the flat directions
178  for (const auto& d: flatDirections)
179  {
180  n[d] = NumberTraits<Integer>::ZERO;
181  }
182 
183  if (n == RealVector::zero)
184  {
185  return trivial;
186  }
187 
188  return n;
189 }
190 
191 // ------------------------------------------------------------------------
192 template < typename TSurface >
193 template < typename SurfelConstIterator, typename OutputIterator >
194 inline
195 OutputIterator
196 DGtal::MaximalSegmentSliceEstimation<TSurface>::
197 eval (SurfelConstIterator itb, SurfelConstIterator ite, OutputIterator out) const
198 {
199  for (auto it = itb; it != ite; ++it)
200  {
201  *out++ = eval(it);
202  }
203 
204  return out;
205 }
206 
207 // ------------------------------------------------------------------------
208 template < typename TSurface >
209 inline
210 typename DGtal::MaximalSegmentSliceEstimation<TSurface>::Scalar
211 DGtal::MaximalSegmentSliceEstimation<TSurface>::
212 h () const
213 {
214  return myH;
215 }
216 
217 // --------------- model of CDigitalSurfaceLocalEstimator ------------------
218 
219 // ------------------------------------------------------------------------
220 template < typename TSurface >
221 inline
222 void DGtal::MaximalSegmentSliceEstimation<TSurface>::
223 attach (ConstAlias<Surface> aSurface)
224 {
225  mySurface = aSurface;
226 }
227 
228 ///////////////////////////////////////////////////////////////////////////////
229 // Interface - public :
230 
231 /**
232  * Writes/Displays the object on an output stream.
233  * @param out the output stream where the object is written.
234  */
235 template <typename TSurface>
236 inline
237 void
238 DGtal::MaximalSegmentSliceEstimation<TSurface>::selfDisplay ( std::ostream & out ) const
239 {
240  out << "[MaximalSegmentSliceEstimation]";
241 }
242 
243 /**
244  * Checks the validity/consistency of the object.
245  * @return 'true' if the object is valid, 'false' otherwise.
246  */
247 template <typename TSurface>
248 inline
249 bool
250 DGtal::MaximalSegmentSliceEstimation<TSurface>::isValid() const
251 {
252  return true;
253 }
254 
255 ///////////////////////////////////////////////////////////////////////////////
256 // ------------------------- Internals ------------------------------------
257 
258 // ------------------------------------------------------------------------
259 template < typename TSurface >
260 inline
261 typename DGtal::MaximalSegmentSliceEstimation<TSurface>::KSpace const&
262 DGtal::MaximalSegmentSliceEstimation<TSurface>::
263 space () const
264 {
265  return mySurface->container().space();
266 }
267 
268 // ------------------------------------------------------------------------
269 template < typename TSurface >
270 inline
271 std::pair<typename DGtal::MaximalSegmentSliceEstimation<TSurface>::Point2,
272  typename DGtal::MaximalSegmentSliceEstimation<TSurface>::Point2>
273 DGtal::MaximalSegmentSliceEstimation<TSurface>::
274 getExtremalPoints (SurfaceSlice const& aSlice, Dimension const& aProjectDir1, Dimension const& aProjectDir2) const
275 {
276  KSpace const& K = space();
277 
278  DSSComputer computer_left(K, aProjectDir1, aProjectDir2);
279  computer_left.init(aSlice.c());
280  while (computer_left.extendBack()) {}
281  while (computer_left.extendFront()) {}
282  Point2 n_left(computer_left.primitive().dsl().b(), computer_left.primitive().dsl().a());
283 
284  DSSComputer computer_right(K, aProjectDir1, aProjectDir2);
285  computer_right.init(aSlice.c());
286  while (computer_right.extendFront()) {}
287  while (computer_right.extendBack()) {}
288  Point2 n_right(computer_right.primitive().dsl().b(), computer_right.primitive().dsl().a());
289 
290  return { n_left, n_right };
291 }
292 
293 // ------------------------------------------------------------------------
294 template < typename TSurface >
295 inline
296 typename DGtal::MaximalSegmentSliceEstimation<TSurface>::Vector
297 DGtal::MaximalSegmentSliceEstimation<TSurface>::
298 trivialNormal (Surfel const& aSurfel) const
299 {
300  const KSpace& K = space();
301 
302  DGtal::Dimension q = K.sOrthDir(aSurfel);
303  bool direct = K.sDirect(aSurfel, q);
304 
305  Vector trivial_normal = Vector::zero;
306  trivial_normal[q] = direct ? -1 : 1;
307 
308  return trivial_normal;
309 }
310 
311 ///////////////////////////////////////////////////////////////////////////////
312 // Implementation of inline functions //
313 
314 template <typename TSurface>
315 inline
316 std::ostream&
317 DGtal::operator<< ( std::ostream & out,
318  const MaximalSegmentSliceEstimation<TSurface> & object )
319 {
320  object.selfDisplay( out );
321  return out;
322 }
323 
324 // //
325 ///////////////////////////////////////////////////////////////////////////////
326 
327