DGtal  1.5.beta
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 
33 template <DGtal::Dimension dim, typename TInteger>
34 size_t
35 DGtal::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 
43 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
44 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::DiscreteExteriorCalculus()
45  : myKSpace(), myCachedOperatorsNeedUpdate(true), myIndexesNeedUpdate(false)
46 {
47 }
48 
49 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
50 template <typename TDomain>
51 void
52 DGtal::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 
66 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
67 bool
68 DGtal::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 
82 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
83 bool
84 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::insertSCell(const SCell& signed_cell)
85 {
86  return insertSCell(signed_cell, 1, 1);
87 }
88 
89 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
90 bool
91 DGtal::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 
122 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
123 void
124 DGtal::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 
135 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
136 std::string
137 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::className() const
138 {
139  return "Calculus";
140 }
141 
142 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
143 void
144 DGtal::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 
152 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
153 template <DGtal::Order order, DGtal::Duality duality, typename TConstIterator>
154 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order, duality>
155 DGtal::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 
186 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
187 template <DGtal::Order order, DGtal::Duality duality>
188 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order, duality>
189 DGtal::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 
199 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
200 template <DGtal::Duality duality>
201 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 0, duality, 0, duality>
202 DGtal::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 
211 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
212 template <DGtal::Duality duality>
213 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 0, duality, 0, duality>
214 DGtal::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 
273 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
274 template <DGtal::Order order, DGtal::Duality duality>
275 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order-1, duality>
276 DGtal::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 
291 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
292 template <DGtal::Order order, DGtal::Duality duality>
293 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order+1, duality>
294 DGtal::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 
345 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
346 template <DGtal::Order order, DGtal::Duality duality>
347 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, dimEmbedded-order, DGtal::OppositeDuality<duality>::duality>
348 DGtal::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 
383 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
384 template <DGtal::Duality duality>
385 DGtal::VectorField<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, duality>
386 DGtal::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 
405 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
406 template <DGtal::Duality duality>
407 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 1, duality, 0, duality>
408 DGtal::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 
420 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
421 template <DGtal::Duality duality>
422 void
423 DGtal::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 
498 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
499 template <DGtal::Duality duality>
500 DGtal::KForm<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 1, duality>
501 DGtal::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 
518 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
519 template <DGtal::Duality duality>
520 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 0, duality, 1, duality>
521 DGtal::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 
532 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
533 template <DGtal::Duality duality>
534 void
535 DGtal::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 
607 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
608 void
609 DGtal::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 
633 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
634 void
635 DGtal::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 
646 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
647 const typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Properties&
648 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getProperties() const
649 {
650  return myCellProperties;
651 }
652 
653 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
654 template <DGtal::Order order, DGtal::Duality duality>
655 const typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::SCells&
656 DGtal::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 
662 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
663 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::SCell
664 DGtal::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 
673 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
674 bool
675 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::containsCell(const Cell& cell) const
676 {
677  return myCellProperties.find(cell) != myCellProperties.end();
678 }
679 
680 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
681 bool
682 DGtal::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 
689 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
690 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Index
691 DGtal::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 
699 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
700 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::ConstIterator
701 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::begin() const
702 {
703  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
704  return myCellProperties.begin();
705 }
706 
707 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
708 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::ConstIterator
709 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::end() const
710 {
711  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
712  return myCellProperties.end();
713 }
714 
715 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
716 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Iterator
717 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::begin()
718 {
719  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
720  return myCellProperties.begin();
721 }
722 
723 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
724 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Iterator
725 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::end()
726 {
727  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
728  return myCellProperties.end();
729 }
730 
731 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
732 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Index
733 DGtal::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 
739 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
740 DGtal::Order
741 DGtal::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 
746 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
747 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Scalar
748 DGtal::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 
756 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
757 typename DGtal::Dimension
758 DGtal::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 
774 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
775 bool
776 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::isValid() const
777 {
778  return true;
779 }
780 
781 ///////////////////////////////////////////////////////////////////////////////
782 // Implementation of inline functions //
783 
784 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
785 std::ostream&
786 DGtal::operator<<(std::ostream & out, const DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>& object)
787 {
788  object.selfDisplay(out);
789  return out;
790 }
791 
792 // //
793 ///////////////////////////////////////////////////////////////////////////////