DGtal 1.4.0
Loading...
Searching...
No Matches
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
34template < typename TDigitalSurfaceContainer >
35double
36DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
37energyArea()
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
59template < typename TDigitalSurfaceContainer >
60double
61DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
62energySnake()
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
93template < typename TDigitalSurfaceContainer >
94double
95DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
96energySquaredCurvature()
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
127template < typename TDigitalSurfaceContainer >
128std::pair<double,double>
129DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
130oneStepAreaMinimization( 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
177template < typename TDigitalSurfaceContainer >
178std::pair<double,double>
179DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
180oneStepSnakeMinimization
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
240template < typename TDigitalSurfaceContainer >
241std::pair<double,double>
242DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
243oneStepSquaredCurvatureMinimization
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
372template < typename TDigitalSurfaceContainer >
373void
374DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
375enforceBounds()
376{
377 for ( Vertex v = 0; v < myT.size(); ++v )
378 myT[ v ] = std::max( myEpsilon, std::min( 1.0 - myEpsilon, myT[ v ] ) );
379}