DGtal  0.9.2
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 
126 template <typename TLinearAlgebraBackend, typename TInteger>
127 template <DGtal::Dimension dimEmbedded, typename TNSCellConstIterator>
128 DGtal::DiscreteExteriorCalculus<dimEmbedded, TNSCellConstIterator::value_type::Point::dimension, TLinearAlgebraBackend, TInteger>
129 DGtal::DiscreteExteriorCalculusFactory<TLinearAlgebraBackend, TInteger>::createFromNSCells(const TNSCellConstIterator& begin, const TNSCellConstIterator& end, const bool add_border)
130 {
131  BOOST_CONCEPT_ASSERT(( boost_concepts::SinglePassIteratorConcept<TNSCellConstIterator> ));
132  BOOST_CONCEPT_ASSERT(( boost_concepts::ReadableIteratorConcept<TNSCellConstIterator> ));
133 
134  typedef DGtal::DiscreteExteriorCalculus<dimEmbedded, TNSCellConstIterator::value_type::Point::dimension, TLinearAlgebraBackend, TInteger> Calculus;
135  typedef typename Calculus::Cell Cell;
136  typedef typename Calculus::SCell SCell;
137  typedef typename Calculus::Scalar Scalar;
138  typedef typename Calculus::KSpace KSpace;
139  typedef typename KSpace::Cells Cells;
140 
141  BOOST_STATIC_ASSERT(( boost::is_convertible<typename TNSCellConstIterator::value_type, const SCell>::value ));
142 
143  Calculus calculus;
144 
145  // compute dimEmbedded-1 cells border
146  typedef std::map<Cell, int> CellsAccum;
147  CellsAccum border_accum;
148  CellsAccum lower_accum;
149  for (TNSCellConstIterator ci=begin; ci!=end; ++ci)
150  {
151  const SCell cell_signed = *ci;
152  const Dimension cell_dim = calculus.myKSpace.sDim(cell_signed);
153  ASSERT_MSG( cell_dim == dimEmbedded, "wrong n-cell dimension" );
154  boost::ignore_unused_variable_warning(cell_dim);
155 
156  calculus.insertSCell(cell_signed);
157 
158  const Cell cell = calculus.myKSpace.unsigns(cell_signed);
159  const Cells border = calculus.myKSpace.uLowerIncident(cell);
160  for (typename Cells::ConstIterator bi=border.begin(), be=border.end(); bi!=be; bi++)
161  {
162  const Cell cell_border = *bi;
163  if (!add_border)
164  {
165  if (border_accum.find(cell_border) == border_accum.end()) border_accum[cell_border] = 0;
166  border_accum[cell_border]++;
167  }
168  accumulateAllLowerIncidentCells(calculus.myKSpace, cell_border, lower_accum);
169  }
170  }
171  ASSERT( !add_border || border_accum.empty() );
172 
173  typedef std::set<Cell> CellsSet;
174  CellsSet border;
175  for (typename CellsAccum::const_iterator bai=border_accum.begin(), bae=border_accum.end(); bai!=bae; bai++)
176  {
177  ASSERT( bai->second > 0 );
178  //ASSERT( bai->second < 3 );
179  if (bai->second >= 2) continue;
180  insertAllLowerIncidentCells(calculus.myKSpace, bai->first, border);
181  }
182  ASSERT( !add_border || border.empty() );
183 
184  // normalize cell size and set flipped flag
185  for (typename CellsAccum::const_iterator lai=lower_accum.begin(), lae=lower_accum.end(); lai!=lae; ++lai)
186  {
187  const Cell cell = lai->first;
188  if (border.find(cell) != border.end()) continue;
189 
190  const DGtal::Dimension dual_dim = Calculus::dimensionEmbedded-calculus.myKSpace.uDim(cell);
191  ASSERT( dual_dim > 0 );
192  ASSERT( dual_dim <= Calculus::dimensionEmbedded );
193 
194  const Scalar factor = pow(.5, 2*dual_dim-1);
195  const Scalar dual_size = lai->second * factor;
196  ASSERT( dual_size > 0 );
197 
198  const SCell cell_signed = calculus.myKSpace.signs(cell, KSpace::POS);
199  calculus.insertSCell(cell_signed, 1, dual_size);
200  }
201 
202  calculus.updateIndexes();
203 
204  return calculus;
205 }
206 
207 // //
208 ///////////////////////////////////////////////////////////////////////////////
209 
210