DGtal  1.5.beta
CorrectedNormalCurrentComputer.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 CorrectedNormalCurrentComputer.ih
19  * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20  * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
21  *
22  * @date 2020/02/18
23  *
24  * Implementation of inline methods defined in CorrectedNormalCurrentComputer.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 //////////////////////////////////////////////////////////////////////////////
33 
34 ///////////////////////////////////////////////////////////////////////////////
35 // IMPLEMENTATION of inline methods.
36 ///////////////////////////////////////////////////////////////////////////////
37 
38 ///////////////////////////////////////////////////////////////////////////////
39 // ----------------------- Standard services ------------------------------
40 
41 //-----------------------------------------------------------------------------
42 template <typename TRealPoint, typename TRealVector>
43 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
44 CorrectedNormalCurrentComputer( ConstAlias< SurfaceMesh > aMesh,
45  bool unit_u )
46  : myMesh( aMesh ), myUnitU( unit_u )
47 {}
48 
49 //-----------------------------------------------------------------------------
50 template <typename TRealPoint, typename TRealVector>
51 typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
52 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
53 computeMu0() const
54 {
55  if ( ! myMesh.vertexNormals().empty() ) return computeMu0InterpolatedU();
56  if ( ! myMesh.faceNormals().empty() ) return computeMu0ConstantU();
57  trace.warning() << "[CorrectedNormalCurrentComputer::computeInterpolatedMu0]"
58  << " Unable to compute measures without vertex or face normals."
59  << std::endl;
60  return ScalarMeasure( &myMesh, 0.0 );
61 }
62 
63 //-----------------------------------------------------------------------------
64 template <typename TRealPoint, typename TRealVector>
65 typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
66 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
67 computeMu1() const
68 {
69  if ( ! myMesh.vertexNormals().empty() ) return computeMu1InterpolatedU();
70  if ( ! myMesh.faceNormals().empty() ) return computeMu1ConstantU();
71  trace.warning() << "[CorrectedNormalCurrentComputer::computeInterpolatedMu1]"
72  << " Unable to compute measures without vertex or face normals."
73  << std::endl;
74  return ScalarMeasure( &myMesh, 0.0 );
75 }
76 
77 //-----------------------------------------------------------------------------
78 template <typename TRealPoint, typename TRealVector>
79 typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
80 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
81 computeMu2() const
82 {
83  if ( ! myMesh.vertexNormals().empty() ) return computeMu2InterpolatedU();
84  if ( ! myMesh.faceNormals().empty() ) return computeMu2ConstantU();
85  trace.warning() << "[CorrectedNormalCurrentComputer::computeInterpolatedMu2]"
86  << " Unable to compute measures without vertex or face normals."
87  << std::endl;
88  return ScalarMeasure( &myMesh, 0.0 );
89 }
90 
91 //-----------------------------------------------------------------------------
92 template <typename TRealPoint, typename TRealVector>
93 typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::TensorMeasure
94 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
95 computeMuXY() const
96 {
97  const RealTensor zeroT { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
98  if ( ! myMesh.vertexNormals().empty() ) return computeMuXYInterpolatedU();
99  if ( ! myMesh.faceNormals().empty() ) return computeMuXYConstantU();
100  trace.warning() << "[CorrectedNormalCurrentComputer::computeInterpolatedMuXY]"
101  << " Unable to compute measures without vertex or face normals."
102  << std::endl;
103  return TensorMeasure( &myMesh, zeroT );
104 }
105 
106 //-----------------------------------------------------------------------------
107 template <typename TRealPoint, typename TRealVector>
108 typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
109 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
110 computeMu0InterpolatedU() const
111 {
112  ScalarMeasure mu0( &myMesh, 0.0 );
113  ASSERT( ! myMesh.vertexNormals().empty() );
114  auto& face_mu0 = mu0.kMeasures( 2 );
115  face_mu0.resize( myMesh.nbFaces() );
116  Index idx_f = 0;
117  for ( const auto& f : myMesh.allIncidentVertices() )
118  {
119  RealPoints p( f.size() );
120  RealVectors u( f.size() );
121  for ( Index idx_v = 0; idx_v < f.size(); ++idx_v )
122  {
123  p[ idx_v ] = myMesh.positions() [ f[ idx_v ] ];
124  u[ idx_v ] = myMesh.vertexNormals()[ f[ idx_v ] ];
125  }
126  face_mu0[ idx_f++ ] = Formula::mu0InterpolatedU( p, u, myUnitU );
127  }
128  return mu0;
129 }
130 
131 //-----------------------------------------------------------------------------
132 template <typename TRealPoint, typename TRealVector>
133 typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
134 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
135 computeMu1InterpolatedU() const
136 {
137  ScalarMeasure mu1( &myMesh, 0.0 );
138  ASSERT( ! myMesh.vertexNormals().empty() );
139  auto& face_mu1 = mu1.kMeasures( 2 );
140  face_mu1.resize( myMesh.nbFaces() );
141  Index idx_f = 0;
142  for ( const auto& f : myMesh.allIncidentVertices() )
143  {
144  RealPoints p( f.size() );
145  RealVectors u( f.size() );
146  for ( Index idx_v = 0; idx_v < f.size(); ++idx_v )
147  {
148  p[ idx_v ] = myMesh.positions() [ f[ idx_v ] ];
149  u[ idx_v ] = myMesh.vertexNormals()[ f[ idx_v ] ];
150  }
151  face_mu1[ idx_f++ ] = Formula::mu1InterpolatedU( p, u, myUnitU );
152  }
153  return mu1;
154 }
155 
156 //-----------------------------------------------------------------------------
157 template <typename TRealPoint, typename TRealVector>
158 typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
159 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
160 computeMu2InterpolatedU() const
161 {
162  ScalarMeasure mu2( &myMesh, 0.0 );
163  ASSERT( ! myMesh.vertexNormals().empty() );
164  auto& face_mu2 = mu2.kMeasures( 2 );
165  face_mu2.resize( myMesh.nbFaces() );
166  Index idx_f = 0;
167  for ( const auto& f : myMesh.allIncidentVertices() )
168  {
169  RealPoints p( f.size() );
170  RealVectors u( f.size() );
171  for ( Index idx_v = 0; idx_v < f.size(); ++idx_v )
172  {
173  p[ idx_v ] = myMesh.positions() [ f[ idx_v ] ];
174  u[ idx_v ] = myMesh.vertexNormals()[ f[ idx_v ] ];
175  }
176  face_mu2[ idx_f++ ] = Formula::mu2InterpolatedU( p, u, myUnitU );
177  }
178  return mu2;
179 }
180 
181 //-----------------------------------------------------------------------------
182 template <typename TRealPoint, typename TRealVector>
183 typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::TensorMeasure
184 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
185 computeMuXYInterpolatedU() const
186 {
187  const RealTensor zeroT { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
188  TensorMeasure muXY( &myMesh, zeroT );
189  ASSERT( ! myMesh.vertexNormals().empty() );
190  auto& face_muXY = muXY.kMeasures( 2 );
191  face_muXY.resize( myMesh.nbFaces() );
192  Index idx_f = 0;
193  for ( const auto& f : myMesh.allIncidentVertices() )
194  {
195  RealPoints p( f.size() );
196  RealVectors u( f.size() );
197  for ( Index idx_v = 0; idx_v < f.size(); ++idx_v )
198  {
199  p[ idx_v ] = myMesh.positions() [ f[ idx_v ] ];
200  u[ idx_v ] = myMesh.vertexNormals()[ f[ idx_v ] ];
201  }
202  face_muXY[ idx_f++ ] = Formula::muXYInterpolatedU( p, u, myUnitU );
203  }
204  return muXY;
205 }
206 
207 //-----------------------------------------------------------------------------
208 template <typename TRealPoint, typename TRealVector>
209 typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
210 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
211 computeMu0ConstantU() const
212 {
213  ScalarMeasure mu0( &myMesh, 0.0 );
214  ASSERT( ! myMesh.faceNormals().empty() );
215  auto& face_mu0 = mu0.kMeasures( 2 );
216  face_mu0.resize( myMesh.nbFaces() );
217  Index idx_f = 0;
218  for ( const auto& f : myMesh.allIncidentVertices() )
219  {
220  RealPoints p( f.size() );
221  const RealVector& u = myMesh.faceNormal( idx_f );
222  for ( Index idx_v = 0; idx_v < f.size(); ++idx_v )
223  p[ idx_v ] = myMesh.positions() [ f[ idx_v ] ];
224  face_mu0[ idx_f ] = Formula::mu0ConstantU( p, u );
225  ++idx_f;
226  }
227  return mu0;
228 }
229 
230 //-----------------------------------------------------------------------------
231 template <typename TRealPoint, typename TRealVector>
232 typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
233 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
234 computeMu1ConstantU() const
235 {
236  ScalarMeasure mu1( &myMesh, 0.0 );
237  ASSERT( ! myMesh.faceNormals().empty() );
238  auto& edge_mu1 = mu1.kMeasures( 1 );
239  edge_mu1.resize( myMesh.nbEdges() );
240 
241  for ( Index idx_e = 0; idx_e < myMesh.nbEdges(); ++idx_e )
242  {
243  const auto& right_f = myMesh.edgeRightFaces( idx_e );
244  const auto& left_f = myMesh.edgeLeftFaces ( idx_e );
245  if ( right_f.size() == 1 && left_f.size() == 1 )
246  {
247  const auto ab = myMesh.edgeVertices( idx_e );
248  const RealPoint& xa = myMesh.position( ab.first );
249  const RealPoint& xb = myMesh.position( ab.second );
250  const RealVector& ur = myMesh.faceNormal( right_f.front() );
251  const RealVector& ul = myMesh.faceNormal( left_f.front() );
252  edge_mu1[ idx_e ] = Formula::mu1ConstantUAtEdge( xa, xb, ur, ul );
253  }
254  else
255  edge_mu1[ idx_e ] = 0.0;
256  }
257  return mu1;
258 }
259 
260 //-----------------------------------------------------------------------------
261 template <typename TRealPoint, typename TRealVector>
262 typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
263 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
264 computeMu2ConstantU() const
265 {
266  ScalarMeasure mu2( &myMesh, 0.0 );
267  ASSERT( ! myMesh.faceNormals().empty() );
268  auto& vertex_mu2 = mu2.kMeasures( 0 );
269  vertex_mu2.resize( myMesh.nbVertices() );
270  Index idx_v = 0;
271  for ( const auto& faces_v : myMesh.allIncidentFaces() )
272  {
273  const RealPoint a = myMesh.positions()[ idx_v ];
274  std::vector< Index > faces;
275  std::vector< Index > prev;
276  std::vector< Index > next;
277  for ( auto f : faces_v )
278  {
279  const auto & vtcs = myMesh.allIncidentVertices()[ f ];
280  const auto nbv = vtcs.size();
281  Index j = std::find( vtcs.cbegin(), vtcs.cend(), idx_v ) - vtcs.cbegin();
282  if ( j == nbv ) continue;
283  faces.push_back( f );
284  prev.push_back( vtcs[ ( j + nbv - 1 ) % nbv ] );
285  next.push_back( vtcs[ ( j + nbv + 1 ) % nbv ] );
286  }
287  // Try to reorder faces as an umbrella. If this is not possible,
288  // the vertex is not a manifold point and its measure is set to
289  // 0.
290  bool manifold = true;
291  const Index nb = faces.size();
292  for ( Index i = 1; i < nb && manifold; i++ )
293  {
294  Index j = std::find( next.cbegin() + i, next.cend(), prev[ i - 1 ] )
295  - next.cbegin();
296  if ( j == nb ) manifold = false;
297  else if ( j > i )
298  {
299  std::swap( faces[ i ], faces[ j ] );
300  std::swap( next [ i ], next [ j ] );
301  std::swap( prev [ i ], prev [ j ] );
302  }
303  }
304  vertex_mu2[ idx_v ] = 0.0;
305  if ( manifold )
306  {
307  RealVectors vu( nb );
308  for ( Index i = 0; i < nb; ++i )
309  vu[ i ] = myMesh.faceNormal( faces[ i ] );
310  vertex_mu2[ idx_v ] = Formula::mu2ConstantUAtVertex( a, vu );
311  }
312  idx_v++;
313  }
314  return mu2;
315 }
316 
317 //-----------------------------------------------------------------------------
318 template <typename TRealPoint, typename TRealVector>
319 typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::TensorMeasure
320 DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
321 computeMuXYConstantU() const
322 {
323  ASSERT( ! myMesh.faceNormals().empty() );
324  const RealTensor zeroT { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
325  TensorMeasure muXY( &myMesh, zeroT );
326  ASSERT( ! myMesh.faceNormals().empty() );
327  auto& edge_muXY = muXY.kMeasures( 1 );
328  edge_muXY.resize( myMesh.nbEdges() );
329 
330  for ( Index idx_e = 0; idx_e < myMesh.nbEdges(); ++idx_e )
331  {
332  const auto& right_f = myMesh.edgeRightFaces( idx_e );
333  const auto& left_f = myMesh.edgeLeftFaces ( idx_e );
334  if ( right_f.size() == 1 && left_f.size() == 1 )
335  {
336  const auto ab = myMesh.edgeVertices( idx_e );
337  const RealPoint& xa = myMesh.position( ab.first );
338  const RealPoint& xb = myMesh.position( ab.second );
339  const RealVector& ur = myMesh.faceNormal( right_f.front() );
340  const RealVector& ul = myMesh.faceNormal( left_f.front() );
341  edge_muXY[ idx_e ] = Formula::muXYConstantUAtEdge( xa, xb, ur, ul );
342  }
343  else
344  edge_muXY[ idx_e ] = zeroT;
345  }
346  return muXY;
347 }
348 
349 ///////////////////////////////////////////////////////////////////////////////
350 ///////////////////////////////////////////////////////////////////////////////