DGtal  1.0.0
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  RealPoint art = rt(t);
106  RealPoint arp = rp(t);
107  RealPoint artt = rtt(t);
108  RealPoint arpp = rpp(t);
109  RealPoint artp = rtp(t);
110 
111  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]);
112  double norme= sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
113  n[0]=n[0]/norme;
114  n[1]=n[1]/norme;
115  n[2]=n[2]/norme;
116  double E = art[0] * art[0]+ art[1] * art[1]+ art[2] * art[2];
117  double F= art[0] * arp[0]+ art[1] * arp[1]+ art[2] * arp[2];
118  double G = arp[0] * arp[0]+ arp[1] * arp[1]+ arp[2] * arp[2];
119  double L = artt[0] * n[0]+ artt[1] * n[1]+ artt[2] * n[2];
120  double M = artp[0] * n[0]+ artp[1] * n[1]+ artp[2] * n[2];
121  double N = arpp[0] * n[0]+ arpp[1] * n[1]+ arpp[2] * n[2];
122  double H = (E*N-2.0f*F*M+G*L)/(2.0f*(E*G-F*F));
123 
124 
125  return H;
126 }
127 
128 
129 
130 
131 /**
132  * @param t is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi] .
133  *
134  * @return the gaussian curvature at point (x(t),y(t),z(t)), positive
135  * is convex, negative is concave when shape is to the left and
136  * the shape boundary is followed counterclockwise.
137  */
138 template<typename TSpace>
139 inline
140 double
141 DGtal::StarShaped3D<TSpace>::gaussianCurvature( AngularCoordinates t ) const
142 {
143  RealPoint art = rt(t);
144  RealPoint arp = rp(t);
145  RealPoint artt = rtt(t);
146  RealPoint arpp = rpp(t);
147  RealPoint artp = rtp(t);
148 
149 
150  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]);
151  double norme= sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
152  n[0]=n[0]/norme;
153  n[1]=n[1]/norme;
154  n[2]=n[2]/norme;
155  double E = art[0] * art[0]+ art[1] * art[1]+ art[2] * art[2];
156  double F= art[0] * arp[0]+ art[1] * arp[1]+ art[2] * arp[2];
157  double G = arp[0] * arp[0]+ arp[1] * arp[1]+ arp[2] * arp[2];
158  double L = artt[0] * n[0]+ artt[1] * n[1]+ artt[2] * n[2];
159  double M = artp[0] * n[0]+ artp[1] * n[1]+ artp[2] * n[2];
160  double N = arpp[0] * n[0]+ arpp[1] * n[1]+ arpp[2] * n[2];
161 
162  double K = (L*N-M*M)/(E*G-F*F);
163 
164  return K;
165 }
166 
167 
168 
169 /**
170  * @param t1 is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi] .
171  * @param t2 is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi] further from [t1].
172  * @param nb the number of points used to estimate the arclength between x(Teta1,Phi1) and x(Teta2,Phi2).
173  * @return the estimated arclength.
174  */
175 template<typename TSpace>
176 inline
177 double
178 DGtal::StarShaped3D<TSpace>::arclength( AngularCoordinates t1, AngularCoordinates t2, unsigned int nb ) const
179 {
180  while ( t2.first < t1.first ) t2.first += 2.0*M_PI;
181  while ( t2.second < t1.second ) t2.second += 2.0*M_PI;
182 
183  RealPoint x0( x( t1 ) );
184  double l = 0.0;
185  for ( unsigned int i = 1; i <= nb; ++i )
186  {
187  AngularCoordinates t;
188  t.first = ( ( t2.first - t1.first ) * i ) / nb;
189  t.second = ( ( t2.second - t1.second ) * i ) / nb;
190  AngularCoordinates h;
191  h.first=t1.first + t.first;
192  h.second=t1.second + t.second;
193  RealPoint x1( x( h ) );
194  l += sqrt( ( x1[0] - x0[0] )*( x1[0] - x0[0] )
195  + ( x1[1] - x0[1] )*( x1[1] - x0[1] ) + ( x1[2] - x0[2] )*( x1[2] - x0[2] ) );
196  x0 = x1;
197  }
198  return l;
199 }
200 
201 
202 /**
203  * @param t1 is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi].
204  * @param t2 is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi], further from [t1].
205  * @param nb the number of points used to estimate the surfacelength between x(Teta1,Phi1) and x(Teta2,Phi2).
206  * @return the estimated surfacelength.
207  */
208 template<typename TSpace>
209 inline
210 double
211 DGtal::StarShaped3D<TSpace>::surfacelength( AngularCoordinates t1, AngularCoordinates t2, unsigned int nb ) const
212 {
213  while ( t2.first < t1.first ) t2.first += 2.0*M_PI;
214  while ( t2.second < t1.second ) t2.second += 2.0*M_PI;
215 
216  double step1 = ( t2.first - t1.first ) / nb;
217  double step2 = ( t2.second - t1.second ) / nb;
218  AngularCoordinates t( std::make_pair(0.0,0.0));
219  double l = 0.0;
220  for ( unsigned int i = 0; i < nb; ++i )
221  {
222  t.first = step1 *i;
223  for ( unsigned int j = 0; j < nb; ++j )
224  {
225 
226  t.second = step2 * j;
227  RealPoint xtp (x( std::make_pair( t1.first + t.first - step1 ,t1.second + t.second - step2 ) ));
228  RealPoint xt1p1( x( std::make_pair(t1.first + t.first,t1.second + t.second) ) );
229  RealPoint xt1p( x( std::make_pair(t1.first + t.first,t1.second + t.second- step2 ) ) );
230  double D = sqrt( ( xt1p[0] - xtp[0] )*( xt1p[0] - xtp[0] ) +
231  ( xt1p[1] - xtp[1] )*( xt1p[1] - xtp[1] ) +
232  ( xt1p[2] - xtp[2] )*( xt1p[2] - xtp[2] ));
233  double d = sqrt( ( xt1p1[0] - xt1p[0] )*( xt1p1[0] - xt1p[0] ) +
234  ( xt1p1[1] - xt1p[1] )*( xt1p1[1] - xt1p[1] ) +
235  ( xt1p1[2] - xt1p[2] )*( xt1p1[2] - xt1p[2] ));
236 
237 
238  l += (d*D);
239  }
240  }
241  return l;
242 }
243 
244 
245 
246 /**
247  * Destructor.
248  */
249 template <typename T>
250 inline
251 DGtal::StarShaped3D<T>::~StarShaped3D()
252 {
253 }
254 
255 ///////////////////////////////////////////////////////////////////////////////
256 // Interface - public :
257 
258 /**
259  * Writes/Displays the object on an output stream.
260  * @param out the output stream where the object is written.
261  */
262 template <typename T>
263 inline
264 void
265 DGtal::StarShaped3D<T>::selfDisplay ( std::ostream & out ) const
266 {
267  out << "[StarShaped2D]";
268 }
269 
270 /**
271  * Checks the validity/consistency of the object.
272  * @return 'true' if the object is valid, 'false' otherwise.
273  */
274 template <typename T>
275 inline
276 bool
277 DGtal::StarShaped3D<T>::isValid() const
278 {
279  return true;
280 }
281 
282 
283 
284 ///////////////////////////////////////////////////////////////////////////////
285 // Implementation of inline functions //
286 
287 template <typename T>
288 inline
289 std::ostream&
290 DGtal::operator<< ( std::ostream & out,
291  const StarShaped3D<T> & object )
292 {
293  object.selfDisplay( out );
294  return out;
295 }
296 
297 // //
298 ///////////////////////////////////////////////////////////////////////////////