Changeset 1528 in Sophya for trunk/SophyaProg/PrgMap/extrap2sph.cc

Timestamp:

Jun 14, 2001, 5:16:20 PM (24 years ago)

Author:

cmv

Message:

optimisation cmv 14/6/01

File:

• trunk/SophyaProg/PrgMap/extrap2sph.cc (modified) (8 diffs)

trunk/SophyaProg/PrgMap/extrap2sph.cc

@@ -1 +1 @@
 // Pour boucher les trous d'une sphere HEALPIX utilisant une autre sphere
 //             cmv 13/6/01
+// extrap2sph -w sph143k05_e.fits sph143k05.fits sphred.fits sphout_2.fits
+// extrap2sph -m -1.e-30 -M +1.e-30 sph143k05.fits sphred.fits sphout_22.fits
 #include "machdefs.h"
 #include <unistd.h>
@@ -14 +16 @@
 void usage()
 {
-cout<<"extrap2sph [-w sphinw.fits -0 nulval]"<<endl
+cout<<"extrap2sph [-w sphinw.fits -m valmin -M valmax]"<<endl
     <<"           sphin.fits sphred.fits sphout.fits"<<endl
     <<" -w sphinw.fits : input sphere  weight (def=NULL)"<<endl
-    <<" -0 nulval : value to identify empty pixel in sphin (def=0.)"<<endl
-    <<"             (used if sphinw.fits==NULL)"<<endl
-    <<" sphin.fits : input sphere"<<endl
+    <<" -m minval : value to identify EMPTY pixel (def=-1.e30.)"<<endl
+    <<" -M maxval : value to identify EMPTY pixel (def=+1.e30)"<<endl
+    <<"    condition for EMPTY pixel is (minval<val<maxval)"<<endl
+    <<"    (bounds are excluded, used only if sphinw.fits==NULL)"<<endl
+    <<" sphin.fits   : input sphere"<<endl
     <<" sphred.fits  : input reducted sphere"<<endl
     <<" sphout.fits  : output sphere"<<endl;
@@ -26 +30 @@
 int main(int narg, char* arg[])
 {
-double nulval=0.;
+bool tstmin=false, tstmax=false;
+double valmin=-1.e30, valmax=+1.e30;
 char * fsphinw = NULL;
 int c;
-while((c = getopt(narg,arg,"h0:w:")) != -1) {
+while((c = getopt(narg,arg,"hm:M:w:")) != -1) {
   switch (c) {
-  case '0' :
-    sscanf(optarg,"%lf",&nulval);
+  case 'm' :
+    sscanf(optarg,"%lf",&valmin);
+    tstmin=true;
+    break;
+  case 'M' :
+    sscanf(optarg,"%lf",&valmax);
+    tstmax=true;
     break;
   case 'w' :
@@ -48 +58 @@
 char * fsphout  = arg[optind+2];
 
-cout<<"Sphere Input : "<<fsphin<<endl
-    <<"Weight Sphere Input : "<<fsphinw<<endl
-    <<"Reducted Sphere Input : "<<fsphred<<endl
-    <<"Sphere Output : "<<fsphout<<endl
-    <<"nulval : "<<nulval<<endl;
+cout<<"Input Sphere          : "<<fsphin<<endl
+    <<"Input Weight Sphere   : "<<fsphinw<<endl
+    <<"Input Reducted Sphere : "<<fsphred<<endl
+    <<"Output Sphere         : "<<fsphout<<endl;
+if(!fsphinw)
+  cout<<"- test min("<<tstmin<<") : "<<valmin<<endl
+      <<"  test max("<<tstmax<<") : "<<valmax<<endl
+      <<"Condition for EMPTY pixel is :\n    ("<<valmin
+      <<" < sphere value < "<<valmax<<") bounds are excluded!"<<endl;
 
 // Lecture des spheres
@@ -58 +72 @@
 {FitsInFile sfits(fsphin); sfits >> sphin;}
 int nlat = sphin.SizeIndex();
-cout<<"Opening Sphere Input :"<<endl
+cout<<"Opening Input Sphere :"<<endl
     <<"          Type of map : "<<sphin.TypeOfMap()<<endl
     <<"     Number of pixels : "<<sphin.NbPixels()<<endl
@@ -67 +81 @@
   FitsInFile sfits(fsphinw);
   sfits >> sphinw;
-  cout<<"Opening Weight Sphere Input :"<<endl
+  cout<<"Opening Input Weight Sphere :"<<endl
       <<"          Type of map : "<<sphinw.TypeOfMap()<<endl
       <<"     Number of pixels : "<<sphinw.NbPixels()<<endl
@@ -79 +93 @@
 SphereHEALPix<r_8> sphred;
 {FitsInFile sfits(fsphred); sfits >> sphred;}
-cout<<"Opening Reducted Sphere Input :"<<endl
+cout<<"Opening Input Reducted Sphere :"<<endl
     <<"          Type of map : "<<sphred.TypeOfMap()<<endl
     <<"     Number of pixels : "<<sphred.NbPixels()<<endl
     <<"                 Nlat : "<<sphred.SizeIndex()<<endl;
 
-// Creation de la sphere output
-cout<<"Creating Sphere for output"<<endl;
-SphereHEALPix<r_8> sphout(sphin,false);
-
 // Filling hole for Output Sphere
 cout<<"...Filling hole for Output Sphere"<<endl;
-uint_4 n=0;
-for(int_4 i=0;i<sphout.NbPixels();i++) {
-  if(fsphinw) {if(sphinw(i)>0.) continue;}
-     else     {if(sphout(i)!=nulval) continue;}
+uint_4 n=0, n0=0;
+for(int_4 i=0;i<sphin.NbPixels();i++) {
+  if(fsphinw) {
+    if(sphinw(i)>0.) continue;
+  } else {
+    // Not Filled pixels ]valmin,valmax[
+    // Filled pixels ]-Infinity,valmin] or [valmax,+Infinity[
+    if(tstmin && sphin(i)<=valmin) continue;
+    if(tstmax && sphin(i)>=valmax) continue;
+  }
   double theta,phi;
-  sphout.PixThetaPhi(i,theta,phi);
+  sphin.PixThetaPhi(i,theta,phi);
   int_4 ir = sphred.PixIndexSph(theta,phi);
-  sphout(i) = sphred(ir);
-  n++;
+  sphin(i) = sphred(ir);
+  n++; if(sphred(ir)!=0.) n0++;
 }
-cout<<"      Number of filled pixels "<<n<<endl;
+cout<<"      Number of pixels filled "
+    <<n<<" (not null values "<<n0<<")"<<endl;
 
 // Ecriture de la sphere
@@ -107 +124 @@
 FitsOutFile sfits(fsphout);
 cout<<"Writing Output Sphere Fits file"<<endl;
-sfits<<sphout;
+sfits<<sphin;
 }