DGtal  1.5.beta
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  */
49 inline
50 DGtal::MeasureOfStraightLines::MeasureOfStraightLines()
51 {
52  ///Default value
53  myEpsilon = 0.0005;
54 }
55 
56 
57 /**
58  * Destructor.
59  */
60 inline
61 DGtal::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  */
69 inline
70 void
71 DGtal::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  */
80 inline
81 bool
82 DGtal::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 ///
90 inline double
91 DGtal::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  */
123 inline void
124 DGtal::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  */
139 inline double
140 DGtal::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  */
163 inline double
164 DGtal::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  */
190 inline double
191 DGtal::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 ///
216 inline
217 double
218 DGtal::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 
228 inline
229 int
230 DGtal::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 ///
243 inline
244 double
245 DGtal::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 ///
283 inline
284 double
285 DGtal::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 ///
325 inline
326 double
327 DGtal::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 ///
354 inline
355 double
356 DGtal::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  */
390 inline
391 std::ostream&
392 DGtal::operator<<( std::ostream & out,
393  const DGtal::MeasureOfStraightLines & object )
394 {
395  object.selfDisplay( out );
396  return out;
397 }
398 
399 // //
400 ///////////////////////////////////////////////////////////////////////////////
401 
402