51 myAccumulator.resize(myNphi*myNtheta,0);
52 myAccumulatorDir.resize(myNphi*myNtheta, Vector::zero);
59 for(
Size posPhi=0; posPhi < myNphi; posPhi++)
60 for(
Size posTheta=0; posTheta < myNtheta; posTheta++)
62 double dphi =
M_PI/((double)myNphi-1);
63 double Ntheta_i = floor(2.0*((
double)myNphi)*sin((
double)posPhi*dphi)) ;
65 if (((posPhi == 0) || (posPhi == (myNphi-1))) && (posTheta==0))
68 if ((posPhi < myNphi) && (posTheta<Ntheta_i) && (posTheta< myNtheta))
88 double norm = aDir.norm();
94 double dphi =
M_PI/(double)(myNphi-1);
96 posPhi =
static_cast<Size>(floor( (phi+dphi/2.) *(myNphi-1)/
M_PI));
97 if(posPhi == 0 || posPhi== (myNphi-1))
109 Nthetai = floor(2.0*(myNphi)*sin(posPhi*dphi));
110 double dtheta = 2.0*
M_PI/(Nthetai);
111 posTheta =
static_cast<Size>(floor( (theta+dtheta/2.0)/dtheta));
113 if (posTheta >= Nthetai)
117 ASSERT(posPhi < myNphi);
118 ASSERT( isValidBin(posPhi,posTheta) );
121 template <
typename T>
125 Size posPhi,posTheta;
127 binCoordinates(aDir , posPhi, posTheta);
128 myAccumulator[posTheta + posPhi*myNtheta] += 1;
129 myAccumulatorDir[posTheta + posPhi*myNtheta] += aDir;
133 if ( myAccumulator[posTheta + posPhi*myNtheta] >
134 myAccumulator[ myMaxBinTheta + myMaxBinPhi*myNtheta])
136 myMaxBinTheta = posTheta;
137 myMaxBinPhi = posPhi;
141 template <
typename T>
149 template <
typename T>
155 theta = myMaxBinTheta;
158 template <
typename T>
162 const Size &posTheta)
const
164 ASSERT( isValidBin(posPhi,posTheta) );
165 return myAccumulator[posTheta + posPhi*myNtheta];
168 template <
typename T>
172 const Size &posTheta)
const
174 ASSERT( isValidBin(posPhi,posTheta) );
175 Vector p=myAccumulatorDir[posTheta + posPhi*myNtheta] ;
179 template <
typename T>
184 Size dist = it - myAccumulator.begin();
185 Vector p = *(myAccumulatorDir.begin() + dist);
190 template <
typename T>
195 Size &posTheta)
const
197 Size dist = it - myAccumulator.begin();
198 posPhi = dist/ myNtheta;
199 posTheta = dist % myNtheta;
202 template <
typename T>
206 const Size &posTheta)
const
208 double dphi =
M_PI/((double)myNphi-1.0);
209 double Ntheta_i = floor(2.0*((
double)myNphi)*sin((
double)posPhi*dphi));
211 if ((posPhi == 0) || (posPhi == (myNphi-1)))
212 return (posTheta==0);
214 return (posPhi < myNphi) && (posTheta<Ntheta_i) && (posTheta< myNtheta);
217 template <
typename T>
221 const Size &posTheta,
227 ASSERT( isValidBin(posPhi,posTheta) );
228 double dphi =
M_PI/(double)((
double)myNphi-1);
229 double phi= (double)posPhi*dphi;
231 double Nthetai = floor(2.0*(myNphi)*sin(phi) );
233 dtheta = 2.0*
M_PI/(Nthetai);
235 double theta=posTheta*dtheta;
237 a[0]=cos(theta-dtheta/2.0)*sin(phi -dphi/2.0);
238 a[1]=sin(theta-dtheta/2.0)*sin(phi -dphi/2.0);
239 a[2]=cos(phi -dphi/2.0);
242 if ((
double)posTheta -1 >= Nthetai)
244 b[0]=cos(-dtheta/2.0)*sin(phi -dphi/2.0);
245 b[1]=sin(-dtheta/2.0)*sin(phi -dphi/2.0);
246 b[2]=cos(phi -dphi/2.0);
248 c[0]=cos(- dtheta/2.0)*sin(phi+dphi/2.0);
249 c[1]=sin( - dtheta/2.0)*sin(phi+dphi/2.0);
250 c[2]=cos(phi+dphi/2.0);
254 b[0]=cos(theta+dtheta/2.0)*sin(phi -dphi/2.0);
255 b[1]=sin(theta+dtheta/2.0)*sin(phi -dphi/2.0);
256 b[2]=cos(phi -dphi/2.0);
258 c[0]=cos(theta+dtheta/2.0)*sin(phi+dphi/2.0);
259 c[1]=sin(theta+dtheta/2.0)*sin(phi+dphi/2.0);
260 c[2]=cos(phi+dphi/2.0);
263 d[0]=cos(theta-dtheta/2.0)*sin(phi+dphi/2.0);
264 d[1]=sin(theta-dtheta/2.0)*sin(phi+dphi/2.0);
265 d[2]=cos(phi+dphi/2.0);
268 template <
typename T>
274 for(
Size i=0; i < myNphi; i++)
275 for(
Size j=0; j < myNtheta; j++)
277 myAccumulator[j + i*myNtheta]=0;
279 myAccumulator[j+ i*myNtheta]=-1;
281 std::fill(myAccumulatorDir.begin(), myAccumulatorDir.end(), Vector::zero);
292 template <
typename T>
297 out <<
"[SphericalAccumulator] Nphi="<<myNphi<<
" Ntheta=" <<myNtheta
298 <<
" Number of samples="<<myTotal
299 <<
" Number of bins="<<myBinNumber;
306 template <
typename T>
319 template <
typename T>
325 object.selfDisplay( out );