DGtal 1.3.0
Loading...
Searching...
No Matches
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
32template <typename TLinearAlgebraBackend, typename TInteger>
33template <typename TDigitalSet>
34DGtal::DiscreteExteriorCalculus<TDigitalSet::Point::dimension, TDigitalSet::Point::dimension, TLinearAlgebraBackend, TInteger>
35DGtal::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
96template <typename TLinearAlgebraBackend, typename TInteger>
97template <typename KSpace, typename CellsSet>
98void
99DGtal::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
110template <typename TLinearAlgebraBackend, typename TInteger>
111template <typename KSpace, typename CellsAccum>
112void
113DGtal::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
125template <typename TLinearAlgebraBackend, typename TInteger>
126template <typename KSpace, typename CellsAccum, typename MeasureAccum>
127void
128DGtal::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
148template <typename TLinearAlgebraBackend, typename TInteger>
149template <DGtal::Dimension dimEmbedded, typename TNSCellConstIterator>
150DGtal::DiscreteExteriorCalculus<dimEmbedded, TNSCellConstIterator::value_type::Point::dimension, TLinearAlgebraBackend, TInteger>
151DGtal::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
229template <typename TLinearAlgebraBackend, typename TInteger>
230template <DGtal::Dimension dimEmbedded, typename TNSCellConstIterator, typename TSCellMeasureFunctor>
231DGtal::DiscreteExteriorCalculus<dimEmbedded, TNSCellConstIterator::value_type::Point::dimension, TLinearAlgebraBackend, TInteger>
232DGtal::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 double primal_measure_accum = 0., dual_measure_accum = 0.;
257
258 typedef std::map<Cell, int> CellsAccum;
259 CellsAccum border_accum;
260 CellsAccum lower_accum;
261 for (TNSCellConstIterator ci=begin; ci!=end; ++ci)
262 {
263 const SCell cell_signed = *ci;
264 const Dimension cell_dim = calculus.myKSpace.sDim(cell_signed);
265 ASSERT_MSG( cell_dim == dimEmbedded, "wrong n-cell dimension" );
266 boost::ignore_unused_variable_warning(cell_dim);
267
268 const int orth_dir = calculus.myKSpace.sOrthDir( cell_signed );
269 const double measure = pow(h, dimEmbedded) * fabs( (*measure_it)[orth_dir] );
270 primal_measure_accum += measure;
271 ++measure_it;
272
273 calculus.insertSCell(cell_signed, measure, 1.);
274
275 const Cell cell = calculus.myKSpace.unsigns(cell_signed);
276 const Cells border = calculus.myKSpace.uLowerIncident(cell);
277 CellsAccum local_accum;
278 for (typename Cells::ConstIterator bi=border.begin(), be=border.end(); bi!=be; bi++)
279 {
280 const Cell cell_border = *bi;
281 if (!add_border)
282 {
283 if (border_accum.find(cell_border) == border_accum.end()) border_accum[cell_border] = 0;
284 border_accum[cell_border]++;
285 }
286 accumulateAllLowerIncidentCells(calculus.myKSpace, cell_border, lower_accum, local_accum, cell_to_measure, measure / 4.);
287 }
288 }
289 ASSERT( !add_border || border_accum.empty() );
290
291 typedef std::set<Cell> CellsSet;
292 CellsSet border;
293 for (typename CellsAccum::const_iterator bai=border_accum.begin(), bae=border_accum.end(); bai!=bae; bai++)
294 {
295 ASSERT( bai->second > 0 );
296 //ASSERT( bai->second < 3 );
297 if (bai->second >= 2) continue;
298 insertAllLowerIncidentCells(calculus.myKSpace, bai->first, border);
299 }
300 ASSERT( !add_border || border.empty() );
301
302 // normalize cell size and set flipped flag
303 for (typename CellsAccum::const_iterator lai=lower_accum.begin(), lae=lower_accum.end(); lai!=lae; ++lai)
304 {
305 const Cell cell = lai->first;
306 if (border.find(cell) != border.end()) continue;
307 const SCell cell_signed = calculus.myKSpace.signs(cell, KSpace::POS);
308
309 if( cell_to_measure.find( cell ) != cell_to_measure.end() )
310 {
311 const double dual_measure = cell_to_measure.find( cell )->second;
312 if(cell_to_measure.find( cell ) == cell_to_measure.end()) trace.warning() << "HOUUU" << std::endl;
313 dual_measure_accum += dual_measure;
314 calculus.insertSCell(cell_signed, 1, dual_measure);
315 continue;
316 }
317
318 const DGtal::Dimension dual_dim = Calculus::dimensionEmbedded-calculus.myKSpace.uDim(cell);
319 ASSERT( dual_dim > 0 );
320 ASSERT( dual_dim <= Calculus::dimensionEmbedded );
321
322 const Scalar factor = pow(.5, 2*dual_dim-1);
323 const Scalar dual_size = lai->second * factor;
324 ASSERT( dual_size > 0 );
325
326 calculus.insertSCell(cell_signed, 1, dual_size);
327 }
328
329 //ASSERT( primal_measure_accum == dual_measure_accum );
330
331 calculus.updateIndexes();
332
333 return calculus;
334}
335
336// //
337///////////////////////////////////////////////////////////////////////////////
338
339