DGtal 1.4.0
Loading...
Searching...
No Matches
ImplicitPolynomial3Shape.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 ImplicitPolynomial3Shape.ih
19 * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20 * Laboratory of Mathematics (CNRS, UMR 5807), University of Savoie, France
21 *
22 * @date 2012/02/14
23 *
24 * Implementation of inline methods defined in ImplicitPolynomial3Shape.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//-----------------------------------------------------------------------------
42template <typename TSpace>
43inline
44DGtal::ImplicitPolynomial3Shape<TSpace>::~ImplicitPolynomial3Shape()
45{
46}
47//-----------------------------------------------------------------------------
48template <typename TSpace>
49inline
50DGtal::ImplicitPolynomial3Shape<TSpace>::
51ImplicitPolynomial3Shape( const Polynomial3 & poly )
52{
53 init( poly );
54}
55//-----------------------------------------------------------------------------
56template <typename TSpace>
57inline
58DGtal::ImplicitPolynomial3Shape<TSpace> &
59DGtal::ImplicitPolynomial3Shape<TSpace>::
60operator=( const ImplicitPolynomial3Shape & other )
61{
62 if ( this != &other )
63 {
64 myPolynomial = other.myPolynomial;
65
66 myFx= other.myFx;
67 myFy= other.myFy;
68 myFz= other.myFz;
69
70 myFxx= other.myFxx;
71 myFxy= other.myFxy;
72 myFxz= other.myFxz;
73
74 myFyx= other.myFyx;
75 myFyy= other.myFyy;
76 myFyz= other.myFyz;
77
78 myFzx= other.myFzx;
79 myFzy= other.myFzy;
80 myFzz= other.myFzz;
81
82 myUpPolynome = other.myUpPolynome;
83 myLowPolynome = other.myLowPolynome;
84 }
85 return *this;
86}
87//-----------------------------------------------------------------------------
88template <typename TSpace>
89inline
90void
91DGtal::ImplicitPolynomial3Shape<TSpace>::
92init( const Polynomial3 & poly )
93{
94 myPolynomial = poly;
95
96 myFx= derivative<0>( poly );
97 myFy= derivative<1>( poly );
98 myFz= derivative<2>( poly );
99
100 myFxx= derivative<0>( myFx );
101 myFxy= derivative<1>( myFx );
102 myFxz= derivative<2>( myFx);
103
104 myFyx= derivative<0>( myFy );
105 myFyy= derivative<1>( myFy );
106 myFyz= derivative<2>( myFy );
107
108 myFzx= derivative<0>( myFz );
109 myFzy= derivative<1>( myFz );
110 myFzz= derivative<2>( myFz );
111
112 // These two polynomials are used for mean curvature estimation.
113 myUpPolynome = myFx*(myFx*myFxx+myFy*myFyx+myFz*myFzx)+
114 myFy*(myFx*myFxy+myFy*myFyy+myFz*myFzy)+
115 myFz*(myFx*myFxz+myFy*myFyz+myFz*myFzz)-
116 ( myFx*myFx +myFy*myFy+myFz*myFz )*(myFxx+myFyy+myFzz);
117
118 myLowPolynome = myFx*myFx +myFy*myFy+myFz*myFz;
119}
120//-----------------------------------------------------------------------------
121template <typename TSpace>
122inline
123double
124DGtal::ImplicitPolynomial3Shape<TSpace>::
125operator()(const RealPoint &aPoint) const
126{
127 return myPolynomial( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
128}
129//-----------------------------------------------------------------------------
130template <typename TSpace>
131inline
132bool
133DGtal::ImplicitPolynomial3Shape<TSpace>::
134isInside(const RealPoint &aPoint) const
135{
136 return orientation( aPoint ) == INSIDE;
137}
138//-----------------------------------------------------------------------------
139template <typename TSpace>
140inline
141DGtal::Orientation
142DGtal::ImplicitPolynomial3Shape<TSpace>::
143orientation(const RealPoint &aPoint) const
144{
145 Ring v = this->operator()(aPoint);
146 if ( v < (Ring)0 )
147 return INSIDE;
148 else if ( v > (Ring)0 )
149 return OUTSIDE;
150 else
151 return ON;
152}
153//-----------------------------------------------------------------------------
154template <typename TSpace>
155inline
156typename DGtal::ImplicitPolynomial3Shape<TSpace>::RealVector
157DGtal::ImplicitPolynomial3Shape<TSpace>::
158gradient( const RealPoint &aPoint ) const
159{
160 // ISO C++ tells that an object created at return time will not be
161 // copied into the caller context, but will be already defined in
162 // the correct context.
163 return RealVector
164 ( myFx ( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] ),
165 myFy ( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] ),
166 myFz ( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] ) );
167
168}
169
170
171// ------------------------------------------------------------ Added by Anis Benyoub
172//-----------------------------------------------------------------------------
173
174/**
175 * @param aPoint any point in the Euclidean space.
176 * This computation is based on the hessian formula of the mean curvature
177 * k=-(∇F ∗ H (F ) ∗ ∇F T − |∇F |^2 *Trace(H (F ))/2|∇F |^3
178 * we define it as positive for a sphere
179 * @return the mean curvature value of the polynomial at \a aPoint.
180 *
181*/
182template <typename TSpace>
183inline
184double
185DGtal::ImplicitPolynomial3Shape<TSpace>::
186meanCurvature( const RealPoint &aPoint ) const
187{
188 double temp= myLowPolynome( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
189 temp = sqrt(temp);
190 double downValue = 2.0*(temp*temp*temp);
191 double upValue = myUpPolynome( aPoint[ 0 ] )( aPoint[ 1 ] )( aPoint[ 2 ] );
192
193
194 return -(upValue/downValue);
195}
196
197
198
199//-----------------------------------------------------------------------------
200template <typename TSpace>
201inline
202double
203DGtal::ImplicitPolynomial3Shape<TSpace>::
204gaussianCurvature( const RealPoint &aPoint ) const
205{
206 /*
207 JOL: new Gaussian curvature formula (in sage)
208 var('Fx','Fy','Fz','Fxx','Fxy','Fxz','Fyy','Fyz','Fzz')
209 M=Matrix(4,4,[[Fxx,Fxy,Fxz,Fx],[Fxy,Fyy,Fyz,Fy],[Fxz,Fyz,Fzz,Fz],[Fx,Fy,Fz,0]])
210 det(M)
211# Fxz^2*Fy^2 - 2*Fx*Fxz*Fy*Fyz + Fx^2*Fyz^2 - 2*Fxy*Fxz*Fy*Fz + 2*Fx*Fxz*Fyy*Fz - 2*Fx*Fxy*Fyz*Fz + 2*Fxx*Fy*Fyz*Fz + Fxy^2*Fz^2 - Fxx*Fyy*Fz^2 + 2*Fx*Fxy*Fy*Fzz - Fxx*Fy^2*Fzz - Fx^2*Fyy*Fzz
212 G = -det(M) / ( Fx^2 + Fy^2 + Fz^2 )^2
213 */
214 const double x = aPoint[ 0 ];
215 const double y = aPoint[ 1 ];
216 const double z = aPoint[ 2 ];
217 const double Fx = myFx( x )( y )( z );
218 const double Fy = myFy( x )( y )( z );
219 const double Fz = myFz( x )( y )( z );
220 const double Fx2 = Fx * Fx;
221 const double Fy2 = Fy * Fy;
222 const double Fz2 = Fz * Fz;
223 const double G2 = Fx2 + Fy2 + Fz2;
224 const double Fxx = myFxx( x )( y )( z );
225 const double Fxy = myFxy( x )( y )( z );
226 const double Fxz = myFxz( x )( y )( z );
227 const double Fyy = myFyy( x )( y )( z );
228 const double Fyz = myFyz( x )( y )( z );
229 const double Fzz = myFzz( x )( y )( z );
230 const double Ax2 = ( Fyz * Fyz - Fyy * Fzz ) * Fx2;
231 const double Ay2 = ( Fxz * Fxz - Fxx * Fzz ) * Fy2;
232 const double Az2 = ( Fxy * Fxy - Fxx * Fyy ) * Fz2;
233 const double Axy = ( Fxy * Fzz - Fxz * Fyz ) * Fx * Fy;
234 const double Axz = ( Fxz * Fyy - Fxy * Fyz ) * Fx * Fz;
235 const double Ayz = ( Fxx * Fyz - Fxy * Fxz ) * Fy * Fz;
236 const double det = Ax2 + Ay2 + Az2 + 2 * ( Axy + Axz + Ayz );
237 return - det / ( G2*G2 );
238}
239
240template< typename TSpace >
241inline
242void
243DGtal::ImplicitPolynomial3Shape<TSpace>::principalCurvatures
244( const RealPoint & aPoint,
245 double & k1,
246 double & k2 ) const
247{
248 double H = meanCurvature( aPoint );
249 double G = gaussianCurvature( aPoint );
250 double tmp = std::sqrt( fabs( H * H - G ));
251 k2 = H + tmp;
252 k1 = H - tmp;
253}
254
255template< typename TSpace >
256inline
257void
258DGtal::ImplicitPolynomial3Shape<TSpace>::principalDirections
259( const RealPoint & aPoint,
260 RealVector & d1,
261 RealVector & d2 ) const
262{
263 const RealVector grad_F = gradient( aPoint );
264 const auto Fn = grad_F.norm();
265 if ( Fn < 1e-8 )
266 {
267 d1 = d2 = RealVector();
268 return;
269 }
270 RealVector u, v;
271 const RealVector n = grad_F / Fn;
272 u = RealVector( 1.0, 0.0, 0.0 ).crossProduct( n );
273 auto u_norm = u.norm();
274 if ( u_norm < 1e-8 )
275 {
276 u = RealVector( 0.0, 1.0, 0.0 ).crossProduct( n );
277 u_norm = u.norm();
278 }
279 u /= u_norm;
280 v = n.crossProduct( u );
281 double k_min, k_max;
282 principalCurvatures( aPoint, k_min, k_max );
283 const double x = aPoint[ 0 ];
284 const double y = aPoint[ 1 ];
285 const double z = aPoint[ 2 ];
286 // Computing Hessian matrix
287 const double Fxx = myFxx( x )( y )( z );
288 const double Fxy = myFxy( x )( y )( z );
289 const double Fxz = myFxz( x )( y )( z );
290 const double Fyy = myFyy( x )( y )( z );
291 const double Fyz = myFyz( x )( y )( z );
292 const double Fzz = myFzz( x )( y )( z );
293 const RealVector HessF_u = { Fxx * u[ 0 ] + Fxy * u[ 1 ] + Fxz * u[ 2 ],
294 Fxy * u[ 0 ] + Fyy * u[ 1 ] + Fyz * u[ 2 ],
295 Fxz * u[ 0 ] + Fyz * u[ 1 ] + Fzz * u[ 2 ] };
296 const RealVector HessF_v = { Fxx * v[ 0 ] + Fxy * v[ 1 ] + Fxz * v[ 2 ],
297 Fxy * v[ 0 ] + Fyy * v[ 1 ] + Fyz * v[ 2 ],
298 Fxz * v[ 0 ] + Fyz * v[ 1 ] + Fzz * v[ 2 ] };
299 const double Fuu = u.dot( HessF_u );
300 const double Fuv = u.dot( HessF_v );
301 const double Fvv = v.dot( HessF_v );
302 if ( fabs( k_min * Fn - Fuu ) >= fabs( k_min * Fn - Fvv ) )
303 {
304 // Choose k1 = k_min and k2 = k_max,
305 // to avoid null k1*Fn - Fuu = -(k2*Fn - Fvv) = 0
306 double k1 = k_min;
307 double k2 = k_max;
308 d1 = RealVector( ( k1 * Fn - Fuu ) * v[ 0 ] + Fuv * u[ 0 ],
309 ( k1 * Fn - Fuu ) * v[ 1 ] + Fuv * u[ 1 ],
310 ( k1 * Fn - Fuu ) * v[ 2 ] + Fuv * u[ 2 ] );
311 d2 = -1.0 * RealVector( ( k2 * Fn - Fvv ) * u[ 0 ] + Fuv * v[ 0 ],
312 ( k2 * Fn - Fvv ) * u[ 1 ] + Fuv * v[ 1 ],
313 ( k2 * Fn - Fvv ) * u[ 2 ] + Fuv * v[ 2 ] );
314 }
315 else
316 {
317 // Choose k2 = k_min and k1 = k_max,
318 // then | k_max*Fn - Fuu | >= | k_max*Fn - Fvv | >= 0
319 double k1 = k_max;
320 double k2 = k_min;
321 d2 = RealVector( ( k1 * Fn - Fuu ) * v[ 0 ] + Fuv * u[ 0 ],
322 ( k1 * Fn - Fuu ) * v[ 1 ] + Fuv * u[ 1 ],
323 ( k1 * Fn - Fuu ) * v[ 2 ] + Fuv * u[ 2 ] );
324 d1 = -1.0 * RealVector( ( k2 * Fn - Fvv ) * u[ 0 ] + Fuv * v[ 0 ],
325 ( k2 * Fn - Fvv ) * u[ 1 ] + Fuv * v[ 1 ],
326 ( k2 * Fn - Fvv ) * u[ 2 ] + Fuv * v[ 2 ] );
327 }
328 d1 /= d1.norm();
329 d2 /= d2.norm();
330}
331
332/**
333 *@param aPoint any point in the Euclidean space.
334 *@param accuracy refers to the precision
335 *@param maxIter refers to the maximum iterations the fonction user authorises
336 *@param gamma refers to the step
337 *@return the nearest point on the surface to the one given in parameter.
338 */
339template <typename TSpace>
340inline
341typename DGtal::ImplicitPolynomial3Shape<TSpace>::RealPoint
342DGtal::ImplicitPolynomial3Shape<TSpace>::nearestPoint
343( const RealPoint &aPoint, const double accuracy,
344 const int maxIter, const double gamma ) const
345{
346 RealPoint X = aPoint;
347 for ( int numberIter = 0; numberIter < maxIter; numberIter++ )
348 {
349 double val_X = (*this)( X );
350 if ( fabs( val_X ) < accuracy ) break;
351 RealVector grad_X = (*this).gradient( X );
352 double n2_grad_X = grad_X.dot( grad_X );
353 if ( n2_grad_X > 0.000001 ) grad_X /= n2_grad_X;
354 X -= val_X * gamma * grad_X ;
355 }
356 return X;
357}
358
359///////////////////////////////////////////////////////////////////////////////
360// Interface - public :
361
362/**
363 * Writes/Displays the object on an output stream.
364 * @param out the output stream where the object is written.
365 */
366template <typename TSpace>
367inline
368void
369DGtal::ImplicitPolynomial3Shape<TSpace>::selfDisplay ( std::ostream & out ) const
370{
371 out << "[ImplicitPolynomial3Shape] P(x,y,z) = " << myPolynomial;
372}
373
374/**
375 * Checks the validity/consistency of the object.
376 * @return 'true' if the object is valid, 'false' otherwise.
377 */
378template <typename TSpace>
379inline
380bool
381DGtal::ImplicitPolynomial3Shape<TSpace>::isValid() const
382{
383 return true;
384}
385
386
387
388///////////////////////////////////////////////////////////////////////////////
389// Implementation of inline functions //
390
391template <typename TSpace>
392inline
393std::ostream&
394DGtal::operator<< ( std::ostream & out,
395 const ImplicitPolynomial3Shape<TSpace> & object )
396{
397 object.selfDisplay( out );
398 return out;
399}
400
401// //
402///////////////////////////////////////////////////////////////////////////////
403
404