DGtal  1.5.beta
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  const RealPoint x_rel( x(parameter(p)) - center() );
59  const RealPoint p_rel( p - center() );
60  const double diff = p_rel.norm() - x_rel.norm();
61 
62  return isAlmostEqual(diff,0.) ? ON : (diff > 0. ? OUTSIDE : INSIDE);
63 }
64 
65 
66 /**
67  * @param t is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi] a.
68  *
69  * @return the vector (gradient(f(x,y,z))) made unitary which is the unit
70  * normal to the shape boundary looking inside the shape.
71  */
72 template<typename TSpace>
73 inline
74 typename DGtal::StarShaped3D<TSpace>::RealPoint
75 DGtal::StarShaped3D<TSpace>::normal( const AngularCoordinates& t ) const
76 {
77  const RealPoint aNormal( gradient(t) );
78  return aNormal / aNormal.norm();
79 }
80 
81 
82 /**
83  * @param t is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi]
84  *
85  * @return the mean curvature at point (x(t),y(t),z(t)), positive
86  * is convex, negative is concave when shape is to the left and
87  * the shape boundary is followed counterclockwise.
88  */
89 template<typename TSpace>
90 inline
91 double
92 DGtal::StarShaped3D<TSpace>::meanCurvature( const AngularCoordinates& t ) const
93 {
94  const RealPoint art = rt(t);
95  const RealPoint arp = rp(t);
96  const RealPoint artt = rtt(t);
97  const RealPoint arpp = rpp(t);
98  const RealPoint artp = rtp(t);
99 
100  RealPoint n(
101  art[1] * arp[2] - art[2] * arp[1],
102  art[2] * arp[0] - art[0] * arp[2],
103  art[0] * arp[1] - art[1] * arp[0]
104  );
105  n /= n.norm();
106 
107  const double E = art[0] * art[0]+ art[1] * art[1]+ art[2] * art[2];
108  const double F = art[0] * arp[0]+ art[1] * arp[1]+ art[2] * arp[2];
109  const double G = arp[0] * arp[0]+ arp[1] * arp[1]+ arp[2] * arp[2];
110  const double L = artt[0] * n[0]+ artt[1] * n[1]+ artt[2] * n[2];
111  const double M = artp[0] * n[0]+ artp[1] * n[1]+ artp[2] * n[2];
112  const double N = arpp[0] * n[0]+ arpp[1] * n[1]+ arpp[2] * n[2];
113 
114  return ( E*N - 2.*F*M + G*L ) / ( 2.* (E*G - F*F) );
115 }
116 
117 
118 
119 
120 /**
121  * @param t is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi] .
122  *
123  * @return the gaussian curvature at point (x(t),y(t),z(t)), positive
124  * is convex, negative is concave when shape is to the left and
125  * the shape boundary is followed counterclockwise.
126  */
127 template<typename TSpace>
128 inline
129 double
130 DGtal::StarShaped3D<TSpace>::gaussianCurvature( const AngularCoordinates& t ) const
131 {
132  const RealPoint art = rt(t);
133  const RealPoint arp = rp(t);
134  const RealPoint artt = rtt(t);
135  const RealPoint arpp = rpp(t);
136  const RealPoint artp = rtp(t);
137 
138 
139  RealPoint n(
140  art[1] * arp[2] - art[2] * arp[1],
141  art[2] * arp[0] - art[0] * arp[2],
142  art[0] * arp[1] - art[1] * arp[0]
143  );
144  n /= n.norm();
145 
146  const double E = art[0] * art[0]+ art[1] * art[1]+ art[2] * art[2];
147  const double F= art[0] * arp[0]+ art[1] * arp[1]+ art[2] * arp[2];
148  const double G = arp[0] * arp[0]+ arp[1] * arp[1]+ arp[2] * arp[2];
149  const double L = artt[0] * n[0]+ artt[1] * n[1]+ artt[2] * n[2];
150  const double M = artp[0] * n[0]+ artp[1] * n[1]+ artp[2] * n[2];
151  const double N = arpp[0] * n[0]+ arpp[1] * n[1]+ arpp[2] * n[2];
152 
153  return ( L*N - M*M ) / ( E*G - F*F );
154 }
155 
156 
157 
158 /**
159  * @param t1 is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi] .
160  * @param t2 is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi] further from [t1].
161  * @param nb the number of points used to estimate the arclength between x(Teta1,Phi1) and x(Teta2,Phi2).
162  * @return the estimated arclength.
163  */
164 template<typename TSpace>
165 inline
166 double
167 DGtal::StarShaped3D<TSpace>::arclength( const AngularCoordinates& t1, AngularCoordinates t2, unsigned int nb ) const
168 {
169  while ( t2.first < t1.first ) t2.first += 2.0*M_PI;
170  while ( t2.second < t1.second ) t2.second += 2.0*M_PI;
171 
172  RealPoint x0( x(t1) );
173  double l = 0.;
174 
175  for ( unsigned int i = 1; i <= nb; ++i )
176  {
177  const AngularCoordinates h(
178  t1.first + ((t2.first - t1.first) * i) / nb,
179  t1.second + ((t2.second - t1.second) * i) / nb
180  );
181  const RealPoint x1( x(h) );
182  l += sqrt( (x1[0] - x0[0]) * (x1[0] - x0[0])
183  + (x1[1] - x0[1]) * (x1[1] - x0[1])
184  + (x1[2] - x0[2]) * (x1[2] - x0[2]) );
185  x0 = x1;
186  }
187  return l;
188 }
189 
190 
191 /**
192  * @param t1 is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi].
193  * @param t2 is a couple of Teta && Phi wich are 2 angles respectivly between [0,2PI] and [0,Pi], further from [t1].
194  * @param nb the number of points used to estimate the surfacelength between x(Teta1,Phi1) and x(Teta2,Phi2).
195  * @return the estimated surfacelength.
196  */
197 template<typename TSpace>
198 inline
199 double
200 DGtal::StarShaped3D<TSpace>::surfacelength( const AngularCoordinates& t1, AngularCoordinates t2, unsigned int nb ) const
201 {
202  while ( t2.first < t1.first ) t2.first += 2.0*M_PI;
203  while ( t2.second < t1.second ) t2.second += 2.0*M_PI;
204 
205  const double step1 = ( t2.first - t1.first ) / nb;
206  const double step2 = ( t2.second - t1.second ) / nb;
207 
208  double l = 0.0;
209 
210  for ( unsigned int i = 0; i < nb; ++i )
211  {
212  const double tFirst = step1 * i;
213  for ( unsigned int j = 0; j < nb; ++j )
214  {
215  const double tSecond = step2 * j;
216  const RealPoint xtp ( x( std::make_pair( t1.first + tFirst - step1 ,t1.second + tSecond - step2 )));
217  const RealPoint xt1p1( x( std::make_pair( t1.first + tFirst, t1.second + tSecond )));
218  const RealPoint xt1p( x( std::make_pair( t1.first + tFirst, t1.second + tSecond - step2 )));
219 
220  l += sqrt( (xt1p[0] - xtp[0]) * (xt1p[0] - xtp[0]) +
221  (xt1p[1] - xtp[1]) * (xt1p[1] - xtp[1]) +
222  (xt1p[2] - xtp[2]) * (xt1p[2] - xtp[2]) )
223  * sqrt( (xt1p1[0] - xt1p[0]) * (xt1p1[0] - xt1p[0]) +
224  (xt1p1[1] - xt1p[1]) * (xt1p1[1] - xt1p[1]) +
225  (xt1p1[2] - xt1p[2]) * (xt1p1[2] - xt1p[2]) );
226  }
227  }
228  return l;
229 }
230 
231 ///////////////////////////////////////////////////////////////////////////////
232 // Interface - public :
233 
234 /**
235  * Writes/Displays the object on an output stream.
236  * @param out the output stream where the object is written.
237  */
238 template <typename T>
239 inline
240 void
241 DGtal::StarShaped3D<T>::selfDisplay ( std::ostream & out ) const
242 {
243  out << "[StarShaped2D]";
244 }
245 
246 /**
247  * Checks the validity/consistency of the object.
248  * @return 'true' if the object is valid, 'false' otherwise.
249  */
250 template <typename T>
251 inline
252 bool
253 DGtal::StarShaped3D<T>::isValid() const
254 {
255  return true;
256 }
257 
258 
259 
260 ///////////////////////////////////////////////////////////////////////////////
261 // Implementation of inline functions //
262 
263 template <typename T>
264 inline
265 std::ostream&
266 DGtal::operator<< ( std::ostream & out,
267  const StarShaped3D<T> & object )
268 {
269  object.selfDisplay( out );
270  return out;
271 }
272 
273 // //
274 ///////////////////////////////////////////////////////////////////////////////