DGtal 1.3.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 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
128template < typename TDigitalSurfaceContainer >
129std::pair<double,double>
130DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
131oneStepAreaMinimization( 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
178template < typename TDigitalSurfaceContainer >
179std::pair<double,double>
180DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
181oneStepSnakeMinimization
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
241template < typename TDigitalSurfaceContainer >
242std::pair<double,double>
243DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
244oneStepSquaredCurvatureMinimization
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
371template < typename TDigitalSurfaceContainer >
372void
373DGtal::ShroudsRegularization< TDigitalSurfaceContainer >::
374enforceBounds()
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 Sun Nov 27 2022 15:15:39 for DGtal by doxygen 1.9.5