DGtal  0.9.2
StarShaped3D.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 Anis Benyoub (\c anis.benyoub@liris.cnrs.fr )
20  * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
21  *
22  * @date 2012/06/05
23  *
24  * Header file for module StarShaped3D.cpp
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 #include <iostream>
33 //////////////////////////////////////////////////////////////////////////////
34 
35 ///////////////////////////////////////////////////////////////////////////////
36 // IMPLEMENTATION of inline methods.
37 ///////////////////////////////////////////////////////////////////////////////
38 
39 ///////////////////////////////////////////////////////////////////////////////
40 // ----------------------- Standard services ------------------------------
41 
42 typedef std::pair<double,double> AngularCoordinates;
43 /////////////////////////////////////////////////////////////////////////////
44 // ------------------------- star-shaped services -------------------------
45 
46 
47 /**
48  * @param p any point in the plane.
49  *
50  * @return 'true' if the point is inside the shape, 'false' if it
51  * is strictly outside.
52  */
53 template<typename TSpace>
54 inline
55 DGtal::Orientation
56 DGtal::StarShaped3D<TSpace>::orientation( const RealPoint & p ) const
57 {
58  RealPoint x_rel = x( parameter( p ));
59  x_rel -= center();
60  RealPoint p_rel( p );
61  p_rel -= center();
62 
63  double d_x = x_rel.norm();
64  double d_p = p_rel.norm();
65  double diff = d_p - d_x;
66 
67  if ( diff > 0.0 )
68  return OUTSIDE;
69  else if ( diff < 0.0 )
70  return INSIDE;
71  else
72  return ON;
73 }
74 
75 
76 /**
77  * @param t is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi] a.
78  *
79  * @return the vector (gradient(f(x,y,z))) made unitary which is the unit
80  * normal to the shape boundary looking inside the shape.
81  */
82 template<typename TSpace>
83 inline
84 typename DGtal::StarShaped3D<TSpace>::RealPoint
85 DGtal::StarShaped3D<TSpace>::normal( AngularCoordinates t ) const
86 {
87  RealPoint aNormal( gradient(t));
88  double norm = sqrt(aNormal[0]*aNormal[0] + aNormal[1]*aNormal[1] +aNormal[2]*aNormal[2]);
89  return (RealPoint(aNormal[0]/norm,aNormal[1]/norm,aNormal[2]/norm));
90 }
91 
92 
93 /**
94  * @param t is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi]
95  *
96  * @return the mean curvature at point (x(t),y(t),z(t)), positive
97  * is convex, negative is concave when shape is to the left and
98  * the shape boundary is followed counterclockwise.
99  */
100 template<typename TSpace>
101 inline
102 double
103 DGtal::StarShaped3D<TSpace>::meanCurvature( AngularCoordinates t) const
104 {
105  typedef typename Space::RealPoint RealPoint;
106 
107 
108  RealPoint art = rt(t);
109  RealPoint arp = rp(t);
110  RealPoint artt = rtt(t);
111  RealPoint arpp = rpp(t);
112  RealPoint artp = rtp(t);
113 
114  RealPoint n(art[1]*arp[2]-art[2]*arp[1],art[2]*arp[0]-art[0]*arp[2],art[0]*arp[1]-art[1]*arp[0]);
115  double norme= sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
116  n[0]=n[0]/norme;
117  n[1]=n[1]/norme;
118  n[2]=n[2]/norme;
119  double E = art[0] * art[0]+ art[1] * art[1]+ art[2] * art[2];
120  double F= art[0] * arp[0]+ art[1] * arp[1]+ art[2] * arp[2];
121  double G = arp[0] * arp[0]+ arp[1] * arp[1]+ arp[2] * arp[2];
122  double L = artt[0] * n[0]+ artt[1] * n[1]+ artt[2] * n[2];
123  double M = artp[0] * n[0]+ artp[1] * n[1]+ artp[2] * n[2];
124  double N = arpp[0] * n[0]+ arpp[1] * n[1]+ arpp[2] * n[2];
125  double H = (E*N-2.0f*F*M+G*L)/(2.0f*(E*G-F*F));
126 
127 
128  return H;
129 }
130 
131 
132 
133 
134 /**
135  * @param t is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi] .
136  *
137  * @return the gaussian curvature at point (x(t),y(t),z(t)), positive
138  * is convex, negative is concave when shape is to the left and
139  * the shape boundary is followed counterclockwise.
140  */
141 template<typename TSpace>
142 inline
143 double
144 DGtal::StarShaped3D<TSpace>::gaussianCurvature( AngularCoordinates t ) const
145 {
146 
147  typedef typename Space::RealPoint RealPoint;
148 
149  RealPoint art = rt(t);
150  RealPoint arp = rp(t);
151  RealPoint artt = rtt(t);
152  RealPoint arpp = rpp(t);
153  RealPoint artp = rtp(t);
154 
155 
156  RealPoint n(art[1]*arp[2]-art[2]*arp[1],art[2]*arp[0]-art[0]*arp[2],art[0]*arp[1]-art[1]*arp[0]);
157  double norme= sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
158  n[0]=n[0]/norme;
159  n[1]=n[1]/norme;
160  n[2]=n[2]/norme;
161  double E = art[0] * art[0]+ art[1] * art[1]+ art[2] * art[2];
162  double F= art[0] * arp[0]+ art[1] * arp[1]+ art[2] * arp[2];
163  double G = arp[0] * arp[0]+ arp[1] * arp[1]+ arp[2] * arp[2];
164  double L = artt[0] * n[0]+ artt[1] * n[1]+ artt[2] * n[2];
165  double M = artp[0] * n[0]+ artp[1] * n[1]+ artp[2] * n[2];
166  double N = arpp[0] * n[0]+ arpp[1] * n[1]+ arpp[2] * n[2];
167 
168  double K = (L*N-M*M)/(E*G-F*F);
169 
170  return K;
171 }
172 
173 
174 
175 /**
176  * @param t1 is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi] .
177  * @param t2 is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi] further from [t1].
178  * @param nb the number of points used to estimate the arclength between x(Teta1,Phi1) and x(Teta2,Phi2).
179  * @return the estimated arclength.
180  */
181 template<typename TSpace>
182 inline
183 double
184 DGtal::StarShaped3D<TSpace>::arclength( AngularCoordinates t1, AngularCoordinates t2, unsigned int nb ) const
185 {
186  while ( t2.first < t1.first ) t2.first += 2.0*M_PI;
187  while ( t2.second < t1.second ) t2.second += 2.0*M_PI;
188 
189  RealPoint x0( x( t1 ) );
190  double l = 0.0;
191  for ( unsigned int i = 1; i <= nb; ++i )
192  {
193  AngularCoordinates t;
194  t.first = ( ( t2.first - t1.first ) * i ) / nb;
195  t.second = ( ( t2.second - t1.second ) * i ) / nb;
196  AngularCoordinates h;
197  h.first=t1.first + t.first;
198  h.second=t1.second + t.second;
199  RealPoint x1( x( h ) );
200  l += sqrt( ( x1[0] - x0[0] )*( x1[0] - x0[0] )
201  + ( x1[1] - x0[1] )*( x1[1] - x0[1] ) + ( x1[2] - x0[2] )*( x1[2] - x0[2] ) );
202  x0 = x1;
203  }
204  return l;
205 }
206 
207 
208 /**
209  * @param t1 is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi].
210  * @param t2 is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi], further from [t1].
211  * @param nb the number of points used to estimate the surfacelength between x(Teta1,Phi1) and x(Teta2,Phi2).
212  * @return the estimated surfacelength.
213  */
214 template<typename TSpace>
215 inline
216 double
217 DGtal::StarShaped3D<TSpace>::surfacelength( AngularCoordinates t1, AngularCoordinates t2, unsigned int nb ) const
218 {
219  while ( t2.first < t1.first ) t2.first += 2.0*M_PI;
220  while ( t2.second < t1.second ) t2.second += 2.0*M_PI;
221 
222  double step1 = ( t2.first - t1.first ) / nb;
223  double step2 = ( t2.second - t1.second ) / nb;
224  AngularCoordinates t( std::make_pair(0.0,0.0));
225  double l = 0.0;
226  for ( unsigned int i = 0; i < nb; ++i )
227  {
228  t.first = step1 *i;
229  for ( unsigned int j = 0; j < nb; ++j )
230  {
231 
232  t.second = step2 * j;
233  RealPoint xtp (x( std::make_pair( t1.first + t.first - step1 ,t1.second + t.second - step2 ) ));
234  RealPoint xt1p1( x( std::make_pair(t1.first + t.first,t1.second + t.second) ) );
235  RealPoint xt1p( x( std::make_pair(t1.first + t.first,t1.second + t.second- step2 ) ) );
236  double D = sqrt( ( xt1p[0] - xtp[0] )*( xt1p[0] - xtp[0] ) +
237  ( xt1p[1] - xtp[1] )*( xt1p[1] - xtp[1] ) +
238  ( xt1p[2] - xtp[2] )*( xt1p[2] - xtp[2] ));
239  double d = sqrt( ( xt1p1[0] - xt1p[0] )*( xt1p1[0] - xt1p[0] ) +
240  ( xt1p1[1] - xt1p[1] )*( xt1p1[1] - xt1p[1] ) +
241  ( xt1p1[2] - xt1p[2] )*( xt1p1[2] - xt1p[2] ));
242 
243 
244  l += (d*D);
245  }
246  }
247  return l;
248 }
249 
250 
251 
252 /**
253  * Destructor.
254  */
255 template <typename T>
256 inline
257 DGtal::StarShaped3D<T>::~StarShaped3D()
258 {
259 }
260 
261 ///////////////////////////////////////////////////////////////////////////////
262 // Interface - public :
263 
264 /**
265  * Writes/Displays the object on an output stream.
266  * @param out the output stream where the object is written.
267  */
268 template <typename T>
269 inline
270 void
271 DGtal::StarShaped3D<T>::selfDisplay ( std::ostream & out ) const
272 {
273  out << "[StarShaped2D]";
274 }
275 
276 /**
277  * Checks the validity/consistency of the object.
278  * @return 'true' if the object is valid, 'false' otherwise.
279  */
280 template <typename T>
281 inline
282 bool
283 DGtal::StarShaped3D<T>::isValid() const
284 {
285  return true;
286 }
287 
288 
289 
290 ///////////////////////////////////////////////////////////////////////////////
291 // Implementation of inline functions //
292 
293 template <typename T>
294 inline
295 std::ostream&
296 DGtal::operator<< ( std::ostream & out,
297  const StarShaped3D<T> & object )
298 {
299  object.selfDisplay( out );
300  return out;
301 }
302 
303 // //
304 ///////////////////////////////////////////////////////////////////////////////