Return to jrebinsmooth.c CVS log

Up to [Development] / JSOC / proj / globalhs / apps

no need to create a temporary output record to determine the name of the history series - use the template record for the output series

/*
this module (optionally) rebins its input data and then (optionally) convolves with a gaussian and subsamples
*/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <sys/time.h>
#include <sys/resource.h>
#include "jsoc_main.h"
#include "drms_dsdsapi.h"
#include "fresize.h"
#include "fstats.h"

#define ARRLENGTH(ARR)  (sizeof(ARR)/sizeof(ARR[0]))
#define QUAL_NODATA	(0x80000000)
#define kNOTSPECIFIED   "not specified"

char *module_name = "jrebinsmooth";
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 $";

ModuleArgs_t module_args[] = 
{
   {ARG_STRING,  "in", "", "input data records"},
   {ARG_STRING,  "out", "", "output data series"},
   {ARG_STRING,  "interout", kNOTSPECIFIED, "output data series for intermediate data (binned only)"},
   {ARG_STRING,  "segin", kNOTSPECIFIED, "name of input segment if not using segment 0"},
   {ARG_STRING,  "segout", kNOTSPECIFIED, "name of output segment if not using segment 0"},
   {ARG_STRING,  "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"},
   {ARG_STRING,  "srclink",  "SRCDATA", "name of link to source data"},
   {ARG_INT,     "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
   {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 */

   {ARG_INT,   "IBIN",        "0",    "flag for binning"},
   {ARG_INT,   "NBIN",        "-1",   "factor for binning"},
   {ARG_INT,   "BIN_XOFF",    "0",    "offset in pixels in x direction to apply before binning"},
   {ARG_INT,   "BIN_YOFF",    "0",    "offset in pixels in y direction to apply before binning"},
   {ARG_FLOAT, "BIN_FILL",    "0.0",  "value to use outside valid area"},

   {ARG_INT,   "IGAUSS",      "0",    "flag for convolving with gaussian"},
   {ARG_FLOAT, "GAUSS_SIG",   "-1.0", "width of gaussian"},
   {ARG_INT,   "GAUSS_WID",   "-1",   "half width of kernel"},
   {ARG_INT,   "GAUSS_NSUB",  "-1",   "distance between sampled points in pixels"},
   {ARG_INT,   "GAUSS_XOFF",  "0",    "offset in pixels to add to x index of input image"},
   {ARG_INT,   "GAUSS_YOFF",  "0",    "offset in pixels to add to y index of input image"},
   {ARG_FLOAT, "GAUSS_FILL",  "0.0",  "value to use outside valid area"},

   {ARG_END}
};

#include "saveparm.c"
#include "timing.c"
#include "set_history.c"

int DoIt(void)
{

  int newstat = 0;
  int status = DRMS_SUCCESS;
  int quality;
  TIME trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
  char trecstr[100];
  struct fresize_struct fress;
  int xpixin, ypixin, xpix1, ypix1, xpix2, ypix2;
  float x0in, y0in, xscalein, yscalein, x01, y01, xscale1, yscale1, x02, y02, xscale2, yscale2;
  float *inptr, *outptr;

  DRMS_RecordSet_t *inrecset = NULL;
  DRMS_Segment_t *segin = NULL;
  DRMS_Segment_t *segout = NULL;
  DRMS_Array_t *inarr = NULL;
  DRMS_Array_t *outarr = NULL;
  DRMS_Record_t *inrec = NULL;
  DRMS_Record_t *outrec = NULL;

  int length[2];
  DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
  DRMS_RecLifetime_t lifetime;

  int nrecs, irec, error;
  long long histrecnum=-1;

  int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
  int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);

  char *inrecquery = (char *)cmdparams_save_str(&cmdparams, "in", &newstat);
  char *outseries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat);
  char *segnamein = (char *)cmdparams_save_str(&cmdparams, "segin", &newstat);
  char *segnameout = (char *)cmdparams_save_str(&cmdparams, "segout", &newstat);
  int seginflag = strcmp(kNOTSPECIFIED, segnamein);
  int segoutflag = strcmp(kNOTSPECIFIED, segnameout);
  char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);
  char *srclinkname = (char *)cmdparams_save_str(&cmdparams, "srclink", &newstat);
  char *interout = (char *)cmdparams_save_str(&cmdparams, "interout", &newstat);
  int interflag = strcmp(kNOTSPECIFIED, interout);
  float *interptr=NULL;

  int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
  int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
  if (permflag)
    lifetime = DRMS_PERMANENT;
  else
    lifetime = DRMS_TRANSIENT;

  int ibin = cmdparams_save_int(&cmdparams, "IBIN", &newstat);
  int nbin = cmdparams_save_int(&cmdparams, "NBIN", &newstat);
  int bxoff = cmdparams_save_int(&cmdparams, "BIN_XOFF", &newstat);
  int byoff = cmdparams_save_int(&cmdparams, "BIN_YOFF", &newstat);
  float bfill = cmdparams_save_float(&cmdparams, "BIN_FILL", &newstat);

  int igauss = cmdparams_save_int(&cmdparams, "IGAUSS", &newstat);
  float sigma = cmdparams_save_float(&cmdparams, "GAUSS_SIG", &newstat);
  int hwidth = cmdparams_save_int(&cmdparams, "GAUSS_WID", &newstat);
  int nsub = cmdparams_save_int(&cmdparams, "GAUSS_NSUB", &newstat);
  int gxoff = cmdparams_save_int(&cmdparams, "GAUSS_XOFF", &newstat);
  int gyoff = cmdparams_save_int(&cmdparams, "GAUSS_YOFF", &newstat);
  float gfill = cmdparams_save_float(&cmdparams, "GAUSS_FILL", &newstat);

  if (newstat) 
  {
    fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
    cpsave_decode_error(newstat);
    return 1;
  }  
  else if (savestrlen != strlen(savestr)) 
  {
    fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr));
    return 1;
  }

  if (ibin == 0 && igauss == 0)
  {
    fprintf(stderr,"ERROR: must set at least one of IBIN and IGAUSS\n");
    return 1;
  }

  if (interflag != 0 && (ibin == 0 || igauss == 0))
  {
    fprintf(stderr,"ERROR: cannot specify interout unless both IBIN and IGAUSS are set\n");
    return 1;
  }

  if (ibin != 0 && nbin < 0)
  {
    fprintf(stderr,"ERROR: must specify NBIN > 0 when IBIN is set\n");
    return 1;
  }

  if (igauss != 0 && (sigma < 0.0 || hwidth < 0 || nsub < 0))
  {
    fprintf(stderr,"ERROR: must specify nonnegative value for each of GAUSS_SIG, GAUSS_WID, and GAUSS_NSUB when IGAUSS is set\n");
    return 1;
  }

  if (ibin == 0)
    nbin=1;
  if (igauss == 0)
    nsub=1;


  DRMS_Record_t *templateOut = drms_template_record(drms_env, outseries, &status);

  if (status != DRMS_SUCCESS || templateOut == NULL) 
  {
    fprintf(stderr,"ERROR: couldn't open output dataseries %s template record, status = %d\n", outseries, status);
    return 1;
  }

  DRMS_Link_t *histlink = hcon_lookup_lower(&templateOut->links, histlinkname);
  DRMS_Link_t *srclink = hcon_lookup_lower(&templateOut->links, srclinkname);

//  char *cvsinfo = strdup("$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jrebinsmooth.c,v 1.7 2018/02/02 19:10:00 arta Exp $");
  if (histlink != NULL) 
  {
    histrecnum=set_history(histlink);
    if (histrecnum < 0)
    {
      return 1;
    }
  }
  else
  {
    fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
  }

  inrecset = drms_open_records(drms_env, inrecquery, &status);
  nrecs = inrecset->n;
  if (status != DRMS_SUCCESS || inrecset == NULL)
  {
    fprintf(stderr, "ERROR: problem opening input recordset: status = %d\n", status);
    return 1;
  }

  if (nrecs == 0)
  {
    printf("ERROR: input recordset contains no records\n");
    drms_close_records(inrecset, DRMS_FREE_RECORD);
    return 1;
  }

  if (verbflag) 
    printf("input recordset opened, nrecs = %d\n",nrecs);

  char *inchecklist[] = {"T_REC", "QUALITY"};
  inrec=inrecset->records[0];
  int itest;
  for (itest=0; itest < ARRLENGTH(inchecklist); itest++)
  {
    DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]);
    if (!inkeytest)
    {
      fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
      drms_close_records(inrecset, DRMS_FREE_RECORD);
      return 1;
    }
  }


  init_fresize_gaussian(&fress, sigma, hwidth, nsub);
//  init_fresize_gaussian_fft(&fress, sigma, hwidth, nsub, nxin, nyin);

  DRMS_RecordSet_t *outrecset = drms_create_records(drms_env, nrecs, outseries, lifetime, &status);
  if (status != DRMS_SUCCESS || outrecset == NULL)
  {
    fprintf(stderr,"ERROR: unable to create record in output dataseries, status = %d\n", status);
    return 1;
  }

  DRMS_RecordSet_t *interoutrecset;
  if (interflag)
  {
    interoutrecset = drms_create_records(drms_env, nrecs, interout, lifetime, &status);
    if (status != DRMS_SUCCESS || interoutrecset == NULL)
    {
      fprintf(stderr,"ERROR: unable to create record in intermediate output dataseries, status = %d\n", status);
      return 1;
    }
  }

  error=0;
  int nodata=0;
  for (irec=0; irec<nrecs; irec++)
  {

    inrec = inrecset->records[irec];
    outrec = outrecset->records[irec];

    if (histlink)
      drms_setlink_static(outrec, histlinkname,  histrecnum);
    if (srclink)
      drms_setlink_static(outrec, srclinkname,  inrec->recnum);

    trec = drms_getkey_time(inrec, "T_REC", &status);
    if (status != DRMS_SUCCESS)
    {
      fprintf(stderr, "SKIP: error getting primekey T_REC: irec = %d, status = %d, recnum = %lld\n", irec, status, inrec->recnum);
      continue;
    }
    sprint_time(trecstr, trec, "TAI", 0); 

    quality = drms_getkey_int(inrec, "QUALITY", &status);
    if (status != DRMS_SUCCESS)
    {
      fprintf(stderr, "SKIP: error getting keyword QUALITY: irec = %d, status = %d, T_REC = %s, recnum = %lld\n", irec, status, trecstr, inrec->recnum);
      continue;
    }

    if (quality & QUAL_NODATA)
    {
      nodata=1;
      if (verbflag > 1)
        fprintf(stderr,"SKIP: record rejected based on quality = %d: irec = %d, T_REC = %s, recnum = %lld\n", quality, irec, trecstr, inrec->recnum);
//      continue;
    }
    else
    {
      nodata=0;

      if (seginflag)
        segin = drms_segment_lookup(inrec, segnamein);
      else
        segin = drms_segment_lookupnum(inrec, 0);

      if (segin)
        inarr = drms_segment_read(segin, usetype, &status);

      if (segin  == NULL || inarr == NULL || status != DRMS_SUCCESS)
      {
        fprintf(stderr, "ERROR: problem with input segment or array: irec = %d, status = %d, T_REC = %s, recnum = %lld \n", irec, status, trecstr, inrec->recnum);
        if (inarr)
          drms_free_array(inarr);
        error=1;
        break;
      }

      xpixin=inarr->axis[0];
      ypixin=inarr->axis[1];
      xpix1=xpixin/nbin;
      ypix1=ypixin/nbin;
      xpix2=xpix1/nsub;
      ypix2=ypix1/nsub;

      x0in = drms_getkey_float(inrec, "CRPIX1", &status) - 1;
      y0in = drms_getkey_float(inrec, "CRPIX2", &status) - 1;
      xscalein = drms_getkey_float(inrec, "CDELT1", &status);
      yscalein = drms_getkey_float(inrec, "CDELT2", &status);

      x01 = (x0in - bxoff + 0.5)/nbin - 0.5;
      y01 = (y0in - byoff + 0.5)/nbin - 0.5;
      xscale1=xscalein*nbin;
      yscale1=yscalein*nbin;

      if (ibin)
      {
        interptr=(float *)malloc(xpix1*ypix1*sizeof(float));
        fbin(
           (float *)inarr->data,
           interptr,
           xpixin,
           ypixin,
           xpixin,
           xpix1,
           ypix1,
           xpix1,
           nbin,
           bxoff,
           byoff,
           bfill);
      }
    }

    if (interflag && nodata)
    {
      DRMS_Record_t *outrec = interoutrecset->records[irec];
      DRMS_Link_t *histlink = hcon_lookup_lower(&outrec->links, histlinkname);
      DRMS_Link_t *srclink = hcon_lookup_lower(&outrec->links, srclinkname);
      if (histlink)
        drms_setlink_static(outrec, histlinkname,  histrecnum);
      if (srclink)
        drms_setlink_static(outrec, srclinkname,  inrec->recnum);
      drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
      tnow = (double)time(NULL);
      tnow += UNIX_epoch;
      drms_setkey_time(outrec, "DATE", tnow);

    }
    else if (interflag)
    {
//      DRMS_Record_t *outrec = drms_create_record(drms_env, interout, lifetime, &status);
      DRMS_Record_t *outrec = interoutrecset->records[irec];
      DRMS_Link_t *histlink = hcon_lookup_lower(&outrec->links, histlinkname);
      DRMS_Link_t *srclink = hcon_lookup_lower(&outrec->links, srclinkname);
      if (histlink)
        drms_setlink_static(outrec, histlinkname,  histrecnum);
      if (srclink)
        drms_setlink_static(outrec, srclinkname,  inrec->recnum);
      if (segoutflag)
        segout = drms_segment_lookup(outrec, segnameout);
      else
        segout = drms_segment_lookupnum(outrec, 0);
      length[0]=xpix1;
      length[1]=ypix1;
      outarr = drms_array_create(usetype, 2, length, interptr, &status);
      outarr->bzero=segout->bzero;
      outarr->bscale=segout->bscale;
      status=drms_segment_write(segout, outarr, 0);
      drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
      drms_setkey_int(outrec, "TOTVALS", xpix1*ypix1);
      set_statistics(segout, outarr, 1);
      drms_setkey_time(outrec, "T_REC", trec);
      drms_setkey_int(outrec, "QUALITY", quality);
      drms_setkey_float(outrec, "CRPIX1", x01+1);
      drms_setkey_float(outrec, "CRPIX2", y01+1);
      drms_setkey_float(outrec, "CDELT1", xscale1);
      drms_setkey_float(outrec, "CDELT2", yscale1);
      tnow = (double)time(NULL);
      tnow += UNIX_epoch;
      drms_setkey_time(outrec, "DATE", tnow);
//      drms_close_record(outrec, DRMS_INSERT_RECORD);
    }

    if (nodata)
    {
      drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
      tnow = (double)time(NULL);
      tnow += UNIX_epoch;
      drms_setkey_time(outrec, "DATE", tnow);
      continue;
    }

    length[0]=xpix2;
    length[1]=ypix2;

    if (segoutflag)
      segout = drms_segment_lookup(outrec, segnameout);
    else
      segout = drms_segment_lookupnum(outrec, 0);

    if (!igauss)
      outarr = drms_array_create(usetype, 2, length, interptr, &status);
    else
      outarr = drms_array_create(usetype, 2, length, NULL, &status);

    if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
    {
      fprintf(stderr, "ERROR: problem with output segment or array: irec = %d, status = %d\n", irec, status);
      drms_free_array(inarr);
      if (outarr)
        drms_free_array(outarr);
      error=1;
      break;
    }

    if (!ibin)
      inptr=(float *)inarr->data;
    else
      inptr=interptr;

    if (igauss)
    {
      fresize(
              &fress, // Must have been initialized by init_fresize_XXX
              inptr,
              (float *)(outarr->data),
              xpix1, // Size of input image
              ypix1, // Size of input image
              xpix1, // Leading dimension of input image. nleadin>=nxin
              xpix2, // Size of xin, yin and image_out
              ypix2, // Size of xin, yin and image_out
              xpix2, // Leading dimension. nlead>=nx
              gxoff, // Offset in x direction
              gyoff, // Offset in y direction
              gfill // Value to use if outside area
             );
    }

    x02 = (x01 - gxoff)/nsub;
    y02 = (y01 - gyoff)/nsub;
    xscale2=xscale1*nsub;
    yscale2=yscale1*nsub;

    outarr->bzero=segout->bzero;
    outarr->bscale=segout->bscale;
    status=drms_segment_write(segout, outarr, 0);
    drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
    drms_setkey_int(outrec, "TOTVALS", xpix2*ypix2);
    set_statistics(segout, outarr, 1);
    drms_setkey_time(outrec, "T_REC", trec);
    drms_setkey_int(outrec, "QUALITY", quality);
    drms_setkey_float(outrec, "CRPIX1", x02+1);
    drms_setkey_float(outrec, "CRPIX2", y02+1);
    drms_setkey_float(outrec, "CDELT1", xscale2);
    drms_setkey_float(outrec, "CDELT2", yscale2);
    tnow = (double)time(NULL);
    tnow += UNIX_epoch;
    drms_setkey_time(outrec, "DATE", tnow);

//    drms_close_record(outrec, DRMS_INSERT_RECORD);
    drms_free_array(inarr);
    drms_free_array(outarr);
    if (interptr != NULL)
      free(interptr);

  }

  drms_close_records(inrecset, DRMS_FREE_RECORD);
  if (error)
  {
    drms_close_records(outrecset, DRMS_FREE_RECORD);
    if (interflag)
      drms_close_records(interoutrecset, DRMS_FREE_RECORD);
  }
  else
  {
    drms_close_records(outrecset, DRMS_INSERT_RECORD);
    if (interflag)
      drms_close_records(interoutrecset, DRMS_INSERT_RECORD);
  }

  free_fresize(&fress);

  if (verbflag) 
  {
    if (!error)
      printf("module %s successful completion\n", cmdparams.argv[0]);
  }

  return 0;
}