Changeset 2991 in Sophya for trunk/SophyaLib/Samba/sphericaltransformserver.cc

Timestamp:

Jun 23, 2006, 12:47:49 PM (19 years ago)

Author:

ansari

Message:

Corrections dans SphericalTransformServer:

• gestion taille tableau sortie + imag(f_nyquist) dans RfourierSynthesisFromB
• Utilisation RfourierSynthesisFromB (TF cmplx->reel) pour synthese

depuis alm, pour tout type de carte

• Optimisation acces aux pixels avec la methode GetThetaSliceDataPtr()

Reza , 23 Juin 2006

File:

• trunk/SophyaLib/Samba/sphericaltransformserver.cc (modified) (2 diffs)

trunk/SophyaLib/Samba/sphericaltransformserver.cc

@@ -180 +180 @@
   Bm<complex<T> > b_m_theta(nmmax);
 
-
-
   // pour chaque tranche en theta
-  for (int_4 ith = premiereTranche; ith <= derniereTranche;ith++)
-    {
-      int_4 nph;
-      r_8 phi0;
-      r_8 theta;
-      TVector<int_4> pixNumber; 
-      TVector<T> datan;
-      
-      map.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
-      nph = pixNumber.NElts();
-      if (nph < 2) continue;  // On laisse tomber les tranches avec un point
-      //       -----------------------------------------------------
-      //              for each theta, and each m, computes
-      //              b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m) 
-      //              ------------------------------------------------------
-      // ===> Optimisation Reza, Mai 2006 
-      /*---  Le bout de code suivant est remplace par l'appel a la nouvelle fonction
-        qui calcule la somme au vol 
+  for (int_4 ith = premiereTranche; ith <= derniereTranche;ith++)  {
+    int_4 nph;
+    r_8 phi0;
+    r_8 theta;
+    TVector<int_4> pixNumber; 
+    TVector<T> datan;
+    
+    map.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
+    nph = pixNumber.NElts();
+    if (nph < 2) continue;  // On laisse tomber les tranches avec un point
+    //       -----------------------------------------------------
+    //        for each theta, and each m, computes
+    //        b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m) 
+    //        ------------------------------------------------------
+    // ===> Optimisation Reza, Mai 2006 
+    /*---  Le bout de code suivant est remplace par l'appel a la nouvelle fonction
+      qui calcule la somme au vol 
       LambdaLMBuilder lb(theta,nlmax,nmmax);
       //  somme sur m de 0 a l'infini
-      for (int_4 m = 0; m <= nmmax; m++)
-        {
-          b_m_theta(m) = (T)( lb.lamlm(m,m) ) * alm(m,m);
-          for (int l = m+1; l<= nlmax; l++)
-            {
-               b_m_theta(m) += (T)( lb.lamlm(l,m) ) * alm(l,m);
-            }
-        }
+      for (int_4 m = 0; m <= nmmax; m++) {
+      b_m_theta(m) = (T)( lb.lamlm(m,m) ) * alm(m,m);
+      for (int l = m+1; l<= nlmax; l++)
+        b_m_theta(m) += (T)( lb.lamlm(l,m) ) * alm(l,m);
+      }
       ------- Fin version PRE-Mai2006 */
-      LambdaLMBuilder::ComputeBmFrAlm(theta,nlmax,nmmax, alm, b_m_theta);
-      //Fin Optimisation Reza, Mai 2006 <==== 
+    LambdaLMBuilder::ComputeBmFrAlm(theta,nlmax,nmmax, alm, b_m_theta);
+    //Fin Optimisation Reza, Mai 2006 <==== 
 
     //        obtains the negative m of b(m,theta) (= complex conjugate)
-
-      for (int_4 m=1;m<=nmmax;m++)
-        {
-          b_m_theta(-m) = conj(b_m_theta(m));
-        }
-      // ---------------------------------------------------------------
-      //    sum_m  b(m,theta)*exp(i*m*phi)   -> f(phi,theta)
-      // ---------------------------------------------------------------*/
-
-      /* ----- Reza, Juin 2006 : 
-         En verifiant la difference entre deux cartes 
-         cl -> map -> alm -> map2 et mapdiff = map-map2 
-         je me suis apercu qu'il y avait des differences importantes - dans les
-         deux zones 'polar cap' de HEALPix - apres avoir pas mal chercher, il semble
-         que la routine RfourierSynthesisFromB en soit responsable 
-         Je fais donc tout passer dans fourierSynthesisFromB 
-      if ( (sphere_type == "RING") || (sphere_type == "NESTED") ) {
-          TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
-          for (int i=0;i< nph;i++) map(pixNumber(i))=Temp(i); 
-        }
-      else 
-      */
-        // pour des pixelisations quelconques (autres que HEALPix
-        //  nph n'est pas toujours pair
-        // ca fait des problemes pour les transformees de Fourier
-        // car le server de TF ajuste la longueur du vecteur reel 
-        // en sortie de TF, bref, la securite veut qu'on prenne une 
-        // TF complexe
-        {
-          TVector<complex<T> > Temp = fourierSynthesisFromB(b_m_theta,nph,phi0);
-          for (int i=0;i< nph;i++) map(pixNumber(i))=Temp(i).real();
-        }
-    }
+    for (int_4 m=1;m<=nmmax;m++)
+      b_m_theta(-m) = conj(b_m_theta(m));
+    // ---------------------------------------------------------------
+    //    sum_m  b(m,theta)*exp(i*m*phi)   -> f(phi,theta)
+    // ---------------------------------------------------------------*/
+
+    /* ----- Reza, Juin 2006 : 
+       En verifiant la difference entre deux cartes 
+       cl -> map -> alm -> map2 et mapdiff = map-map2 
+       je me suis apercu qu'il y avait des differences importantes - dans les
+       deux zones 'polar cap' de HEALPix - qui utilisait RfourierSynthesisFromB
+       TF complex -> reel . Le probleme venant de l'ambiguite de taille, lie 
+       a la partie imaginaire de la composante a f_nyquist , j'ai corrige et 
+       tout mis en TF complexe -> reel 
+    */
+    TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
+    // Si on peut acceder directement les pixels d'un tranche, on le fait
+    T* pix = map.GetThetaSliceDataPtr(ith);
+    if (pix != NULL) 
+      for (int_4 i=0;i< nph;i++) pix[i] = Temp(i);
+    else 
+      for (int_4 i=0;i< nph;i++) map(pixNumber(i))=Temp(i); 
+  }
 }
 
@@ -380 +364 @@
   //**********************************************************************
   TVector< complex<T> > bw(nph);
-  TVector< complex<T> > dataout(nph);
   TVector< complex<T> > data(nph/2+1);
-
 
   for (int kk=0; kk<bw.NElts(); kk++) bw(kk)=(T)0.;
   int m;
-  for (m=-b_m.Mmax();m<=-1;m++)
-    {
-      int maux=m;
-      while (maux<0) maux+=nph;
-      int iw=maux%nph;
-      double aux=(m-iw)*phi0;
-      bw(iw) += b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) )  ;
-    }
-  for (m=0;m<=b_m.Mmax();m++)
-    {
-      //      int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
-      int iw=m%nph;
-      double aux=(m-iw)*phi0;
-      bw(iw)+=b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) );
-    }
+  for (m=-b_m.Mmax();m<=-1;m++)  {
+    int maux=m;
+    while (maux<0) maux+=nph;
+    int iw=maux%nph;
+    double aux=(m-iw)*phi0;
+    bw(iw) += b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) )  ;
+  }
+  for (m=0;m<=b_m.Mmax();m++) {
+    //      int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
+    int iw=m%nph;
+    double aux=(m-iw)*phi0;
+    bw(iw)+=b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) );
+  }
 
   //     applies the shift in position <-> phase factor in Fourier space
-  // cout << " TF : nph= " << nph << " vec. entree " << data.Size() << endl;
-  for (int mprime=0; mprime <= nph/2; mprime++)
-    {
-      complex<double> aux(cos(mprime*phi0),sin(mprime*phi0));
-      data(mprime)=bw(mprime)*
-                       (complex<T>)(complex<double>(cos(mprime*phi0),sin(mprime*phi0)));
-    }
-
-  TVector<T> sortie;
-  fftIntfPtr_-> FFTBackward(data, sortie);
-  
+  for (int mprime=0; mprime <= nph/2; mprime++)  
+    data(mprime)=bw(mprime)*complex<T>((T)cos(mprime*phi0),(T)sin(mprime*phi0));    
+  TVector<T> sortie(nph);
+// On met la partie imaginaire du dernier element du data a zero pour nph pair
+  if (nph%2 == 0) data(nph/2) = complex<T>(data(nph/2).real(), (T)0.);
+// et on impose l'utilisation de la taille en sortie pour FFTBack (..., ..., true)
+  fftIntfPtr_-> FFTBackward(data, sortie, true); 
   return sortie;
 }