DGtal 1.4.0
Loading...
Searching...
No Matches
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//-----------------------------------------------------------------------------
42template <typename TRealPoint, typename TRealVector>
43DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
44CorrectedNormalCurrentComputer( ConstAlias< SurfaceMesh > aMesh,
45 bool unit_u )
46 : myMesh( aMesh ), myUnitU( unit_u )
47{}
48
49//-----------------------------------------------------------------------------
50template <typename TRealPoint, typename TRealVector>
51typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
52DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
53computeMu0() 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//-----------------------------------------------------------------------------
64template <typename TRealPoint, typename TRealVector>
65typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
66DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
67computeMu1() 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//-----------------------------------------------------------------------------
78template <typename TRealPoint, typename TRealVector>
79typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
80DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
81computeMu2() 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//-----------------------------------------------------------------------------
92template <typename TRealPoint, typename TRealVector>
93typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::TensorMeasure
94DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
95computeMuXY() 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//-----------------------------------------------------------------------------
107template <typename TRealPoint, typename TRealVector>
108typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
109DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
110computeMu0InterpolatedU() 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//-----------------------------------------------------------------------------
132template <typename TRealPoint, typename TRealVector>
133typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
134DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
135computeMu1InterpolatedU() 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//-----------------------------------------------------------------------------
157template <typename TRealPoint, typename TRealVector>
158typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
159DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
160computeMu2InterpolatedU() 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//-----------------------------------------------------------------------------
182template <typename TRealPoint, typename TRealVector>
183typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::TensorMeasure
184DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
185computeMuXYInterpolatedU() 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//-----------------------------------------------------------------------------
208template <typename TRealPoint, typename TRealVector>
209typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
210DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
211computeMu0ConstantU() 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//-----------------------------------------------------------------------------
231template <typename TRealPoint, typename TRealVector>
232typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
233DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
234computeMu1ConstantU() 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//-----------------------------------------------------------------------------
261template <typename TRealPoint, typename TRealVector>
262typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::ScalarMeasure
263DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
264computeMu2ConstantU() 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//-----------------------------------------------------------------------------
318template <typename TRealPoint, typename TRealVector>
319typename DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::TensorMeasure
320DGtal::CorrectedNormalCurrentComputer<TRealPoint, TRealVector>::
321computeMuXYConstantU() 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///////////////////////////////////////////////////////////////////////////////