// Pour boucher les trous d'une sphere HEALPIX en utilisant
// une autre sphere HEALPIX
//             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 "sopnamsp.h"
#include "machdefs.h"
#include <unistd.h>
#include <stdexcept>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include "skymapinit.h"
#include "skymap.h"
#include "fitsspherehealpix.h" 

/*!
  \ingroup PrgMap
  \file extrap2sph.cc
  \brief \b extrap2sph: Fill holes in a sphere of datas by using
  an other sphere of datas.
  A mask to define the place to be tested in the sphere might be given.
  \verbatim
csh> extrap2sph -h
extrap2sph [-m valmin -M valmax -r rvalmin -R rvalmax]
           [-w sphinw.fits] sphin.fits sphred.fits sphout.fits
 -w sphinw.fits : input sphere  weight, if sphinw_value<=0. (def=NULL)
                  if sphinw_value<=0. pixel has not to be filled
 -m minval : min value to identify EMPTY pixel of sphin (def=no_test)
 -M maxval : max value to identify EMPTY pixel of sphin (def=no_test)
    condition for EMPTY pixel is (minval<=val<=maxval)
 -n : negate the previous condition, condition for EMPTY pixel
      becomes: (val<minval || maxval<val)
 -r rminval : min value to identify GOOD pixel of sphred (def=no_test)
 -R rmaxval : max value to identify GOOD pixel of sphred (def=no_test)
    condition for GOOD pixel is (minval<=val<=maxval)
 -n : negate the previous condition, condition for GOOD pixel
      becomes: (val<minval || maxval<val)
 sphin.fits   : input sphere
 sphred.fits  : input reducted sphere
 sphout.fits  : output sphere
  \endverbatim
*/

void usage();
void usage()
{
cout<<"extrap2sph [-m valmin -M valmax -r rvalmin -R rvalmax]"<<endl
    <<"           [-w sphinw.fits] sphin.fits sphred.fits sphout.fits"<<endl
    <<" -w sphinw.fits : input sphere  weight, if sphinw_value<=0. (def=NULL)"<<endl
    <<"                  if sphinw_value<=0. pixel has not to be filled"<<endl
    <<" -m minval : min value to identify EMPTY pixel of sphin (def=no_test)"<<endl
    <<" -M maxval : max value to identify EMPTY pixel of sphin (def=no_test)"<<endl
    <<"    condition for EMPTY pixel is (minval<=val<=maxval)"<<endl
    <<" -n : negate the previous condition, condition for EMPTY pixel"<<endl
    <<"      becomes: (val<minval || maxval<val)"<<endl
    <<" -r rminval : min value to identify GOOD pixel of sphred (def=no_test)"<<endl
    <<" -R rmaxval : max value to identify GOOD pixel of sphred (def=no_test)"<<endl
    <<"    condition for GOOD pixel is (minval<=val<=maxval)"<<endl
    <<" -n : negate the previous condition, condition for GOOD pixel"<<endl
    <<"      becomes: (val<minval || maxval<val)"<<endl
    <<" sphin.fits   : input sphere"<<endl
    <<" sphred.fits  : input reducted sphere"<<endl
    <<" sphout.fits  : output sphere"<<endl;
}

int main(int narg, char* arg[])
{
bool tstmin=false,  tstmax=false,  neg=false;
bool rtstmin=false, rtstmax=false, rneg=false;
double valmin=0., valmax=0., rvalmin=0., rvalmax=0.;
char * fsphinw = NULL;
string dum;
int c;
while((c = getopt(narg,arg,"hnNm:M:r:R:w:")) != -1) {
  switch (c) {
  case 'n' :
    neg = true;
    break;
  case 'N' :
    rneg = true;
    break;
  case 'm' :
    sscanf(optarg,"%lf",&valmin);
    tstmin=true;
    break;
  case 'M' :
    sscanf(optarg,"%lf",&valmax);
    tstmax=true;
    break;
  case 'r' :
    sscanf(optarg,"%lf",&rvalmin);
    rtstmin=true;
    break;
  case 'R' :
    sscanf(optarg,"%lf",&rvalmax);
    rtstmax=true;
    break;
  case 'w' :
    fsphinw = optarg;
    break;
  case 'h' :
  default:
    usage(); exit(1);
  }
}

if(optind+2>=narg) {usage(); exit(1);}
char * fsphin  = arg[optind];
char * fsphred = arg[optind+1];
char * fsphout = arg[optind+2];

cout<<"Input Sphere          : "<<fsphin<<endl
    <<"Input Weight Sphere   : "<<fsphinw<<endl
    <<"Input Reducted Sphere : "<<fsphred<<endl
    <<"Output Sphere         : "<<fsphout<<endl;
if(neg) dum = ".NOT."; else dum="";
cout<<"- Sphere - test min("<<tstmin<<") : "<<valmin<<endl
    <<"           test max("<<tstmax<<") : "<<valmax<<endl
    <<"           negate("<<neg<<")"<<endl
    <<"  - Condition for EMPTY pixel in Sphere is :"<<endl
    <<"    "<<dum<<"( "<<valmin<<" <= V <= "<<valmax<<" )"<<endl;
if(rneg) dum = ".NOT."; else dum="";
cout<<"- Reducted Sphere - test min("<<rtstmin<<") : "<<rvalmin<<endl
    <<"                    test max("<<rtstmax<<") : "<<rvalmax<<endl
    <<"                    negate("<<rneg<<")"<<endl
    <<"  - Condition for GOOD pixel in Reducted Sphere is :"<<endl
    <<"    "<<dum<<"( "<<rvalmin<<" <= V <= "<<rvalmax<<" )"<<endl;

// Lecture des spheres
SphereHEALPix<r_8> sphin;
{FitsInFile sfits(fsphin); sfits >> sphin;}
int nlat = sphin.SizeIndex();
cout<<"Opening Input Sphere :"<<endl
    <<"          Type of map : "<<sphin.TypeOfMap()<<endl
    <<"     Number of pixels : "<<sphin.NbPixels()<<endl
    <<"                 Nlat : "<<sphin.SizeIndex()<<endl;

SphereHEALPix<r_8> sphinw;
if(fsphinw) {
  FitsInFile sfits(fsphinw);
  sfits >> sphinw;
  cout<<"Opening Input Weight Sphere :"<<endl
      <<"          Type of map : "<<sphinw.TypeOfMap()<<endl
      <<"     Number of pixels : "<<sphinw.NbPixels()<<endl
      <<"                 Nlat : "<<sphinw.SizeIndex()<<endl;
  if(sphinw.SizeIndex()!=nlat) {
    cout<<"Incompatible sphin sphinw"<<endl;
    exit(2);
  }
}

SphereHEALPix<r_8> sphred;
{FitsInFile sfits(fsphred); sfits >> sphred;}
cout<<"Opening Input Reducted Sphere :"<<endl
    <<"          Type of map : "<<sphred.TypeOfMap()<<endl
    <<"     Number of pixels : "<<sphred.NbPixels()<<endl
    <<"                 Nlat : "<<sphred.SizeIndex()<<endl;

// Filling hole for Output Sphere
cout<<"...Filling hole for Output Sphere"<<endl;
uint_4 n=0, n0=0;
for(int_4 i=0;i<sphin.NbPixels();i++) {
  bool skp;
  if(fsphinw) if(sphinw(i)<=0.) continue; // Pixel out of acceptance!

  // Test if pixel of Sphere has to be extrapolated
  skp = (tstmin || tstmax) ? neg : false;
  if((tstmin && sphin(i)<valmin) || (tstmax && sphin(i)>valmax)) skp = !neg;
  if(skp) continue;   // Do nothing

  double theta,phi;
  sphin.PixThetaPhi(i,theta,phi);
  int_4 ir = sphred.PixIndexSph(theta,phi);

  // Test if pixel of Reducted Sphere is good
  skp = (rtstmin || rtstmax) ? rneg : false;
  if((rtstmin && sphred(i)<rvalmin) || (rtstmax && sphred(i)>rvalmax)) skp = !rneg;
  if(skp) continue;   // Do nothing

  sphin(i) = sphred(ir);
  n++; if(sphred(ir)!=0.) n0++;
}
cout<<"      Number of pixels filled "
    <<n<<" (not null values "<<n0<<")"<<endl;

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

exit(0);
}