JSOC_Documentation: proj/globalhs/apps/jrebinsmooth.c Source File

00001 /*
00002 this module (optionally) rebins its input data and then (optionally) convolves with a gaussian and subsamples
00003 */
00004 
00005 #include <stdio.h>
00006 #include <stdlib.h>
00007 #include <time.h>
00008 #include <math.h>
00009 #include <sys/time.h>
00010 #include <sys/resource.h>
00011 #include "jsoc_main.h"
00012 #include "drms_dsdsapi.h"
00013 #include "fresize.h"
00014 #include "fstats.h"
00015 
00016 #define ARRLENGTH(ARR)  (sizeof(ARR)/sizeof(ARR[0]))
00017 #define QUAL_NODATA (0x80000000)
00018 #define kNOTSPECIFIED   "not specified"
00019 
00020 char *module_name = "jrebinsmooth";
00021 char *cvsinfo_jrebinsmooth = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jrebinsmooth.c,v 1.7 2018/02/02 19:10:00 arta Exp $";
00022 
00023 ModuleArgs_t module_args[] = 
00024 {
00025    {ARG_STRING,  "in", "", "input data records"},
00026    {ARG_STRING,  "out", "", "output data series"},
00027    {ARG_STRING,  "interout", kNOTSPECIFIED, "output data series for intermediate data (binned only)"},
00028    {ARG_STRING,  "segin", kNOTSPECIFIED, "name of input segment if not using segment 0"},
00029    {ARG_STRING,  "segout", kNOTSPECIFIED, "name of output segment if not using segment 0"},
00030    {ARG_STRING,  "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"},
00031    {ARG_STRING,  "srclink",  "SRCDATA", "name of link to source data"},
00032    {ARG_INT,     "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
00033    {ARG_INT,     "VERB", "1", "option for level of verbosity: 0=only error and warning messages; >0=print messages outside of loop; >1=print messages inside loop; >2=debugging output", ""},  /* debug messages */
00034 
00035    {ARG_INT,   "IBIN",        "0",    "flag for binning"},
00036    {ARG_INT,   "NBIN",        "-1",   "factor for binning"},
00037    {ARG_INT,   "BIN_XOFF",    "0",    "offset in pixels in x direction to apply before binning"},
00038    {ARG_INT,   "BIN_YOFF",    "0",    "offset in pixels in y direction to apply before binning"},
00039    {ARG_FLOAT, "BIN_FILL",    "0.0",  "value to use outside valid area"},
00040 
00041    {ARG_INT,   "IGAUSS",      "0",    "flag for convolving with gaussian"},
00042    {ARG_FLOAT, "GAUSS_SIG",   "-1.0", "width of gaussian"},
00043    {ARG_INT,   "GAUSS_WID",   "-1",   "half width of kernel"},
00044    {ARG_INT,   "GAUSS_NSUB",  "-1",   "distance between sampled points in pixels"},
00045    {ARG_INT,   "GAUSS_XOFF",  "0",    "offset in pixels to add to x index of input image"},
00046    {ARG_INT,   "GAUSS_YOFF",  "0",    "offset in pixels to add to y index of input image"},
00047    {ARG_FLOAT, "GAUSS_FILL",  "0.0",  "value to use outside valid area"},
00048 
00049    {ARG_END}
00050 };
00051 
00052 #include "saveparm.c"
00053 #include "timing.c"
00054 #include "set_history.c"
00055 
00056 int DoIt(void)
00057 {
00058 
00059   int newstat = 0;
00060   int status = DRMS_SUCCESS;
00061   int quality;
00062   TIME trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
00063   char trecstr[100];
00064   struct fresize_struct fress;
00065   int xpixin, ypixin, xpix1, ypix1, xpix2, ypix2;
00066   float x0in, y0in, xscalein, yscalein, x01, y01, xscale1, yscale1, x02, y02, xscale2, yscale2;
00067   float *inptr, *outptr;
00068 
00069   DRMS_RecordSet_t *inrecset = NULL;
00070   DRMS_Segment_t *segin = NULL;
00071   DRMS_Segment_t *segout = NULL;
00072   DRMS_Array_t *inarr = NULL;
00073   DRMS_Array_t *outarr = NULL;
00074   DRMS_Record_t *inrec = NULL;
00075   DRMS_Record_t *outrec = NULL;
00076 
00077   int length[2];
00078   DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
00079   DRMS_RecLifetime_t lifetime;
00080 
00081   int nrecs, irec, error;
00082   long long histrecnum=-1;
00083 
00084   int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
00085   int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);
00086 
00087   char *inrecquery = (char *)cmdparams_save_str(&cmdparams, "in", &newstat);
00088   char *outseries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat);
00089   char *segnamein = (char *)cmdparams_save_str(&cmdparams, "segin", &newstat);
00090   char *segnameout = (char *)cmdparams_save_str(&cmdparams, "segout", &newstat);
00091   int seginflag = strcmp(kNOTSPECIFIED, segnamein);
00092   int segoutflag = strcmp(kNOTSPECIFIED, segnameout);
00093   char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);
00094   char *srclinkname = (char *)cmdparams_save_str(&cmdparams, "srclink", &newstat);
00095   char *interout = (char *)cmdparams_save_str(&cmdparams, "interout", &newstat);
00096   int interflag = strcmp(kNOTSPECIFIED, interout);
00097   float *interptr=NULL;
00098 
00099   int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
00100   int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
00101   if (permflag)
00102     lifetime = DRMS_PERMANENT;
00103   else
00104     lifetime = DRMS_TRANSIENT;
00105 
00106   int ibin = cmdparams_save_int(&cmdparams, "IBIN", &newstat);
00107   int nbin = cmdparams_save_int(&cmdparams, "NBIN", &newstat);
00108   int bxoff = cmdparams_save_int(&cmdparams, "BIN_XOFF", &newstat);
00109   int byoff = cmdparams_save_int(&cmdparams, "BIN_YOFF", &newstat);
00110   float bfill = cmdparams_save_float(&cmdparams, "BIN_FILL", &newstat);
00111 
00112   int igauss = cmdparams_save_int(&cmdparams, "IGAUSS", &newstat);
00113   float sigma = cmdparams_save_float(&cmdparams, "GAUSS_SIG", &newstat);
00114   int hwidth = cmdparams_save_int(&cmdparams, "GAUSS_WID", &newstat);
00115   int nsub = cmdparams_save_int(&cmdparams, "GAUSS_NSUB", &newstat);
00116   int gxoff = cmdparams_save_int(&cmdparams, "GAUSS_XOFF", &newstat);
00117   int gyoff = cmdparams_save_int(&cmdparams, "GAUSS_YOFF", &newstat);
00118   float gfill = cmdparams_save_float(&cmdparams, "GAUSS_FILL", &newstat);
00119 
00120   if (newstat) 
00121   {
00122     fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
00123     cpsave_decode_error(newstat);
00124     return 1;
00125   }  
00126   else if (savestrlen != strlen(savestr)) 
00127   {
00128     fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr));
00129     return 1;
00130   }
00131 
00132   if (ibin == 0 && igauss == 0)
00133   {
00134     fprintf(stderr,"ERROR: must set at least one of IBIN and IGAUSS\n");
00135     return 1;
00136   }
00137 
00138   if (interflag != 0 && (ibin == 0 || igauss == 0))
00139   {
00140     fprintf(stderr,"ERROR: cannot specify interout unless both IBIN and IGAUSS are set\n");
00141     return 1;
00142   }
00143 
00144   if (ibin != 0 && nbin < 0)
00145   {
00146     fprintf(stderr,"ERROR: must specify NBIN > 0 when IBIN is set\n");
00147     return 1;
00148   }
00149 
00150   if (igauss != 0 && (sigma < 0.0 || hwidth < 0 || nsub < 0))
00151   {
00152     fprintf(stderr,"ERROR: must specify nonnegative value for each of GAUSS_SIG, GAUSS_WID, and GAUSS_NSUB when IGAUSS is set\n");
00153     return 1;
00154   }
00155 
00156   if (ibin == 0)
00157     nbin=1;
00158   if (igauss == 0)
00159     nsub=1;
00160 
00161 
00162   DRMS_Record_t *templateOut = drms_template_record(drms_env, outseries, &status);
00163 
00164   if (status != DRMS_SUCCESS || templateOut == NULL) 
00165   {
00166     fprintf(stderr,"ERROR: couldn't open output dataseries %s template record, status = %d\n", outseries, status);
00167     return 1;
00168   }
00169 
00170   DRMS_Link_t *histlink = hcon_lookup_lower(&templateOut->links, histlinkname);
00171   DRMS_Link_t *srclink = hcon_lookup_lower(&templateOut->links, srclinkname);
00172 
00173 //  char *cvsinfo = strdup("$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jrebinsmooth.c,v 1.7 2018/02/02 19:10:00 arta Exp $");
00174   if (histlink != NULL) 
00175   {
00176     histrecnum=set_history(histlink);
00177     if (histrecnum < 0)
00178     {
00179       return 1;
00180     }
00181   }
00182   else
00183   {
00184     fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
00185   }
00186 
00187   inrecset = drms_open_records(drms_env, inrecquery, &status);
00188   nrecs = inrecset->n;
00189   if (status != DRMS_SUCCESS || inrecset == NULL)
00190   {
00191     fprintf(stderr, "ERROR: problem opening input recordset: status = %d\n", status);
00192     return 1;
00193   }
00194 
00195   if (nrecs == 0)
00196   {
00197     printf("ERROR: input recordset contains no records\n");
00198     drms_close_records(inrecset, DRMS_FREE_RECORD);
00199     return 1;
00200   }
00201 
00202   if (verbflag) 
00203     printf("input recordset opened, nrecs = %d\n",nrecs);
00204 
00205   char *inchecklist[] = {"T_REC", "QUALITY"};
00206   inrec=inrecset->records[0];
00207   int itest;
00208   for (itest=0; itest < ARRLENGTH(inchecklist); itest++)
00209   {
00210     DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]);
00211     if (!inkeytest)
00212     {
00213       fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
00214       drms_close_records(inrecset, DRMS_FREE_RECORD);
00215       return 1;
00216     }
00217   }
00218 
00219 
00220   init_fresize_gaussian(&fress, sigma, hwidth, nsub);
00221 //  init_fresize_gaussian_fft(&fress, sigma, hwidth, nsub, nxin, nyin);
00222 
00223   DRMS_RecordSet_t *outrecset = drms_create_records(drms_env, nrecs, outseries, lifetime, &status);
00224   if (status != DRMS_SUCCESS || outrecset == NULL)
00225   {
00226     fprintf(stderr,"ERROR: unable to create record in output dataseries, status = %d\n", status);
00227     return 1;
00228   }
00229 
00230   DRMS_RecordSet_t *interoutrecset;
00231   if (interflag)
00232   {
00233     interoutrecset = drms_create_records(drms_env, nrecs, interout, lifetime, &status);
00234     if (status != DRMS_SUCCESS || interoutrecset == NULL)
00235     {
00236       fprintf(stderr,"ERROR: unable to create record in intermediate output dataseries, status = %d\n", status);
00237       return 1;
00238     }
00239   }
00240 
00241   error=0;
00242   int nodata=0;
00243   for (irec=0; irec<nrecs; irec++)
00244   {
00245 
00246     inrec = inrecset->records[irec];
00247     outrec = outrecset->records[irec];
00248 
00249     if (histlink)
00250       drms_setlink_static(outrec, histlinkname,  histrecnum);
00251     if (srclink)
00252       drms_setlink_static(outrec, srclinkname,  inrec->recnum);
00253 
00254     trec = drms_getkey_time(inrec, "T_REC", &status);
00255     if (status != DRMS_SUCCESS)
00256     {
00257       fprintf(stderr, "SKIP: error getting primekey T_REC: irec = %d, status = %d, recnum = %lld\n", irec, status, inrec->recnum);
00258       continue;
00259     }
00260     sprint_time(trecstr, trec, "TAI", 0); 
00261 
00262     quality = drms_getkey_int(inrec, "QUALITY", &status);
00263     if (status != DRMS_SUCCESS)
00264     {
00265       fprintf(stderr, "SKIP: error getting keyword QUALITY: irec = %d, status = %d, T_REC = %s, recnum = %lld\n", irec, status, trecstr, inrec->recnum);
00266       continue;
00267     }
00268 
00269     if (quality & QUAL_NODATA)
00270     {
00271       nodata=1;
00272       if (verbflag > 1)
00273         fprintf(stderr,"SKIP: record rejected based on quality = %d: irec = %d, T_REC = %s, recnum = %lld\n", quality, irec, trecstr, inrec->recnum);
00274 //      continue;
00275     }
00276     else
00277     {
00278       nodata=0;
00279 
00280       if (seginflag)
00281         segin = drms_segment_lookup(inrec, segnamein);
00282       else
00283         segin = drms_segment_lookupnum(inrec, 0);
00284 
00285       if (segin)
00286         inarr = drms_segment_read(segin, usetype, &status);
00287 
00288       if (segin  == NULL || inarr == NULL || status != DRMS_SUCCESS)
00289       {
00290         fprintf(stderr, "ERROR: problem with input segment or array: irec = %d, status = %d, T_REC = %s, recnum = %lld \n", irec, status, trecstr, inrec->recnum);
00291         if (inarr)
00292           drms_free_array(inarr);
00293         error=1;
00294         break;
00295       }
00296 
00297       xpixin=inarr->axis[0];
00298       ypixin=inarr->axis[1];
00299       xpix1=xpixin/nbin;
00300       ypix1=ypixin/nbin;
00301       xpix2=xpix1/nsub;
00302       ypix2=ypix1/nsub;
00303 
00304       x0in = drms_getkey_float(inrec, "CRPIX1", &status) - 1;
00305       y0in = drms_getkey_float(inrec, "CRPIX2", &status) - 1;
00306       xscalein = drms_getkey_float(inrec, "CDELT1", &status);
00307       yscalein = drms_getkey_float(inrec, "CDELT2", &status);
00308 
00309       x01 = (x0in - bxoff + 0.5)/nbin - 0.5;
00310       y01 = (y0in - byoff + 0.5)/nbin - 0.5;
00311       xscale1=xscalein*nbin;
00312       yscale1=yscalein*nbin;
00313 
00314       if (ibin)
00315       {
00316         interptr=(float *)malloc(xpix1*ypix1*sizeof(float));
00317         fbin(
00318            (float *)inarr->data,
00319            interptr,
00320            xpixin,
00321            ypixin,
00322            xpixin,
00323            xpix1,
00324            ypix1,
00325            xpix1,
00326            nbin,
00327            bxoff,
00328            byoff,
00329            bfill);
00330       }
00331     }
00332 
00333     if (interflag && nodata)
00334     {
00335       DRMS_Record_t *outrec = interoutrecset->records[irec];
00336       DRMS_Link_t *histlink = hcon_lookup_lower(&outrec->links, histlinkname);
00337       DRMS_Link_t *srclink = hcon_lookup_lower(&outrec->links, srclinkname);
00338       if (histlink)
00339         drms_setlink_static(outrec, histlinkname,  histrecnum);
00340       if (srclink)
00341         drms_setlink_static(outrec, srclinkname,  inrec->recnum);
00342       drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
00343       tnow = (double)time(NULL);
00344       tnow += UNIX_epoch;
00345       drms_setkey_time(outrec, "DATE", tnow);
00346 
00347     }
00348     else if (interflag)
00349     {
00350 //      DRMS_Record_t *outrec = drms_create_record(drms_env, interout, lifetime, &status);
00351       DRMS_Record_t *outrec = interoutrecset->records[irec];
00352       DRMS_Link_t *histlink = hcon_lookup_lower(&outrec->links, histlinkname);
00353       DRMS_Link_t *srclink = hcon_lookup_lower(&outrec->links, srclinkname);
00354       if (histlink)
00355         drms_setlink_static(outrec, histlinkname,  histrecnum);
00356       if (srclink)
00357         drms_setlink_static(outrec, srclinkname,  inrec->recnum);
00358       if (segoutflag)
00359         segout = drms_segment_lookup(outrec, segnameout);
00360       else
00361         segout = drms_segment_lookupnum(outrec, 0);
00362       length[0]=xpix1;
00363       length[1]=ypix1;
00364       outarr = drms_array_create(usetype, 2, length, interptr, &status);
00365       outarr->bzero=segout->bzero;
00366       outarr->bscale=segout->bscale;
00367       status=drms_segment_write(segout, outarr, 0);
00368       drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
00369       drms_setkey_int(outrec, "TOTVALS", xpix1*ypix1);
00370       set_statistics(segout, outarr, 1);
00371       drms_setkey_time(outrec, "T_REC", trec);
00372       drms_setkey_int(outrec, "QUALITY", quality);
00373       drms_setkey_float(outrec, "CRPIX1", x01+1);
00374       drms_setkey_float(outrec, "CRPIX2", y01+1);
00375       drms_setkey_float(outrec, "CDELT1", xscale1);
00376       drms_setkey_float(outrec, "CDELT2", yscale1);
00377       tnow = (double)time(NULL);
00378       tnow += UNIX_epoch;
00379       drms_setkey_time(outrec, "DATE", tnow);
00380 //      drms_close_record(outrec, DRMS_INSERT_RECORD);
00381     }
00382 
00383     if (nodata)
00384     {
00385       drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
00386       tnow = (double)time(NULL);
00387       tnow += UNIX_epoch;
00388       drms_setkey_time(outrec, "DATE", tnow);
00389       continue;
00390     }
00391 
00392     length[0]=xpix2;
00393     length[1]=ypix2;
00394 
00395     if (segoutflag)
00396       segout = drms_segment_lookup(outrec, segnameout);
00397     else
00398       segout = drms_segment_lookupnum(outrec, 0);
00399 
00400     if (!igauss)
00401       outarr = drms_array_create(usetype, 2, length, interptr, &status);
00402     else
00403       outarr = drms_array_create(usetype, 2, length, NULL, &status);
00404 
00405     if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
00406     {
00407       fprintf(stderr, "ERROR: problem with output segment or array: irec = %d, status = %d\n", irec, status);
00408       drms_free_array(inarr);
00409       if (outarr)
00410         drms_free_array(outarr);
00411       error=1;
00412       break;
00413     }
00414 
00415     if (!ibin)
00416       inptr=(float *)inarr->data;
00417     else
00418       inptr=interptr;
00419 
00420     if (igauss)
00421     {
00422       fresize(
00423               &fress, // Must have been initialized by init_fresize_XXX
00424               inptr,
00425               (float *)(outarr->data),
00426               xpix1, // Size of input image
00427               ypix1, // Size of input image
00428               xpix1, // Leading dimension of input image. nleadin>=nxin
00429               xpix2, // Size of xin, yin and image_out
00430               ypix2, // Size of xin, yin and image_out
00431               xpix2, // Leading dimension. nlead>=nx
00432               gxoff, // Offset in x direction
00433               gyoff, // Offset in y direction
00434               gfill // Value to use if outside area
00435              );
00436     }
00437 
00438     x02 = (x01 - gxoff)/nsub;
00439     y02 = (y01 - gyoff)/nsub;
00440     xscale2=xscale1*nsub;
00441     yscale2=yscale1*nsub;
00442 
00443     outarr->bzero=segout->bzero;
00444     outarr->bscale=segout->bscale;
00445     status=drms_segment_write(segout, outarr, 0);
00446     drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
00447     drms_setkey_int(outrec, "TOTVALS", xpix2*ypix2);
00448     set_statistics(segout, outarr, 1);
00449     drms_setkey_time(outrec, "T_REC", trec);
00450     drms_setkey_int(outrec, "QUALITY", quality);
00451     drms_setkey_float(outrec, "CRPIX1", x02+1);
00452     drms_setkey_float(outrec, "CRPIX2", y02+1);
00453     drms_setkey_float(outrec, "CDELT1", xscale2);
00454     drms_setkey_float(outrec, "CDELT2", yscale2);
00455     tnow = (double)time(NULL);
00456     tnow += UNIX_epoch;
00457     drms_setkey_time(outrec, "DATE", tnow);
00458 
00459 //    drms_close_record(outrec, DRMS_INSERT_RECORD);
00460     drms_free_array(inarr);
00461     drms_free_array(outarr);
00462     if (interptr != NULL)
00463       free(interptr);
00464 
00465   }
00466 
00467   drms_close_records(inrecset, DRMS_FREE_RECORD);
00468   if (error)
00469   {
00470     drms_close_records(outrecset, DRMS_FREE_RECORD);
00471     if (interflag)
00472       drms_close_records(interoutrecset, DRMS_FREE_RECORD);
00473   }
00474   else
00475   {
00476     drms_close_records(outrecset, DRMS_INSERT_RECORD);
00477     if (interflag)
00478       drms_close_records(interoutrecset, DRMS_INSERT_RECORD);
00479   }
00480 
00481   free_fresize(&fress);
00482 
00483   if (verbflag) 
00484   {
00485     if (!error)
00486       printf("module %s successful completion\n", cmdparams.argv[0]);
00487   }
00488 
00489   return 0;
00490 }