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 <sstream>
00033 #include <iomanip>
00034 #include <cstdlib>
00035 #include "healpix_map.h"
00036 #include "healpix_map_fitsio.h"
00037 #include "rotmatrix.h"
00038 #include "pointing.h"
00039 #include "tga_image.h"
00040 #include "paramfile.h"
00041
00042 using namespace std;
00043
00044 void histo_eq (arr2<float> &img, float &minv, float &maxv, arr<double> &newpos)
00045 {
00046 const int nbins=100;
00047 arr<int> bincnt (nbins);
00048 bincnt.fill(0);
00049 int pixels=0;
00050
00051 double fact = 1./(maxv-minv);
00052 for (int i=0; i<img.size1(); ++i)
00053 for (int j=0; j<img.size2(); ++j)
00054 if (img[i][j]>-1e30)
00055 {
00056 img[i][j] = (img[i][j]-minv)*fact;
00057 int idx = int(img[i][j]*nbins);
00058 idx=max(0,min(idx,nbins-1));
00059 ++bincnt[idx];
00060 ++pixels;
00061 }
00062
00063 newpos.alloc(nbins+1);
00064 int accu=0;
00065 for (int m=0; m<nbins; ++m)
00066 {
00067 newpos[m] = double(accu)/pixels;
00068 accu += bincnt[m];
00069 }
00070 newpos[nbins]=1.;
00071
00072 for (int i=0; i<img.size1(); ++i)
00073 for (int j=0; j<img.size2(); ++j)
00074 if (img[i][j]>-1e30)
00075 {
00076 int idx = int(img[i][j]*nbins);
00077 idx=max(0,min(idx,nbins-1));
00078 double frac = nbins*img[i][j] - idx;
00079 img[i][j] = (1-frac)*newpos[idx] + frac*newpos[idx+1];
00080 img[i][j] = minv+(maxv-minv)*img[i][j];
00081 }
00082 }
00083
00084 void pro_mollw (const Healpix_Map<float> &map, double lon0, double lat0,
00085 int xsize, arr2<float> &img, float &minv, float &maxv, bool smooth)
00086 {
00087 int ysize=xsize/2;
00088 img.alloc(xsize,ysize);
00089 img.fill(-1e35);
00090 double xc=(xsize-1)/2., yc=(ysize-1)/2.;
00091 double lon0rad = lon0*degr2rad;
00092 double lat0rad = lat0*degr2rad;
00093
00094 rotmatrix rot;
00095 rot.Make_CPAC_Euler_Matrix(0,-lat0rad,-lon0rad);
00096
00097 minv=1e30;
00098 maxv=-1e30;
00099 for (int i=0; i<img.size1(); ++i)
00100 for (int j=0; j<img.size2(); ++j)
00101 {
00102 double u = 2*(i-xc)/(xc/1.02);
00103 double v = (j-yc)/(yc/1.02);
00104 bool mask = ((u*u/4 + v*v) <= 1);
00105 if (mask)
00106 {
00107 pointing ptg (halfpi-(asin(2/pi*(asin(v) + v*sqrt((1-v)*(1+v))))),
00108 -halfpi*u/max(sqrt((1-v)*(1+v)),1e-6));
00109 vec3 pnt = rot.Transform(ptg.to_vec3());
00110 if (smooth)
00111 img[i][j] = map.interpolated_value(pnt);
00112 else
00113 img[i][j]=map[map.ang2pix(pnt)];
00114 if (!approx<double>(img[i][j],Healpix_undef))
00115 {
00116 if (img[i][j]<minv) minv=img[i][j];
00117 if (img[i][j]>maxv) maxv=img[i][j];
00118 }
00119 }
00120 }
00121 }
00122
00123 void pro_gno (const Healpix_Map<float> &map, double lon0, double lat0,
00124 int xsize, double resgrid, arr2<float> &img, float &minv, float &maxv,
00125 bool smooth)
00126 {
00127 double lon0rad = lon0*degr2rad;
00128 double lat0rad = lat0*degr2rad;
00129
00130 rotmatrix rot;
00131 rot.Make_CPAC_Euler_Matrix(lon0rad,-lat0rad,0);
00132
00133 double delta=resgrid*degr2rad/60.;
00134 double start=-(xsize/2.)*delta;
00135 img.alloc(xsize,xsize);
00136 minv=1e30;
00137 maxv=-1e30;
00138 for (int i=0; i<img.size1(); ++i)
00139 for (int j=0; j<img.size2(); ++j)
00140 {
00141 vec3 pnt (1,-(start+i*delta), start+j*delta);
00142 pnt = rot.Transform(pnt);
00143 if (smooth)
00144 img[i][j] = map.interpolated_value(pnt);
00145 else
00146 img[i][j]=map[map.ang2pix(pnt)];
00147 if (!approx<double>(img[i][j],Healpix_undef))
00148 {
00149 if (img[i][j]<minv) minv=img[i][j];
00150 if (img[i][j]>maxv) maxv=img[i][j];
00151 }
00152 }
00153 }
00154
00155 void colorbar (TGA_Image &img, const Palette &pal, int xmin, int xmax,
00156 int ymin, int ymax, bool flippal, const arr<double> &newpos)
00157 {
00158 int nbins = newpos.size()-1;
00159 for (int i=xmin; i<=xmax; ++i)
00160 {
00161 double val = (double(i)-xmin)/(xmax-xmin);
00162 if (nbins>0)
00163 {
00164 int idx = int(val*nbins);
00165 idx=max(0,min(idx,nbins-1));
00166 double frac = nbins*val - idx;
00167 val = (1-frac)*newpos[idx] + frac*newpos[idx+1];
00168 }
00169 if (flippal) val=1-val;
00170 Colour c = pal.Get_Colour(val);
00171 for (int j=ymin; j<=ymax; ++j)
00172 img.put_pixel(i,j,c);
00173 }
00174 }
00175
00176 string conv (double val)
00177 {
00178 ostringstream os;
00179 if (abs(val)>100 || abs(val)<0.01)
00180 {
00181 os << setw(10) << setprecision(3) << scientific << val;
00182 return os.str();
00183 }
00184 os << setw(10) << setprecision(6) << fixed << val;
00185 return trim(os.str());
00186 }
00187
00188 void usage()
00189 {
00190 cout <<
00191 "\nUsage:\n"
00192 " map2tga <parameter file>\n\n"
00193 "or:\n"
00194 " map2tga <input file> <output file> [-sig <int>] [-pal <int>]\n"
00195 " [-xsz <int>] [-bar] [-log] [-asinh] [-lon <float>] [-lat <float>]\n"
00196 " [-mul <float>] [-add <float>] [-min <float>] [-max <float>]\n"
00197 " [-res <float>] [-title <string>] [-flippal] [-gnomonic]\n"
00198 " [-interpol] [-equalize] [-viewer <viewer>]\n\n";
00199 throw Message_error();
00200 }
00201
00202 int main (int argc, const char **argv)
00203 {
00204 PLANCK_DIAGNOSIS_BEGIN
00205 announce ("map2tga");
00206 if (argc<2) usage();
00207
00208 string infile = "";
00209 string outfile = "";
00210 string title = "";
00211 int xres = 1024;
00212 bool bar = false, mollpro=true;
00213 double lat0=0, lon0=0, res=1;
00214 double usermin=-1e30, usermax=1e30, offset=0, factor=1;
00215 bool min_supplied=false, max_supplied=false;
00216 bool logflag=false, asinhflag=false, eqflag=false;
00217 int palnr = 4;
00218 int colnum=1;
00219 bool flippal = false;
00220 bool interpol = false;
00221 string viewer = "";
00222
00223 if (argc>2)
00224 {
00225 infile = argv[1];
00226 outfile = argv[2];
00227 int m=3;
00228 while (m<argc)
00229 {
00230 int mstart = m;
00231 string arg = argv[m];
00232 if (arg == "-sig") { stringToData(argv[m+1],colnum); m+=2; }
00233 if (arg == "-pal") { stringToData(argv[m+1],palnr); m+=2; }
00234 if (arg == "-xsz") { stringToData(argv[m+1],xres); m+=2; }
00235 if (arg == "-bar") { bar=true; ++m; }
00236 if (arg == "-log") { logflag=true; ++m; }
00237 if (arg == "-equalize") { eqflag=true; ++m; }
00238 if (arg == "-asinh") { asinhflag=true; ++m; }
00239 if (arg == "-lon") { stringToData(argv[m+1],lon0); m+=2; }
00240 if (arg == "-lat") { stringToData(argv[m+1],lat0); m+=2; }
00241 if (arg == "-mul") { stringToData(argv[m+1],factor); m+=2; }
00242 if (arg == "-add") { stringToData(argv[m+1],offset); m+=2; }
00243 if (arg == "-min")
00244 {
00245 stringToData(argv[m+1],usermin);
00246 min_supplied=true;
00247 m+=2;
00248 }
00249 if (arg == "-max")
00250 {
00251 stringToData(argv[m+1],usermax);
00252 max_supplied=true;
00253 m+=2;
00254 }
00255 if (arg == "-res") { stringToData(argv[m+1],res); m+=2; }
00256 if (arg == "-title") { title = argv[m+1]; m+=2; }
00257 if (arg == "-viewer") { viewer = argv[m+1]; m+=2; }
00258 if (arg == "-flippal") { flippal=true; ++m; }
00259 if (arg == "-gnomonic") { mollpro=false; ++m; }
00260 if (arg == "-interpol") { interpol=true; ++m; }
00261
00262 if (mstart==m)
00263 {
00264 cout << "unrecognized option: " + arg << endl;
00265 usage();
00266 }
00267 }
00268 }
00269 else
00270 {
00271 paramfile params (argv[1]);
00272 infile = params.find<string>("infile");
00273 outfile = params.find<string>("outfile");
00274 colnum = params.find<int>("sig",colnum);
00275 palnr = params.find<int>("pal",palnr);
00276 flippal = params.find<bool>("flippal",flippal);
00277 xres = params.find<int>("xsz",xres);
00278 bar = params.find<bool>("bar",bar);
00279 logflag = params.find<bool>("log",logflag);
00280 eqflag = params.find<bool>("equalize",eqflag);
00281 asinhflag = params.find<bool>("asinh",asinhflag);
00282 lon0 = params.find<double>("lon",lon0);
00283 lat0 = params.find<double>("lat",lat0);
00284 factor = params.find<double>("mul",factor);
00285 offset = params.find<double>("add",offset);
00286 usermin = params.find<double>("min", usermin);
00287 if (usermin>-1e29) min_supplied = true;
00288 usermax = params.find<double>("max", usermax);
00289 if (usermax<1e29) max_supplied = true;
00290 res = params.find<double>("res",res);
00291 title = params.find<string>("title","");
00292 viewer = params.find<string>("viewer","");
00293 string tmp = params.find<string>("pro","");
00294 if (tmp == "gno") mollpro=false;
00295 interpol = params.find<bool>("interpol");
00296 }
00297
00298 Healpix_Map<float> map(0,RING);
00299 read_Healpix_map_from_fits(infile,map,colnum,2);
00300 for (int m=0; m<map.Npix(); ++m)
00301 {
00302 if (!approx<double>(map[m],Healpix_undef))
00303 {
00304 map[m] = (map[m]+offset)*factor;
00305 if (logflag)
00306 {
00307 if (map[m]<=0)
00308 map[m] = Healpix_undef;
00309 else
00310 map[m] = log(double(map[m]))/ln10;
00311 }
00312 if (asinhflag)
00313 {
00314 if (map[m]>=0)
00315 map[m] = log(double(map[m]+sqrt(map[m]*map[m]+1)));
00316 else
00317 map[m] = -log(double(-map[m]+sqrt(map[m]*map[m]+1)));
00318 }
00319 if (min_supplied) if (map[m] < usermin) map[m] = usermin;
00320 if (max_supplied) if (map[m] > usermax) map[m] = usermax;
00321 }
00322 }
00323
00324 arr2<float> imgarr;
00325 float minv, maxv;
00326 if (mollpro)
00327 pro_mollw (map,lon0,lat0,xres,imgarr,minv,maxv,interpol);
00328 else
00329 pro_gno (map,lon0,lat0,xres,res,imgarr,minv,maxv,interpol);
00330
00331 arr<double> newpos;
00332 if (eqflag) histo_eq(imgarr,minv,maxv,newpos);
00333
00334 if (min_supplied) minv = usermin;
00335 if (max_supplied) maxv = usermax;
00336 if (maxv==minv) maxv=minv+1e-10;
00337
00338 int xsz=imgarr.size1();
00339 int ysz=imgarr.size2();
00340 int yofs=ysz;
00341 int scale = max(1,xsz/800);
00342 if (bar) ysz+=60*scale;
00343 if (title != "") { ysz+=50*scale; yofs+=50*scale; }
00344 TGA_Image img(xsz,ysz);
00345 img.fill(Colour(1,1,1));
00346 img.set_font (giant_font);
00347 Palette pal;
00348 pal.setPredefined(palnr);
00349 if (title != "")
00350 img.annotate_centered(xsz/2,25*scale,Colour(0,0,0),title,scale);
00351 for (int i=0; i<imgarr.size1(); ++i)
00352 for (int j=0; j<imgarr.size2(); ++j)
00353 {
00354 if (imgarr[i][j]>-1e32)
00355 {
00356 Colour c(0.5,0.5,0.5);
00357 if (!approx<double>(imgarr[i][j],Healpix_undef))
00358 {
00359 int col = int((imgarr[i][j]-minv)/(maxv-minv)*256);
00360 col = min(255,max(col,0));
00361 float colfrac = (imgarr[i][j]-minv)/(maxv-minv);
00362 if (flippal) colfrac = 1-colfrac;
00363 c = pal.Get_Colour(colfrac);
00364 }
00365 img.put_pixel(i,yofs-j-1,c);
00366 }
00367 }
00368
00369 if (bar)
00370 {
00371 colorbar (img,pal,xsz/10,(xsz*9)/10,ysz-40*scale,ysz-20*scale,flippal,
00372 newpos);
00373 img.set_font (medium_bold_font);
00374 img.annotate_centered (xsz/20,ysz-30*scale,Colour(0,0,0),conv(minv),scale);
00375 img.annotate_centered ((xsz*19)/20,ysz-30*scale,Colour(0,0,0),conv(maxv),
00376 scale);
00377 }
00378 img.write(outfile);
00379
00380 if (viewer!="")
00381 system ((viewer+" "+outfile).c_str());
00382 PLANCK_DIAGNOSIS_END
00383 }
