(JEC) 20/4/15 add pixelization aids

lagSphericTransform.cc

@@ -67,9 +67,9 @@ void BaseLagSphTransform::SetThetaPhiMap(string choice, int nrings, int nphi){
   Nalm_ = get_nalms(sharpjob_.get_alm_info());
 
   //fill the ANN-kdTree to prepare the filling of the 2D spheres
-  geom_t vecPix;
-  vecPix = GetPixelPositions();
-  ann.SetkdTree(vecPix);
+  //  geom_t vecPix;
+  SetPixelPositions();
+  ann.SetkdTree(vecPix_); 
 
   geomDone_ = true;
     
@@ -82,8 +82,7 @@ void BaseLagSphTransform::SetThetaPhiMap(string choice, int nrings, int nphi){
 
 }//SetThetaPhiMap
 
-geom_t BaseLagSphTransform::GetPixelPositions() {
-  geom_t pixels;
+void BaseLagSphTransform::SetPixelPositions() {
   sharp_geom_info *ginfo; 
   ginfo = sharpjob_.get_geom_info();
   for (int i=0; i<ginfo->npairs; ++i) {
@@ -93,8 +92,8 @@ geom_t BaseLagSphTransform::GetPixelPositions() {
       double phi0 = ginfo->pair[i].r1.phi0;
       for(int j=0; j<nph1; j++) {
 	double phi1 = phi0 + j*(2.0*M_PI/nph1);
-	//	cout << "pixels ( "<< theta1 << ", " << phi1 << ")" << endl;
-	pixels.push_back(make_pair(theta1,phi1));
+	//	cout << "vecPix_s ( "<< theta1 << ", " << phi1 << ")" << endl;
+	vecPix_.push_back(make_pair(theta1,phi1));
       }
       if (ginfo->pair[i].r2.nph>0) {// the second ring in this pair exists
 	double theta2 = ginfo->pair[i].r2.theta;
@@ -102,16 +101,16 @@ geom_t BaseLagSphTransform::GetPixelPositions() {
 	double phi0 = ginfo->pair[i].r2.phi0;
 	for(int j=0; j<nph2; j++) {
 	  double phi2 = phi0 + j*(2.0*M_PI/nph2);
-	  //	  cout << "pixels ( "<< theta2 << ", " << phi2 << ")" << endl;
-	  pixels.push_back(make_pair(theta2,phi2));
+	  //	  cout << "vecPix_s ( "<< theta2 << ", " << phi2 << ")" << endl;
+	  vecPix_.push_back(make_pair(theta2,phi2));
 	}
       }
   }
 
-  sort(pixels.begin(),pixels.end());
+
+  sort(vecPix_.begin(),vecPix_.end()); 
   
-  return pixels;
-}//GetPixelPositions
+}//SetPixelPositions
 
 
 

lagSphericTransform.h

@@ -35,6 +35,8 @@ namespace LagSHT {
     ANNServer(): nNear_(1), dim_(2), eps_(0), 
 		 npts_(0), dataPts_(NULL), kdTree_(NULL), nnIdx_(NULL), dists_(NULL) {
       
+      nnIdx_ = new ANNidx[nNear_];  // allocate near neigh indices
+      dists_ = new ANNdist[nNear_]; // allocate near neighbor dists      
     }
 
     //! Constructor
@@ -78,9 +80,9 @@ namespace LagSHT {
       queryPt[0] = theta;
       queryPt[1] = phi;
       
-      if(kdTree_ == NULL) throw LagSHTError("Get Failure kdTree alloc");
-      if(nnIdx_ == NULL) throw LagSHTError("Get Failure nnIdx alloc");
-      if(dists_ == NULL) throw LagSHTError("Get Failure dists alloc");
+      if (kdTree_ == NULL) throw LagSHTError("Get Failure kdTree alloc");
+      if (nnIdx_ == NULL) throw LagSHTError("Get Failure nnIdx alloc");
+      if (dists_ == NULL) throw LagSHTError("Get Failure dists alloc");
       
       kdTree_->annkSearch(queryPt,nNear_,nnIdx_,dists_,eps_);
       annDeallocPt(queryPt);
@@ -128,8 +130,10 @@ class BaseLagSphTransform {
   //! Total number of Alm coefficients
   int Nalm()   { return Nalm_; }
   
+
   //! Sorted List of pixel center (theta, phi) positions
-  geom_t GetPixelPositions();
+  geom_t GetPixelPositions() { return vecPix_; }
+
     
   /*! Choose the 2D map representation & initialization of libsharp
     \input ECP/Fejer1, Gauss, and Healpix
@@ -143,10 +147,20 @@ class BaseLagSphTransform {
     \input phi   in [0, 2pi]
   */
   int PixIndexSph(r_8 theta, r_8 phi){ return ann.GetClosestId(theta, phi); }
+  /*! Get Theta Phi position of a given pixel index
+    \input k in [0, Npix-1] no boundchckeing done
+   */
 
+  void PixThetaPhi(int k, r_8& theta, r_8& phi) { 
+    pixpos_t pix = vecPix_[k];
+    theta = pix.first;
+    phi   = pix.second;
+  }
 
 protected:
 
+  void SetPixelPositions();
+
   //from sharp test suite
   static ptrdiff_t get_npix(const sharp_geom_info *ginfo) {
     ptrdiff_t res=0;
@@ -175,6 +189,9 @@ protected:
   int Nphi_; //!< Number of phi values (pixels) per colatitude rings
   int Npix_; //!< Total number of 2D pixels
   int Nalm_; //!< Total number of Alm coefficients
+
+  geom_t vecPix_; //!< vector of pixel positions
+
   bool geomDone_; //!< Set to true once the geometry is properly set (See SetThetaPhiMap(string) )
 
   sharp_cxxjob<r_8> sharpjob_; //!< libsharp job // voir si r_8 -> r_16 extend the validity of the computation

laguerreTransform.cc

@@ -73,10 +73,9 @@ LaguerreTransform::LaguerreTransform(int N, r_8 R, int alpha) :
 
 }//Ctor
 
-int LaguerreTransform::RadIndex(r_8 r) const {
+int LaguerreTransform::R2Index(r_8 r) const {
   using namespace std;
   r_8 rscaled = r*nodes_[N_-1]/R_;
-  
   vector<r_8> dist(N_);
   transform(nodes_.begin(),nodes_.end(),dist.begin(), distance(rscaled));
   return min_element(dist.begin(),dist.end()) - dist.begin();

laguerreTransform.h

@@ -42,8 +42,14 @@ class LaguerreTransform : public BaseLaguerreTransform {
   vector<r_8> GetWeights() const { return weights_; }
 
 
-  //! Find the closest nodes
-  int RadIndex(r_8 r) const;
+  //! Find the closest nodes: no boundchecking
+  int R2Index(r_8 r) const;
+  //! return the nodes position: no boundchecking
+  r_8 Index2R(int k) const { 
+    return nodes_[k];
+  }
+  
+  
 
 
   //! Multi Analysis at once (Values -> Coefficients)