DGtal  1.4.beta
DiscreteExteriorCalculusFactory.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 DiscreteExteriorCalculusFactory.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 2015/05/04
23  *
24  * Implementation of inline methods for DiscreteExteriorCalculusFactory
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 ///////////////////////////////////////////////////////////////////////////////
30 // Implementation of inline methods //
31 
32 template <typename TLinearAlgebraBackend, typename TInteger>
33 template <typename TDigitalSet>
34 DGtal::DiscreteExteriorCalculus<TDigitalSet::Point::dimension, TDigitalSet::Point::dimension, TLinearAlgebraBackend, TInteger>
35 DGtal::DiscreteExteriorCalculusFactory<TLinearAlgebraBackend, TInteger>::createFromDigitalSet(const TDigitalSet& _set, const bool add_border)
36 {
37  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDigitalSet<TDigitalSet> ));
38 
39  typedef DGtal::DiscreteExteriorCalculus<TDigitalSet::Point::dimension, TDigitalSet::Point::dimension, TLinearAlgebraBackend, TInteger> Calculus;
40  typedef typename Calculus::Cell Cell;
41  typedef typename Calculus::Scalar Scalar;
42  typedef typename Calculus::Index Index;
43  typedef typename Calculus::Point Point;
44  typedef typename Calculus::Property Property;
45 
46  Calculus calculus;
47  calculus.template initKSpace<typename TDigitalSet::Domain>(_set.domain());
48 
49  // compute raw cell size
50  typedef std::map<Cell, Scalar> Accum;
51  Accum cell_size_accum;
52  for (typename TDigitalSet::ConstIterator ri=_set.begin(), rie=_set.end(); ri!=rie; ri++)
53  {
54  const Point& point = *ri;
55  const Cell cell_point = calculus.myKSpace.uSpel(point);
56 
57  typedef DGtal::SpaceND<Calculus::dimensionAmbient, TInteger> Space;
58  typedef DGtal::HyperRectDomain<Space> Neighborbood;
59  const Point cell_coords = calculus.myKSpace.uKCoords(cell_point);
60  const Neighborbood neighborhood(cell_coords-Point::diagonal(1), cell_coords+Point::diagonal(1));
61  for (typename Neighborbood::ConstIterator pi=neighborhood.begin(), pie=neighborhood.end(); pi!=pie; pi++)
62  {
63  const Cell cell = calculus.myKSpace.uCell(*pi);
64  if (cell_size_accum.find(cell) == cell_size_accum.end()) cell_size_accum[cell] = 0;
65  cell_size_accum[cell] += 1;
66  }
67  }
68 
69  // normalize cell size and set flipped flag
70  for (typename Accum::const_iterator csi=cell_size_accum.begin(), csie=cell_size_accum.end(); csie!=csi; csi++)
71  {
72  const Cell& cell = csi->first;
73  const DGtal::Dimension dual_dim = Calculus::dimensionEmbedded-calculus.myKSpace.uDim(cell);
74  const Scalar factor = pow(.5, dual_dim);
75 
76  const Scalar normalized_size = csi->second * factor;
77  ASSERT(normalized_size > 0 && normalized_size <= 1);
78 
79  if (!add_border && normalized_size < 1) continue;
80 
81  Property property;
82  property.primal_size = 1;
83  property.dual_size = normalized_size;
84  property.index = std::numeric_limits<Index>::max();
85  property.flipped = false;
86 
87  calculus.myCellProperties[cell] = property;
88  }
89 
90  calculus.myIndexesNeedUpdate = true;
91  calculus.updateIndexes();
92 
93  return calculus;
94 }
95 
96 template <typename TLinearAlgebraBackend, typename TInteger>
97 template <typename KSpace, typename CellsSet>
98 void
99 DGtal::DiscreteExteriorCalculusFactory<TLinearAlgebraBackend, TInteger>::insertAllLowerIncidentCells(const KSpace& kspace, const typename CellsSet::value_type& cell, CellsSet& cells_set)
100 {
101  typedef typename KSpace::Cells Cells;
102 
103  cells_set.insert(cell);
104 
105  const Cells border = kspace.uLowerIncident(cell);
106  for (typename Cells::ConstIterator bi=border.begin(), be=border.end(); bi!=be; bi++)
107  insertAllLowerIncidentCells(kspace, *bi, cells_set);
108 }
109 
110 template <typename TLinearAlgebraBackend, typename TInteger>
111 template <typename KSpace, typename CellsAccum>
112 void
113 DGtal::DiscreteExteriorCalculusFactory<TLinearAlgebraBackend, TInteger>::accumulateAllLowerIncidentCells(const KSpace& kspace, const typename CellsAccum::key_type& cell, CellsAccum& cells_accum)
114 {
115  typedef typename KSpace::Cells Cells;
116 
117  if (cells_accum.find(cell) == cells_accum.end()) cells_accum[cell] = 0;
118  cells_accum[cell]++;
119 
120  const Cells border = kspace.uLowerIncident(cell);
121  for (typename Cells::ConstIterator bi=border.begin(), be=border.end(); bi!=be; bi++)
122  accumulateAllLowerIncidentCells(kspace, *bi, cells_accum);
123 }
124 
125 template <typename TLinearAlgebraBackend, typename TInteger>
126 template <typename KSpace, typename CellsAccum, typename MeasureAccum>
127 void
128 DGtal::DiscreteExteriorCalculusFactory<TLinearAlgebraBackend, TInteger>::accumulateAllLowerIncidentCells(const KSpace& kspace, const typename CellsAccum::key_type& cell, CellsAccum& cells_accum,
129  CellsAccum& local_accum, MeasureAccum& cell_to_measure, const double measure)
130 {
131  typedef typename KSpace::Cells Cells;
132 
133  if (cells_accum.find(cell) == cells_accum.end()) cells_accum[cell] = 0;
134  cells_accum[cell]++;
135 
136  if( kspace.uDim( cell ) == 0 && local_accum.find( cell ) == local_accum.end())
137  {
138  local_accum[cell] = 0;
139  if( cell_to_measure.find(cell) == cell_to_measure.end() ) cell_to_measure[cell] = 0.;
140  cell_to_measure[cell] += measure;
141  }
142 
143  const Cells border = kspace.uLowerIncident(cell);
144  for (typename Cells::ConstIterator bi=border.begin(), be=border.end(); bi!=be; bi++)
145  accumulateAllLowerIncidentCells(kspace, *bi, cells_accum, local_accum, cell_to_measure, measure);
146 }
147 
148 template <typename TLinearAlgebraBackend, typename TInteger>
149 template <DGtal::Dimension dimEmbedded, typename TNSCellConstIterator>
150 DGtal::DiscreteExteriorCalculus<dimEmbedded, TNSCellConstIterator::value_type::Point::dimension, TLinearAlgebraBackend, TInteger>
151 DGtal::DiscreteExteriorCalculusFactory<TLinearAlgebraBackend, TInteger>::createFromNSCells(const TNSCellConstIterator& begin, const TNSCellConstIterator& end, const bool add_border)
152 {
153  BOOST_CONCEPT_ASSERT(( boost_concepts::SinglePassIteratorConcept<TNSCellConstIterator> ));
154  BOOST_CONCEPT_ASSERT(( boost_concepts::ReadableIteratorConcept<TNSCellConstIterator> ));
155 
156  typedef DGtal::DiscreteExteriorCalculus<dimEmbedded, TNSCellConstIterator::value_type::Point::dimension, TLinearAlgebraBackend, TInteger> Calculus;
157  typedef typename Calculus::Cell Cell;
158  typedef typename Calculus::SCell SCell;
159  typedef typename Calculus::Scalar Scalar;
160  typedef typename Calculus::KSpace KSpace;
161  typedef typename KSpace::Cells Cells;
162 
163  BOOST_STATIC_ASSERT(( boost::is_convertible<typename TNSCellConstIterator::value_type, const SCell>::value ));
164 
165  Calculus calculus;
166 
167  // compute dimEmbedded-1 cells border
168  typedef std::map<Cell, int> CellsAccum;
169  CellsAccum border_accum;
170  CellsAccum lower_accum;
171  for (TNSCellConstIterator ci=begin; ci!=end; ++ci)
172  {
173  const SCell cell_signed = *ci;
174  const Dimension cell_dim = calculus.myKSpace.sDim(cell_signed);
175  ASSERT_MSG( cell_dim == dimEmbedded, "wrong n-cell dimension" );
176  boost::ignore_unused_variable_warning(cell_dim);
177 
178  calculus.insertSCell(cell_signed);
179 
180  const Cell cell = calculus.myKSpace.unsigns(cell_signed);
181  const Cells border = calculus.myKSpace.uLowerIncident(cell);
182  for (typename Cells::ConstIterator bi=border.begin(), be=border.end(); bi!=be; bi++)
183  {
184  const Cell cell_border = *bi;
185  if (!add_border)
186  {
187  if (border_accum.find(cell_border) == border_accum.end()) border_accum[cell_border] = 0;
188  border_accum[cell_border]++;
189  }
190  accumulateAllLowerIncidentCells(calculus.myKSpace, cell_border, lower_accum);
191  }
192  }
193  ASSERT( !add_border || border_accum.empty() );
194 
195  typedef std::set<Cell> CellsSet;
196  CellsSet border;
197  for (typename CellsAccum::const_iterator bai=border_accum.begin(), bae=border_accum.end(); bai!=bae; bai++)
198  {
199  ASSERT( bai->second > 0 );
200  //ASSERT( bai->second < 3 );
201  if (bai->second >= 2) continue;
202  insertAllLowerIncidentCells(calculus.myKSpace, bai->first, border);
203  }
204  ASSERT( !add_border || border.empty() );
205 
206  // normalize cell size and set flipped flag
207  for (typename CellsAccum::const_iterator lai=lower_accum.begin(), lae=lower_accum.end(); lai!=lae; ++lai)
208  {
209  const Cell cell = lai->first;
210  if (border.find(cell) != border.end()) continue;
211 
212  const DGtal::Dimension dual_dim = Calculus::dimensionEmbedded-calculus.myKSpace.uDim(cell);
213  ASSERT( dual_dim > 0 );
214  ASSERT( dual_dim <= Calculus::dimensionEmbedded );
215 
216  const Scalar factor = pow(.5, 2*dual_dim-1);
217  const Scalar dual_size = lai->second * factor;
218  ASSERT( dual_size > 0 );
219 
220  const SCell cell_signed = calculus.myKSpace.signs(cell, KSpace::POS);
221  calculus.insertSCell(cell_signed, 1, dual_size);
222  }
223 
224  calculus.updateIndexes();
225 
226  return calculus;
227 }
228 
229 template <typename TLinearAlgebraBackend, typename TInteger>
230 template <DGtal::Dimension dimEmbedded, typename TNSCellConstIterator, typename TSCellMeasureFunctor>
231 DGtal::DiscreteExteriorCalculus<dimEmbedded, TNSCellConstIterator::value_type::Point::dimension, TLinearAlgebraBackend, TInteger>
232 DGtal::DiscreteExteriorCalculusFactory<TLinearAlgebraBackend, TInteger>::createFromNSCells(const TNSCellConstIterator& begin, const TNSCellConstIterator& end,
233  const TSCellMeasureFunctor& normalFunctor, const double h, const bool add_border)
234 {
235  BOOST_CONCEPT_ASSERT(( boost_concepts::SinglePassIteratorConcept<TNSCellConstIterator> ));
236  BOOST_CONCEPT_ASSERT(( boost_concepts::ReadableIteratorConcept<TNSCellConstIterator> ));
237 
238  typedef DGtal::DiscreteExteriorCalculus<dimEmbedded, TNSCellConstIterator::value_type::Point::dimension, TLinearAlgebraBackend, TInteger> Calculus;
239  typedef typename Calculus::Cell Cell;
240  typedef typename Calculus::SCell SCell;
241  typedef typename Calculus::Scalar Scalar;
242  typedef typename Calculus::KSpace KSpace;
243  typedef typename KSpace::Cells Cells;
244 
245  BOOST_STATIC_ASSERT(( boost::is_convertible<typename TNSCellConstIterator::value_type, const SCell>::value ));
246 
247  typedef typename TSCellMeasureFunctor::Quantity Quantity;
248  std::vector<Quantity> measures;
249  std::back_insert_iterator< std::vector<Quantity> > back_insert( measures );
250  normalFunctor.eval(begin, end, back_insert);
251  typename std::vector<Quantity>::const_iterator measure_it = measures.begin();
252 
253  Calculus calculus;
254  // compute dimEmbedded-1 cells border
255  std::map<Cell, double> cell_to_measure;
256 
257  typedef std::map<Cell, int> CellsAccum;
258  CellsAccum border_accum;
259  CellsAccum lower_accum;
260  for (TNSCellConstIterator ci=begin; ci!=end; ++ci)
261  {
262  const SCell cell_signed = *ci;
263  const Dimension cell_dim = calculus.myKSpace.sDim(cell_signed);
264  ASSERT_MSG( cell_dim == dimEmbedded, "wrong n-cell dimension" );
265  boost::ignore_unused_variable_warning(cell_dim);
266 
267  const int orth_dir = calculus.myKSpace.sOrthDir( cell_signed );
268  const double measure = pow(h, dimEmbedded) * fabs( (*measure_it)[orth_dir] );
269  ++measure_it;
270 
271  calculus.insertSCell(cell_signed, measure, 1.);
272 
273  const Cell cell = calculus.myKSpace.unsigns(cell_signed);
274  const Cells border = calculus.myKSpace.uLowerIncident(cell);
275  CellsAccum local_accum;
276  for (typename Cells::ConstIterator bi=border.begin(), be=border.end(); bi!=be; bi++)
277  {
278  const Cell cell_border = *bi;
279  if (!add_border)
280  {
281  if (border_accum.find(cell_border) == border_accum.end()) border_accum[cell_border] = 0;
282  border_accum[cell_border]++;
283  }
284  accumulateAllLowerIncidentCells(calculus.myKSpace, cell_border, lower_accum, local_accum, cell_to_measure, measure / 4.);
285  }
286  }
287  ASSERT( !add_border || border_accum.empty() );
288 
289  typedef std::set<Cell> CellsSet;
290  CellsSet border;
291  for (typename CellsAccum::const_iterator bai=border_accum.begin(), bae=border_accum.end(); bai!=bae; bai++)
292  {
293  ASSERT( bai->second > 0 );
294  //ASSERT( bai->second < 3 );
295  if (bai->second >= 2) continue;
296  insertAllLowerIncidentCells(calculus.myKSpace, bai->first, border);
297  }
298  ASSERT( !add_border || border.empty() );
299 
300  // normalize cell size and set flipped flag
301  for (typename CellsAccum::const_iterator lai=lower_accum.begin(), lae=lower_accum.end(); lai!=lae; ++lai)
302  {
303  const Cell cell = lai->first;
304  if (border.find(cell) != border.end()) continue;
305  const SCell cell_signed = calculus.myKSpace.signs(cell, KSpace::POS);
306 
307  if( cell_to_measure.find( cell ) != cell_to_measure.end() )
308  {
309  const double dual_measure = cell_to_measure.find( cell )->second;
310  if(cell_to_measure.find( cell ) == cell_to_measure.end()) trace.warning() << "HOUUU" << std::endl;
311  calculus.insertSCell(cell_signed, 1, dual_measure);
312  continue;
313  }
314 
315  const DGtal::Dimension dual_dim = Calculus::dimensionEmbedded-calculus.myKSpace.uDim(cell);
316  ASSERT( dual_dim > 0 );
317  ASSERT( dual_dim <= Calculus::dimensionEmbedded );
318 
319  const Scalar factor = pow(.5, 2*dual_dim-1);
320  const Scalar dual_size = lai->second * factor;
321  ASSERT( dual_size > 0 );
322 
323  calculus.insertSCell(cell_signed, 1, dual_size);
324  }
325 
326  calculus.updateIndexes();
327 
328  return calculus;
329 }
330 
331 // //
332 ///////////////////////////////////////////////////////////////////////////////
333 
334