(file) Return to sharp.c CVS log (file) (dir) Up to [Development] / JSOC / proj / sharp / apps

Diff for /JSOC/proj/sharp/apps/sharp.c between version 1.16 and 1.17

version 1.16, 2013/07/04 02:14:14 version 1.17, 2013/08/13 01:35:07
Line 22 
Line 22 
  *              v0.2    Sep 04 2012  *              v0.2    Sep 04 2012
  *              v0.3    Dec 18 2012  *              v0.3    Dec 18 2012
  *              v0.4    Jan 02 2013  *              v0.4    Jan 02 2013
  *      v0.5    Jan 23 2012   *    v0.5  Jan 23 2013
    *              v0.6    Aug 12 2013
  *  *
  *      Notes:  *      Notes:
  *              v0.0  *              v0.0
Line 42 
Line 43 
  *              Took out convert_inplace(). Was causing all the images to be int  *              Took out convert_inplace(). Was causing all the images to be int
  *      v0.5  *      v0.5
  *      Corrected ephemeris keywords, added argument mInfo for setKeys()  *      Corrected ephemeris keywords, added argument mInfo for setKeys()
    *              v0.6
    *              Changes in remapping of bitmap and conf_disambig, now near neighbor without anti-aliasing
    *
  *  *
  *      Example:  *      Example:
  *      sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \  *      sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \
Line 239  int getEphemeris(DRMS_Record_t *inRec, s
Line 243  int getEphemeris(DRMS_Record_t *inRec, s
 void findCoord(struct mapInfo *mInfo); void findCoord(struct mapInfo *mInfo);
  
 /* Mapping function */ /* Mapping function */
 int performSampling(float *outData, float *inData, struct mapInfo *mInfo);  int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt);
  
 /* Performing local vector transformation */ /* Performing local vector transformation */
 void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo); void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo);
Line 322  char *CEASegs[] = {"magnetogram", "bitma
Line 326  char *CEASegs[] = {"magnetogram", "bitma
 /* ========================================================================================================== */ /* ========================================================================================================== */
  
 char *module_name = "sharp"; char *module_name = "sharp";
 char *version_id = "2013 Jun 26";  /* Version number */  
 int seed; int seed;
  
 ModuleArgs_t module_args[] = ModuleArgs_t module_args[] =
Line 599  int createCeaRecord(DRMS_Record_t *mharp
Line 602  int createCeaRecord(DRMS_Record_t *mharp
         mInfo.proj = (enum projection) cyleqa;          // projection method         mInfo.proj = (enum projection) cyleqa;          // projection method
         mInfo.xscale = XSCALE;         mInfo.xscale = XSCALE;
         mInfo.yscale = YSCALE;         mInfo.yscale = YSCALE;
         mInfo.nbin = NBIN;  
     int ncol0, nrow0;             // oversampled map size
  
         // Get ephemeris         // Get ephemeris
  
Line 615  int createCeaRecord(DRMS_Record_t *mharp
Line 619  int createCeaRecord(DRMS_Record_t *mharp
                 return 1;                 return 1;
         }         }
  
           // ========================================
           // Do this for all bitmaps, Aug 12 2013 XS
           // ========================================
   
     mInfo.nbin = 1;                       // for bitmaps. suppress anti-aliasing
           ncol0 = mInfo.ncol;
           nrow0 = mInfo.nrow;
   
           mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
           mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
   
           findCoord(&mInfo);              // compute it here so it could be shared by the following 4 functions
   
           if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
                   SHOW("CEA: mapping bitmap error\n");
                   return 1;
           }
           printf("Bitmap mapping done.\n");
   
     if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {
                   SHOW("CEA: mapping conf_disambig error\n");
                   return 1;
           }
           printf("Conf disambig mapping done.\n");
   
     free(mInfo.xi_out);
           free(mInfo.zeta_out);
   
           // ========================================
           // Do this again for floats, Aug 12 2013 XS
           // ========================================
         // Create xi_out, zeta_out array in mInfo:         // Create xi_out, zeta_out array in mInfo:
         // Coordinates to sample in original full disk image         // Coordinates to sample in original full disk image
  
         int ncol0, nrow0;               // oversampled map size          mInfo.nbin = NBIN;
         ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN         ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
         nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;         nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;
  
Line 633  int createCeaRecord(DRMS_Record_t *mharp
Line 668  int createCeaRecord(DRMS_Record_t *mharp
                 SHOW("CEA: mapping magnetogram error\n");                 SHOW("CEA: mapping magnetogram error\n");
                 return 1;                 return 1;
         }         }
         if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {  
                 SHOW("CEA: mapping bitmap error\n");  
                 return 1;  
         }  
         printf("Magnetogram mapping done.\n");         printf("Magnetogram mapping done.\n");
   
         if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {         if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {
                 SHOW("CEA: mapping dopplergram error\n");                 SHOW("CEA: mapping dopplergram error\n");
                 return 1;                 return 1;
Line 650  int createCeaRecord(DRMS_Record_t *mharp
Line 682  int createCeaRecord(DRMS_Record_t *mharp
         }         }
         printf("Intensitygram mapping done.\n");         printf("Intensitygram mapping done.\n");
  
         if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {  
                 SHOW("CEA: mapping conf_disambig error\n");  
                 return 1;  
         }  
         printf("Conf disambig mapping done.\n");  
   
         // Mapping vector B         // Mapping vector B
  
         if (mapVectorB(sharpRec, bharpRec, &mInfo)) {         if (mapVectorB(sharpRec, bharpRec, &mInfo)) {
Line 724  int mapScaler(DRMS_Record_t *sharpRec, D
Line 750  int mapScaler(DRMS_Record_t *sharpRec, D
         int status = 0;         int status = 0;
         int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;         int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
         int dims[2] = {nx, ny};         int dims[2] = {nx, ny};
           int interpOpt = INTERP;         // Aug 12 XS, default, overridden below for bitmaps and conf_disambig
  
         // Input full disk array         // Input full disk array
  
Line 738  int mapScaler(DRMS_Record_t *sharpRec, D
Line 765  int mapScaler(DRMS_Record_t *sharpRec, D
         float *inData;         float *inData;
         int xsz = inArray->axis[0], ysz = inArray->axis[1];         int xsz = inArray->axis[0], ysz = inArray->axis[1];
         if ((xsz != FOURK) || (ysz != FOURK)) {         // for bitmap, make tmp full disk         if ((xsz != FOURK) || (ysz != FOURK)) {         // for bitmap, make tmp full disk
                   interpOpt = 3;          // Aug 12 XS, near neighbor
                 float *inData0 = (float *) inArray->data;                 float *inData0 = (float *) inArray->data;
                 inData = (float *) (calloc(FOURK2, sizeof(float)));                 inData = (float *) (calloc(FOURK2, sizeof(float)));
                 int x0 = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1;                 int x0 = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1;
Line 757  int mapScaler(DRMS_Record_t *sharpRec, D
Line 785  int mapScaler(DRMS_Record_t *sharpRec, D
         // Mapping         // Mapping
  
         float *map = (float *) (malloc(nxny * sizeof(float)));         float *map = (float *) (malloc(nxny * sizeof(float)));
         if (performSampling(map, inData, mInfo))          if (performSampling(map, inData, mInfo, interpOpt))             // Add interpOpt for different types, Aug 12 XS
         {if (inArray) drms_free_array(inArray); free(map); return 1;}         {if (inArray) drms_free_array(inArray); free(map); return 1;}
  
         // Write out         // Write out
Line 820  int mapVectorB(DRMS_Record_t *sharpRec,
Line 848  int mapVectorB(DRMS_Record_t *sharpRec,
         float *bx_map = NULL, *by_map = NULL, *bz_map = NULL;   // intermediate maps, in CCD bxyz representation         float *bx_map = NULL, *by_map = NULL, *bz_map = NULL;   // intermediate maps, in CCD bxyz representation
  
         bx_map = (float *) (malloc(nxny * sizeof(float)));         bx_map = (float *) (malloc(nxny * sizeof(float)));
         if (performSampling(bx_map, bx_img, mInfo))          if (performSampling(bx_map, bx_img, mInfo, INTERP))
         {free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}         {free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
  
         by_map = (float *) (malloc(nxny * sizeof(float)));         by_map = (float *) (malloc(nxny * sizeof(float)));
         if (performSampling(by_map, by_img, mInfo))          if (performSampling(by_map, by_img, mInfo, INTERP))
         {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}         {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
  
         bz_map = (float *) (malloc(nxny * sizeof(float)));         bz_map = (float *) (malloc(nxny * sizeof(float)));
         if (performSampling(bz_map, bz_img, mInfo))          if (performSampling(bz_map, bz_img, mInfo, INTERP))
         {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}         {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
  
         free(bx_img); free(by_img); free(bz_img);         free(bx_img); free(by_img); free(bz_img);
Line 1107  void findCoord(struct mapInfo *mInfo)
Line 1135  void findCoord(struct mapInfo *mInfo)
  *  *
  */  */
  
 int performSampling(float *outData, float *inData, struct mapInfo *mInfo)  int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt)
 { {
  
         int status = 0;         int status = 0;
           int ind_map;
  
         int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2;  // pad with nbin/2 on edge to avoid NAN         int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2;  // pad with nbin/2 on edge to avoid NAN
         int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;         int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
  
         float *outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));          // Changed Aug 12 2013, XS, for bitmaps
           float *outData0;
           if (interpOpt == 3 && mInfo->nbin == 1) {
             outData0 = outData;
           } else {
             outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
           }
  
         float *xi_out = mInfo->xi_out;         float *xi_out = mInfo->xi_out;
         float *zeta_out = mInfo->zeta_out;         float *zeta_out = mInfo->zeta_out;
Line 1123  int performSampling(float *outData, floa
Line 1158  int performSampling(float *outData, floa
         // Interpolation         // Interpolation
  
         struct fint_struct pars;         struct fint_struct pars;
         int interpOpt = INTERP;         // Use Wiener by default, 6 order, 1 constraint          // Aug 12 2013, passed in as argument now
  
         switch (interpOpt) {         switch (interpOpt) {
                 case 0:                 // Wiener, 6 order, 1 constraint                 case 0:                 // Wiener, 6 order, 1 constraint
Line 1135  int performSampling(float *outData, floa
Line 1170  int performSampling(float *outData, floa
                 case 2:                 // Bilinear                 case 2:                 // Bilinear
                         init_finterpolate_linear(&pars, 1.);                         init_finterpolate_linear(&pars, 1.);
                         break;                         break;
                   case 3:                 // Near neighbor
                     break;
                 default:                 default:
                         return 1;                         return 1;
         }         }
  
           printf("interpOpt = %d, nbin = %d ", interpOpt, mInfo->nbin);
           if (interpOpt == 3) {                   // Aug 6 2013, Xudong
                   for (int row0 = 0; row0 < nrow0; row0++) {
                                   for (int col0 = 0; col0 < ncol0; col0++) {
                                           ind_map = row0 * ncol0 + col0;
                                           outData0[ind_map] = nnb(inData, FOURK, FOURK, xi_out[ind_map], zeta_out[ind_map]);
                                   }
                           }
           } else {
         finterpolate(&pars, inData, xi_out, zeta_out, outData0,         finterpolate(&pars, inData, xi_out, zeta_out, outData0,
                                  FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);                                  FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
           }
  
         // Rebinning, smoothing         // Rebinning, smoothing
  
           if (interpOpt == 3 && mInfo->nbin == 1) {
             return 0;
           } else {
         frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian         frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1);                // Gaussian
         free(outData0);         // Dec 18 2012         free(outData0);         // Dec 18 2012
           }
  
         //         //
  
Line 2021  void frebin (float *image_in, float *ima
Line 2072  void frebin (float *image_in, float *ima
         fresize(&fresizes, image_in, image_out, nx, ny, nlead, nxout, nyout, nxout, xoff, yoff, DRMS_MISSING_FLOAT);         fresize(&fresizes, image_in, image_out, nx, ny, nlead, nxout, nyout, nxout, xoff, yoff, DRMS_MISSING_FLOAT);
  
 } }
   

Legend:
Removed from v.1.16  
changed lines
  Added in v.1.17
Karen Tian
 Powered by
ViewCVS 0.9.4