JSOC_Documentation: proj/mag/remapmags/apps/resizemappingmag.c Source File

00001 /*
00002  * Module name:     resizemappingmag.c
00003  *                   copy from Xudong's code testresize.c
00004  *
00005  * Example: resizemappingmag in='su_yang.fd_Mremap_large[2010.08.12/1h]' out='su_yang.resizetest' nbin=3
00006  *
00007  * Note: resize remapped mags to 1801x1440 using Jesper's rebin code. The average
00008  *       was done with a gaussian function.
00009  */
00010 
00011 
00012 #include <jsoc_main.h>
00013 #include <stdio.h>
00014 #include <stdlib.h>
00015 #include <math.h>
00016 #include <sys/time.h>
00017 #include <sys/resource.h>
00018   
00019 #include <mkl_blas.h>
00020 #include <mkl_service.h>
00021 #include <mkl_lapack.h>
00022 #include <mkl_vml_functions.h>
00023 #include <omp.h>
00024 
00025 #include "fresize.h"
00026 #include "fstats.h"
00027 
00028 #define DEG2RAD (M_PI / 180.)
00029 #define RAD2ARCSEC  (648000. / M_PI)
00030 
00031 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
00032 #define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status=%d\n", msg, status); return(status);}
00033 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
00034 
00035 
00036 /* ################## Timing from T. Larson ################## */
00037 
00038 double getwalltime(void)
00039 {
00040   struct timeval tv;
00041   gettimeofday(&tv, NULL);
00042   return tv.tv_sec * 1000.0 + tv.tv_usec/1000.0;
00043 }
00044 
00045 double getcputime(double *utime, double *stime)
00046 {
00047 
00048   struct rusage ru;
00049   getrusage(RUSAGE_SELF, &ru);
00050   *utime = ru.ru_utime.tv_sec * 1000.0 + ru.ru_utime.tv_usec / 1000.0;
00051   *stime = ru.ru_stime.tv_sec * 1000.0 + ru.ru_stime.tv_usec / 1000.0;
00052   return *utime + *stime;
00053 }
00054 
00055 
00056 /* ################## Wrapper for Jesper's code ################## */
00057 
00058 void frebin(float *image_in, float *image_out, int nx_in, int ny_in, int nx_out, int ny_out, int nbin)
00059 {
00060   struct fresize_struct fresizes;
00061 //  int nxout, nyout;
00062   int nlead_in = nx_in, nlead_out = nx_out;
00063   
00064 //  nxout = nx / nbin; nyout = ny / nbin;
00065   init_fresize_gaussian2(&fresizes, nbin, (3 * nbin)/2, (3 * nbin)/2, nbin);
00066 
00067 //  fresize(&fresizes, image_in, image_out, nx, ny, nlead, nxout, nyout, nxout, (nbin / 2) * 2, (nbin / 2) * 2, DRMS_MISSING_FLOAT);
00068 
00069   fresize(&fresizes, image_in, image_out, nx_in, ny_in, nlead_in, nx_out, ny_out, nlead_out, 4*nbin + 1, 4*nbin + 1, DRMS_MISSING_FLOAT);
00070   free_fresize(&fresizes);
00071 
00072 }
00073 
00074 /* ################## Main Module ################## */
00075 
00076 char *module_name = "resizemappingmag"; /* Module name */
00077 
00078 ModuleArgs_t module_args[] =
00079 {
00080     {ARG_STRING, "in", NULL, "Input magnetogram series."},
00081     {ARG_STRING, "out", NULL, "Output data series."},
00082     {ARG_INT, "nbin", "3", "nbin"},
00083     {ARG_INT, "VERB", "1", "Level of verbosity: 0=errors/warnings; 1=all messages"},
00084     {ARG_END}
00085 };
00086 
00087 
00088 int DoIt(void)
00089 {
00090     int status = DRMS_SUCCESS;
00091     char *inQuery, *outQuery;
00092     DRMS_RecordSet_t *inRS = NULL;
00093     DRMS_RecordSet_t *outRS = NULL;
00094     int irec, nrecs;
00095 //    DRMS_Record_t *outRec = NULL;
00096 //    DRMS_Segment_t *outSeg;
00097 //    DRMS_Array_t *outArray;
00098 //    float *inData, *outData;
00099     float fltemp;
00100     int intemp;
00101     int counting = 0;
00102 
00103     int i, j, itest;
00104     int nbin;
00105     int verbflag;
00106     int outDims[2];
00107     int n0, n1, xout, yout;
00108     TIME t_rec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */;
00109 
00110     double statMin, statMax, statMedn, statMean, statSig, statSkew, statKurt;
00111     int statNgood;
00112 
00113     // Time measuring
00114     double wt0, wt1, wt;
00115     double ut0, ut1, ut;
00116     double st0, st1, st;
00117     double ct0, ct1, ct;
00118     wt0 = getwalltime();
00119     ct0 = getcputime(&ut0, &st0);
00120 
00121     char historyofthemodule[2048]; // put history info into the data
00122     sprintf(historyofthemodule,"Smoothing edge-pixel bug corrected -- July 2013");
00123 
00124     /* Get parameters */
00125     inQuery = (char *)params_get_str(&cmdparams, "in");
00126     outQuery = (char *)params_get_str(&cmdparams, "out");
00127     verbflag = params_get_int(&cmdparams, "VERB");
00128     nbin = params_get_int(&cmdparams, "nbin");
00129 
00130     /* Open input */
00131     inRS = drms_open_records(drms_env, inQuery, &status);
00132     if (status || inRS->n == 0) DIE("No input data found");
00133     nrecs = inRS->n;
00134 
00135     /* Create output */
00136 //    outRS = drms_create_records(drms_env, nrecs, outQuery, DRMS_PERMANENT, &status);
00137 //    if (status) DIE("Output recordset not created");
00138 
00139     /* Do this for each record */
00140     for (irec = 0; irec < nrecs; irec++)
00141     {   
00142         DRMS_Record_t *inRec = NULL;
00143         DRMS_Array_t *inArray;
00144         DRMS_Segment_t *inSeg = NULL;
00145         DRMS_Record_t *outRec = NULL;
00146         DRMS_Segment_t *outSeg = NULL;
00147         DRMS_Array_t *outArray;
00148         float *inData = NULL;
00149         float *outData = NULL;
00150  
00151         if (verbflag) {
00152             wt1 = getwalltime();
00153             ct1 = getcputime(&ut1, &st1);
00154             printf("processing record %d...\n", irec);
00155         }
00156 
00157         /* Input record and data */
00158         inRec = inRS->records[irec];        
00159         t_rec = drms_getkey_time(inRec, "T_REC", &status);
00160 
00161         if (!(inSeg = drms_segment_lookupnum(inRec, 0))) {
00162           printf("cannot open file\n");
00163           continue;
00164         }
00165 
00166         inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
00167         if (status) continue;
00168         inData = inArray->data;
00169         n0 = inArray->axis[0]; n1 = inArray->axis[1];
00170 // create an array larger than the original array. The extra data at edge are 
00171 // mirrowed with the data at the original edge.
00172         int xNewdim = n0 + 8 * nbin, yNewdim = n1 + 8 * nbin;
00173         int xStart = 4 * nbin, yStart = 4 * nbin;
00174         long long iDataOld, yOffOld, iDataNew, yOffNew;
00175         int ii, jj;
00176         float *tmpData;
00177 
00178         xout = n0 / nbin; yout = n1 / nbin;
00179 //        xout = xNewdim; yout = yNewdim;
00180         if (verbflag) printf("n0=%d, n1=%d, xout=%d, yout=%d\n", n0, n1, xout, yout); //fflush(stdout);
00181         
00182         /* Output array */
00183         outDims[0] = xout; outDims[1] = yout;       // Odd pixels
00184         outData = (float *) malloc(xout * yout * sizeof(float));
00185     outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, outDims, outData, &status);
00186 
00187         tmpData = (float *) malloc(xNewdim * yNewdim * sizeof(float));
00188 
00189         for (jj = 0; jj < 4 * nbin; jj++) {
00190             yOffNew = jj * xNewdim;
00191             yOffOld = (4 * nbin - jj) * n0;
00192             for (ii = 0; ii < n0; ii++) {
00193                 iDataNew = yOffNew + xStart + ii;
00194                 iDataOld = yOffOld + ii;
00195                 tmpData[iDataNew] = inData[iDataOld];
00196             }
00197         }
00198       
00199         for (jj = 0; jj < 4 * nbin; jj++) {
00200             yOffNew = (jj + n1 + 4 * nbin) * xNewdim;
00201             yOffOld = (n1 - 1 - jj) * n0;
00202             for (ii = 0; ii < n0; ii++) {
00203                 iDataNew = yOffNew + xStart + ii;
00204                 iDataOld = yOffOld + ii;
00205                 tmpData[iDataNew] = inData[iDataOld];
00206             }
00207         }
00208 
00209         for (jj = 0; jj < n1; jj++) {
00210             yOffNew = (jj + 4 * nbin) * xNewdim;
00211             yOffOld = jj * n0;
00212             for (ii = 0; ii < n0; ii++) {
00213                 iDataNew = yOffNew + xStart + ii;
00214                 iDataOld = yOffOld + ii;
00215                 tmpData[iDataNew] = inData[iDataOld];
00216             }
00217         }
00218 
00219         for (jj = 0; jj < yNewdim; jj++) {
00220             yOffNew = jj * xNewdim;
00221             for (ii = 0; ii < xStart; ii++) {
00222                 iDataNew = yOffNew + (xStart - 1) - ii;
00223                 iDataOld = yOffNew + xStart + ii;
00224                 tmpData[iDataNew] = tmpData[iDataOld];
00225             }
00226         }
00227 
00228         for (jj = 0; jj < yNewdim; jj++) {
00229             yOffNew = jj * xNewdim;
00230             for (ii = 0; ii < xStart; ii++) {
00231                 iDataNew = yOffNew + (xNewdim - xStart) + ii;
00232                 iDataOld = yOffNew + (xNewdim - xStart) - 1 - ii;
00233                 tmpData[iDataNew] = tmpData[iDataOld];
00234             }
00235         }
00236  
00237         /************************************************************/
00238         /************************ Rebin here ************************/
00239         /************************************************************/
00240         frebin(tmpData, outData, xNewdim, yNewdim, xout, yout, nbin);
00241         counting++;
00242 printf("processing files number = %d\n", counting);
00243 
00244 // compute statistics
00245         if (fstats(outDims[0]*outDims[1], outData, &statMin, &statMax, &statMedn, &statMean, &statSig,
00246             &statSkew, &statKurt, &statNgood)) printf("\n Statistics computation failed\n");
00247 
00248         /* Output record */
00249 
00250         outRec = drms_create_record(drms_env, outQuery, DRMS_PERMANENT, &status);
00251         outSeg = drms_segment_lookupnum(outRec, 0);
00252         status = drms_segment_write(outSeg, outArray, 0);
00253         if (status) DIE("Problem writing file");
00254 
00255         // Prime key
00256         drms_copykey(outRec, inRec, "T_REC");
00257 
00258         // other keywords
00259         drms_setkey_string(outRec, "HISTORY", historyofthemodule);
00260         tnow = (double)time(NULL);
00261         tnow += UNIX_epoch;
00262         drms_setkey_time(outRec, "DATE", tnow);
00263         drms_copykey(outRec, inRec, "INSTRUME");
00264         drms_copykey(outRec, inRec, "CAMERA");
00265         drms_copykey(outRec, inRec, "COMMENT");
00266         drms_copykey(outRec, inRec, "BLD_VERS");
00267 
00268         drms_copykey(outRec, inRec, "QUALITY");
00269         drms_copykey(outRec, inRec, "CADENCE");
00270         drms_setkey_double(outRec, "DATAMIN", statMin);
00271         drms_setkey_double(outRec, "DATAMAX", statMax);
00272         drms_setkey_double(outRec, "DATAMEDN", statMedn);
00273         drms_setkey_double(outRec, "DATAMEAN", statMean);
00274         drms_setkey_double(outRec, "DATARMS", statSig);
00275         drms_setkey_double(outRec, "DATASKEW", statSkew);
00276         drms_setkey_double(outRec, "DATAKURT", statKurt);
00277         i = outDims[0]*outDims[1];
00278         drms_setkey_int(outRec, "TOTVALS", i);
00279         drms_setkey_int(outRec, "DATAVALS", statNgood);
00280         i = outDims[0]*outDims[1]-statNgood;
00281         drms_setkey_int(outRec, "MISSVALS", i);
00282 
00283         fltemp = drms_getkey_float(inRec, "CRPIX1", &status);
00284         fltemp = (outDims[0]+1.0)/2.0;
00285         drms_setkey_float(outRec, "CRPIX1", fltemp);
00286         fltemp = drms_getkey_float(inRec, "CRPIX2", &status);        
00287         fltemp = (outDims[1]+1.0)/2.0;
00288         drms_setkey_float(outRec, "CRPIX2", fltemp);
00289         drms_copykey(outRec, inRec, "CRVAL1");
00290         drms_copykey(outRec, inRec, "CRVAL2");
00291         fltemp = drms_getkey_float(inRec, "CDELT1", &status);
00292         drms_setkey_float(outRec, "CDELT1", fltemp*nbin);
00293         fltemp = drms_getkey_float(inRec, "CDELT2", &status);
00294         drms_setkey_float(outRec, "CDELT2", fltemp*nbin);
00295         drms_copykey(outRec, inRec, "CROTA2");
00296         drms_copykey(outRec, inRec, "CRDER1");
00297         drms_copykey(outRec, inRec, "CRDER2");
00298         drms_copykey(outRec, inRec, "CSYSER1");
00299         drms_copykey(outRec, inRec, "CSYSER2");
00300         drms_copykey(outRec, inRec, "WCSNAME");
00301       
00302         drms_copykey(outRec, inRec, "DSUN_OBS");
00303         drms_copykey(outRec, inRec, "CRLN_OBS");
00304         drms_copykey(outRec, inRec, "CRLT_OBS");
00305         drms_copykey(outRec, inRec, "HGLN_OBS");
00306         drms_copykey(outRec, inRec, "CAR_ROT");
00307         drms_copykey(outRec, inRec, "OBS_VR");
00308         drms_copykey(outRec, inRec, "OBS_VW");
00309         drms_copykey(outRec, inRec, "OBS_VN");
00310 
00311         drms_copykey(outRec, inRec, "T_OBS");
00312         drms_copykey(outRec, inRec, "DATASIGN");
00313         drms_copykey(outRec, inRec, "MAPLGMAX");
00314         drms_copykey(outRec, inRec, "MAPLGMIN");
00315         intemp = drms_getkey_int(inRec, "MAPMMAX", &status);
00316         drms_setkey_float(outRec, "MAPMMAX", intemp/nbin);
00317         intemp = drms_getkey_int(inRec, "SINBDIVS", &status);
00318         drms_setkey_float(outRec, "SINBDIVS", intemp/nbin);
00319 
00320         drms_copykey(outRec, inRec, "MAPBMAX");
00321         drms_copykey(outRec, inRec, "MAPRMAX");
00322         drms_copykey(outRec, inRec, "LGSHIFT");
00323         drms_copykey(outRec, inRec, "INTERPO");
00324         drms_copykey(outRec, inRec, "MCORLEV");
00325         drms_copykey(outRec, inRec, "MOFFSET");
00326         drms_copykey(outRec, inRec, "CARSTRCH");
00327         drms_copykey(outRec, inRec, "DIFROT_A");
00328         drms_copykey(outRec, inRec, "DIFROT_B");
00329         drms_copykey(outRec, inRec, "DIFROT_C");
00330         drms_copykey(outRec, inRec, "CALVER64");
00331 
00332         /* Time measure */
00333         if (verbflag) {
00334             wt = getwalltime();
00335             ct = getcputime(&ut, &st);
00336             printf("record %d done, %.2f ms wall time, %.2f ms cpu time\n", 
00337                      irec, wt - wt1, ct - ct1);
00338         }
00339 
00340         /* Clean up */
00341         free(tmpData);
00342         drms_free_array(inArray);
00343         drms_free_array(outArray);
00344         drms_close_record(outRec, DRMS_INSERT_RECORD);
00345     } // end of loop
00346      
00347     drms_close_records(inRS, DRMS_FREE_RECORD);
00348     if (verbflag) {
00349         wt = getwalltime();
00350         ct = getcputime(&ut, &st);
00351         printf("total time spent: %.2f ms wall time, %.2f ms cpu time\n", 
00352                 wt - wt0, ct - ct0);
00353     }
00354 
00355     return(DRMS_SUCCESS);
00356 
00357 }  // End of module