source: Sophya/trunk/SophyaProg/PrgMap/msksphere.cc@ 1841

// Creation d'une sphere de masquee a partir d'une sphere de masques
//         cmv 15/6/01
#include "machdefs.h"
#include <unistd.h>
#include <stdexcept>
#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>
#include "skymapinit.h"
#include "skymap.h"
#include "fitsspherehealpix.h" 

/*!
  \ingroup PrgMap
  \file msksphere.cc
  \brief \b msksphere: mask a sphere of datas with a mask sphere.
  \verbatim
csh> msksphere -h
msksphere [-m min -M max -v valmsk] sphval.fits sphmask.fits sphout.fits
 sphval.fits : inpout sphere of datas where the pixels have to be masked
 sphmask.fits : inpout masked sphere used to mask pixels
 sphout.fits : output masked sphere of datas
 -m min : min mask_sphere_value for MASKING the sphere pixel
 -M max : max mask_sphere_value for MASKING the sphere pixel
    condition for MASKING sphere pixel is (min<=mask_sphere_value<=max)
 -n : negate the previous condition, condition for MASKING sphere pixel
      becomes: (mask_sphere_value<minval || maxval<mask_sphere_value)
 Default (no -m and -M): mask pixel if mask_sphere_value<0.
 -v valmsk : set sphere value for masked pixels (def=0.)
  \endverbatim
*/

void usage();
void usage()
{
cout<<"msksphere [-m min -M max -v valmsk]"
    <<" sphval.fits sphmask.fits sphout.fits"<<endl
    <<" sphval.fits : inpout sphere of datas where the pixels have to be masked"<<endl
    <<" sphmask.fits : inpout masked sphere used to mask pixels"<<endl
    <<" sphout.fits : output masked sphere of datas"<<endl
    <<" -m min : min mask_sphere_value for MASKING the sphere pixel"<<endl
    <<" -M max : max mask_sphere_value for MASKING the sphere pixel"<<endl
    <<"    condition for MASKING sphere pixel is (min<=mask_sphere_value<=max)"<<endl
    <<" -n : negate the previous condition, condition for MASKING sphere pixel"<<endl
    <<"      becomes: (mask_sphere_value<minval || maxval<mask_sphere_value)"<<endl
    <<" Default (no -m and -M): mask pixel if mask_sphere_value<0."<<endl
    <<" -v valmsk : set sphere value for masked pixels (def=0.)"<<endl;
}

int main(int narg, char* arg[])
{
double vmask=0.,vmin=-1.e30,vmax=1.e30;
bool tstmin=false,tstmax=false,negate=false;
string dum;
int c;
while((c = getopt(narg,arg,"hnm:M:v:")) != -1) {
  switch (c) {
  case 'n' :
    negate = true;
    break;
  case 'm' :
    sscanf(optarg,"%lf",&vmin);
    tstmin = true;
    break;
  case 'M' :
    sscanf(optarg,"%lf",&vmax);
    tstmax = true;
    break;
  case 'v' :
    sscanf(optarg,"%lf",&vmask);
    break;
  case 'h' :
  default:
    usage(); exit(1);
  }
}

if(optind+2>=narg) {usage(); exit(1);}
char * sphval = arg[optind];
char * sphmsk = arg[optind+1];
char * sphout = arg[optind+2];

if(!tstmin && !tstmax) {tstmin=true; vmin=0.; negate=true;}

if(negate) dum = ".NOT."; else dum="";
cout<<"Sphere values : "<<sphval<<endl
    <<"Sphere mask   : "<<sphmsk<<endl
    <<"Sphere out    : "<<sphout<<endl
    <<"  ...min("<<tstmin<<") "<<vmin<<endl
    <<"     max("<<tstmax<<") "<<vmax<<endl
    <<"  ...mask set value "<<vmask<<endl
    <<"  Condition for masking pixel is :"<<endl
    <<"    "<<dum<<"( "<<vmin<<" <= V <= "<<vmax<<" )"<<endl;

// Lecture de la sphere Healpix des valeurs
SphereHEALPix<r_8> sph;        
{FitsInFile sfits(sphval); sfits >> sph;}
cout<<"Opening Sphere HEALPix for testing values :"<<endl
    <<"          Type of map : "<<sph.TypeOfMap()<<endl
    <<"     Number of pixels : "<<sph.NbPixels()<<endl
    <<"                 Nlat : "<<sph.SizeIndex()<<endl;

// Lecture de la sphere Healpix des masques
SphereHEALPix<r_8> msph;        
{FitsInFile sfits(sphmsk); sfits >> msph;}
cout<<"Opening Sphere HEALPix for mask :"<<endl
    <<"          Type of map : "<<msph.TypeOfMap()<<endl
    <<"     Number of pixels : "<<msph.NbPixels()<<endl
    <<"                 Nlat : "<<msph.SizeIndex()<<endl;
if(msph.SizeIndex()!=sph.SizeIndex()) {
  cout<<"Sphere of values and sphere of mask must have same NLat"<<endl;
  exit(1);
}

// Filling Mask Sphere
cout<<"Testing For Masked Pixels"<<endl;
uint_4 nmask = 0;
for(int_4 i=0;i<sph.NbPixels();i++) {

  bool skp = (tstmin || tstmax) ? negate : false;
  if((tstmin && msph(i)<vmin) || (tstmax && msph(i)>vmax)) skp = !negate;
  if(skp) continue;   // Do nothing

  sph(i) = vmask;
  nmask++;
}
cout<<"    .... Number of values masked   : "<<nmask<<endl;
cout<<"    .... Fraction of values masked : "
    <<(r_8)nmask/sph.NbPixels()*100.<<" %"<<endl;

// Ecriture de la sphere Healpix sur fits
{
dum = "rm -f "; dum += sphout; system(dum.c_str());
FitsOutFile sfits(sphout, FitsFile::clear);
cout<<"Writing Output Masked Sphere Fits file"<<endl;
sfits<<sph;
}

exit(0);
}