(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.15 and 1.18

version 1.15, 2013/06/27 01:58:39 version 1.18, 2013/11/02 20:31:08
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 92 
Line 96 
 // Nyqvist rate at disk center is 0.03 degree. Oversample above 0.015 degree // Nyqvist rate at disk center is 0.03 degree. Oversample above 0.015 degree
 #define NYQVIST         (0.015) #define NYQVIST         (0.015)
  
   // Some other things
 #ifndef MIN #ifndef MIN
 #define MIN(a,b) (((a)<(b)) ? (a) : (b)) #define MIN(a,b) (((a)<(b)) ? (a) : (b))
 #endif #endif
Line 119 
Line 124 
 #define INTERP                  0 #define INTERP                  0
 #define dpath    "/home/jsoc/cvs/Development/JSOC" #define dpath    "/home/jsoc/cvs/Development/JSOC"
  
   
 /* ========================================================================================================== */ /* ========================================================================================================== */
  
 // Space weather keywords // Space weather keywords
Line 197  struct mapInfo {
Line 203  struct mapInfo {
         float *xi_out, *zeta_out;       // coordinate on full disk image to sample at         float *xi_out, *zeta_out;       // coordinate on full disk image to sample at
 }; };
  
   
 /* ========================================================================================================== */ /* ========================================================================================================== */
  
 /* Get all input data series */ /* Get all input data series */
Line 238  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 321  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 598  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 614  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 632  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 649  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 723  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 737  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 756  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 819  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 1003  int getEphemeris(DRMS_Record_t *inRec, s
Line 1032  int getEphemeris(DRMS_Record_t *inRec, s
         float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);         float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
         float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);         float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
         float cdelt = drms_getkey_float(inRec, "CDELT1", &status);  // in arcsec, assumimg dx=dy         float cdelt = drms_getkey_float(inRec, "CDELT1", &status);  // in arcsec, assumimg dx=dy
         printf("cdelt=%f\n",cdelt);  
         ephem->disk_xc = PIX_X(0.0,0.0) - 1.0;          // Center of disk in pixel, starting at 0         ephem->disk_xc = PIX_X(0.0,0.0) - 1.0;          // Center of disk in pixel, starting at 0
         ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;         ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
  
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 1722  void computeSWIndex(struct swIndex *swKe
Line 1773  void computeSWIndex(struct swIndex *swKe
         float *bx_err = (float *) bx_errArray->data;            // bx_err         float *bx_err = (float *) bx_errArray->data;            // bx_err
  
         // Get emphemeris         // Get emphemeris
           float  cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);
         //float cdelt1_orig = drms_getkey_float(inRec, "CDELT1",   &status);  
         float  cdelt1      = drms_getkey_float(inRec, "CDELT1",   &status);  
         float  dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);         float  dsun_obs    = drms_getkey_float(inRec, "DSUN_OBS",   &status);
         double rsun_ref   = drms_getkey_double(inRec, "RSUN_REF", &status);         double rsun_ref   = drms_getkey_double(inRec, "RSUN_REF", &status);
         double rsun_obs   = drms_getkey_double(inRec, "RSUN_OBS", &status);         double rsun_obs   = drms_getkey_double(inRec, "RSUN_OBS", &status);
Line 1733  void computeSWIndex(struct swIndex *swKe
Line 1782  void computeSWIndex(struct swIndex *swKe
         float crpix1      = drms_getkey_float(inRec, "CRPIX1", &status);         float crpix1      = drms_getkey_float(inRec, "CRPIX1", &status);
         float crpix2      = drms_getkey_float(inRec, "CRPIX2", &status);         float crpix2      = drms_getkey_float(inRec, "CRPIX2", &status);
  
         // Temp arrays          // convert cdelt1_orig from degrees to arcsec
           float cdelt1       = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
  
           // Temp arrays
         float *bh      = (float *) (malloc(nxny * sizeof(float)));         float *bh      = (float *) (malloc(nxny * sizeof(float)));
         float *bt      = (float *) (malloc(nxny * sizeof(float)));         float *bt      = (float *) (malloc(nxny * sizeof(float)));
         float *jz      = (float *) (malloc(nxny * sizeof(float)));         float *jz      = (float *) (malloc(nxny * sizeof(float)));
Line 1847  void computeSWIndex(struct swIndex *swKe
Line 1898  void computeSWIndex(struct swIndex *swKe
         swKeys_ptr->totpot_err         = DRMS_MISSING_FLOAT;         swKeys_ptr->totpot_err         = DRMS_MISSING_FLOAT;
         }         }
  
         if (computeShearAngle(bx_err, by_err, bh_err, bx, by, bz, bpx, bpy, bpz, dims,          if (computeShearAngle(bx_err, by_err, bz_err, bx, by, bz, bpx, bpy, bpz, dims,
                                                   &(swKeys_ptr->meanshear_angle), &(swKeys_ptr->meanshear_angle_err), &(swKeys_ptr->area_w_shear_gt_45),                                                   &(swKeys_ptr->meanshear_angle), &(swKeys_ptr->meanshear_angle_err), &(swKeys_ptr->area_w_shear_gt_45),
                                                   mask, bitmask)) {                                                   mask, bitmask)) {
                 swKeys_ptr->meanshear_angle    = DRMS_MISSING_FLOAT; // If fail, fill in NaN                 swKeys_ptr->meanshear_angle    = DRMS_MISSING_FLOAT; // If fail, fill in NaN
Line 1975  void setKeys(DRMS_Record_t *outRec, DRMS
Line 2026  void setKeys(DRMS_Record_t *outRec, DRMS
  
         // set cvs commit version into keyword HEADER         // set cvs commit version into keyword HEADER
         char *cvsinfo = strdup("$Id$");         char *cvsinfo = strdup("$Id$");
         //   status = drms_setkey_string(outRec, "HEADER", cvsinfo);          char *cvsinfo2 = sw_functions_version();
         status = drms_setkey_string(outRec, "CODEVER7", cvsinfo);          char cvsinfoall[2048];
           strcat(cvsinfoall,cvsinfo);
           strcat(cvsinfoall,"\n");
           strcat(cvsinfoall,cvsinfo2);
           status = drms_setkey_string(outRec, "CODEVER7", cvsinfoall);
  
 }; };
  
Line 2017  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.15  
changed lines
  Added in v.1.18

Karen Tian
Powered by
ViewCVS 0.9.4