hotspots_cxx.cc

00001 /*
00002  *  This file is part of Healpix_cxx.
00003  *
00004  *  Healpix_cxx is free software; you can redistribute it and/or modify
00005  *  it under the terms of the GNU General Public License as published by
00006  *  the Free Software Foundation; either version 2 of the License, or
00007  *  (at your option) any later version.
00008  *
00009  *  Healpix_cxx is distributed in the hope that it will be useful,
00010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *  GNU General Public License for more details.
00013  *
00014  *  You should have received a copy of the GNU General Public License
00015  *  along with Healpix_cxx; if not, write to the Free Software
00016  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017  *
00018  *  For more information about HEALPix, see http://healpix.jpl.nasa.gov
00019  */
00020 
00021 /*
00022  *  Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
00023  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
00024  *  (DLR).
00025  */
00026 
00027 /*
00028  *  Copyright (C) 2003, 2004 Max-Planck-Society
00029  *  Author: Martin Reinecke
00030  */
00031 
00032 #include <fstream>
00033 #include "cxxutils.h"
00034 #include "paramfile.h"
00035 #include "simparams.h"
00036 #include "healpix_map.h"
00037 #include "healpix_map_fitsio.h"
00038 #include "fitshandle.h"
00039 
00040 using namespace std;
00041 
00042 int main (int argc, const char **argv)
00043   {
00044 PLANCK_DIAGNOSIS_BEGIN
00045   module_startup ("hotspots_cxx", argc, argv, 2, "<parameter file>");
00046   paramfile params (argv[1]);
00047   simparams par;
00048 
00049   string infile = params.find<string>("infile");
00050   string mapfile = params.find<string>("outmap","");
00051   bool have_mapfile = mapfile!="";
00052   string minfile = params.find<string>("minfile","");
00053   bool have_minfile = minfile!="";
00054   string maxfile = params.find<string>("maxfile","");
00055   bool have_maxfile = maxfile!="";
00056   planck_assert (have_mapfile || have_minfile || have_maxfile,
00057     "no output file specified");
00058 
00059   Healpix_Map<float> inmap;
00060   read_Healpix_map_from_fits(infile,inmap,1,2);
00061   Healpix_Map<float> outmap;
00062   if (have_mapfile) outmap.Set(inmap.Order(),inmap.Scheme());
00063 
00064   ofstream minout, maxout;
00065   if (have_minfile) minout.open(minfile.c_str());
00066   if (have_maxfile) maxout.open(maxfile.c_str());
00067 
00068   fix_arr<int,8> nb;
00069 // FIXME: This should be parallelized
00070   for (int m=0; m<inmap.Npix(); ++m)
00071     {
00072     float value = inmap[m];
00073     if (!approx<double>(value, Healpix_undef))
00074       {
00075       inmap.neighbors(m,nb);
00076       bool ismax=true, ismin=true;
00077       for (int n=0; n<nb.size(); ++n)
00078         {
00079         if (nb[n] >=0)
00080           {
00081           float nbval = inmap[nb[n]];
00082           if (!approx<double>(nbval, Healpix_undef))
00083             {
00084             if (nbval>=value) ismax=false;
00085             if (nbval<=value) ismin=false;
00086             }
00087           }
00088         }
00089       if (have_mapfile) outmap[m] = (ismax||ismin) ? value : Healpix_undef;
00090       if (have_minfile && ismin) minout << m << " " << value << endl;
00091       if (have_maxfile && ismax) maxout << m << " " << value << endl;
00092       }
00093     }
00094   if (have_mapfile)
00095     {
00096     fitshandle out;
00097     out.create (mapfile);
00098     write_Healpix_map_to_fits (out,outmap,TFLOAT);
00099     }
00100 PLANCK_DIAGNOSIS_END
00101   }