DGtal  1.2.0
ShroudsRegularization.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 ShroudsRegularization.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 2020/06/09
23  *
24  * Implementation of inline methods defined in ShroudsRegularization.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 //////////////////////////////////////////////////////////////////////////////
33 
34 template < typename TDigitalSurfaceContainer >
35 double
36 DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
37 energyArea()
38 {
39  parameterize();
40  double E = 0.0;
41  for ( Vertex v = 0; v < myT.size(); ++v )
42  {
43  double area = 1.0;
44  const auto s = myPtrIdxSurface->surfel( v );
45  const auto k = myPtrK->sOrthDir( s );
46  for ( Dimension i = 0; i < 3; ++i )
47  {
48  if ( i == k ) continue; // not a valid slice
49  const Scalar dn = myNextD[ i ][ v ];
50  const Scalar dp = myPrevD[ i ][ v ];
51  const Scalar l = 0.5 * ( dn + dp ); // local length
52  area *= l;
53  }
54  E += area;
55  }
56  return E;
57 }
58 
59 template < typename TDigitalSurfaceContainer >
60 double
61 DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
62 energySnake()
63 {
64  parameterize();
65  double E = 0.0;
66  for ( Vertex v = 0; v < myT.size(); ++v )
67  {
68  const auto s = myPtrIdxSurface->surfel( v );
69  const auto k = myPtrK->sOrthDir( s );
70  for ( Dimension i = 0; i < 3; ++i )
71  {
72  if ( i == k ) continue; // not a valid slice
73  const Scalar dn = myNextD[ i ][ v ];
74  const Scalar dp = myPrevD[ i ][ v ];
75  const Scalar l = 0.5 * ( dn + dp ); // local length
76  const Scalar cn = 2.0 / ( dn*dn+dn*dp );
77  const Scalar cp = 2.0 / ( dp*dn+dp*dp );
78  const Scalar ci = 2.0 / ( dp*dn );
79  const RealPoint vn = position( myNext[ i ][ v ] );
80  const RealPoint vp = position( myPrev[ i ][ v ] );
81  const RealPoint vi = position( v );
82  const Scalar xp = ( vn[ i ] - vp[ i ] ) / ( dn + dp );
83  const Scalar yp = ( vn[ k ] - vp[ k ] ) / ( dn + dp );
84  const Scalar xpp = cn * vn[ i ] - ci * vi[ i ] + cp * vp[ i ];
85  const Scalar ypp = cn * vn[ k ] - ci * vi[ k ] + cp * vp[ k ];
86  E += l * ( myAlpha * ( xp * xp + yp * yp )
87  + myBeta * ( xpp * xpp + ypp * ypp ) );
88  }
89  }
90  return E;
91 }
92 
93 template < typename TDigitalSurfaceContainer >
94 double
95 DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
96 energySquaredCurvature()
97 {
98  parameterize();
99  double E = 0.0;
100  for ( Vertex v = 0; v < myT.size(); ++v )
101  {
102  double area = 1.0;
103  const auto s = myPtrIdxSurface->surfel( v );
104  const auto k = myPtrK->sOrthDir( s );
105  for ( Dimension i = 0; i < 3; ++i )
106  {
107  if ( i == k ) continue; // not a valid slice
108  const Scalar dn = myNextD[ i ][ v ];
109  const Scalar dp = myPrevD[ i ][ v ];
110  const Scalar l = 0.5 * ( dn + dp ); // local length
111  const Scalar cn = 2.0 / ( dn*dn+dn*dp );
112  const Scalar cp = 2.0 / ( dp*dn+dp*dp );
113  const Scalar ci = 2.0 / ( dp*dn );
114  const RealPoint vn = position( myNext[ i ][ v ] );
115  const RealPoint vp = position( myPrev[ i ][ v ] );
116  const RealPoint vi = position( v );
117  const Scalar xp = ( vn[ i ] - vp[ i ] ) / ( dn + dp );
118  const Scalar yp = ( vn[ k ] - vp[ k ] ) / ( dn + dp );
119  const Scalar xpp = cn * vn[ i ] - ci * vi[ i ] + cp * vp[ i ];
120  const Scalar ypp = cn * vn[ k ] - ci * vi[ k ] + cp * vp[ k ];
121  E += l * ( pow( xp * ypp - yp * xpp, 2.0 )
122  / pow( xp * xp + yp * yp, 3.0 ) );
123  }
124  }
125  return E;
126 }
127 
128 template < typename TDigitalSurfaceContainer >
129 std::pair<double,double>
130 DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
131 oneStepAreaMinimization( const double randomization )
132 {
133  parameterize();
134  Scalars newT = myT;
135  for ( Vertex v = 0; v < myT.size(); ++v )
136  {
137  double right = 0.0;
138  double left = 0.0;
139  double coef = 0.0;
140  const auto s = myPtrIdxSurface->surfel( v );
141  const auto k = myPtrK->sOrthDir( s );
142  for ( Dimension i = 0; i < 3; ++i )
143  {
144  if ( i == k ) continue; // not a valid slice
145  const Scalar dn = myNextD[ i ][ v ];
146  const Scalar dp = myPrevD[ i ][ v ];
147  const Scalar cn = 2.0 / ( dn*dn+dn*dp );
148  const Scalar cp = 2.0 / ( dp*dn+dp*dp );
149  const Scalar ci = 2.0 / ( dp*dn );
150  const RealPoint vn = position( myNext[ i ][ v ] );
151  const RealPoint vp = position( myPrev[ i ][ v ] );
152  const RealPoint vi = position( v );
153  const Scalar yp = ( vn[ k ] - vp[ k ] ) / ( dn + dp );
154  const Scalar xp = ( vn[ i ] - vp[ i ] ) / ( dn + dp );
155  const Scalar xpp = cn * vn[ i ] - ci * vi[ i ] + cp * vp[ i ];
156  right += yp * xpp / xp;
157  left += cn * vn[ k ] + cp * vp[ k ] - ci * myInsV[ v ][ k ];
158  coef += ci * ( myInsV[ v ][ k ] - myOutV[ v ][ k ] );
159  }
160  newT[ v ] = ( right - left ) / coef
161  + ( (double) rand() / (double) RAND_MAX - 0.49 ) * randomization;
162  }
163  // Weak damping since problem is convex.
164  auto X = positions();
165  for ( Vertex v = 0; v < myT.size(); ++v )
166  myT[ v ] = 0.9 * newT[ v ] + 0.1 * myT[ v ];
167  enforceBounds();
168  auto Xnext = positions();
169  Scalar l2 = 0.0;
170  Scalar loo = 0.0;
171  for ( Vertex v = 0; v < myT.size(); ++v ) {
172  loo = std::max( loo, ( Xnext[ v ] - X[ v ] ).norm() );
173  l2 += ( Xnext[ v ] - X[ v ] ).squaredNorm();
174  }
175  return std::make_pair( loo, sqrt( l2 / myT.size() ) );
176 }
177 
178 template < typename TDigitalSurfaceContainer >
179 std::pair<double,double>
180 DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
181 oneStepSnakeMinimization
182 ( const double alpha, const double beta, const double randomization )
183 {
184  parameterize();
185  Scalars newT = myT;
186  for ( Vertex v = 0; v < myT.size(); ++v )
187  {
188  double right = 0.0;
189  double left = 0.0;
190  double coef = 0.0;
191  const auto s = myPtrIdxSurface->surfel( v );
192  const auto k = myPtrK->sOrthDir( s );
193  for ( Dimension i = 0; i < 3; ++i )
194  {
195  if ( i == k ) continue; // not a valid slice
196  const auto v_i = std::make_pair( v, i );
197  const auto vn_i = next( v_i );
198  const auto vnn_i = next( vn_i );
199  const auto vp_i = prev( v_i );
200  const auto vpp_i = prev( vp_i );
201  const auto c = c2_all( v_i );
202  const auto cn = c2_all( vn_i );
203  const auto cp = c2_all( vp_i );
204  const RealPoint Xnn = position( vnn_i.first );
205  const RealPoint Xn = position( vn_i.first );
206  const RealPoint X = position( v );
207  const RealPoint Xpp = position( vpp_i.first );
208  const RealPoint Xp = position( vp_i.first );
209  right += beta * ( - get<0>( c ) * get<0>( cn ) * Xnn[ k ]
210  + ( get<0>( c ) * get<1>( cn )
211  + get<1>( c ) * get<0>( c ) ) * Xn[ k ]
212  + ( get<2>( c ) * get<1>( cp )
213  + get<1>( c ) * get<2>( c ) ) * Xp[ k ]
214  - get<2>( c ) * get<2>( cp ) * Xpp[ k ] )
215  + alpha * ( get<0>( c ) * Xn[ k ] + get<2>( c ) * Xp[ k ] );
216  Scalar a = get<0>( c ) * get<2>( cn ) + get<1>( c ) * get<1>( c )
217  + get<2>( c ) * get<0>( cp );
218  left += ( beta * a + alpha * get<1>( c ) ) * myInsV[ v ][ k ];
219  coef += ( beta * a + alpha * get<1>( c ) )
220  * ( myOutV[ v ][ k ] - myInsV[ v ][ k ] );
221  }
222  // Possibly randomization to avoid local minima.
223  newT[ v ] = ( right - left ) / coef
224  + ( (double) rand() / (double) RAND_MAX - 0.49 ) * randomization;
225  }
226  auto X = positions();
227  // Damping between old and new positions.
228  for ( Vertex v = 0; v < myT.size(); ++v )
229  myT[ v ] = 0.5 * newT[ v ] + 0.5 * myT[ v ];
230  enforceBounds();
231  auto Xnext = positions();
232  Scalar l2 = 0.0;
233  Scalar loo = 0.0;
234  for ( Vertex v = 0; v < myT.size(); ++v ) {
235  loo = std::max( loo, ( Xnext[ v ] - X[ v ] ).norm() );
236  l2 += ( Xnext[ v ] - X[ v ] ).squaredNorm();
237  }
238  return std::make_pair( loo, sqrt( l2 / myT.size() ) );
239 }
240 
241 template < typename TDigitalSurfaceContainer >
242 std::pair<double,double>
243 DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
244 oneStepSquaredCurvatureMinimization
245 ( const double randomization )
246 {
247  parameterize();
248  Scalars newT = myT;
249  for ( Vertex v = 0; v < myT.size(); ++v )
250  {
251  double right = 0.0;
252  double left = 0.0;
253  double coef = 0.0;
254  const auto s = myPtrIdxSurface->surfel( v );
255  const auto k = myPtrK->sOrthDir( s );
256  for ( Dimension i = 0; i < 3; ++i )
257  {
258  if ( i == k ) continue; // not a valid slice
259  const auto v_i = std::make_pair( v, i );
260  const auto vn_i = next( v_i );
261  const auto vnn_i = next( vn_i );
262  const auto vp_i = prev( v_i );
263  const auto vpp_i = prev( vp_i );
264  const auto c_1 = c1( v_i );
265  const auto c = c2_all( v_i );
266  const auto cn = c2_all( vn_i );
267  const auto cp = c2_all( vp_i );
268  const RealPoint xnn = position( vnn_i.first );
269  const RealPoint xn = position( vn_i.first );
270  const RealPoint x = position( v );
271  const RealPoint xpp = position( vpp_i.first );
272  const RealPoint xp = position( vp_i.first );
273  const Scalar x_1 = c_1 * ( xn[ i ] - xp[ i ] );
274  const Scalar y_1 = c_1 * ( xn[ k ] - xp[ k ] );
275  const Scalar x_2 = get<0>( c ) * xn[ i ]
276  - get<1>( c ) * x[ i ] + get<2>( c ) * xp[ i ];
277  const Scalar xn_2 = get<0>( cn ) * xnn[ i ]
278  - get<1>( cn ) * xn[ i ] + get<2>( cn ) * x[ i ];
279  const Scalar xp_2 = get<0>( cp ) * x[ i ]
280  - get<1>( cp ) * xp[ i ] + get<2>( cp ) * xpp[ i ];
281  const Scalar x_3 = c_1 * ( xn_2 - xp_2 );
282  const Scalar x_4 = get<0>( c ) * xn_2
283  - get<1>( c ) * x_2 + get<2>( c ) * xp_2;
284  const Scalar y_2 = get<0>( c ) * xn[ k ]
285  - get<1>( c ) * x[ k ] + get<2>( c ) * xp[ k ];
286  const Scalar yn_2 = get<0>( cn ) * xnn[ k ]
287  - get<1>( cn ) * xn[ k ] + get<2>( cn ) * x[ k ];
288  const Scalar yp_2 = get<0>( cp ) * x[ k ]
289  - get<1>( cp ) * xp[ k ] + get<2>( cp ) * xpp[ k ];
290  const Scalar y_3 = c_1 * ( yn_2 - yp_2 );
291  const Scalar y_4 = get<0>( c ) * yn_2
292  - get<1>( c ) * y_2 + get<2>( c ) * yp_2;
293 
294  // We linearize Euler-Lagrange
295  // EL = 24*x'^3*x''^3*y'
296  // + x'''*( - 13*x'^4*x''*y' - 14*x'^2*x''*y'^3 - x''*y'^5 )
297  // + y'''*( 8*x'^5*x'' + 4*x'^3*x''*y'^2 - 4*x'*x''*y'^4 )
298  // + x''''*( x'^5*y' + 2*x'^3*y'^3 + x'*y'^5 )
299  // + y''*( 12*x'^4*y'*y''' + 12*x'^2*y'^3*y''' - 24*x'^4*x''^2 + 5*x'^5*x''' + 51*x'^2*x''^2*y'^2 - 2*x'^3*x'''*y'^2 + 3*x''^2*y'^4 - 7*x'*x'''*y'^4 )
300  // + y''^2*( -54*x'^3*x''*y' + 18*x'*x''*y'^3 )
301  // + y''^3*( 3*x'^4 - 21*x'^2*y'^2 )
302  // + y''''*( - x'^6 - 2*x'^4*y'^2 - x'^2*y'^4 )
303 
304  const Scalar x_1_pow2 = x_1 * x_1;
305  const Scalar x_1_pow3 = x_1_pow2 * x_1;
306  const Scalar x_1_pow4 = x_1_pow2 * x_1_pow2;
307  const Scalar x_1_pow5 = x_1_pow4 * x_1;
308  const Scalar x_1_pow6 = x_1_pow3 * x_1_pow3;
309  const Scalar x_2_pow2 = x_2 * x_2;
310  const Scalar x_2_pow3 = x_2_pow2 * x_1;
311  const Scalar y_1_pow2 = y_1 * y_1;
312  const Scalar y_1_pow3 = y_1_pow2 * y_1;
313  const Scalar y_1_pow4 = y_1_pow2 * y_1_pow2;
314  const Scalar y_1_pow5 = y_1_pow4 * y_1;
315  // EL = 24*x'^3*x''^3*y'
316  right += - 24. * x_1_pow3 * x_2_pow3 * y_1;
317  // + x'''*( - 13*x'^4*x''*y' - 14*x'^2*x''*y'^3 - x''*y'^5 )
318  right += - x_3*( -13. * x_1_pow4 * x_2 * y_1 - 14. * x_1_pow2 * x_2 * y_1_pow3
319  - x_2 * y_1_pow5 );
320  // + y'''*( 8*x'^5*x'' + 4*x'^3*x''*y'^2 - 4*x'*x''*y'^4 )
321  right += - y_3*( 8. * x_1_pow5 * x_2 + 4. * x_1_pow3 * x_2 * y_1_pow2
322  - 4. * x_1 * x_2 * y_1_pow4 );
323  // + x''''*( x'^5*y' + 2*x'^3*y'^3 + x'*y'^5 )
324  right += - x_4 * ( x_1_pow5 * y_1 + 2. * x_1_pow3 * y_1_pow3 + x_1 * y_1_pow5 );
325  // + y''*( 12*x'^4*y'*y''' + 12*x'^2*y'^3*y''' - 24*x'^4*x''^2 + 5*x'^5*x''' + 51*x'^2*x''^2*y'^2 - 2*x'^3*x'''*y'^2 + 3*x''^2*y'^4 - 7*x'*x'''*y'^4 )
326  const Scalar y_2_fact =
327  12. * x_1_pow4 * y_1 * y_3
328  + 12. * x_1_pow2 * y_1_pow3 * y_3 - 24. * x_1_pow4 * x_2_pow2
329  + 5. * x_1_pow5 * x_3 + 51. * x_1_pow2 * x_2_pow2 * y_1_pow2
330  - 2. * x_1_pow3 * x_3 * y_1_pow2 + 3. * x_2_pow2 * y_1_pow4
331  - 7. * x_1 * x_3 * y_1_pow4;
332  // + y''^2*( -54*x'^3*x''*y' + 18*x'*x''*y'^3 )
333  const Scalar y_2_pow2_fact =
334  - 54. * x_1_pow3 * x_2 * y_1 + 18. * x_1 * x_2 * y_1_pow3;
335  // + y''^3*( 3*x'^4 - 21*x'^2*y'^2 )
336  const Scalar y_2_pow3_fact =
337  3. * x_1_pow4 - 21. * x_1_pow2 * y_1_pow2;
338  // + y''''*( - x'^6 - 2*x'^4*y'^2 - x'^2*y'^4 )
339  const Scalar y_4_fact =
340  - x_1_pow6 - 2. * x_1_pow4 * y_1_pow2 - x_1_pow2 * y_1_pow4;
341  const Scalar gbl_y_2_fact =
342  y_2_fact + y_2_pow2_fact * y_2 + y_2_pow3_fact * y_2 * y_2
343  - get<1>( c ) * y_4_fact;
344  right += - y_4_fact * ( get<0>( c ) * yn_2 + get<2>( c ) * yp_2 );
345  right += - gbl_y_2_fact * ( get<0>( c ) * xn[ k ] + get<2>( c ) * xp[ k ] );
346  left += -get<1>( c ) * gbl_y_2_fact * myInsV[ v ][ k ];
347  coef += -get<1>( c ) * gbl_y_2_fact * ( myOutV[ v ][ k ] - myInsV[ v ][ k ] );
348 
349  }
350  // Possible randomization to avoid local minima.
351  newT[ v ] = ( right - left ) / coef
352  + ( (double) rand() / (double) RAND_MAX - 0.49 ) * randomization;
353  }
354  // Damping between old and new positions.
355  // Move vertices slightly toward optimal solution (since the
356  // problem has been linearized).
357  auto X = positions();
358  for ( Vertex v = 0; v < myT.size(); ++v )
359  myT[ v ] = 0.2 * newT[ v ] + 0.8 * myT[ v ];
360  enforceBounds();
361  auto Xnext = positions();
362  Scalar l2 = 0.0;
363  Scalar loo = 0.0;
364  for ( Vertex v = 0; v < myT.size(); ++v ) {
365  loo = std::max( loo, ( Xnext[ v ] - X[ v ] ).norm() );
366  l2 += ( Xnext[ v ] - X[ v ] ).squaredNorm();
367  }
368  return std::make_pair( loo, sqrt( l2 / myT.size() ) );
369 }
370 
371 template < typename TDigitalSurfaceContainer >
372 void
373 DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
374 enforceBounds()
375 {
376  for ( Vertex v = 0; v < myT.size(); ++v )
377  myT[ v ] = std::max( myEpsilon, std::min( 1.0 - myEpsilon, myT[ v ] ) );
378 }
Generated on Tue Jun 1 2021 17:05:42 for DGtal by doxygen 1.9.1