DGtal 1.4.2
Loading...
Searching...
No Matches
PlaneProbingLNeighborhood.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
19 * @author Tristan Roussillon (\c tristan.roussillon@liris.cnrs.fr )
20 * Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
21 *
22 * @date 2024/09/16
23 *
24 * Implementation of inline methods defined in PlaneProbingLNeighborhood.h
25 *
26 * This file is part of the DGtal library.
27 */
28
29
30//////////////////////////////////////////////////////////////////////////////
31#include <cstdlib>
32//////////////////////////////////////////////////////////////////////////////
33
34
35///////////////////////////////////////////////////////////////////////////////
36// IMPLEMENTATION of inline methods.
37///////////////////////////////////////////////////////////////////////////////
38
39///////////////////////////////////////////////////////////////////////////////
40// ----------------------- Standard services ------------------------------
41
42// ------------------------------------------------------------------------
43template < typename TPredicate >
44inline
45DGtal::PlaneProbingLNeighborhood<TPredicate>::
46PlaneProbingLNeighborhood(Predicate const& aPredicate, Point const& aQ, Triangle const& aM)
47 : DGtal::PlaneProbingRNeighborhood<TPredicate>(aPredicate, aQ, aM)
48{
49 for (int k = 0; k < 3; k++) {
50 myGrids.push_back(closestInGrid(k));
51 }
52}
53
54// ------------------------------------------------------------------------
55template < typename TPredicate >
56inline
57DGtal::PlaneProbingLNeighborhood<TPredicate>::~PlaneProbingLNeighborhood()
58{}
59
60///////////////////////////////////////////////////////////////////////////////
61// ----------------------- Plane Probing services ------------------------------
62
63// ------------------------------------------------------------------------
64template < typename TPredicate >
65inline
66typename DGtal::PlaneProbingLNeighborhood<TPredicate>::HexagonState
67DGtal::PlaneProbingLNeighborhood<TPredicate>::hexagonState ()
68{
69 for (int k = 0; k < 3; k++) {
70 updateGrid(k);
71 }
72
73 std::array<bool, 6> state({ false, false, false, false, false, false });
74 for (int k = 0; k < 3; k++) {
75 state[2*k] = myGrids[k].myPair.first;
76 state[2*k+1] = myGrids[k].myPair.second;
77 }
78
79 return this->classify(state);
80}
81
82// ------------------------------------------------------------------------
83template < typename TPredicate >
84inline
85typename DGtal::PlaneProbingLNeighborhood<TPredicate>::UpdateOperation
86DGtal::PlaneProbingLNeighborhood<TPredicate>::
87closestCandidate ()
88{
89
90 std::vector<GridPoint> validGridPoints;
91 for (int k = 0; k < 3; k++) {
92 GridPoint gp = myGrids[k].myGridPoint;
93 if (gp.isValid())
94 validGridPoints.push_back(gp);
95 }
96
97 if (validGridPoints.size() == 1) {
98
99 return getOperationFromGridPoint(validGridPoints[0]);
100
101 } else {
102 // One should call hexagonState before closestCandidate, and check the return value
103 // to ensure that there is at least one point in the plane in the H-neighbhorhood
104 ASSERT(validGridPoints.size() == 2);
105
106 if (this->isSmallest(this->relativePoint(validGridPoints[0]),
107 this->relativePoint(validGridPoints[1]))) {
108 return getOperationFromGridPoint(validGridPoints[1]);
109 } else {
110 return getOperationFromGridPoint(validGridPoints[0]);
111 }
112 }
113
114}
115
116///////////////////////////////////////////////////////////////////////////////
117// ----------------------- Helpers to find a closest point ----------------
118
119// ------------------------------------------------------------------------
120template < typename TPredicate >
121inline
122typename DGtal::PlaneProbingLNeighborhood<TPredicate>::ClosestGridPoint
123DGtal::PlaneProbingLNeighborhood<TPredicate>::
124closestInGrid (const Index& aIdx) const
125{
126 std::vector<GridPoint> candidates;
127
128 Index k = aIdx;
129 GridPoint y1(1,0,k); //y1 = vk + m1
130 GridPoint y2(0,1,k); //y2 = vk + m2
131
132 if (this->myPredicate(this->absolutePoint(y1))) {
133 if (this->myPredicate(this->absolutePoint(y2))) {
134 //at least 2 candidates
135 if (this->isSmallest(this->relativePoint(y1), this->relativePoint(y2))) {
136 candidates.push_back(y2);
137 } else {
138 candidates.push_back(y1);
139 }
140
141 ASSERT(candidates.size() == 1);
142 candidatesInGrid(y1, y2, std::back_inserter(candidates));
143 GridPoint closest = this->closestPointInList(candidates);
144 return ClosestGridPoint( closest, true, true );
145
146 } else {
147 //1 candidate
148 return ClosestGridPoint( y1, true, false );
149 }
150 } else {
151 if (this->myPredicate(this->absolutePoint(y2))) {
152 //1 candidate
153 return ClosestGridPoint( y2, false, true );
154 } else {
155 //no candidate
156 return ClosestGridPoint( GridPoint(0, 0, k), false, false );
157 }
158 }
159}
160
161// ------------------------------------------------------------------------
162template < typename TPredicate >
163inline
164void
165DGtal::PlaneProbingLNeighborhood<TPredicate>::
166candidatesInGrid (const GridPoint& y1, const GridPoint& y2,
167 std::back_insert_iterator<std::vector<GridPoint> > out) const
168{
169 using DGtal::detail::squaredNorm;
170
171 ASSERT( (this->myPredicate(this->absolutePoint(y1)) &&
172 (this->myPredicate(this->absolutePoint(y2)))) );
173
174 Integer a = direction(y1).dot(direction(y2));
175 if (a < 0) {
176
177 GridPoint y = y1 + y2;
178 Integer a1 = direction(y1).dot(direction(y));
179 Integer a2 = direction(y2).dot(direction(y));
180 if ( (a1 < 0)||(a2 < 0) ) {
181
182 //if a2 < 0
183 GridPoint u = y1;
184 GridPoint w = y2;
185 if (a1 < 0) {
186 std::swap(u, w);
187 }
188 ASSERT(squaredNorm(direction(u)) > squaredNorm(direction(w)));
189
190 Integer gamma = (-a)/squaredNorm(direction(w));
191 GridPoint w2 = u + w*gamma;
192 GridPoint w2Next = u + w*(gamma+1);
193
194 if (direction(w2).dot(direction(w2Next)) < 0) {
195
196 if (this->myPredicate(this->absolutePoint(w2))) {
197 candidatesInGrid(w, w2, out);
198 } else {
199 //we look for a closest point on the ray u +cw, c<gamma
200 GridPointOnProbingRay gr = closestOnBoundedRayLogWithPredicate(GridPointOnProbingRay(u,w.direction()),gamma);
201 ASSERT( gr == closestOnBoundedRayLinearWithPredicate(GridPointOnProbingRay(u,w.direction()),gamma) );
202 *out++ = gr.gridPoint();
203 }
204 } else {
205 //excepted the first one, the other angles along
206 // the ray are all acute (or right),
207 ASSERT( direction(w2).dot(direction(w2Next)) >= 0 );
208 ASSERT( direction(w).dot(direction(w2Next)) >= 0 );
209
210 //whatever w2Next is in plane or not,
211 //we look for a closest point on the ray u +cw, c<=gamma+1
212 GridPointOnProbingRay gr = closestOnBoundedRayLogWithPredicate(GridPointOnProbingRay(u,w.direction()),(gamma+2));
213 ASSERT( gr == closestOnBoundedRayLinearWithPredicate(GridPointOnProbingRay(u,w.direction()),(gamma+2)) );
214 *out++ = gr.gridPoint();
215 }
216
217 } else { //a1 and a2 are both acute (or right),
218 //only y has yet to be considered (in addition to y1, y2)
219 //(because the angle between y1 and y2 is obtuse)
220 if (this->myPredicate(this->absolutePoint(y))) {
221 *out++ = y;
222 }
223 }
224 } //if a >= 0, we have an acute or right angle
225 //and no other point has to be considered.
226 //(if right angle, then case of cospherical points)
227}
228
229// ------------------------------------------------------------------------
230template < typename TPredicate >
231inline
232typename DGtal::PlaneProbingLNeighborhood<TPredicate>::GridPointOnProbingRay
233DGtal::PlaneProbingLNeighborhood<TPredicate>::closestOnBoundedRayLogWithPredicate (GridPointOnProbingRay const& aRay, Integer const& aBound) const
234{
235 GridPointOnProbingRay Xk = aRay, Xl = aRay.next(aRay.index()+aBound);
236 ASSERT(this->myPredicate(this->absolutePoint(Xk)));
237
238 // Binary search
239 Integer d = Xl.index() - Xk.index();
240 while (d > 4) {
241 ASSERT(this->myPredicate(this->absolutePoint(Xk)));
242
243 GridPointOnProbingRay Xalpha = Xk.next(Integer(d / 4)),
244 Xbeta = Xk.next(Integer(d / 2)),
245 Xgamma = Xk.next(Integer(3*d/4));
246
247 ASSERT(Xk.index() < Xalpha.index() && Xalpha.index() < Xbeta.index() &&
248 Xbeta.index() < Xgamma.index() && Xgamma.index() < Xl.index());
249
250 if (this->myPredicate(this->absolutePoint(Xbeta)) &&
251 this->isSmallest(this->relativePoint(Xbeta), this->relativePoint(Xgamma))) {
252 Xk = Xbeta;
253 } else if (! this->myPredicate(this->absolutePoint(Xalpha)) ||
254 this->isSmallest(this->relativePoint(Xbeta), this->relativePoint(Xalpha))) {
255 Xl = Xbeta;
256 } else {
257 Xk = Xalpha;
258 Xl = Xgamma;
259 }
260
261 d = Xl.index() - Xk.index();
262 }
263
264 return closestOnBoundedRayLinearWithPredicate(Xk, d);
265}
266
267// ------------------------------------------------------------------------
268template < typename TPredicate >
269inline
270typename DGtal::PlaneProbingLNeighborhood<TPredicate>::GridPointOnProbingRay
271DGtal::PlaneProbingLNeighborhood<TPredicate>::closestOnBoundedRayLinearWithPredicate (GridPointOnProbingRay const& aRay, const Integer& aBound) const
272{
273 ASSERT(this->myPredicate(this->absolutePoint(aRay)));
274
275 Integer counter = 0;
276 GridPointOnProbingRay previousX = aRay, currentX = previousX.next(1);
277 while (this->myPredicate(this->absolutePoint(currentX)) &&
278 this->isSmallest(this->relativePoint(previousX), this->relativePoint(currentX)) &&
279 counter < aBound) {
280 previousX = currentX;
281 currentX = previousX.next(1);
282 }
283
284 ASSERT(this->myPredicate(this->absolutePoint(previousX)));
285
286 return previousX;
287}
288
289// ------------------------------------------------------------------------
290template < typename TPredicate >
291inline
292typename DGtal::PlaneProbingLNeighborhood<TPredicate>::UpdateOperation
293DGtal::PlaneProbingLNeighborhood<TPredicate>::
294getOperationFromGridPoint (GridPoint const& aP) const
295{
296 auto k = aP.k();
297 auto d = aP.direction();
298 return { { k, (k+1)%3, (k+2)%3 }, //permutation
299 { 1, -d.first, -d.second } //coefficients
300 };
301}
302
303// ------------------------------------------------------------------------
304template < typename TPredicate >
305inline
306void
307DGtal::PlaneProbingLNeighborhood<TPredicate>::
308updateGrid (const Index& aIdx)
309{
310 //Let r1 and r2 be the two rays bound to the vertex of index 'aIdx'
311 PointOnProbingRay r1 = PointOnProbingRay({aIdx, (aIdx+1)%3, (aIdx+2)%3});
312 PointOnProbingRay r2 = PointOnProbingRay({aIdx, (aIdx+2)%3, (aIdx+1)%3});
313 //Check: do they are in the allowed rays?
314 if (this->isNeighbor(r1)) {
315 if (this->isNeighbor(r2)) {
316 //if both,
317 myGrids[aIdx] = closestInGrid(aIdx);
318 } else {
319 //if only r1,
320 myGrids[aIdx] = ClosestGridPoint( GridPoint(1, 0, aIdx), true, false );
321 }
322 } else {
323 if (this->isNeighbor(r2)) {
324 //if only r2,
325 myGrids[aIdx] = ClosestGridPoint( GridPoint(0, 1, aIdx), false, true );
326 } else {
327 //if none of them, no candidate
328 myGrids[aIdx] = ClosestGridPoint( GridPoint(0, 0, aIdx), false, false );
329 }
330 }
331}
332
333// ------------------------------------------------------------------------
334template < typename TPredicate >
335inline
336typename DGtal::PlaneProbingLNeighborhood<TPredicate>::Point
337DGtal::PlaneProbingLNeighborhood<TPredicate>::
338direction (GridPoint const& aP) const
339{
340 return aP.directionVector(this->myM);
341}
342
343///////////////////////////////////////////////////////////////////////////////
344// Interface - public :
345
346/**
347 * Writes/Displays the object on an output stream.
348 * @param out the output stream where the object is written.
349 */
350template <typename TPredicate>
351inline
352void
353DGtal::PlaneProbingLNeighborhood<TPredicate>::selfDisplay ( std::ostream & out ) const
354{
355 out << "[PlaneProbingLNeighborhood]";
356}
357
358/**
359 * Checks the validity/consistency of the object.
360 * @return 'true' if the object is valid, 'false' otherwise.
361 */
362template <typename TPredicate>
363inline bool DGtal::PlaneProbingLNeighborhood<TPredicate>::isValid() const
364{
365 return true;
366}
367
368
369
370///////////////////////////////////////////////////////////////////////////////
371// Implementation of inline functions //
372
373template <typename TPredicate>
374inline
375std::ostream&
376DGtal::operator<< ( std::ostream & out,
377 const DGtal::PlaneProbingLNeighborhood<TPredicate> & object )
378{
379 object.selfDisplay( out );
380 return out;
381}
382
383// //
384///////////////////////////////////////////////////////////////////////////////
385
386