DGtal 1.4.0
Loading...
Searching...
No Matches
DiscreteExteriorCalculus.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 DiscreteExteriorCalculus.ih
19 * @author Pierre Gueth (\c pierre.gueth@liris.cnrs.fr )
20 * Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
21 *
22 * @date 2014/03/27
23 *
24 * Implementation of inline methods defined in DiscreteExteriorCalculus.h
25 *
26 * This file is part of the DGtal library.
27 */
28
29///////////////////////////////////////////////////////////////////////////////
30// IMPLEMENTATION of inline methods.
31///////////////////////////////////////////////////////////////////////////////
32
33template <DGtal::Dimension dim, typename TInteger>
34size_t
35DGtal::hash_value(const DGtal::KhalimskyCell<dim, TInteger>& cell)
36{
37 return std::hash<size_t>()( boost::hash_range( cell.preCell().coordinates.begin(), cell.preCell().coordinates.end()) );
38}
39
40///////////////////////////////////////////////////////////////////////////////
41// ----------------------- Standard services ------------------------------
42
43template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
44DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::DiscreteExteriorCalculus()
45 : myKSpace(), myCachedOperatorsNeedUpdate(true), myIndexesNeedUpdate(false)
46{
47}
48
49template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
50template <typename TDomain>
51void
52DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::initKSpace(DGtal::ConstAlias<TDomain> _domain)
53{
54 BOOST_CONCEPT_ASSERT(( concepts::CDomain<TDomain> ));
55
56 // FIXME borders are removed from set => better not initialize kspace
57 // FIXME should be open or closed? => closed = true
58 const bool kspace_init_ok = const_cast<KSpace&>(myKSpace).init(_domain->lowerBound(), _domain->upperBound(), true);
59 ASSERT(kspace_init_ok);
60 boost::ignore_unused_variable_warning(kspace_init_ok);
61}
62
63///////////////////////////////////////////////////////////////////////////////
64// Interface - public :
65
66template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
67bool
68DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::eraseCell(const Cell& _cell)
69{
70 typename Properties::iterator iter_property = myCellProperties.find(_cell);
71 if (iter_property == myCellProperties.end())
72 return false;
73
74 myCellProperties.erase(iter_property);
75
76 myIndexesNeedUpdate = true;
77 myCachedOperatorsNeedUpdate = true;
78
79 return true;
80}
81
82template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
83bool
84DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::insertSCell(const SCell& signed_cell)
85{
86 return insertSCell(signed_cell, 1, 1);
87}
88
89template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
90bool
91DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::insertSCell(const SCell& signed_cell, const Scalar& primal_size, const Scalar& dual_size)
92{
93 const Cell cell = myKSpace.unsigns(signed_cell);
94 const DGtal::Dimension cell_dim = myKSpace.uDim(cell);
95 Property property;
96 property.primal_size = primal_size;
97 property.dual_size = dual_size;
98 property.index = std::numeric_limits<Index>::max();
99 property.flipped = ( myKSpace.sSign(signed_cell) == KSpace::NEG );
100
101 if (cell_dim == 0) ASSERT_MSG( property.primal_size == 1, "0-cell primal size must be equal to 1" );
102 if (cell_dim == dimEmbedded) ASSERT_MSG( property.dual_size == 1, "n-cell dual size must be equal to 1" );
103 ASSERT_MSG( cell_dim <= dimEmbedded, "wrong cell dimension" );
104 ASSERT_MSG( cell_dim != 0 || !property.flipped , "can't insert negative 0-cells" );
105 ASSERT_MSG( cell_dim != dimAmbient || !property.flipped , "can't insert negative n-cells" );
106
107 std::pair<typename Properties::iterator, bool> insert_pair = myCellProperties.insert(std::make_pair(cell, property));
108 if (!insert_pair.second) insert_pair.first->second = property;
109
110 ASSERT( insert_pair.first->first == cell );
111 ASSERT( insert_pair.first->second.dual_size == property.dual_size );
112 ASSERT( insert_pair.first->second.primal_size == property.primal_size );
113 ASSERT( insert_pair.first->second.index == property.index );
114 ASSERT( insert_pair.first->second.flipped == property.flipped );
115
116 myIndexesNeedUpdate = true;
117 myCachedOperatorsNeedUpdate = true;
118
119 return insert_pair.second;
120}
121
122template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
123void
124DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::resetSizes()
125{
126 for (typename Properties::iterator pi=myCellProperties.begin(), pe=myCellProperties.end(); pi!=pe; pi++)
127 {
128 pi->second.primal_size = 1;
129 pi->second.dual_size = 1;
130 }
131
132 myCachedOperatorsNeedUpdate = true;
133}
134
135template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
136std::string
137DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::className() const
138{
139 return "Calculus";
140}
141
142template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
143void
144DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::selfDisplay(std::ostream & os) const
145{
146 os << "[dec";
147 for (DGtal::Order order=0; order<=dimEmbedded; order++)
148 os << " | primal " << order << "-cells <-> dual " << dimEmbedded-order << "-cells (" << kFormLength(order, PRIMAL) << ")";
149 os << "]";
150}
151
152template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
153template <DGtal::Order order, DGtal::Duality duality, typename TConstIterator>
154DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order, duality>
155DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::reorder(const TConstIterator& begin_range, const TConstIterator& end_range) const
156{
157 BOOST_STATIC_ASSERT(( boost::is_convertible<typename TConstIterator::value_type, const SCell>::value ));
158
159 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
160
161 typedef typename TLinearAlgebraBackend::Triplet Triplet;
162 typedef std::vector<Triplet> Triplets;
163 Triplets triplets;
164
165 Index original_index = 0;
166 for (TConstIterator ci=begin_range; ci!=end_range; ci++)
167 {
168 const SCell signed_cell = *ci;
169 const Dimension dim = myKSpace.sDim(signed_cell);
170 if (dim != actualOrder(order, duality)) continue;
171 const Cell cell = myKSpace.unsigns(signed_cell);
172 const Index calculus_index = getCellIndex(cell);
173 triplets.push_back( Triplet(original_index, calculus_index, 1) );
174 original_index++;
175 }
176
177 const Index length = kFormLength(order, duality);
178 ASSERT( triplets.size() == static_cast<typename Triplets::size_type>(length) );
179 SparseMatrix reorder_matrix(length, length);
180 reorder_matrix.setFromTriplets(triplets.begin(), triplets.end());
181
182 typedef LinearOperator<Self, order, duality, order, duality> ReorderOperator;
183 return ReorderOperator(*this, reorder_matrix);
184}
185
186template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
187template <DGtal::Order order, DGtal::Duality duality>
188DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order, duality>
189DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::identity() const
190{
191 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
192
193 typedef LinearOperator<Self, order, duality, order, duality> Operator;
194 Operator id(*this);
195 id.myContainer.setIdentity();
196 return id;
197}
198
199template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
200template <DGtal::Duality duality>
201DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 0, duality, 0, duality>
202DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::laplace() const
203{
204 typedef DGtal::LinearOperator<Self, 0, duality, 1, duality> Derivative;
205 typedef DGtal::LinearOperator<Self, 1, duality, 0, duality> Antiderivative;
206 const Derivative d = derivative<0, duality>();
207 const Antiderivative ad = antiderivative<1, duality>();
208 return ad * d;
209}
210
211template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
212template <DGtal::Duality duality>
213DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 0, duality, 0, duality>
214DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::heatLaplace(const typename DenseVector::Scalar& h,
215 const typename DenseVector::Scalar& t, const typename DenseVector::Scalar& K) const
216{
217 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
218
219 typedef typename TLinearAlgebraBackend::Triplet Triplet;
220 typedef LinearOperator<Self, 0, duality, 0, duality> Operator;
221
222 const typename DenseVector::Scalar cut = K * sqrt(2. * t);
223
224 DGtal::CanonicSCellEmbedder<KSpace> canonicSCellEmbedder(myKSpace);
225
226 typedef std::vector<Triplet> Triplets;
227 Triplets triplets;
228
229 for (Index i = 0; i < kFormLength(0, duality); i++)
230 {
231 trace.progressBar( i, kFormLength( 0, duality ) - 1 );
232
233 const SCell signed_cell_i = myIndexSignedCells[ actualOrder(0, duality) ][i];
234 const auto p_i = double( h ) * canonicSCellEmbedder( signed_cell_i );
235
236 typename DenseVector::Scalar total_weight = 0.;
237
238 for(Index j = 0; j < kFormLength(0, duality); j++)
239 {
240 if(i == j) continue;
241
242 const SCell signed_cell_j = myIndexSignedCells[ actualOrder(0, duality) ][j];
243 const auto p_j = double( h ) * canonicSCellEmbedder( signed_cell_j );
244
245 const typename DenseVector::Scalar l2_distance = (p_i - p_j).norm();
246 if(l2_distance < cut)
247 {
248 const typename Properties::const_iterator iter_property = myCellProperties.find( myKSpace.unsigns( signed_cell_j ) );
249 if ( iter_property == myCellProperties.end() )
250 continue;
251
252 const typename DenseVector::Scalar measure = (duality == DUAL) ? iter_property->second.primal_size : iter_property->second.dual_size;
253 const typename DenseVector::Scalar laplace_value = measure * exp(- l2_distance * l2_distance / (4. * t)) * ( 1. / (t * pow(4. * M_PI * t, dimEmbedded / 2.)) );
254
255 triplets.push_back( Triplet(i, j, laplace_value) );
256 total_weight -= laplace_value;
257 }
258 }
259
260 //_operator.myContainer.insert( i, i ) = total_weight;
261 triplets.push_back( Triplet( i, i, total_weight ) );
262 }
263
264 Operator _operator( *this );
265 ASSERT( _operator.myContainer.rows() == kFormLength(0, duality) );
266 ASSERT( _operator.myContainer.cols() == kFormLength(0, duality) );
267 _operator.myContainer.setFromTriplets(triplets.begin(), triplets.end());
268 _operator.myContainer.makeCompressed();
269
270 return _operator;
271}
272
273template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
274template <DGtal::Order order, DGtal::Duality duality>
275DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order-1, duality>
276DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::antiderivative() const
277{
278 BOOST_STATIC_ASSERT(( order > 0 ));
279 BOOST_STATIC_ASSERT(( order <= dimEmbedded ));
280
281 typedef DGtal::LinearOperator<Self, order, duality, dimEmbedded-order, OppositeDuality<duality>::duality> FirstHodge;
282 typedef DGtal::LinearOperator<Self, dimEmbedded-order, OppositeDuality<duality>::duality, dimEmbedded-order+1, OppositeDuality<duality>::duality> Derivative;
283 typedef DGtal::LinearOperator<Self, dimEmbedded-order+1, OppositeDuality<duality>::duality, order-1, duality> SecondHodge;
284 const FirstHodge h_first = hodge<order, duality>();
285 const Derivative d = derivative<dimEmbedded-order, OppositeDuality<duality>::duality>();
286 const SecondHodge h_second = hodge<dimEmbedded-order+1, OppositeDuality<duality>::duality>();
287 const Scalar sign = ( order*(dimEmbedded-order)%2 == 0 ? 1 : -1 );
288 return sign * h_second * d * h_first;
289}
290
291template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
292template <DGtal::Order order, DGtal::Duality duality>
293DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order+1, duality>
294DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::derivative() const
295{
296 BOOST_STATIC_ASSERT(( order >= 0 ));
297 BOOST_STATIC_ASSERT(( order < dimEmbedded ));
298
299 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
300
301 typedef typename TLinearAlgebraBackend::Triplet Triplet;
302 typedef std::vector<Triplet> Triplets;
303 Triplets triplets;
304
305 // iterate over output form values
306 for (Index index_output=0; index_output<kFormLength(order+1, duality); index_output++)
307 {
308 const SCell signed_cell = myIndexSignedCells[actualOrder(order+1, duality)][index_output];
309
310 // find cell border
311 typedef typename KSpace::SCells Border;
312 const Border border = ( duality == PRIMAL ? myKSpace.sLowerIncident(signed_cell) : myKSpace.sUpperIncident(signed_cell) );
313
314 // iterate over cell border
315 for (typename Border::const_iterator bi=border.begin(), bie=border.end(); bi!=bie; bi++)
316 {
317 const SCell signed_cell_border = *bi;
318 ASSERT( myKSpace.sDim(signed_cell_border) == actualOrder(order, duality) );
319
320 const typename Properties::const_iterator iter_property = myCellProperties.find(myKSpace.unsigns(signed_cell_border));
321 if ( iter_property == myCellProperties.end() )
322 continue;
323
324 const Index index_input = iter_property->second.index;
325 ASSERT( index_input < kFormLength(order, duality) );
326
327 const bool flipped_border = ( myKSpace.sSign(signed_cell_border) == KSpace::NEG );
328 const Scalar orientation = ( flipped_border == iter_property->second.flipped ? 1 : -1 );
329
330 triplets.push_back( Triplet(index_output, index_input, orientation) );
331
332 }
333 }
334
335 typedef LinearOperator<Self, order, duality, order+1, duality> Derivative;
336 Derivative _derivative(*this);
337 ASSERT( _derivative.myContainer.rows() == kFormLength(order+1, duality) );
338 ASSERT( _derivative.myContainer.cols() == kFormLength(order, duality) );
339 _derivative.myContainer.setFromTriplets(triplets.begin(), triplets.end());
340
341 if ( duality == DUAL && order*(dimEmbedded-order)%2 != 0 ) return -1 * _derivative;
342 return _derivative;
343}
344
345template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
346template <DGtal::Order order, DGtal::Duality duality>
347DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, dimEmbedded-order, DGtal::OppositeDuality<duality>::duality>
348DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::hodge() const
349{
350 BOOST_STATIC_ASSERT(( order >= 0 ));
351 BOOST_STATIC_ASSERT(( order <= dimEmbedded ));
352
353 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
354
355 typedef typename TLinearAlgebraBackend::Triplet Triplet;
356 typedef std::vector<Triplet> Triplets;
357 Triplets triplets;
358
359 // iterate over output form values
360 for (Index index=0; index<kFormLength(order, duality); index++)
361 {
362 const Cell cell = myKSpace.unsigns(myIndexSignedCells[actualOrder(order, duality)][index]);
363
364 const typename Properties::const_iterator iter_property = myCellProperties.find(cell);
365 ASSERT( iter_property != myCellProperties.end() );
366 ASSERT( iter_property->second.index == index );
367
368 const Scalar size_ratio = ( duality == DGtal::PRIMAL ?
369 iter_property->second.dual_size/iter_property->second.primal_size :
370 iter_property->second.primal_size/iter_property->second.dual_size );
371 triplets.push_back( Triplet(index, index, hodgeSign(cell, duality) * size_ratio) );
372 }
373
374 typedef LinearOperator<Self, order, duality, dimEmbedded-order, OppositeDuality<duality>::duality> Hodge;
375 Hodge _hodge(*this);
376 ASSERT( _hodge.myContainer.rows() == _hodge.myContainer.cols() );
377 ASSERT( _hodge.myContainer.rows() == kFormLength(order, duality) );
378 _hodge.myContainer.setFromTriplets(triplets.begin(), triplets.end());
379
380 return _hodge;
381}
382
383template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
384template <DGtal::Duality duality>
385DGtal::VectorField<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, duality>
386DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::sharp(const DGtal::KForm<Self, 1, duality>& one_form) const
387{
388 ASSERT( one_form.myCalculus == this );
389
390 const_cast<Self*>(this)->updateCachedOperators();
391 ASSERT( !myCachedOperatorsNeedUpdate );
392
393 const boost::array<SparseMatrix, dimAmbient>& sharp_operator_matrix = mySharpOperatorMatrixes[static_cast<int>(duality)];
394
395 typedef VectorField<Self, duality> Field;
396 Field field(*this);
397 ASSERT( field.myCoordinates.cols() == dimAmbient);
398 ASSERT( field.myCoordinates.rows() == kFormLength(0, duality));
399 for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
400 field.myCoordinates.col(direction) = sharp_operator_matrix[direction]*one_form.myContainer;
401
402 return field;
403}
404
405template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
406template <DGtal::Duality duality>
407DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 1, duality, 0, duality>
408DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::sharpDirectional(const DGtal::Dimension& direction) const
409{
410 const_cast<Self*>(this)->updateCachedOperators();
411 ASSERT( !myCachedOperatorsNeedUpdate );
412
413 ASSERT( direction < dimAmbient );
414
415 typedef LinearOperator<Self, 1, duality, 0, duality> Operator;
416 return Operator(*this, mySharpOperatorMatrixes[static_cast<int>(duality)][direction]);
417}
418
419
420template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
421template <DGtal::Duality duality>
422void
423DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::updateSharpOperator()
424{
425 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
426 ASSERT( myCachedOperatorsNeedUpdate );
427
428 typedef typename TLinearAlgebraBackend::Triplet Triplet;
429 typedef std::vector<Triplet> Triplets;
430 typedef typename Properties::const_iterator PropertiesConstIterator;
431
432 boost::array<Triplets, dimAmbient> triplets;
433
434 // iterate over points
435 for (Index point_index=0; point_index<kFormLength(0, duality); point_index++)
436 {
437 const SCell signed_point = myIndexSignedCells[actualOrder(0, duality)][point_index];
438 ASSERT( myKSpace.sDim(signed_point) == actualOrder(0, duality) );
439 const Scalar point_orientation = ( myKSpace.sSign(signed_point) == KSpace::POS ? 1 : -1 );
440 const Cell point = myKSpace.unsigns(signed_point);
441
442 typedef typename KSpace::Cells Edges;
443 typedef typename Edges::const_iterator EdgesConstIterator;
444 const Edges edges = ( duality == PRIMAL ? myKSpace.uUpperIncident(point) : myKSpace.uLowerIncident(point) );
445 ASSERT( edges.size() <= 2*dimAmbient );
446
447 // collect 1-form values over neighboring edges
448 typedef boost::array< std::pair<Scalar, std::list< std::pair<Index, Scalar> > >, dimAmbient > EdgeIndexes;
449 EdgeIndexes edge_indexes;
450 for (EdgesConstIterator ei=edges.begin(), eie=edges.end(); ei!=eie; ei++)
451 {
452 const Cell edge = *ei;
453 ASSERT( myKSpace.uDim(edge) == actualOrder(1, duality) );
454
455 const PropertiesConstIterator edge_property_iter = myCellProperties.find(edge);
456 if (edge_property_iter == myCellProperties.end())
457 continue;
458
459 const Index edge_index = edge_property_iter->second.index;
460 const Scalar edge_length = ( duality == PRIMAL ? edge_property_iter->second.primal_size : edge_property_iter->second.dual_size );
461
462 const Scalar edge_orientation = ( edge_property_iter->second.flipped ? 1 : -1 );
463 const DGtal::Dimension edge_direction = edgeDirection(edge, duality); //FIXME iterate over direction
464
465 edge_indexes[edge_direction].second.push_back(std::make_pair(edge_index, edge_orientation));
466 edge_indexes[edge_direction].first += edge_length;
467 }
468
469 for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
470 {
471 const Scalar edge_sign = ( duality == DUAL && (direction*(dimAmbient-direction))%2 == 0 ? -1 : 1 );
472 const Scalar edge_length_sum = edge_indexes[direction].first;
473
474 for (typename std::list< std::pair<Index, Scalar> >::const_iterator ei=edge_indexes[direction].second.begin(), ee=edge_indexes[direction].second.end(); ei!=ee; ei++)
475 {
476 const Index edge_index = ei->first;
477 const Scalar edge_orientation = ei->second;
478 ASSERT( edge_index < static_cast<Index>(myIndexSignedCells[actualOrder(1, duality)].size()) );
479 ASSERT( edge_index < kFormLength(1, duality) );
480 ASSERT( edge_length_sum > 0 );
481
482 triplets[direction].push_back( Triplet(point_index, edge_index, point_orientation*edge_sign*edge_orientation/edge_length_sum) );
483 }
484 }
485 }
486
487 boost::array<SparseMatrix, dimAmbient> sharp_operator_matrix;
488
489 for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
490 {
491 sharp_operator_matrix[direction] = SparseMatrix(kFormLength(0, duality), kFormLength(1, duality));
492 sharp_operator_matrix[direction].setFromTriplets(triplets[direction].begin(), triplets[direction].end());
493 }
494
495 mySharpOperatorMatrixes[static_cast<int>(duality)] = sharp_operator_matrix;
496}
497
498template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
499template <DGtal::Duality duality>
500DGtal::KForm<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 1, duality>
501DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::flat(const DGtal::VectorField<Self, duality>& vector_field) const
502{
503 ASSERT( vector_field.myCalculus == this );
504
505 const_cast<Self*>(this)->updateCachedOperators();
506 ASSERT( !myCachedOperatorsNeedUpdate );
507
508 const boost::array<SparseMatrix, dimAmbient>& flat_operator_matrix = myFlatOperatorMatrixes[static_cast<int>(duality)];
509
510 typedef KForm<Self, 1, duality> OneForm;
511 OneForm one_form(*this);
512 for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
513 one_form.myContainer += flat_operator_matrix[direction]*vector_field.myCoordinates.col(direction);
514
515 return one_form;
516}
517
518template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
519template <DGtal::Duality duality>
520DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 0, duality, 1, duality>
521DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::flatDirectional(const DGtal::Dimension& direction) const
522{
523 const_cast<Self*>(this)->updateCachedOperators();
524 ASSERT( !myCachedOperatorsNeedUpdate );
525
526 ASSERT( direction < dimAmbient );
527
528 typedef LinearOperator<Self, 0, duality, 1, duality> Operator;
529 return Operator(*this, myFlatOperatorMatrixes[static_cast<int>(duality)][direction]);
530}
531
532template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
533template <DGtal::Duality duality>
534void
535DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::updateFlatOperator()
536{
537 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
538 ASSERT( myCachedOperatorsNeedUpdate );
539
540 typedef typename TLinearAlgebraBackend::Triplet Triplet;
541 typedef std::vector<Triplet> Triplets;
542 typedef typename Properties::const_iterator PropertiesConstIterator;
543
544 boost::array<Triplets, dimAmbient> triplets;
545
546 // iterate over edges
547 for (Index edge_index=0; edge_index<kFormLength(1, duality); edge_index++)
548 {
549 const SCell signed_edge = myIndexSignedCells[actualOrder(1, duality)][edge_index];
550 ASSERT( myKSpace.sDim(signed_edge) == actualOrder(1, duality) );
551 const Cell edge = myKSpace.unsigns(signed_edge);
552
553 const Scalar edge_orientation = ( myKSpace.sSign(signed_edge) == KSpace::NEG ? 1 : -1 );
554 const DGtal::Dimension& edge_direction = edgeDirection(edge, duality); //FIXME iterate over edge direction
555 const Scalar edge_sign = ( duality == DUAL && (edge_direction*(dimAmbient-edge_direction))%2 == 0 ? -1 : 1 );
556 const PropertiesConstIterator edge_property_iter = myCellProperties.find(edge);
557 ASSERT( edge_property_iter != myCellProperties.end() );
558 const Scalar edge_length = ( duality == PRIMAL ? edge_property_iter->second.primal_size : edge_property_iter->second.dual_size );
559
560 typedef typename KSpace::Cells Points;
561 const Points points = ( duality == PRIMAL ? myKSpace.uLowerIncident(edge) : myKSpace.uUpperIncident(edge) );
562
563 // project vector field along edge from neighboring points
564 typedef std::pair<Index, Scalar> BorderInfo;
565 typedef std::list<BorderInfo> BorderInfos;
566 BorderInfos border_infos;
567 for (typename Points::const_iterator pi=points.begin(), pie=points.end(); pi!=pie; pi++)
568 {
569 const Cell point = *pi;
570 ASSERT( myKSpace.uDim(point) == actualOrder(0, duality) );
571
572 const PropertiesConstIterator point_property_iter = myCellProperties.find(point);
573 if (point_property_iter == myCellProperties.end())
574 continue;
575
576 const Index point_index = point_property_iter->second.index;
577 const Scalar point_orientation = ( point_property_iter->second.flipped ? -1 : 1 );
578
579 border_infos.push_back(std::make_pair(point_index, point_orientation));
580 }
581
582 ASSERT( border_infos.size() <= 2 );
583
584 for (typename BorderInfos::const_iterator bi=border_infos.begin(), bie=border_infos.end(); bi!=bie; bi++)
585 {
586 const Index point_index = bi->first;
587 const Scalar point_orientation = bi->second;
588 ASSERT( point_index < static_cast<Index>(myIndexSignedCells[actualOrder(0, duality)].size()) );
589 ASSERT( point_index < kFormLength(0, duality) );
590
591 triplets[edge_direction].push_back( Triplet(edge_index, point_index, point_orientation*edge_length*edge_sign*edge_orientation/border_infos.size()) );
592 }
593
594 }
595
596 boost::array<SparseMatrix, dimAmbient> flat_operator_matrix;
597
598 for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
599 {
600 flat_operator_matrix[direction] = SparseMatrix(kFormLength(1, duality), kFormLength(0, duality));
601 flat_operator_matrix[direction].setFromTriplets(triplets[direction].begin(), triplets[direction].end());
602 }
603
604 myFlatOperatorMatrixes[static_cast<int>(duality)] = flat_operator_matrix;
605}
606
607template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
608void
609DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::updateIndexes()
610{
611 if (!myIndexesNeedUpdate) return;
612
613 // clear index signed cells
614 for (DGtal::Dimension dim=0; dim<dimEmbedded+1; dim++)
615 myIndexSignedCells[dim].clear();
616
617 // compute cell index
618 for (typename Properties::iterator csi=myCellProperties.begin(), csie=myCellProperties.end(); csie!=csi; csi++)
619 {
620 const Cell& cell = csi->first;
621 const DGtal::Dimension cell_dim = myKSpace.uDim(cell);
622
623 csi->second.index = myIndexSignedCells[cell_dim].size();
624
625 const SCell& signed_cell = myKSpace.signs(cell, csi->second.flipped ? KSpace::NEG : KSpace::POS);
626 myIndexSignedCells[cell_dim].push_back(signed_cell);
627 }
628
629 myIndexesNeedUpdate = false;
630 myCachedOperatorsNeedUpdate = true;
631}
632
633template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
634void
635DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::updateCachedOperators()
636{
637 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
638 if (!myCachedOperatorsNeedUpdate) return;
639 updateFlatOperator<PRIMAL>();
640 updateFlatOperator<DUAL>();
641 updateSharpOperator<PRIMAL>();
642 updateSharpOperator<DUAL>();
643 myCachedOperatorsNeedUpdate = false;
644}
645
646template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
647const typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Properties&
648DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getProperties() const
649{
650 return myCellProperties;
651}
652
653template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
654template <DGtal::Order order, DGtal::Duality duality>
655const typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::SCells&
656DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getIndexedSCells() const
657{
658 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
659 return myIndexSignedCells[actualOrder(order, duality)];
660}
661
662template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
663typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::SCell
664DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getSCell(const Order& order, const Duality& duality, const Index& index) const
665{
666 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
667 const Order& actual_order = actualOrder(order, duality);
668 const SCell& signed_cell = myIndexSignedCells[actual_order][index];
669 ASSERT( myKSpace.sDim(signed_cell) == actual_order );
670 return signed_cell;
671}
672
673template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
674bool
675DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::containsCell(const Cell& cell) const
676{
677 return myCellProperties.find(cell) != myCellProperties.end();
678}
679
680template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
681bool
682DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::isCellFlipped(const Cell& cell) const
683{
684 const typename Properties::const_iterator iter_property = myCellProperties.find(cell);
685 ASSERT( iter_property != myCellProperties.end() );
686 return iter_property->second.flipped;
687}
688
689template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
690typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Index
691DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getCellIndex(const Cell& cell) const
692{
693 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
694 const typename Properties::const_iterator iter_property = myCellProperties.find(cell);
695 ASSERT( iter_property != myCellProperties.end() );
696 return iter_property->second.index;
697}
698
699template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
700typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::ConstIterator
701DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::begin() const
702{
703 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
704 return myCellProperties.begin();
705}
706
707template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
708typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::ConstIterator
709DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::end() const
710{
711 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
712 return myCellProperties.end();
713}
714
715template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
716typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Iterator
717DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::begin()
718{
719 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
720 return myCellProperties.begin();
721}
722
723template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
724typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Iterator
725DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::end()
726{
727 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
728 return myCellProperties.end();
729}
730
731template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
732typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Index
733DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::kFormLength(const DGtal::Order& order, const DGtal::Duality& duality) const
734{
735 ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
736 return myIndexSignedCells[actualOrder(order, duality)].size();
737}
738
739template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
740DGtal::Order
741DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::actualOrder(const DGtal::Order& order, const DGtal::Duality& duality) const
742{
743 return duality == PRIMAL ? order : dimEmbedded-order;
744}
745
746template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
747typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Scalar
748DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::hodgeSign(const Cell& cell, const DGtal::Duality& duality) const
749{
750 if (duality == PRIMAL) return 1;
751
752 const DGtal::Dimension& primal_dim = myKSpace.uDim(cell);
753 return (dimEmbedded-primal_dim)*primal_dim % 2 != 0 ? -1 : 1;
754}
755
756template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
757typename DGtal::Dimension
758DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::edgeDirection(const Cell& cell, const DGtal::Duality& duality) const
759{
760 ASSERT( myKSpace.uDim(cell) == actualOrder(1, duality) );
761
762 typename KSpace::DirIterator di = myKSpace.uDirs(cell);
763 if (duality == DUAL) di = myKSpace.uOrthDirs(cell);
764
765 ASSERT( di != 0 );
766 const DGtal::Dimension direction = *di;
767 ++di;
768 //ASSERT( !(di != 0) ); //FIXME multiple edge direction with embedding
769 ASSERT( direction < dimAmbient );
770
771 return direction;
772}
773
774template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
775bool
776DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::isValid() const
777{
778 return true;
779}
780
781///////////////////////////////////////////////////////////////////////////////
782// Implementation of inline functions //
783
784template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
785std::ostream&
786DGtal::operator<<(std::ostream & out, const DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>& object)
787{
788 object.selfDisplay(out);
789 return out;
790}
791
792// //
793///////////////////////////////////////////////////////////////////////////////