DGtal 1.4.0
Loading...
Searching...
No Matches
SphericalAccumulator.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 SphericalAccumulator.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 2012/09/17
23 *
24 * Implementation of inline methods defined in SphericalAccumulator.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/**
42 * Constructor.
43 */
44template <typename T>
45inline
46DGtal::SphericalAccumulator<T>::SphericalAccumulator(const Size aP)
47{
48 myTotal = 0;
49 myNphi = aP;
50 myNtheta = 2*myNphi;
51 myAccumulator.resize(myNphi*myNtheta,0);
52 myAccumulatorDir.resize(myNphi*myNtheta, Vector::zero);
53 myBinNumber = 0;
54
55 myMaxBinTheta = 0;
56 myMaxBinPhi= 0;
57
58
59 for(Size posPhi=0; posPhi < myNphi; posPhi++)
60 for(Size posTheta=0; posTheta < myNtheta; posTheta++)
61 {
62 double dphi = M_PI/((double)myNphi-1);
63 double Ntheta_i = floor(2.0*((double)myNphi)*sin((double)posPhi*dphi)) ;
64
65 if (((posPhi == 0) || (posPhi == (myNphi-1))) && (posTheta==0))
66 myBinNumber++;
67 else
68 if ((posPhi < myNphi) && (posTheta<Ntheta_i) && (posTheta< myNtheta))
69 myBinNumber ++;
70 }
71}
72/**
73 * Destructor.
74 */
75template <typename T>
76inline
77DGtal::SphericalAccumulator<T>::~SphericalAccumulator()
78{
79}
80// --------------------------------------------------------
81template <typename T>
82inline
83void DGtal::SphericalAccumulator<T>::binCoordinates(const Vector &aDir,
84 Size &posPhi,
85 Size &posTheta) const
86{
87 double theta,phi, theta2;
88 double norm = aDir.norm();
89
90 ASSERT(norm != 0);
91
92 phi = acos(NumberTraits<typename T::Component>::castToDouble(aDir[2])/norm);
93
94 double dphi = M_PI/(double)(myNphi-1);
95 double Nthetai;
96 posPhi = static_cast<Size>(floor( (phi+dphi/2.) *(myNphi-1)/ M_PI));
97 if(posPhi == 0 || posPhi== (myNphi-1))
98 {
99 posTheta =0;
100 }
101 else
102 {
103 theta2=atan2(NumberTraits<typename T::Component>::castToDouble(aDir[1]),
104 NumberTraits<typename T::Component>::castToDouble(aDir[0]));
105
106 if(theta2<0)
107 theta = theta2 + 2.0*M_PI;
108 else
109 theta = theta2;
110 Nthetai = floor(2.0*(myNphi)*sin(posPhi*dphi));
111 double dtheta = 2.0*M_PI/(Nthetai);
112 posTheta = static_cast<Size>(floor( (theta+dtheta/2.0)/dtheta));
113
114 if (posTheta >= Nthetai)
115 posTheta -= Nthetai;
116 }
117
118 ASSERT(posPhi < myNphi);
119 ASSERT( isValidBin(posPhi,posTheta) );
120}
121// --------------------------------------------------------
122template <typename T>
123inline
124void DGtal::SphericalAccumulator<T>::addDirection(const Vector &aDir)
125{
126 Size posPhi,posTheta;
127
128 binCoordinates(aDir , posPhi, posTheta);
129 myAccumulator[posTheta + posPhi*myNtheta] += 1;
130 myAccumulatorDir[posTheta + posPhi*myNtheta] += aDir;
131 myTotal ++;
132
133 //Max bin update
134 if ( myAccumulator[posTheta + posPhi*myNtheta] >
135 myAccumulator[ myMaxBinTheta + myMaxBinPhi*myNtheta])
136 {
137 myMaxBinTheta = posTheta;
138 myMaxBinPhi = posPhi;
139 }
140}
141// --------------------------------------------------------
142template <typename T>
143inline
144typename DGtal::SphericalAccumulator<T>::Quantity
145DGtal::SphericalAccumulator<T>::samples() const
146{
147 return myTotal;
148}
149// --------------------------------------------------------
150template <typename T>
151inline
152void
153DGtal::SphericalAccumulator<T>::maxCountBin(Size &phi, Size &theta) const
154{
155 phi = myMaxBinPhi;
156 theta = myMaxBinTheta;
157}
158// --------------------------------------------------------
159template <typename T>
160inline
161typename DGtal::SphericalAccumulator<T>::Quantity
162DGtal::SphericalAccumulator<T>::count(const Size &posPhi,
163 const Size &posTheta) const
164{
165 ASSERT( isValidBin(posPhi,posTheta) );
166 return myAccumulator[posTheta + posPhi*myNtheta];
167}
168// --------------------------------------------------------
169template <typename T>
170inline
171T
172DGtal::SphericalAccumulator<T>::representativeDirection(const Size &posPhi,
173 const Size &posTheta) const
174{
175 ASSERT( isValidBin(posPhi,posTheta) );
176 Vector p=myAccumulatorDir[posTheta + posPhi*myNtheta] ;
177 return p;
178}
179// --------------------------------------------------------
180template <typename T>
181inline
182T
183DGtal::SphericalAccumulator<T>::representativeDirection(ConstIterator &it) const
184{
185 Size dist = it - myAccumulator.begin();
186 Vector p = *(myAccumulatorDir.begin() + dist);
187 return p;
188}
189// --------------------------------------------------------
190// --------------------------------------------------------
191template <typename T>
192inline
193void
194DGtal::SphericalAccumulator<T>::binCoordinates(ConstIterator &it,
195 Size &posPhi,
196 Size &posTheta) const
197{
198 Size dist = it - myAccumulator.begin();
199 posPhi = dist/ myNtheta;
200 posTheta = dist % myNtheta;
201}
202// --------------------------------------------------------
203template <typename T>
204inline
205bool
206DGtal::SphericalAccumulator<T>::isValidBin(const Size &posPhi,
207 const Size &posTheta) const
208{
209 ASSERT( myNphi != 1 );
210 double dphi = M_PI/((double)myNphi-1.0);
211 double Ntheta_i = floor(2.0*((double)myNphi)*sin((double)posPhi*dphi));
212
213 if ((posPhi == 0) || (posPhi == (myNphi-1)))
214 return (posTheta==0);
215 else
216 return (posPhi < myNphi) && (posTheta<Ntheta_i) && (posTheta< myNtheta);
217}
218// --------------------------------------------------------
219template <typename T>
220inline
221void
222DGtal::SphericalAccumulator<T>::getBinGeometry(const Size &posPhi,
223 const Size &posTheta,
224 RealVector &a,
225 RealVector &b,
226 RealVector &c,
227 RealVector &d) const
228{
229 ASSERT( isValidBin(posPhi,posTheta) );
230 double dphi = M_PI/(double)((double)myNphi-1);
231 double phi= (double)posPhi*dphi;
232 double dtheta;
233 double Nthetai = floor(2.0*(myNphi)*sin(phi) );
234 double theta = 0.0;
235
236 if ((posPhi == 0) ||(posPhi == (myNphi-1)))
237 {
238 a[0]=cos(theta-dphi/2.0)*sin(phi -dphi/2.0);
239 a[1]=sin(theta-dphi/2.0)*sin(phi -dphi/2.0);
240 a[2]=cos(phi -dphi/2.0);
241
242 b[0]=cos(theta+dphi/2.0)*sin(phi -dphi/2.0);
243 b[1]=sin(theta+dphi/2.0)*sin(phi -dphi/2.0);
244 b[2]=cos(phi -dphi/2.0);
245
246 c[0]=cos(theta-dphi/2.0)*sin(phi+dphi/2.0);
247 c[1]=sin(theta-dphi/2.0)*sin(phi+dphi/2.0);
248 c[2]=cos(phi+dphi/2.0);
249
250 d[0]=cos(theta+dphi/2.0)*sin(phi+dphi/2.0);
251 d[1]=sin(theta+dphi/2.0)*sin(phi+dphi/2.0);
252 d[2]=cos(phi+dphi/2.0);
253 }
254
255 else
256 {
257 dtheta = 2.0*M_PI/(Nthetai);
258 theta=posTheta*dtheta;
259
260 a[0]=cos(theta-dtheta/2.0)*sin(phi -dphi/2.0);
261 a[1]=sin(theta-dtheta/2.0)*sin(phi -dphi/2.0);
262 a[2]=cos(phi -dphi/2.0);
263
264
265 if ((double)posTheta -1 >= Nthetai)
266 {
267 b[0]=cos(-dtheta/2.0)*sin(phi -dphi/2.0);
268 b[1]=sin(-dtheta/2.0)*sin(phi -dphi/2.0);
269 b[2]=cos(phi -dphi/2.0);
270
271 c[0]=cos(- dtheta/2.0)*sin(phi+dphi/2.0);
272 c[1]=sin( - dtheta/2.0)*sin(phi+dphi/2.0);
273 c[2]=cos(phi+dphi/2.0);
274 }
275 else
276 {
277 b[0]=cos(theta+dtheta/2.0)*sin(phi -dphi/2.0);
278 b[1]=sin(theta+dtheta/2.0)*sin(phi -dphi/2.0);
279 b[2]=cos(phi -dphi/2.0);
280
281 c[0]=cos(theta+dtheta/2.0)*sin(phi+dphi/2.0);
282 c[1]=sin(theta+dtheta/2.0)*sin(phi+dphi/2.0);
283 c[2]=cos(phi+dphi/2.0);
284 }
285
286 d[0]=cos(theta-dtheta/2.0)*sin(phi+dphi/2.0);
287 d[1]=sin(theta-dtheta/2.0)*sin(phi+dphi/2.0);
288 d[2]=cos(phi+dphi/2.0);
289 }
290}
291// --------------------------------------------------------
292template <typename T>
293inline
294typename DGtal::SphericalAccumulator<T>::RealVector
295DGtal::SphericalAccumulator<T>::getBinDirection(const Size &posPhi,
296 const Size &posTheta) const
297{
298 ASSERT( isValidBin(posPhi,posTheta) );
299 double dphi = M_PI/(double)((double)myNphi-1);
300 double phi= (double)posPhi*dphi;
301 double Nthetai = floor(2.0*(myNphi)*sin(phi) );
302 RealVector midpoint;
303 double theta;
304
305 if ((posPhi == 0) || (posPhi == (myNphi-1)))
306 theta = 0;
307 else
308 theta=posTheta*2.0*M_PI/(Nthetai);
309
310 midpoint[0]=cos(theta)*sin(phi);
311 midpoint[1]=sin(theta)*sin(phi);
312 midpoint[2]=cos(phi);
313
314 return midpoint;
315}
316// --------------------------------------------------------
317template <typename T>
318inline
319void
320DGtal::SphericalAccumulator<T>::clear()
321{
322 myTotal = 0;
323 for(Size i=0; i < myNphi; i++)
324 for(Size j=0; j < myNtheta; j++)
325 if (isValidBin(i,j))
326 myAccumulator[j + i*myNtheta]=0;
327 else
328 myAccumulator[j+ i*myNtheta]=-1;
329
330 std::fill(myAccumulatorDir.begin(), myAccumulatorDir.end(), Vector::zero);
331}
332
333
334///////////////////////////////////////////////////////////////////////////////
335// Interface - public :
336
337/**
338 * Writes/Displays the object on an output stream.
339 * @param out the output stream where the object is written.
340 */
341template <typename T>
342inline
343void
344DGtal::SphericalAccumulator<T>::selfDisplay ( std::ostream & out ) const
345{
346 out << "[SphericalAccumulator] Nphi="<<myNphi<<" Ntheta=" <<myNtheta
347 <<" Number of samples="<<myTotal
348 <<" Number of bins="<<myBinNumber;
349}
350
351/**
352 * Checks the validity/consistency of the object.
353 * @return 'true' if the object is valid, 'false' otherwise.
354 */
355template <typename T>
356inline
357bool
358DGtal::SphericalAccumulator<T>::isValid() const
359{
360 return true;
361}
362
363
364
365///////////////////////////////////////////////////////////////////////////////
366// Implementation of inline functions //
367
368template <typename T>
369inline
370std::ostream&
371DGtal::operator<< ( std::ostream & out,
372 const SphericalAccumulator<T> & object )
373{
374 object.selfDisplay( out );
375 return out;
376}
377
378// //
379///////////////////////////////////////////////////////////////////////////////
380
381