DGtal 2.1.0
Loading...
Searching...
No Matches
BoundedLatticePolytopeCounter.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 BoundedLatticePolytopeCounter.ih
19 * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20 * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
21 *
22 * @date 2022/06/17
23 *
24 * Implementation of inline methods defined in BoundedLatticePolytopeCounter.h
25 *
26 * This file is part of the DGtal library.
27 */
28
29
30//////////////////////////////////////////////////////////////////////////////
31#include <cstdlib>
32//////////////////////////////////////////////////////////////////////////////
33
34///////////////////////////////////////////////////////////////////////////////
35// IMPLEMENTATION of inline methods.
36///////////////////////////////////////////////////////////////////////////////
37
38///////////////////////////////////////////////////////////////////////////////
39// ----------------------- Standard services ------------------------------
40
41//-----------------------------------------------------------------------------
42template <typename TSpace>
43DGtal::BoundedLatticePolytopeCounter<TSpace>::
44BoundedLatticePolytopeCounter
45( const Polytope& P )
46{
47 init( &P );
48}
49
50//-----------------------------------------------------------------------------
51template <typename TSpace>
52void
53DGtal::BoundedLatticePolytopeCounter<TSpace>::
54init
55( const Polytope* ptrP )
56{
57 myPolytope = ptrP;
58 if ( ptrP == nullptr ) return;
59 myLower = ptrP->getDomain().lowerBound();
60 myUpper = ptrP->getDomain().upperBound();
61}
62
63
64//-----------------------------------------------------------------------------
65template <typename TSpace>
66typename DGtal::BoundedLatticePolytopeCounter<TSpace>::Interval
67DGtal::BoundedLatticePolytopeCounter<TSpace>::
68intersectionIntervalAlongAxis( Point p, Dimension a ) const
69{
70 ASSERT( myPolytope != nullptr );
71 const Polytope& P = *myPolytope;
72 const InequalityMatrix& A = P.getA();
73 const InequalityVector& B = P.getB();
74 const std::vector<bool>& I = P.getI();
75 Integer x_min = myLower[ a ];
76 Integer x_max = myUpper[ a ]+1;
77 Integer x = 0;
78 const Integer x_a = x_min;
79 p[ a ] = x_a;
80 bool empty = false;
81 for ( Dimension k = 0; k < A.size(); k++ )
82 {
83 const Integer c = A[ k ].dot( p );
84 const Integer n = A[ k ][ a ];
85 const Integer b = B[ k ];
86 if ( n == 0 )
87 { // constraint is // to the specified axis.
88 empty = ! ( I[ k ] ? ( c <= b ) : c < b );
89 }
90 else if ( n > 0 )
91 {
92 Integer d = b - c;
93 if ( d < 0 ) empty = true;
94 else
95 {
96 x = I[ k ] ? ( d / n + 1 ) : ( (d+n-1) / n ) ;
97 x_max = std::min( x_max, x_a + x );
98 }
99 }
100 else // ( n < 0 )
101 {
102 Integer d = c - b;
103 if ( d >= 0 )
104 {
105 x = I[ k ] ? ( (d-n-1) / -n ) : ( d / -n + 1 );
106 x_min = std::max( x_min, x_a + x );
107 }
108 // otherwise the constraint is true
109 }
110 if ( empty || ( x_max <= x_min ) ) return Interval( 0, 0 );
111 }
112 return Interval( x_min, x_max );
113}
114
115//-----------------------------------------------------------------------------
116template <typename TSpace>
117typename DGtal::BoundedLatticePolytopeCounter<TSpace>::Interval
118DGtal::BoundedLatticePolytopeCounter<TSpace>::
119interiorIntersectionIntervalAlongAxis( Point p, Dimension a ) const
120{
121 ASSERT( myPolytope != nullptr );
122 const Polytope& P = *myPolytope;
123 const InequalityMatrix& A = P.getA();
124 const InequalityVector& B = P.getB();
125 Integer x_min = myLower[ a ];
126 Integer x_max = myUpper[ a ]+1;
127 Integer x = 0;
128 const Integer x_a = x_min;
129 p[ a ] = x_a;
130 bool empty = false;
131 // We must take into account also bounding box constraints for interior points.
132 for ( Dimension k = 0; k < A.size(); k++ )
133 {
134 const Integer c = A[ k ].dot( p );
135 const Integer n = A[ k ][ a ];
136 const Integer b = B[ k ];
137 if ( n == 0 )
138 { // constraint is // to the specified axis.
139 empty = ( b <= c );
140 }
141 else if ( n > 0 )
142 {
143 Integer d = b - c;
144 if ( d < 0 ) empty = true;
145 else
146 {
147 x = (d+n-1) / n;
148 x_max = std::min( x_max, x_a + x );
149 }
150 }
151 else // ( n < 0 )
152 {
153 Integer d = c - b;
154 if ( d >= 0 )
155 {
156 x = d / -n + 1;
157 x_min = std::max( x_min, x_a + x );
158 }
159 // otherwise the constraint is true
160 }
161 if ( empty || ( x_max <= x_min ) ) return Interval( 0, 0 );
162 }
163 return Interval( x_min, x_max );
164}
165
166//-----------------------------------------------------------------------------
167template <typename TSpace>
168typename DGtal::BoundedLatticePolytopeCounter<TSpace>::Integer
169DGtal::BoundedLatticePolytopeCounter<TSpace>::
170countAlongAxis( Dimension a ) const
171{
172 ASSERT( myPolytope != nullptr );
173 Point lo = myLower;
174 Point hi = myUpper;
175 hi[ a ] = lo[ a ];
176 Domain D( lo, hi );
177 Integer nb = 0;
178 for ( auto&& p : D )
179 {
180 auto I = intersectionIntervalAlongAxis( p, a );
181 nb += I.second - I.first;
182 }
183 return nb;
184}
185
186//-----------------------------------------------------------------------------
187template <typename TSpace>
188typename DGtal::BoundedLatticePolytopeCounter<TSpace>::Integer
189DGtal::BoundedLatticePolytopeCounter<TSpace>::
190countInteriorAlongAxis( Dimension a ) const
191{
192 ASSERT( myPolytope != nullptr );
193 Point lo = myLower;
194 Point hi = myUpper;
195 hi[ a ] = lo[ a ];
196 Domain D( lo, hi );
197 Integer nb = 0;
198 for ( auto&& p : D )
199 {
200 auto I = interiorIntersectionIntervalAlongAxis( p, a );
201 nb += I.second - I.first;
202 }
203 return nb;
204}
205
206//-----------------------------------------------------------------------------
207template <typename TSpace>
208void
209DGtal::BoundedLatticePolytopeCounter<TSpace>::
210getPointsAlongAxis( PointRange& pts, Dimension a ) const
211{
212 ASSERT( myPolytope != nullptr );
213 Point lo = myLower;
214 Point hi = myUpper;
215 hi[ a ] = lo[ a ];
216 Domain D( lo, hi );
217 for ( auto&& p : D )
218 {
219 auto I = intersectionIntervalAlongAxis( p, a );
220 Point q = p;
221 for ( Integer x = I.first; x != I.second; x++ )
222 {
223 q[ a ] = x;
224 pts.push_back( q );
225 }
226 }
227}
228
229//-----------------------------------------------------------------------------
230template <typename TSpace>
231void
232DGtal::BoundedLatticePolytopeCounter<TSpace>::
233getInteriorPointsAlongAxis( PointRange& pts, Dimension a ) const
234{
235 ASSERT( myPolytope != nullptr );
236 Point lo = myLower;
237 Point hi = myUpper;
238 hi[ a ] = lo[ a ];
239 Domain D( lo, hi );
240 //Integer nb = 0; not used
241 for ( auto&& p : D )
242 {
243 auto I = interiorIntersectionIntervalAlongAxis( p, a );
244 Point q = p;
245 for ( Integer x = I.first; x != I.second; x++ )
246 {
247 q[ a ] = x;
248 pts.push_back( q );
249 }
250 }
251}
252
253
254//-----------------------------------------------------------------------------
255template <typename TSpace>
256typename DGtal::BoundedLatticePolytopeCounter<TSpace>::LatticeSetByInterval
257DGtal::BoundedLatticePolytopeCounter<TSpace>::
258getLatticeSet( Dimension a ) const
259{
260 ASSERT( myPolytope != nullptr );
261 Point lo = myLower;
262 Point hi = myUpper;
263 hi[ a ] = 0;
264 lo[ a ] = 0;
265 Domain D( lo, hi );
266 LatticeSetByInterval L;
267 for ( auto&& p : D )
268 {
269 auto I = intersectionIntervalAlongAxis( p, a );
270 if ( I.second > I.first )
271 L[ p ] = Interval( I.first, I.second - 1 );
272 }
273 return L;
274}
275
276//-----------------------------------------------------------------------------
277template <typename TSpace>
278typename DGtal::BoundedLatticePolytopeCounter<TSpace>::LatticeSetByInterval
279DGtal::BoundedLatticePolytopeCounter<TSpace>::
280getLatticeCells( Dimension a ) const
281{
282 ASSERT( myPolytope != nullptr );
283 Point lo = myLower;
284 Point hi = myUpper;
285 hi[ a ] = 0;
286 lo[ a ] = 0;
287 Domain D( lo, hi );
288 LatticeSetByInterval L; //< stores the intersected cells
289 const Point One = Point::diagonal( 1 );
290 Point q;
291 for ( auto&& p : D )
292 {
293 q = 2*p - One; q[ a ] = 0;
294 const auto I = intersectionIntervalAlongAxis( p, a );
295 const auto n = I.second - I.first;
296 if ( n != 0 )
297 {
298 // Now the second bound is included
299 L[ q ] = Interval( 2 * I.first - 1, 2 * I.second - 3 );
300 }
301 }
302 // It remains to compute all the k-cells, 0 <= k < d, intersected by Cvxh( Z )
303 for ( Dimension k = 0; k < dimension; k++ )
304 {
305 if ( k == a ) continue;
306 std::vector< Point > q_computed;
307 std::vector< Interval > I_computed;
308 for ( const auto& value : L )
309 {
310 Point p = value.first;
311 Interval I = value.second;
312 Point r = p; r[ k ] += 2;
313 const auto it = L.find( r );
314 if ( it == L.end() ) continue; // neighbor is empty
315 // Otherwise compute common part.
316 Interval J = it->second;
317 auto f = std::max( I.first, J.first );
318 auto s = std::min( I.second, J.second );
319 if ( f <= s )
320 {
321 Point qq = p; qq[ k ] += 1;
322 q_computed.push_back( qq );
323 I_computed.push_back( Interval( f, s ) );
324 }
325 }
326 // Add new columns to map Point -> column
327 for ( typename Point::Index i = 0; i < q_computed.size(); ++i )
328 {
329 L[ q_computed[ i ] ] = I_computed[ i ];
330 }
331 }
332 return L;
333}
334
335//-----------------------------------------------------------------------------
336template <typename TSpace>
337DGtal::Dimension
338DGtal::BoundedLatticePolytopeCounter<TSpace>::
339longestAxis( ) const
340{
341 ASSERT( myPolytope != nullptr );
342 Dimension b = 0;
343 auto b_size = myUpper[ 0 ] - myLower[ 0 ];
344 for ( Dimension a = 1; a < dimension; a++ )
345 {
346 const auto a_size = myUpper[ a ] - myLower[ a ];
347 if ( b_size < a_size ) { b = a; b_size = a_size; }
348 }
349 return b;
350}
351
352
353// //
354///////////////////////////////////////////////////////////////////////////////