DGtal  0.9.2
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 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 >= 0 && 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::Order order, DGtal::Duality duality>
213 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order-1, duality>
214 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::antiderivative() const
215 {
216  BOOST_STATIC_ASSERT(( order > 0 ));
217  BOOST_STATIC_ASSERT(( order <= dimEmbedded ));
218 
219  typedef DGtal::LinearOperator<Self, order, duality, dimEmbedded-order, OppositeDuality<duality>::duality> FirstHodge;
220  typedef DGtal::LinearOperator<Self, dimEmbedded-order, OppositeDuality<duality>::duality, dimEmbedded-order+1, OppositeDuality<duality>::duality> Derivative;
221  typedef DGtal::LinearOperator<Self, dimEmbedded-order+1, OppositeDuality<duality>::duality, order-1, duality> SecondHodge;
222  const FirstHodge h_first = hodge<order, duality>();
223  const Derivative d = derivative<dimEmbedded-order, OppositeDuality<duality>::duality>();
224  const SecondHodge h_second = hodge<dimEmbedded-order+1, OppositeDuality<duality>::duality>();
225  const Scalar sign = ( order*(dimEmbedded-order)%2 == 0 ? 1 : -1 );
226  return sign * h_second * d * h_first;
227 }
228 
229 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
230 template <DGtal::Order order, DGtal::Duality duality>
231 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, order+1, duality>
232 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::derivative() const
233 {
234  BOOST_STATIC_ASSERT(( order >= 0 ));
235  BOOST_STATIC_ASSERT(( order < dimEmbedded ));
236 
237  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
238 
239  typedef typename TLinearAlgebraBackend::Triplet Triplet;
240  typedef std::vector<Triplet> Triplets;
241  Triplets triplets;
242 
243  // iterate over output form values
244  for (Index index_output=0; index_output<kFormLength(order+1, duality); index_output++)
245  {
246  const SCell signed_cell = myIndexSignedCells[actualOrder(order+1, duality)][index_output];
247 
248  // find cell border
249  typedef typename KSpace::SCells Border;
250  const Border border = ( duality == PRIMAL ? myKSpace.sLowerIncident(signed_cell) : myKSpace.sUpperIncident(signed_cell) );
251 
252  // iterate over cell border
253  for (typename Border::const_iterator bi=border.begin(), bie=border.end(); bi!=bie; bi++)
254  {
255  const SCell signed_cell_border = *bi;
256  ASSERT( myKSpace.sDim(signed_cell_border) == actualOrder(order, duality) );
257 
258  const typename Properties::const_iterator iter_property = myCellProperties.find(myKSpace.unsigns(signed_cell_border));
259  if ( iter_property == myCellProperties.end() )
260  continue;
261 
262  const Index index_input = iter_property->second.index;
263  ASSERT( index_input < kFormLength(order, duality) );
264 
265  const bool flipped_border = ( myKSpace.sSign(signed_cell_border) == KSpace::NEG );
266  const Scalar orientation = ( flipped_border == iter_property->second.flipped ? 1 : -1 );
267 
268  triplets.push_back( Triplet(index_output, index_input, orientation) );
269 
270  }
271  }
272 
273  typedef LinearOperator<Self, order, duality, order+1, duality> Derivative;
274  Derivative _derivative(*this);
275  ASSERT( _derivative.myContainer.rows() == kFormLength(order+1, duality) );
276  ASSERT( _derivative.myContainer.cols() == kFormLength(order, duality) );
277  _derivative.myContainer.setFromTriplets(triplets.begin(), triplets.end());
278 
279  if ( duality == DUAL && order*(dimEmbedded-order)%2 != 0 ) return -1 * _derivative;
280  return _derivative;
281 }
282 
283 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
284 template <DGtal::Order order, DGtal::Duality duality>
285 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, order, duality, dimEmbedded-order, DGtal::OppositeDuality<duality>::duality>
286 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::hodge() const
287 {
288  BOOST_STATIC_ASSERT(( order >= 0 ));
289  BOOST_STATIC_ASSERT(( order <= dimEmbedded ));
290 
291  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
292 
293  typedef typename TLinearAlgebraBackend::Triplet Triplet;
294  typedef std::vector<Triplet> Triplets;
295  Triplets triplets;
296 
297  // iterate over output form values
298  for (Index index=0; index<kFormLength(order, duality); index++)
299  {
300  const Cell cell = myKSpace.unsigns(myIndexSignedCells[actualOrder(order, duality)][index]);
301 
302  const typename Properties::const_iterator iter_property = myCellProperties.find(cell);
303  ASSERT( iter_property != myCellProperties.end() );
304  ASSERT( iter_property->second.index == index );
305 
306  const Scalar size_ratio = ( duality == DGtal::PRIMAL ?
307  iter_property->second.dual_size/iter_property->second.primal_size :
308  iter_property->second.primal_size/iter_property->second.dual_size );
309  triplets.push_back( Triplet(index, index, hodgeSign(cell, duality) * size_ratio) );
310  }
311 
312  typedef LinearOperator<Self, order, duality, dimEmbedded-order, OppositeDuality<duality>::duality> Hodge;
313  Hodge _hodge(*this);
314  ASSERT( _hodge.myContainer.rows() == _hodge.myContainer.cols() );
315  ASSERT( _hodge.myContainer.rows() == kFormLength(order, duality) );
316  _hodge.myContainer.setFromTriplets(triplets.begin(), triplets.end());
317 
318  return _hodge;
319 }
320 
321 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
322 template <DGtal::Duality duality>
323 DGtal::VectorField<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, duality>
324 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::sharp(const DGtal::KForm<Self, 1, duality>& one_form) const
325 {
326  ASSERT( one_form.myCalculus == this );
327 
328  const_cast<Self*>(this)->updateCachedOperators();
329  ASSERT( !myCachedOperatorsNeedUpdate );
330 
331  const boost::array<SparseMatrix, dimAmbient>& sharp_operator_matrix = mySharpOperatorMatrixes[static_cast<int>(duality)];
332 
333  typedef VectorField<Self, duality> Field;
334  Field field(*this);
335  ASSERT( field.myCoordinates.cols() == dimAmbient);
336  ASSERT( field.myCoordinates.rows() == kFormLength(0, duality));
337  for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
338  field.myCoordinates.col(direction) = sharp_operator_matrix[direction]*one_form.myContainer;
339 
340  return field;
341 }
342 
343 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
344 template <DGtal::Duality duality>
345 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 1, duality, 0, duality>
346 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::sharpDirectional(const DGtal::Dimension& direction) const
347 {
348  const_cast<Self*>(this)->updateCachedOperators();
349  ASSERT( !myCachedOperatorsNeedUpdate );
350 
351  ASSERT( direction < dimAmbient );
352 
353  typedef LinearOperator<Self, 1, duality, 0, duality> Operator;
354  return Operator(*this, mySharpOperatorMatrixes[static_cast<int>(duality)][direction]);
355 }
356 
357 
358 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
359 template <DGtal::Duality duality>
360 void
361 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::updateSharpOperator()
362 {
363  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
364  ASSERT( myCachedOperatorsNeedUpdate );
365 
366  typedef typename TLinearAlgebraBackend::Triplet Triplet;
367  typedef std::vector<Triplet> Triplets;
368  typedef typename Properties::const_iterator PropertiesConstIterator;
369 
370  boost::array<Triplets, dimAmbient> triplets;
371 
372  // iterate over points
373  for (Index point_index=0; point_index<kFormLength(0, duality); point_index++)
374  {
375  const SCell signed_point = myIndexSignedCells[actualOrder(0, duality)][point_index];
376  ASSERT( myKSpace.sDim(signed_point) == actualOrder(0, duality) );
377  const Scalar point_orientation = ( myKSpace.sSign(signed_point) == KSpace::POS ? 1 : -1 );
378  const Cell point = myKSpace.unsigns(signed_point);
379 
380  typedef typename KSpace::Cells Edges;
381  typedef typename Edges::const_iterator EdgesConstIterator;
382  const Edges edges = ( duality == PRIMAL ? myKSpace.uUpperIncident(point) : myKSpace.uLowerIncident(point) );
383  ASSERT( edges.size() <= 2*dimAmbient );
384 
385  // collect 1-form values over neighboring edges
386  typedef boost::array< std::pair<Scalar, std::list< std::pair<Index, Scalar> > >, dimAmbient > EdgeIndexes;
387  EdgeIndexes edge_indexes;
388  for (EdgesConstIterator ei=edges.begin(), eie=edges.end(); ei!=eie; ei++)
389  {
390  const Cell edge = *ei;
391  ASSERT( myKSpace.uDim(edge) == actualOrder(1, duality) );
392 
393  const PropertiesConstIterator edge_property_iter = myCellProperties.find(edge);
394  if (edge_property_iter == myCellProperties.end())
395  continue;
396 
397  const Index edge_index = edge_property_iter->second.index;
398  const Scalar edge_length = ( duality == PRIMAL ? edge_property_iter->second.primal_size : edge_property_iter->second.dual_size );
399 
400  const Scalar edge_orientation = ( edge_property_iter->second.flipped ? 1 : -1 );
401  const DGtal::Dimension edge_direction = edgeDirection(edge, duality); //FIXME iterate over direction
402 
403  edge_indexes[edge_direction].second.push_back(std::make_pair(edge_index, edge_orientation));
404  edge_indexes[edge_direction].first += edge_length;
405  }
406 
407  for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
408  {
409  const Scalar edge_sign = ( duality == DUAL && (direction*(dimAmbient-direction))%2 == 0 ? -1 : 1 );
410  const Scalar edge_length_sum = edge_indexes[direction].first;
411 
412  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++)
413  {
414  const Index edge_index = ei->first;
415  const Scalar edge_orientation = ei->second;
416  ASSERT( edge_index < static_cast<Index>(myIndexSignedCells[actualOrder(1, duality)].size()) );
417  ASSERT( edge_index < kFormLength(1, duality) );
418  ASSERT( edge_length_sum > 0 );
419 
420  triplets[direction].push_back( Triplet(point_index, edge_index, point_orientation*edge_sign*edge_orientation/edge_length_sum) );
421  }
422  }
423  }
424 
425  boost::array<SparseMatrix, dimAmbient> sharp_operator_matrix;
426 
427  for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
428  {
429  sharp_operator_matrix[direction] = SparseMatrix(kFormLength(0, duality), kFormLength(1, duality));
430  sharp_operator_matrix[direction].setFromTriplets(triplets[direction].begin(), triplets[direction].end());
431  }
432 
433  mySharpOperatorMatrixes[static_cast<int>(duality)] = sharp_operator_matrix;
434 }
435 
436 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
437 template <DGtal::Duality duality>
438 DGtal::KForm<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 1, duality>
439 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::flat(const DGtal::VectorField<Self, duality>& vector_field) const
440 {
441  ASSERT( vector_field.myCalculus == this );
442 
443  const_cast<Self*>(this)->updateCachedOperators();
444  ASSERT( !myCachedOperatorsNeedUpdate );
445 
446  const boost::array<SparseMatrix, dimAmbient>& flat_operator_matrix = myFlatOperatorMatrixes[static_cast<int>(duality)];
447 
448  typedef KForm<Self, 1, duality> OneForm;
449  OneForm one_form(*this);
450  for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
451  one_form.myContainer += flat_operator_matrix[direction]*vector_field.myCoordinates.col(direction);
452 
453  return one_form;
454 }
455 
456 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
457 template <DGtal::Duality duality>
458 DGtal::LinearOperator<DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>, 0, duality, 1, duality>
459 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::flatDirectional(const DGtal::Dimension& direction) const
460 {
461  const_cast<Self*>(this)->updateCachedOperators();
462  ASSERT( !myCachedOperatorsNeedUpdate );
463 
464  ASSERT( direction < dimAmbient );
465 
466  typedef LinearOperator<Self, 0, duality, 1, duality> Operator;
467  return Operator(*this, myFlatOperatorMatrixes[static_cast<int>(duality)][direction]);
468 }
469 
470 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
471 template <DGtal::Duality duality>
472 void
473 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::updateFlatOperator()
474 {
475  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
476  ASSERT( myCachedOperatorsNeedUpdate );
477 
478  typedef typename TLinearAlgebraBackend::Triplet Triplet;
479  typedef std::vector<Triplet> Triplets;
480  typedef typename Properties::const_iterator PropertiesConstIterator;
481 
482  boost::array<Triplets, dimAmbient> triplets;
483 
484  // iterate over edges
485  for (Index edge_index=0; edge_index<kFormLength(1, duality); edge_index++)
486  {
487  const SCell signed_edge = myIndexSignedCells[actualOrder(1, duality)][edge_index];
488  ASSERT( myKSpace.sDim(signed_edge) == actualOrder(1, duality) );
489  const Cell edge = myKSpace.unsigns(signed_edge);
490 
491  const Scalar edge_orientation = ( myKSpace.sSign(signed_edge) == KSpace::NEG ? 1 : -1 );
492  const DGtal::Dimension& edge_direction = edgeDirection(edge, duality); //FIXME iterate over edge direction
493  const Scalar edge_sign = ( duality == DUAL && (edge_direction*(dimAmbient-edge_direction))%2 == 0 ? -1 : 1 );
494  const PropertiesConstIterator edge_property_iter = myCellProperties.find(edge);
495  ASSERT( edge_property_iter != myCellProperties.end() );
496  const Scalar edge_length = ( duality == PRIMAL ? edge_property_iter->second.primal_size : edge_property_iter->second.dual_size );
497 
498  typedef typename KSpace::Cells Points;
499  const Points points = ( duality == PRIMAL ? myKSpace.uLowerIncident(edge) : myKSpace.uUpperIncident(edge) );
500 
501  // project vector field along edge from neighboring points
502  typedef std::pair<Index, Scalar> BorderInfo;
503  typedef std::list<BorderInfo> BorderInfos;
504  BorderInfos border_infos;
505  for (typename Points::const_iterator pi=points.begin(), pie=points.end(); pi!=pie; pi++)
506  {
507  const Cell point = *pi;
508  ASSERT( myKSpace.uDim(point) == actualOrder(0, duality) );
509 
510  const PropertiesConstIterator point_property_iter = myCellProperties.find(point);
511  if (point_property_iter == myCellProperties.end())
512  continue;
513 
514  const Index point_index = point_property_iter->second.index;
515  const Scalar point_orientation = ( point_property_iter->second.flipped ? -1 : 1 );
516 
517  border_infos.push_back(std::make_pair(point_index, point_orientation));
518  }
519 
520  ASSERT( border_infos.size() <= 2 );
521 
522  for (typename BorderInfos::const_iterator bi=border_infos.begin(), bie=border_infos.end(); bi!=bie; bi++)
523  {
524  const Index point_index = bi->first;
525  const Scalar point_orientation = bi->second;
526  ASSERT( point_index < static_cast<Index>(myIndexSignedCells[actualOrder(0, duality)].size()) );
527  ASSERT( point_index < kFormLength(0, duality) );
528 
529  triplets[edge_direction].push_back( Triplet(edge_index, point_index, point_orientation*edge_length*edge_sign*edge_orientation/border_infos.size()) );
530  }
531 
532  }
533 
534  boost::array<SparseMatrix, dimAmbient> flat_operator_matrix;
535 
536  for (DGtal::Dimension direction=0; direction<dimAmbient; direction++)
537  {
538  flat_operator_matrix[direction] = SparseMatrix(kFormLength(1, duality), kFormLength(0, duality));
539  flat_operator_matrix[direction].setFromTriplets(triplets[direction].begin(), triplets[direction].end());
540  }
541 
542  myFlatOperatorMatrixes[static_cast<int>(duality)] = flat_operator_matrix;
543 }
544 
545 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
546 void
547 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::updateIndexes()
548 {
549  if (!myIndexesNeedUpdate) return;
550 
551  // clear index signed cells
552  for (DGtal::Dimension dim=0; dim<dimEmbedded+1; dim++)
553  myIndexSignedCells[dim].clear();
554 
555  // compute cell index
556  for (typename Properties::iterator csi=myCellProperties.begin(), csie=myCellProperties.end(); csie!=csi; csi++)
557  {
558  const Cell& cell = csi->first;
559  const DGtal::Dimension cell_dim = myKSpace.uDim(cell);
560 
561  csi->second.index = myIndexSignedCells[cell_dim].size();
562 
563  const SCell& signed_cell = myKSpace.signs(cell, csi->second.flipped ? KSpace::NEG : KSpace::POS);
564  myIndexSignedCells[cell_dim].push_back(signed_cell);
565  }
566 
567  myIndexesNeedUpdate = false;
568  myCachedOperatorsNeedUpdate = true;
569 }
570 
571 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
572 void
573 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::updateCachedOperators()
574 {
575  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
576  if (!myCachedOperatorsNeedUpdate) return;
577  updateFlatOperator<PRIMAL>();
578  updateFlatOperator<DUAL>();
579  updateSharpOperator<PRIMAL>();
580  updateSharpOperator<DUAL>();
581  myCachedOperatorsNeedUpdate = false;
582 }
583 
584 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
585 const typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Properties&
586 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getProperties() const
587 {
588  return myCellProperties;
589 }
590 
591 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
592 template <DGtal::Order order, DGtal::Duality duality>
593 const typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::SCells&
594 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getIndexedSCells() const
595 {
596  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
597  return myIndexSignedCells[actualOrder(order, duality)];
598 }
599 
600 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
601 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::SCell
602 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getSCell(const Order& order, const Duality& duality, const Index& index) const
603 {
604  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
605  const Order& actual_order = actualOrder(order, duality);
606  const SCell& signed_cell = myIndexSignedCells[actual_order][index];
607  ASSERT( myKSpace.sDim(signed_cell) == actual_order );
608  return signed_cell;
609 }
610 
611 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
612 bool
613 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::containsCell(const Cell& cell) const
614 {
615  return myCellProperties.find(cell) != myCellProperties.end();
616 }
617 
618 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
619 bool
620 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::isCellFlipped(const Cell& cell) const
621 {
622  const typename Properties::const_iterator iter_property = myCellProperties.find(cell);
623  ASSERT( iter_property != myCellProperties.end() );
624  return iter_property->second.flipped;
625 }
626 
627 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
628 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Index
629 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::getCellIndex(const Cell& cell) const
630 {
631  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
632  const typename Properties::const_iterator iter_property = myCellProperties.find(cell);
633  ASSERT( iter_property != myCellProperties.end() );
634  return iter_property->second.index;
635 }
636 
637 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
638 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::ConstIterator
639 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::begin() const
640 {
641  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
642  return myCellProperties.begin();
643 }
644 
645 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
646 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::ConstIterator
647 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::end() const
648 {
649  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
650  return myCellProperties.end();
651 }
652 
653 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
654 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Iterator
655 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::begin()
656 {
657  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
658  return myCellProperties.begin();
659 }
660 
661 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
662 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Iterator
663 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::end()
664 {
665  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
666  return myCellProperties.end();
667 }
668 
669 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
670 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Index
671 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::kFormLength(const DGtal::Order& order, const DGtal::Duality& duality) const
672 {
673  ASSERT_MSG( !myIndexesNeedUpdate, "call updateIndexes() after manual structure modification" );
674  return myIndexSignedCells[actualOrder(order, duality)].size();
675 }
676 
677 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
678 DGtal::Order
679 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::actualOrder(const DGtal::Order& order, const DGtal::Duality& duality) const
680 {
681  return duality == PRIMAL ? order : dimEmbedded-order;
682 }
683 
684 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
685 typename DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::Scalar
686 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::hodgeSign(const Cell& cell, const DGtal::Duality& duality) const
687 {
688  if (duality == PRIMAL) return 1;
689 
690  const DGtal::Dimension& primal_dim = myKSpace.uDim(cell);
691  return (dimEmbedded-primal_dim)*primal_dim % 2 != 0 ? -1 : 1;
692 }
693 
694 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
695 typename DGtal::Dimension
696 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::edgeDirection(const Cell& cell, const DGtal::Duality& duality) const
697 {
698  ASSERT( myKSpace.uDim(cell) == actualOrder(1, duality) );
699 
700  typename KSpace::DirIterator di = myKSpace.uDirs(cell);
701  if (duality == DUAL) di = myKSpace.uOrthDirs(cell);
702 
703  ASSERT( di != 0 );
704  const DGtal::Dimension direction = *di;
705  ++di;
706  //ASSERT( !(di != 0) ); //FIXME multiple edge direction with embedding
707  ASSERT( direction < dimAmbient );
708 
709  return direction;
710 }
711 
712 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
713 bool
714 DGtal::DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>::isValid() const
715 {
716  return true;
717 }
718 
719 ///////////////////////////////////////////////////////////////////////////////
720 // Implementation of inline functions //
721 
722 template <DGtal::Dimension dimEmbedded, DGtal::Dimension dimAmbient, typename TLinearAlgebraBackend, typename TInteger>
723 std::ostream&
724 DGtal::operator<<(std::ostream & out, const DiscreteExteriorCalculus<dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger>& object)
725 {
726  object.selfDisplay(out);
727  return out;
728 }
729 
730 // //
731 ///////////////////////////////////////////////////////////////////////////////