DGtal  1.5.beta
DigitalSurfaceRegularization.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
19  * @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr )
20  * Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
21  *
22  * @date 2019/10/25
23  *
24  * Implementation of inline methods defined in DigitalSurfaceRegularization.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 #include <algorithm>
33 //////////////////////////////////////////////////////////////////////////////
34 
35 ///////////////////////////////////////////////////////////////////////////////
36 // Interface - public :
37 
38 /**
39  * Writes/Displays the object on an output stream.
40  * @param out the output stream where the object is written.
41  */
42 template <typename T>
43 inline
44 void
45 DGtal::DigitalSurfaceRegularization<T>::selfDisplay ( std::ostream & out ) const
46 {
47  out << "[DigitalSurfaceRegularization] alpha= "<<myAlpha
48  <<" beta= "<<myBeta <<" gamma= "<<myGamma
49  << " number of surfels= "<<mySurfels.size();
50 }
51 
52 /**
53  * Checks the validity/consistency of the object.
54  * @return 'true' if the object is valid, 'false' otherwise.
55  */
56 template <typename T>
57 inline
58 bool
59 DGtal::DigitalSurfaceRegularization<T>::isValid() const
60 {
61  return myDigitalSurface.isValid();
62 }
63 
64 ///////////////////////////////////////////////////////////////////////////////
65 // Implementation of inline functions //
66 ///////////////////////////////////////////////////////////////////////////////
67 template <typename T>
68 inline
69 void DGtal::DigitalSurfaceRegularization<T>::attachConvolvedTrivialNormalVectors(const Parameters someParams)
70 {
71  ASSERT_MSG(myInit, "The init() method must be called before setting the normals");
72  myNormals = SHG3::getCTrivialNormalVectors( myDigitalSurface, mySurfels, someParams);
73 }
74 
75 ///////////////////////////////////////////////////////////////////////////////
76 template <typename T>
77 inline
78 void
79 DGtal::DigitalSurfaceRegularization<T>::attachNormalVectors(const std::function<SHG3::RealVector(SH3::SCell&)> &normalFunc)
80 {
81  ASSERT_MSG(myInit, "The init() method must be called before setting the normals");
82  myNormals.resize( mySurfels.size());
83  for(auto i=0u; i< mySurfels.size();++i)
84  myNormals[i] = normalFunc( mySurfels[i] );
85 }
86 ///////////////////////////////////////////////////////////////////////////////
87 //
88 template <typename T>
89 inline
90 void DGtal::DigitalSurfaceRegularization<T>::cacheInit()
91 {
92  //Collecting the surfels range
93  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
94  mySurfels.clear();
95  mySurfels = SH3::getSurfelRange( myDigitalSurface, params );
96 
97  myPointelIndex.clear();
98  //Collecting the input points
99  auto pointels = SH3::getPointelRange(myPointelIndex,myDigitalSurface);
100  auto embedder = SH3::getCellEmbedder(myK); // /!\ no grid step here
101  myOriginalPositions.clear();
102  std::transform(pointels.begin(), pointels.end(), std::back_inserter(myOriginalPositions),
103  [embedder](SH3::Cell &pointel ) {return embedder(pointel); } );
104 
105  //Init the regularized positions.
106  myRegularizedPositions.resize(myOriginalPositions.size());
107  std::copy(myOriginalPositions.begin(), myOriginalPositions.end(), myRegularizedPositions.begin());
108 
109  //Allocating Gradient vector
110  myGradient.clear();
111  myGradientAlign.clear();
112  myGradient.resize(myOriginalPositions.size());
113  myGradientAlign.resize(myOriginalPositions.size());
114 
115  /////
116  ///Cacheing some topological information
117  auto polySurf = SH3::makeDualPolygonalSurface(mySurfelIndex,myDigitalSurface);
118  myFaces = polySurf->allFaces() ; //All faces of the dual digital surface
119  auto dsurf_faces = myDigitalSurface->allClosedFaces(); // Faces of digital surface (umbrellas)
120  // dsurf_pointels list the pointels in the same order as the faces of polySurf.
121  std::vector<SH3::Cell> dsurf_pointels( dsurf_faces.size() );
122  std::transform( dsurf_faces.cbegin(), dsurf_faces.cend(),
123  dsurf_pointels.begin(),
124  [&] ( const SH3::DigitalSurface::Face f ) { return myK.unsigns(myDigitalSurface->pivot( f )); } );
125 
126 
127  myNumberAdjEdgesToPointel.resize(myOriginalPositions.size(),0);
128 
129  // Precompute all relations for align energy
130  myAlignPointelsIdx.resize( mySurfels.size() * 4 );
131  myAlignPointels.resize( mySurfels.size() * 4 );
132 
133  for(size_t i = 0; i < mySurfels.size(); ++i)
134  {
135  auto pointelsAround = myDigitalSurface->facesAroundVertex(mySurfels[i],true);
136  ASSERT(pointelsAround.size() == 4);
137  for(auto j = 0; j < 4; ++j)
138  {
139  auto p = myK.unsigns(myDigitalSurface->pivot(pointelsAround[ j ]));
140  auto cell_p = myPointelIndex[ p ];
141  myAlignPointelsIdx[ 4*i + j ] = cell_p;
142  myAlignPointels[4*i + j] = p;
143  myNumberAdjEdgesToPointel[ cell_p ] ++;
144  }
145  }
146 
147  // Precompute all relations for fairness energy
148  myNbAdjacent.resize( myFaces.size() );
149  for(size_t faceId=0 ; faceId < myFaces.size(); ++faceId)
150  {
151  auto idx = myPointelIndex[ dsurf_pointels[ faceId ] ];
152  myFairnessPointelsIdx.push_back( idx );
153  unsigned char nbAdj = 0;
154  auto arcs = polySurf->arcsAroundFace(faceId);
155  for(auto anArc : arcs)
156  {
157  // /!\ we must check that the opposite exists (aka surface with boundaries)
158  auto op = polySurf->opposite(anArc);
159  auto adjFace = polySurf->faceAroundArc(op);
160  ASSERT(adjFace != faceId);
161  myFairnessPointelsIdx.push_back( myPointelIndex[ dsurf_pointels[ adjFace] ] );
162  nbAdj++;
163  }
164  ASSERT(nbAdj>0);
165  myNbAdjacent[ faceId ] = nbAdj;
166  }
167 }
168 ///////////////////////////////////////////////////////////////////////////////
169 //
170 template <typename T>
171 inline
172 void DGtal::DigitalSurfaceRegularization<T>::init(const double alpha,
173  const double beta,
174  const double gamma)
175 {
176  myAlpha = alpha;
177  myBeta = beta;
178  myGamma = gamma;
179  myInit = true;
180  myConstantCoeffs = true;
181  cacheInit();
182 }
183 ///////////////////////////////////////////////////////////////////////////////
184 //
185 template <typename T>
186 inline
187 void DGtal::DigitalSurfaceRegularization<T>::init(ConstAlias< std::vector<double> > alphas,
188  ConstAlias< std::vector<double> > betas,
189  ConstAlias< std::vector<double> > gammas)
190 {
191  myAlpha = 0.0;
192  myBeta = 0.0;
193  myGamma = 0.0;
194 
195  myAlphas = &alphas;
196  myBetas = &betas;
197  myGammas = &gammas;
198 
199  myInit = true;
200  myConstantCoeffs = false;
201  cacheInit();
202 
203  ASSERT_MSG(alphas->size() == myOriginalPositions.size(), "The alpha vector size must equal the number of pointels");
204  ASSERT_MSG(betas->size() == myOriginalPositions.size(), "The beta vector size must equal the number of pointels");
205  ASSERT_MSG(gammas->size() == myOriginalPositions.size(), "The gamma vector size must equal the number of pointels");
206 }
207 ///////////////////////////////////////////////////////////////////////////////
208 template <typename T>
209 inline
210 double
211 DGtal::DigitalSurfaceRegularization<T>::computeGradient()
212 {
213  double energy= 0.0;
214  auto zero = SH3::RealPoint(0,0,0);
215 
216  ASSERT_MSG(myInit, "The init() method must be called before computing the gradient");
217  ASSERT_MSG(myNormals.size() != 0, "Some normal vectors must be attached to the digital surface before computing the gradient");
218 
219  //data attachment term
220  for(size_t i = 0; i < myOriginalPositions.size(); ++i)
221  {
222  const auto delta_d = myOriginalPositions[i] - myRegularizedPositions[i];
223  energy += myAlpha * delta_d.squaredNorm() ;
224  myGradient[i] = 2.0*myAlpha * delta_d;
225  //reset for align
226  myGradientAlign[i] = zero;
227  }
228 
229  //align
230  auto it = myAlignPointelsIdx.cbegin();
231  for(size_t i = 0; i < mySurfels.size(); ++i)
232  {
233  const auto cell_p0 = *it++;
234  const auto cell_p1 = *it++;
235  const auto cell_p2 = *it++;
236  const auto cell_p3 = *it++;
237  const auto e0 = myRegularizedPositions[ cell_p0 ] - myRegularizedPositions[ cell_p1 ];
238  const auto e1 = myRegularizedPositions[ cell_p1 ] - myRegularizedPositions[ cell_p2 ];
239  const auto e2 = myRegularizedPositions[ cell_p2 ] - myRegularizedPositions[ cell_p3 ];
240  const auto e3 = myRegularizedPositions[ cell_p3 ] - myRegularizedPositions[ cell_p0 ];
241  const auto cos_a0 = e0.dot( myNormals[i] );
242  const auto cos_a1 = e1.dot( myNormals[i] );
243  const auto cos_a2 = e2.dot( myNormals[i] );
244  const auto cos_a3 = e3.dot( myNormals[i] );
245  energy += myBeta * ( cos_a0 * cos_a0 + cos_a1 * cos_a1
246  + cos_a2 * cos_a2 + cos_a3 * cos_a3 );
247  myGradientAlign[ cell_p0 ] += cos_a0 * myNormals[i];
248  myGradientAlign[ cell_p1 ] += cos_a1 * myNormals[i];
249  myGradientAlign[ cell_p2 ] += cos_a2 * myNormals[i];
250  myGradientAlign[ cell_p3 ] += cos_a3 * myNormals[i];
251  }
252  for(auto i=0u; i < myOriginalPositions.size(); ++i)
253  {
254  ASSERT(myNumberAdjEdgesToPointel[i] >0);
255  myGradient[i] += 2.0*myBeta * myGradientAlign[i] / (double)myNumberAdjEdgesToPointel[i];
256  }
257 
258  //fairness
259  auto itP = myFairnessPointelsIdx.cbegin();
260  auto itNb = myNbAdjacent.cbegin();
261  SH3::RealPoint barycenter, phat ;
262  for(size_t faceId=0 ; faceId < myFaces.size(); ++faceId)
263  {
264  auto idx = *itP++;
265  unsigned int nbAdj = *itNb++;
266  barycenter = zero;
267  phat = myRegularizedPositions[ idx ];
268  for ( unsigned int i = 0; i < nbAdj; ++i )
269  barycenter += myRegularizedPositions[ *itP++ ];
270  ASSERT(nbAdj>0);
271  barycenter /= (double)nbAdj;
272  auto delta_f = phat - barycenter;
273  energy += myGamma * delta_f.squaredNorm() ;
274  myGradient[ idx ] += 2.0*myGamma * delta_f;
275  }
276 
277  return energy;
278 }
279 ///////////////////////////////////////////////////////////////////////////////
280 template <typename T>
281 inline
282 double
283 DGtal::DigitalSurfaceRegularization<T>::computeGradientLocalWeights()
284 {
285  double energy= 0.0;
286  auto zero = SH3::RealPoint(0,0,0);
287 
288  ASSERT_MSG(myInit, "The init() method must be called before computing the gradient");
289  ASSERT_MSG(myNormals.size() != 0, "Some normal vectors must be attached to the digital surface before computing the gradient");
290 
291  //data attachment term
292  for(size_t i = 0; i < myOriginalPositions.size(); ++i)
293  {
294  const auto delta_d = myOriginalPositions[i] - myRegularizedPositions[i];
295  energy += (*myAlphas)[i] * delta_d.squaredNorm() ;
296  myGradient[i] = 2.0*(*myAlphas)[i] * delta_d;
297  //reset for align
298  myGradientAlign[i] = zero;
299  }
300 
301  //align
302  auto it = myAlignPointelsIdx.cbegin();
303  for(size_t i = 0; i < mySurfels.size(); ++i)
304  {
305  const auto cell_p0 = *it++;
306  const auto cell_p1 = *it++;
307  const auto cell_p2 = *it++;
308  const auto cell_p3 = *it++;
309  const auto e0 = myRegularizedPositions[ cell_p0 ] - myRegularizedPositions[ cell_p1 ];
310  const auto e1 = myRegularizedPositions[ cell_p1 ] - myRegularizedPositions[ cell_p2 ];
311  const auto e2 = myRegularizedPositions[ cell_p2 ] - myRegularizedPositions[ cell_p3 ];
312  const auto e3 = myRegularizedPositions[ cell_p3 ] - myRegularizedPositions[ cell_p0 ];
313  const auto cos_a0 = e0.dot( myNormals[i] );
314  const auto cos_a1 = e1.dot( myNormals[i] );
315  const auto cos_a2 = e2.dot( myNormals[i] );
316  const auto cos_a3 = e3.dot( myNormals[i] );
317  energy += (*myBetas)[ cell_p0 ] * cos_a0 * cos_a0
318  + (*myBetas)[ cell_p1 ] * cos_a1 * cos_a1
319  + (*myBetas)[ cell_p2 ] * cos_a2 * cos_a2
320  + (*myBetas)[ cell_p3 ] * cos_a3 * cos_a3;
321  myGradientAlign[ cell_p0 ] += cos_a0 * myNormals[i];
322  myGradientAlign[ cell_p1 ] += cos_a1 * myNormals[i];
323  myGradientAlign[ cell_p2 ] += cos_a2 * myNormals[i];
324  myGradientAlign[ cell_p3 ] += cos_a3 * myNormals[i];
325  }
326  for(size_t i=0; i < myOriginalPositions.size(); ++i)
327  {
328  ASSERT(myNumberAdjEdgesToPointel[i] >0);
329  myGradient[i] += 2.0*(*myBetas)[i] * myGradientAlign[i] / (double)myNumberAdjEdgesToPointel[i];
330  }
331 
332  //fairness
333  auto itP = myFairnessPointelsIdx.cbegin();
334  auto itNb = myNbAdjacent.cbegin();
335  SH3::RealPoint barycenter, phat ;
336  for(auto faceId=0u; faceId < myFaces.size(); ++faceId)
337  {
338  const auto idx = *itP++;
339  const unsigned int nbAdj = *itNb++;
340  barycenter = zero;
341  phat = myRegularizedPositions[ idx ];
342  for ( unsigned int i = 0; i < nbAdj; ++i )
343  barycenter += myRegularizedPositions[ *itP++ ];
344  ASSERT(nbAdj>0);
345  barycenter /= (double)nbAdj;
346  auto delta_f = phat - barycenter;
347  energy += (*myGammas)[ idx ] * delta_f.squaredNorm() ;
348  myGradient[ idx ] += 2.0*(*myGammas)[ idx ] * delta_f;
349  }
350 
351  return energy;
352 }
353 
354 ///////////////////////////////////////////////////////////////////////////////
355 template <typename T>
356 template <typename AdvectionFunction>
357 inline
358 double
359 DGtal::DigitalSurfaceRegularization<T>::regularize(const unsigned int nbIters,
360  const double dt,
361  const double epsilon,
362  const AdvectionFunction &advectionFunc)
363 {
364  double energy = 0.0;
365  double last_energy = 0.0;
366  double mydt=dt;
367  bool first_iter = true;
368  for(unsigned int i = 0; i < nbIters; ++i)
369  {
370  if (myConstantCoeffs)
371  energy = computeGradient();
372  else
373  energy = computeGradientLocalWeights();
374 
375  double gradnorm=0.0;
376  for(auto &v: myGradient)
377  gradnorm = std::max(gradnorm, v.norm());
378 
379  if (myVerbose)
380  trace.info()<< "Step " << i
381  << " dt=" << mydt
382  << " energy = " << energy
383  << " gradnorm= " << gradnorm << std::endl;
384 
385 
386  //Naive linesearch by doubling the learning rate
387  if ( ! first_iter )
388  mydt *= ( energy > last_energy ) ? 0.5 : 1.1;
389 
390  //Stopping criterion
391  if ( mydt < epsilon ) return energy;
392 
393  last_energy = energy;
394  first_iter = false;
395 
396  //One step advection
397  SHG3::RealVector v;
398  for(size_t ii=0; ii < myRegularizedPositions.size(); ++ii)
399  {
400  v = - mydt * myGradient[ii] ;
401  advectionFunc( myRegularizedPositions[ii], myOriginalPositions[ii], v );
402  }
403  }
404  return energy;
405 }
406 ///////////////////////////////////////////////////////////////////////////////
407 
408 template <typename T>
409 inline
410 std::ostream&
411 DGtal::operator<< ( std::ostream & out,
412  const DGtal::DigitalSurfaceRegularization<T> & object )
413 {
414  object.selfDisplay( out );
415  return out;
416 }
417 // //
418 ///////////////////////////////////////////////////////////////////////////////
419