JSOC_Documentation: proj/lev0/apps/add_small

00001 #include "mypng.h"
00002 
00003 #define PNGDIE(msg) {fprintf(stderr,"%s\n",msg); \
00004     if (png_ptr) png_destroy_write_struct(&png_ptr, &info_ptr); \
00005     if (row_pointers) free(row_pointers); \
00006     if (fp) fclose(fp); \
00007     return(1);}
00008 
00009 #define GREY    1
00010 #define MINMAX  1
00011 #define MAXMIN  2
00012 #define HISTEQ  3
00013 #define MINMAX99  4
00014 #define MAXMIN99  5
00015 
00016 int make_png(char *filename, unsigned char *data, int height, int width, int pallette, int bytepercolor, int colors)
00017   { // make png
00018   int row;
00019   FILE *fp = NULL;
00020   png_byte **row_pointers = NULL;
00021   png_structp png_ptr = NULL;
00022   png_infop info_ptr = NULL;
00023 
00024   fp = fopen(filename, "wb");
00025   if (!fp)
00026     PNGDIE("cant open output file");
00027   png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
00028   if (!png_ptr)
00029     PNGDIE("cant create png_struct");
00030   info_ptr = png_create_info_struct(png_ptr);
00031   if (!info_ptr)
00032     PNGDIE("cant create png_infop");
00033   if (setjmp(png_jmpbuf(png_ptr)))
00034     PNGDIE("png fatal error");
00035   png_init_io(png_ptr, fp);
00036   png_set_IHDR(png_ptr, info_ptr, width, height, 8*bytepercolor,
00037     PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE,
00038     PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
00039   png_write_info(png_ptr, info_ptr);
00040   row_pointers = (png_byte **)malloc(height * sizeof(png_byte *));
00041   for (row=0; row<height ; row++)
00042     {
00043     if (bytepercolor == 1)
00044       row_pointers[height - row - 1] = (png_byte *)((char *)data + row*width);
00045     else
00046       row_pointers[height - row - 1] = (png_byte *)((short *)data + row*width);
00047     }
00048   png_set_rows(png_ptr, info_ptr, row_pointers);
00049   png_write_png(png_ptr, info_ptr, PNG_TRANSFORM_SWAP_ENDIAN, NULL);
00050   fclose(fp);
00051   png_destroy_write_struct(&png_ptr, &info_ptr);
00052   free(row_pointers);
00053   return(0);
00054   }
00055 
00056 char *set_scaling(DRMS_Array_t *in, double *minp, double *maxp, int *nmissp,
00057         int pallette, int colors, int bytepercolor, int scaling)
00058   {
00059   char *out;
00060   int idata, ndata, nmiss = 0;
00061   float *inData = in->data;
00062   ndata = in->axis[0] * in->axis[1];
00063   out = (char *)malloc(ndata * bytepercolor * sizeof(char));
00064   int missingcolor = (bytepercolor == 1 ? 255 : 65535);
00065   int maxcolor = missingcolor - 1;
00066   if (!out)
00067     return(NULL);
00068 
00069   switch(scaling)
00070     {
00071     case MINMAX:
00072     case MAXMIN:
00073     {
00074     double min = 1.0e20;
00075     double max = -min;
00076     if (bytepercolor != 1 && bytepercolor != 2)
00077         {fprintf(stderr,"wrong bytes per color\n"); return(NULL);}
00078     for (idata=0; idata<ndata; idata++)
00079       {
00080       float val = inData[idata];
00081       if (!isnan(val))
00082         {
00083             if (val > max) max = val;
00084             if (val < min) min = val;
00085             }
00086           }
00087         for (idata=0; idata<ndata; idata++)
00088           {
00089           float val = inData[idata];
00090       unsigned short newval;
00091           if (!isnan(val))
00092         {
00093             newval = maxcolor*(val-min)/(max-min);
00094         if (scaling == MAXMIN)
00095         newval = maxcolor - newval;
00096         }
00097           else
00098         {
00099             newval = missingcolor;
00100         nmiss += 1;
00101         }
00102       if (bytepercolor == 1)
00103         *((unsigned char *)out + idata) = newval;
00104       else
00105             *((unsigned short *)out + idata) = newval;
00106           }
00107     *maxp = max;
00108     *minp = min;
00109     *nmissp = nmiss;
00110     break;
00111     }
00112     case HISTEQ:
00113     {
00114     int ihist, nhist=65536;
00115     int *hist = (int *)malloc(nhist * sizeof(int));
00116     int icolor, ncolor = (bytepercolor == 1 ? 255 : 65535);
00117     int total, want;
00118     unsigned int newval;
00119     double valspercolor;
00120     char *outa;
00121     outa = set_scaling(in, minp, maxp, nmissp, GREY, 1, 2, MINMAX);
00122     valspercolor = (ndata - *nmissp)/(double)ncolor;
00123     for (ihist=0; ihist<nhist; ihist++)
00124       hist[ihist] = 0;
00125     for (idata=0; idata<ndata; idata++)
00126       hist[*((unsigned short *)outa+idata)] += 1;
00127     for (icolor=0, total=0, ihist=0; icolor<=ncolor; icolor++)
00128       {
00129       want = icolor * valspercolor;
00130       for ( ; ihist<nhist; ihist++)
00131         {
00132         if (total >= want)
00133         break;
00134         total += hist[ihist];
00135         hist[ihist] = icolor;
00136         }
00137       }
00138     for ( ; ihist<nhist; ihist++)
00139         hist[ihist] = ncolor;
00140     for (idata=0; idata<ndata; idata++)
00141           if (!isnan(inData[idata]))
00142         {
00143         if (bytepercolor == 1)
00144           *((unsigned char *)out + idata) = hist[*((unsigned short *)outa+idata)];
00145         else
00146               *((unsigned short *)out + idata) = hist[*((unsigned short *)outa+idata)];
00147         }
00148           else
00149             {
00150         if (bytepercolor == 1)
00151           *((unsigned char *)out + idata) = missingcolor;
00152         else
00153               *((unsigned short *)out + idata) = missingcolor;
00154             }
00155     free(outa);  
00156     free(hist);
00157     break;
00158     }
00159     case MINMAX99:
00160     case MAXMIN99:
00161     {
00162     int ihist, nhist=65536;
00163     int *hist = (int *)malloc(nhist * sizeof(int));
00164     int icolor, ncolor = (bytepercolor == 1 ? 255 : 65535);
00165     int total, want;
00166     unsigned int newval;
00167     double newmin, newmax;
00168     char *outa;
00169     outa = set_scaling(in, minp, maxp, nmissp, GREY, 1, 2, MINMAX);
00170 
00171         // create histogram
00172     for (ihist=0; ihist<nhist; ihist++)
00173       hist[ihist] = 0;
00174     for (idata=0; idata<ndata; idata++)
00175       hist[*((unsigned short *)outa+idata)] += 1;
00176 
00177         // find new min and max vals at cut points
00178         want = ndata/200;
00179         total = 0;
00180     for (ihist=0; ihist<nhist; total += hist[ihist++])
00181       if (total >= want)
00182         break;
00183         newmin = *minp + ((double)ihist/(double)nhist) * (*maxp - *minp);
00184         total = 0;
00185     for (ihist=nhist; ihist>0; total += hist[--ihist])
00186       if (total >= want)
00187         break;
00188         newmax = *minp + ((double)ihist/(double)nhist) * (*maxp - *minp);
00189         if (newmax <= newmin)
00190           {
00191           newmax = *maxp;
00192           newmin = *minp;
00193           }
00194 
00195         for (idata=0; idata<ndata; idata++)
00196           {
00197           float val = inData[idata];
00198           unsigned short newval;
00199           if (!isnan(val))
00200             {
00201         if (val <= newmin) val=newmin;
00202         if (val >= newmax) val=newmax;
00203             newval = maxcolor*(val-newmin)/(newmax-newmin);
00204             if (scaling == MAXMIN99)
00205                 newval = maxcolor - newval;
00206             }
00207           else
00208             newval = missingcolor;
00209           if (bytepercolor == 1)
00210             *((unsigned char *)out + idata) = newval;
00211           else
00212             *((unsigned short *)out + idata) = newval;
00213           }
00214     free(outa);  
00215     free(hist);
00216     break;
00217     }
00218     }
00219   return(out);
00220   }
00221 
00222 // expects short in and returns float DRMS_Array_t's
00223 
00224 rebinArraySF(DRMS_Array_t *out, DRMS_Array_t *in)
00225   {
00226   short *inData = in->data;
00227   short *inp;
00228   float *outData = out->data;
00229   float *outp;
00230   int inNx=in->axis[1], inNy=in->axis[0], outNx=out->axis[1], outNy=out->axis[0];
00231   int inI, inJ, outI, outJ, i;
00232   int binsize = inNx/outNx;
00233   double inRow[inNx], outRow[outNx];
00234   int inRowN[inNx], outRowN[outNx];
00235   for (outJ=0; outJ<outNy; outJ++)
00236     {
00237     int inRow0 = outJ * binsize;
00238     for (outI=0; outI<outNx; outI++)
00239       {
00240       outRowN[outI] = 0;
00241       outRow[outI] = 0.0;
00242       }
00243     for (inJ = inRow0; inJ < inRow0+binsize; inJ++)
00244       {
00245       inp = inData + inJ*inNx;
00246       for (outI=0; outI<outNx; outI++)
00247     for (i=0; i<binsize; i++, inp++)
00248       if (*inp != DRMS_MISSING_SHORT)
00249         {
00250             outRow[outI] += *inp;
00251         outRowN[outI]++;
00252         }
00253       }
00254     for (outI=0; outI<outNx; outI++)
00255       outData[outJ*outNx + outI] = (outRowN[outI] ? (outRow[outI]/outRowN[outI]) : DRMS_MISSING_FLOAT);
00256     }
00257   }
00258 
00259 // Given rec and big, create smaller scale segment reduced by doScaleFits in each
00260 // dimension and create a smaller scale image reduced by doScaleImage.
00261 // Setting either doScale{Fits,Image} to 0 will prevent that component
00262 // from being generated.  Typical use, to convert a 4096x4096 image to
00263 // a 512x512 fits file and 256x256 pixel image file call with 
00264 //  add_small_array(rs, array, 8, 16);
00265 // Scaling for the image will be equal spaced by histogram
00266 // Filenames will be put into segments number 1 and 2 respectively if
00267 // the segments are in the .jsd.
00268 // segment 1 should be a vardim FITS segment and segment 2 should
00269 // be GENERIC.  The
00270 // small FITS and image files will be placed in the record directory
00271 // with filenames derived from segment number 0.
00272 //
00273 // return value is 0 if all OK, else error code.
00274 
00275 // Now expects data in "big" to be shorts.
00276 
00277 int add_small_array(DRMS_Record_t *rec, DRMS_Array_t *big, int doScaleFits, int doScaleImage)
00278   {
00279   DRMS_Segment_t *bigSeg = drms_segment_lookupnum(rec, 0);
00280   char recdir[DRMS_MAXPATHLEN];
00281   int bigNx, bigNy;
00282   int status;
00283 
00284   drms_record_directory(rec, recdir, 1);
00285   bigNy = big->axis[0];
00286   bigNx = big->axis[1];
00287   /* assume second segment, if it is present, is vardims FITS for small image */
00288   if (doScaleFits && rec->segments.num_total > 1)
00289     {
00290     DRMS_Segment_t *smallSeg;
00291     DRMS_Array_t *small;
00292     DRMS_SegmentDimInfo_t smallDimInfo;
00293     int smallNx, smallNy;
00294     int smallDims[2];
00295     char smallFileName[DRMS_MAXSEGFILENAME];
00296 
00297     smallNy = bigNy/doScaleFits;
00298     smallNx = bigNx/doScaleFits;
00299     smallDims[0] = smallNy;
00300     smallDims[1] = smallNx;
00301     small = drms_array_create(DRMS_TYPE_FLOAT, 2, smallDims, NULL, &status);
00302 
00303     rebinArraySF(small, big);
00304     strcpy(smallFileName, bigSeg->info->name);
00305     strcat(smallFileName, "_sm.fits");
00306     smallSeg = drms_segment_lookupnum(rec, 1);
00307     strcpy(smallSeg->filename, smallFileName);
00308     smallDimInfo.naxis = 2;
00309     smallDimInfo.axis[0] = smallDims[0];
00310     smallDimInfo.axis[1] = smallDims[1];
00311     drms_segment_setdims(smallSeg, &smallDimInfo);
00312     // For SDO CCD images of 14 bits use FITS shorts as unsigned and keep full range.
00313     small->bzero = -(4.0*8192.0);
00314     small->bscale = 4.0;
00315     small->israw = 0;
00316     drms_array_convert_inplace(DRMS_TYPE_SHORT, small->bzero, small->bscale, small);
00317 
00318     small->bzero = 8192.0;
00319     small->bscale = 0.25;
00320     small->israw = 1;
00321     drms_segment_writewithkeys(smallSeg, small, 0);
00322 
00323     drms_free_array(small);
00324     }
00325 
00326   /* assume third segment, if exists, is for small image, generic segment */
00327   if (doScaleImage && rec->segments.num_total > 2)
00328     {
00329     DRMS_Segment_t *imageSeg;
00330     DRMS_Array_t *imageArray = NULL;
00331     char *imgData = NULL;
00332     int imageDims[2];
00333     int imageNx, imageNy;
00334     char imageFileName[DRMS_MAXSEGFILENAME];
00335     char imagePathName[DRMS_MAXPATHLEN];
00336     int bytepercolor = 1;
00337     int nmissing;
00338     double min, max;
00339     imageNy = bigNy/doScaleImage;
00340     imageNx = bigNx/doScaleImage;
00341     imageDims[0] = imageNy;
00342     imageDims[1] = imageNx;
00343     imageArray = drms_array_create(DRMS_TYPE_FLOAT, 2, imageDims, NULL, &status);
00344 
00345     strcpy(imageFileName, bigSeg->info->name);
00346     strcat(imageFileName, ".png");
00347     strcpy(imagePathName, recdir);
00348     strcat(imagePathName, "/");
00349     strcat(imagePathName, imageFileName);
00350     imageSeg = drms_segment_lookupnum(rec, 2);
00351     strcpy(imageSeg->filename, imageFileName);
00352     rebinArraySF(imageArray, big);
00353     imgData = (unsigned char *)set_scaling(imageArray, &min, &max, &nmissing, GREY, 1, bytepercolor, MINMAX99);
00354     drms_free_array(imageArray);
00355     if (!imgData)
00356       {
00357       fprintf(stderr,"scaling failed, status=%d\n",status);
00358       return(status);
00359       }
00360     // should add text info to .png file with min, max, nmissing, etc.
00361     if (make_png(imagePathName, imgData, imageNy, imageNx, GREY, 1, bytepercolor) != 0)
00362       {
00363       fprintf(stderr,"png failed, status=%d\n",status);
00364       free(imgData);
00365       return(status);
00366       }
00367     free(imgData);
00368     }
00369   return(0);
00370   }
00371
proj/lev0/apps/add_small_image.c
