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