median_filter.cc
00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "cxxutils.h"
00033 #include "healpix_map.h"
00034 #include "healpix_map_fitsio.h"
00035 #include "fitshandle.h"
00036
00037 using namespace std;
00038
00039 template<typename Iterator> typename iterator_traits<Iterator>::value_type
00040 median(Iterator first, Iterator last)
00041 {
00042 Iterator mid = first+(last-first-1)/2;
00043 nth_element(first,mid,last);
00044 if ((last-first)&1) return *mid;
00045 return 0.5*((*mid)+(*min_element(mid+1,last)));
00046 }
00047
00048 int main (int argc, const char **argv)
00049 {
00050 module_startup ("median_filter", argc, argv, 4,
00051 "<input map> <output map> <radius in arcmin>");
00052
00053 double radius = stringToData<double>(argv[3])/60*degr2rad;
00054
00055 Healpix_Map<float> inmap;
00056 read_Healpix_map_from_fits(argv[1],inmap,1,2);
00057 Healpix_Map<float> outmap (inmap.Nside(), inmap.Scheme(), SET_NSIDE);
00058
00059 vector<int> list;
00060 vector<float> list2;
00061 for (int m=0; m<inmap.Npix(); ++m)
00062 {
00063 inmap.query_disc(inmap.pix2ang(m),radius,list);
00064 list2.resize(list.size());
00065 for (unsigned int i=0; i<list.size(); ++i)
00066 list2[i] = inmap[list[i]];
00067 outmap[m]=median(list2.begin(),list2.end());
00068 }
00069
00070 fitshandle out;
00071 out.create(argv[2]);
00072 write_Healpix_map_to_fits (out,outmap,TFLOAT);
00073 }