DGtal 1.4.0
Loading...
Searching...
No Matches
MeasureOfStraightLines.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 MeasureOfStraightLines.ih
19 * @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr )
20 * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
21 *
22 * @date 2010/03/04
23 *
24 * Implementation of inline methods defined in MeasureOfStraightLines.h
25 *
26 * This file is part of the DGtal library.
27 */
28
29///////////////////////////////////////////////////////////////////////////////
30// IMPLEMENTATION of inline methods.
31///////////////////////////////////////////////////////////////////////////////
32
33//////////////////////////////////////////////////////////////////////////////
34#include <cstdlib>
35#include <vector>
36#include <cmath>
37
38#include <boost/math/special_functions/asinh.hpp>
39
40//////////////////////////////////////////////////////////////////////////////
41
42///////////////////////////////////////////////////////////////////////////////
43// Internal functions //
44
45
46/**
47 * Constructor.
48 */
49inline
50DGtal::MeasureOfStraightLines::MeasureOfStraightLines()
51{
52 ///Default value
53 myEpsilon = 0.0005;
54}
55
56
57/**
58 * Destructor.
59 */
60inline
61DGtal::MeasureOfStraightLines:: ~MeasureOfStraightLines()
62{
63}
64
65/**
66 * Writes/Displays the object on an output stream.
67 * @param out the output stream where the object is written.
68 */
69inline
70void
71DGtal::MeasureOfStraightLines::selfDisplay( std::ostream & out ) const
72{
73 out << "[MeasureOfStraightLines]";
74}
75
76/**
77 * Checks the validity/consistency of the object.
78 * @return 'true' if the object is valid, 'false' otherwise.
79 */
80inline
81bool
82DGtal::MeasureOfStraightLines::isValid() const
83{
84 return true;
85}
86
87///
88/// Compute the measure associated to the Edge (a0,b0) -- (a1,b1) in the (a,b)-space
89///
90inline double
91DGtal::MeasureOfStraightLines::computeMeasureEdge ( double a0,double b0, double a1, double b1 )
92{
93 //double measure=0;
94
95 //Aligned with the 'b' axis
96 if ( ( b0 == 0 ) && ( b1 == 0 ) )
97 return 0;
98 //Aligned with the 'b' axis
99 if ( ( a0 == 0 ) && ( a1 == 0 ) )
100 return 0;
101
102 double delta = ( a0*b1 - a1*b0 );
103
104 if ( a0 == 0 )
105 return delta * ( 1.0/ ( 1.0 + sqrt ( 1 + pow ( a1,2 ) ) ) );
106 if ( a1 == 0 )
107 return delta * ( ( -1.0 + sqrt ( 1.0 + pow ( a0,2 ) ) ) /pow ( a0,2 ) );
108
109 if ( a0 != a1 )
110 {
111 return delta * ( ( a0 - a1 + sqrt ( 1 + pow ( a0,2 ) ) *a1 - a0*sqrt ( 1 + pow ( a1,2 ) ) ) / ( pow ( a0,2 ) *a1 - a0*pow ( a1,2 ) ) );
112 }
113 else
114 return delta * a1/ ( a0* ( 1 + pow ( a1,2 ) + sqrt ( 1 + pow ( a1,2 ) ) ) );
115}
116
117
118/**
119 * Set the epsilon threshold in the numerical approximation.
120 *
121 * @param aValue the new epsilon value
122 */
123inline void
124DGtal::MeasureOfStraightLines::setEpsilon(const double aValue)
125{
126 myEpsilon = aValue;
127}
128
129///////////////////////////////////////////////////////////////////////////////
130// Implementation of inline methods //
131
132/**
133 * Compute the measure of the polygon {(a_i,b_i)} in the (a,b)-parameter space
134 * @param a the a-value of polygon vertices
135 * @param b the b-value of polygon vertices
136 *
137 * @return the measure value (positive value)
138 */
139inline double
140DGtal::MeasureOfStraightLines::computeMeasure(const std::vector<double> &a,
141 const std::vector<double> &b)
142{
143 double measure = 0;
144
145 ASSERT(a.size() == b.size());
146
147 for (unsigned int i=0 ; i < a.size() ; i++)
148 {
149 measure += computeMeasureEdge ( a[i], b[i],
150 a[ ( i+1 ) % a.size()],b[ ( i+1 ) %a.size()] );
151 }
152
153 return measure;
154}
155
156/**
157 * Compute the abscissa of the centroid of the polygon {(a_i,b_i)} in the (a,b)-parameter space
158 * @param a the a-value of polygon vertices
159 * @param b the b-value of polygon vertices
160 *
161 * @return the measure value (positive value)
162 */
163inline double
164DGtal::MeasureOfStraightLines::computeCentroidA(const std::vector<double> &a,
165 const std::vector<double> &b)
166{
167 double measure = 0;
168 double C_a = 0;
169
170 ASSERT(a.size() == b.size());
171
172 for (unsigned int i=0 ; i < a.size() ; i++)
173 {
174 measure += computeMeasureEdge( a[i], b[i],
175 a[ ( i+1 ) % a.size()],b[ ( i+1 ) %a.size()] );
176 C_a += computeCentroidEdge_a ( a[i], b[i],
177 a[ ( i+1 ) % a.size()],b[ ( i+1 ) %a.size()] );
178 }
179
180 return C_a/measure;
181}
182
183/**
184 * Compute the abscissa of the centroid of the polygon {(a_i,b_i)} in the (a,b)-parameter space
185 * @param a the a-value of polygon vertices
186 * @param b the b-value of polygon vertices
187 *
188 * @return the measure value (positive value)
189 */
190inline double
191DGtal::MeasureOfStraightLines::computeCentroidB(const std::vector<double> &a,
192 const std::vector<double> &b)
193{
194 double measure = 0;
195 double C_b = 0;
196
197 ASSERT(a.size() == b.size());
198
199 for (unsigned int i=0 ; i < a.size() ; i++)
200 {
201 measure += computeMeasureEdge( a[i], b[i],
202 a[ ( i+1 ) % a.size()],b[ ( i+1 ) %a.size()] );
203 C_b += computeCentroidEdge_b ( a[i], b[i],
204 a[ ( i+1 ) % a.size()],b[ ( i+1 ) %a.size()] );
205 }
206
207 return C_b/measure;
208}
209
210
211///
212/// Compute the centroid on 'b' on the rectangular domain with vertices (x1,,y1) - (x2,y2)
213/// PRECONDITION: y1<y2
214///
215///
216inline
217double
218DGtal::MeasureOfStraightLines::__computeCentroidSquare_b ( double x1, double y1, double x2,double y2 )
219{
220 double val;
221 val = ( ( -1.0* ( sqrt ( 1.0 + pow ( x1,2 ) ) *x2 ) + x1*sqrt ( 1 + pow ( x2,2 ) ) ) * ( pow ( y1,2 ) - pow ( y2,2 ) ) ) / ( 2.0*sqrt ( ( 1.0 + pow ( x1,2 ) ) * ( 1.0 + pow ( x2,2 ) ) ) );
222
223 ASSERT ( val>=0 );
224 return val;
225}
226
227
228inline
229int
230DGtal::MeasureOfStraightLines::sign ( const double a )
231{
232 if ( a>=0 )
233 return 1;
234 else
235 return -1;
236}
237
238
239///
240/// Approximate the centroid on 'b' on the trapezioid (a0,0)-(a0,b0)-(a1,b1)-(a1,0)
241/// (internal function)
242///
243inline
244double
245DGtal::MeasureOfStraightLines::__computeCentroidEdgeApprox_b ( double a0, double b0,double a1,double b1 )
246{
247 double measure,y2;
248 unsigned int nb_step;
249 /// A0 -> A1
250 if ( a1 < a0 )
251 {
252 double tmp = a0;
253 a0 = a1;
254 a1 = tmp;
255 tmp = b0;
256 b0 = b1;
257 b1 = tmp;
258 }
259
260 measure = 0;
261 nb_step = ( unsigned int ) floor ( fabs ( a1-a0 ) /myEpsilon );
262 double slope = ( b1-b0 ) / ( a1-a0 );
263 double decal = b1 - ( b1-b0 ) / ( a1-a0 ) *a1;
264
265 if ( slope == 0 )
266 return __computeCentroidSquare_b ( a0, 0, a1, b0 );
267
268 for ( unsigned int i=0; i < nb_step ; i++ )
269 {
270 y2 = ( b1-b0 ) / ( a1-a0 ) * ( a0+ ( i+1 ) *myEpsilon ) + decal;
271 if ( y2>0 )
272 measure += __computeCentroidSquare_b ( a0+i*myEpsilon, 0, a0+ ( i+1 ) *myEpsilon, y2 );
273 else
274 measure += __computeCentroidSquare_b ( a0+i*myEpsilon, y2, a0+ ( i+1 ) *myEpsilon, 0 );
275 }
276 return measure;
277}
278
279///
280/// Approximate the centroid on 'b' on the triangle (0,0)-(a0,b0)-(a1,b1)
281/// (internal function)
282///
283inline
284double
285DGtal::MeasureOfStraightLines::__computeCentroidTriApprox_b ( double a0, double b0,double a1,double b1 )
286{
287 int signe = sign ( a0*b1 - a1*b0 );
288 double measure = 0;
289
290 double m_A0A1=0, m_0A0=0, m_0A1=0;
291 if ( ( b0 != 0 ) && ( a0 != 0 ) )
292 m_0A0 = __computeCentroidEdgeApprox_b ( 0,0,a0,b0 );
293 if ( ( b1 != 0 ) && ( a1 != 0 ) )
294 m_0A1 = __computeCentroidEdgeApprox_b ( 0,0,a1,b1 );
295 if ( ( a1 != a0 ) )
296 m_A0A1 = __computeCentroidEdgeApprox_b ( a0,b0,a1,b1 );
297
298 ASSERT ( m_0A0>=0 );
299 ASSERT ( m_0A1>=0 );
300 ASSERT ( m_A0A1>=0 );
301
302 if ( a0 < a1 )
303 {
304 double det = ( a1 * ( b0 - b1 ) - b1* ( a0 - a1 ) );
305 if ( det >0 )
306 measure = m_0A0 + m_A0A1 - m_0A1;
307 else
308 measure = m_0A1 - m_0A0 - m_A0A1;
309 }
310 else
311 {
312 double det = ( a0 * ( b1 - b0 ) - b0* ( a1 - a0 ) );
313 if ( det >0 )
314 measure = m_0A1 + m_A0A1 - m_0A0;
315 else
316 measure = m_0A0 - m_0A1 - m_A0A1;
317 }
318 ASSERT ( measure>=0 );
319 return signe*measure;
320}
321
322///
323/// Compute the centroid on 'b' on the triangle (0,0)-(a0,b0)-(a1,b1)
324///
325inline
326double
327DGtal::MeasureOfStraightLines::computeCentroidEdge_b ( double a0,double b0, double a1, double b1 )
328{
329 //double measure=0;
330 double delta= ( a0*b1 - a1*b0 );
331
332 if ( ( b0 == 0 ) && ( b1 == 0 ) )
333 return 0;
334 if ( ( a0 == 0 ) && ( a1 == 0 ) )
335 return 0;
336
337 if ( a0 == 0 )
338 return delta* ( a1* ( ( -2 + sqrt ( 1 + a1*a1 ) ) *b0 + 2*b1 ) + ( b0 - 2*b1 ) *boost::math::asinh ( a1 ) ) / ( 2*a1*a1*a1 );
339 if ( a1 == 0 )
340 return delta* ( a0* ( 2*b0 + ( -2 + sqrt ( 1 + a0*a0 ) ) *b1 ) + ( -2*b0 + b1 ) *boost::math::asinh ( a0 ) ) / ( 2*a0*a0*a0 );
341
342 if ( a0 == a1 )
343 return delta* ( ( - ( ( a1* ( b0 + a1* ( -a1 + sqrt ( 1 + pow ( a1,-2 ) ) *sqrt ( 1 + pow ( a1,2 ) ) ) *
344 b1 ) ) /sqrt ( 1.0 + pow ( a1,2 ) ) ) + ( b0 + b1 ) *boost::math::asinh ( a1 ) ) /
345 ( 2.*pow ( a1,3 ) ) );
346
347 return __computeCentroidTriApprox_b ( a0,b0,a1,b1 );
348
349}
350
351///
352/// Compute the centroid on 'a' on the triangle (0,0)-(a0,b0)-(a1,b1)
353///
354inline
355double
356DGtal::MeasureOfStraightLines::computeCentroidEdge_a ( double a0,double b0, double a1, double b1 )
357{
358 //double measure=0;
359 double delta= ( a0*b1 - a1*b0 );
360
361 if ( ( b0 == 0 ) && ( b1 == 0 ) )
362 return 0;
363 if ( ( a0 == 0 ) && ( a1 == 0 ) )
364 return 0;
365
366 if ( a0 == 0 )
367 return delta* ( a1 - boost::math::asinh ( a1 ) ) / ( a1*a1 );
368 if ( a1 == 0 )
369 return delta* ( a0 - boost::math::asinh ( a0 ) ) / ( a0*a0 );
370
371 if ( a0 == a1 )
372 return delta* ( ( -sqrt ( 1 + pow ( a1,-2 ) ) - 1.0/ ( a1*sqrt ( 1 + pow ( a1,2 ) ) ) + a1/sqrt ( 1 + pow ( a1,2 ) ) +
373 ( 2*boost::math::asinh ( a1 ) ) /pow ( a1,2 ) ) /2. );
374
375 return delta* ( ( - ( a1*boost::math::asinh ( a0 ) ) + a0*boost::math::asinh ( a1 ) ) / ( a0* ( a0 - a1 ) *a1 ) );
376}
377
378
379
380
381///////////////////////////////////////////////////////////////////////////////
382// Implementation of inline functions and external operators //
383
384/**
385 * Overloads 'operator<<' for displaying objects of class 'MeasureOfStraightLines'.
386 * @param out the output stream where the object is written.
387 * @param object the object of class 'MeasureOfStraightLines' to write.
388 * @return the output stream after the writing.
389 */
390inline
391std::ostream&
392DGtal::operator<<( std::ostream & out,
393 const DGtal::MeasureOfStraightLines & object )
394{
395 object.selfDisplay( out );
396 return out;
397}
398
399// //
400///////////////////////////////////////////////////////////////////////////////
401
402