DGtal  1.5.beta
CircleFrom3Points.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 CircleFrom3Points.ih
19  * @author Tristan Roussillon (\c tristan.roussillon@liris.cnrs.fr )
20  * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
21  *
22  * @date 2011/09/22
23  *
24  * @brief Implementation of inline methods defined in CircleFrom3Points.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 template <typename TPoint>
42 inline
43 DGtal::CircleFrom3Points<TPoint>::~CircleFrom3Points()
44 {
45 }
46 
47 template <typename TPoint>
48 inline
49 DGtal::CircleFrom3Points<TPoint>::CircleFrom3Points()
50 {
51 }
52 
53 template <typename TPoint>
54 inline
55 void
56 DGtal::CircleFrom3Points<TPoint>::init(
57  const Point& aFirstPoint,
58  const Point& aSecondPoint,
59  const Point& aThirdPoint)
60 {
61  myP = aFirstPoint;
62  myQ = aSecondPoint;
63  myR = aThirdPoint;
64 }
65 
66 
67 template <typename TPoint>
68 inline
69 DGtal::CircleFrom3Points<TPoint>::CircleFrom3Points(
70  const Point& aFirstPoint,
71  const Point& aSecondPoint,
72  const Point& aThirdPoint):
73  myP(aFirstPoint),
74  myQ(aSecondPoint),
75  myR(aThirdPoint)
76 {
77 }
78 
79 
80 template <typename TPoint>
81 inline
82 DGtal::CircleFrom3Points<TPoint>::CircleFrom3Points(
83  const CircleFrom3Points & other):
84  myP(other.myP),
85  myQ(other.myQ),
86  myR(other.myR)
87 {
88 }
89 
90 
91 template <typename TPoint>
92 inline
93 DGtal::CircleFrom3Points<TPoint>&
94 DGtal::CircleFrom3Points<TPoint>::operator=(
95  const CircleFrom3Points & other)
96 {
97  myP = other.myP;
98  myQ = other.myQ;
99  myR = other.myR;
100  return *this;
101 }
102 
103 
104 
105 template <typename TPoint>
106 inline
107 typename DGtal::CircleFrom3Points<TPoint>::Distance
108 DGtal::CircleFrom3Points<TPoint>::signedDistance(const Point& aP) const
109 {
110  Vector u( (myP[0]-aP[0])*(myR[1]-aP[1])-(myR[0]-aP[0])*(myP[1]-aP[1]),
111  (myR[0]-aP[0])*(myR[0]-myP[0])+(myR[1]-aP[1])*(myR[1]-myP[1]) );
112  Vector v( (myP[0]-aP[0])*(myQ[1]-aP[1])-(myQ[0]-aP[0])*(myP[1]-aP[1]),
113  (myQ[0]-aP[0])*(myQ[0]-myP[0])+(myQ[1]-aP[1])*(myQ[1]-myP[1]) );
114 
115  return -( (u[0] * v[1]) - (u[1] * v[0]) );
116 }
117 
118 template <typename TPoint>
119 inline
120 void
121 DGtal::CircleFrom3Points<TPoint>::getParameters(double& cx, double& cy, double& rr) const
122 {
123  Vector u(myP - myQ);
124  Vector v(myR - myP);
125 
126  Coordinate area = (u[0] * v[1]) - (u[1] * v[0]);
127 
128  if (area == 0) {
129  std::cerr << "Infinite radius detected in method getParameters() of CircleFrom3Points" << std::endl;
130  throw InfiniteNumberException();
131  } else {
132 
133  //center
134  double areaD = NumberTraits<Coordinate>::castToDouble(area)*2.0;
135  double pz, qz, rz;
136  pz = std::pow( NumberTraits<Coordinate>::castToDouble(myP[0]),2.0 )/areaD
137  + std::pow( NumberTraits<Coordinate>::castToDouble(myP[1]),2.0 )/areaD;
138  qz = std::pow( NumberTraits<Coordinate>::castToDouble(myQ[0]),2.0 )/areaD
139  + std::pow( NumberTraits<Coordinate>::castToDouble(myQ[1]),2.0 )/areaD;
140  rz = std::pow( NumberTraits<Coordinate>::castToDouble(myR[0]),2.0 )/areaD
141  + std::pow( NumberTraits<Coordinate>::castToDouble(myR[1]),2.0 )/areaD;
142 
143  cx = NumberTraits<Coordinate>::castToDouble(myP[1])*(qz-rz)
144  -NumberTraits<Coordinate>::castToDouble(myQ[1])*(pz-rz)
145  +NumberTraits<Coordinate>::castToDouble(myR[1])*(pz-qz);
146  cy = -NumberTraits<Coordinate>::castToDouble(myP[0])*(qz-rz)
147  +NumberTraits<Coordinate>::castToDouble(myQ[0])*(pz-rz)
148  -NumberTraits<Coordinate>::castToDouble(myR[0])*(pz-qz);
149 
150  //radius
151  double a, b, c;
152  a = std::pow( NumberTraits<Coordinate>::castToDouble(myQ[0]-myP[0]),2.0 )
153  +std::pow( NumberTraits<Coordinate>::castToDouble(myQ[1]-myP[1]),2.0 );
154  b = std::pow( NumberTraits<Coordinate>::castToDouble(myR[0]-myQ[0]),2.0 )
155  +std::pow( NumberTraits<Coordinate>::castToDouble(myR[1]-myQ[1]),2.0 );
156  c = std::pow( NumberTraits<Coordinate>::castToDouble(myR[0]-myP[0]),2.0 )
157  +std::pow( NumberTraits<Coordinate>::castToDouble(myR[1]-myP[1]),2.0 );
158 
159  rr = ( ( std::sqrt(a)/std::abs(areaD) )*std::sqrt(b)*std::sqrt(c) );
160 
161  }
162 }
163 
164 template <typename TPoint>
165 inline
166 double
167 DGtal::CircleFrom3Points<TPoint>::getCurvature() const
168 {
169  Vector u(myP - myQ);
170  Vector v(myR - myP);
171 
172  //4 times the area of the triangle
173  Coordinate area = (u[0] * v[1]) - (u[1] * v[0]);
174  double areaD = NumberTraits<Coordinate>::castToDouble(area)*2.0;
175 
176  //length of the sides of the triangle
177  double a, b, c;
178  a = std::pow( NumberTraits<Coordinate>::castToDouble(myQ[0]-myP[0]),2.0 )
179  +std::pow( NumberTraits<Coordinate>::castToDouble(myQ[1]-myP[1]),2.0 );
180  b = std::pow( NumberTraits<Coordinate>::castToDouble(myR[0]-myQ[0]),2.0 )
181  +std::pow( NumberTraits<Coordinate>::castToDouble(myR[1]-myQ[1]),2.0 );
182  c = std::pow( NumberTraits<Coordinate>::castToDouble(myR[0]-myP[0]),2.0 )
183  +std::pow( NumberTraits<Coordinate>::castToDouble(myR[1]-myP[1]),2.0 );
184 
185  //curvature
186  return -( areaD/ (std::sqrt(a)*std::sqrt(b)*std::sqrt(c)) );
187 }
188 
189 
190 ///////////////////////////////////////////////////////////////////////////////
191 // Interface - public :
192 
193 template <typename TPoint>
194 inline
195 std::string
196 DGtal::CircleFrom3Points<TPoint>::className() const
197 {
198  return "CircleFrom3Points";
199 }
200 
201 template <typename TPoint>
202 inline
203 void
204 DGtal::CircleFrom3Points<TPoint>::selfDisplay ( std::ostream & out ) const
205 {
206  out << "[CircleFrom3Points] passing through:\n";
207  out << myP << myQ << myR;
208 }
209 
210 template <typename TPoint>
211 inline
212 bool
213 DGtal::CircleFrom3Points<TPoint>::isValid() const
214 {
215  return true;
216 }
217 
218