util/apps/rebin2.c

/**
   @defgroup rebin2 rebin2 reduce/increase resolution to that of MDI/HMI
   @ingroup su_util

   @brief Reduce (increase) the resolution of the input data to that of MDI (HMI) by factor of 2.

   @par Synopsis:
   @code
   rebin2  in=input data out=output data  factor=2*n mode=simple
   where n is a multiple/fraction of 2.
   @endcode

   This is a general purpose module that takes a series of input 
   data and modifies its spatial resolution by a factor (multiples 
   or fractions of 2) as required and gives out a set of output 
   data. The method for avaraging (interpolation) can be specified 
   through the input "mode". The current version handles a simple 
   boxcar average. If 'scale' < 0 then the input is reduced in size to 1/|scale|.

   Make sure you created the appropriate output series before running 
   the program. For example, su_bala.rebin2up.jsd

   @par Flags:
   @c none
   Currently it doesn't have any flags.

   @par GEN_FLAGS:
   Ubiquitous flags present in every module.
   @ref jsoc_main

   @param in  The input data series.
   @param out The output series.

   @par Exit_Status:
   Brief description of abnormal, non-zero, exit values.

   @par Example:
   Takes a series of 1024 X 1024 MDI full disk magnetogram and
   produces images with the resolution of HMI, 4096 X 4096.

   @code
   rebin2 in='mdi.fd_M_96m_lev18[2003.10.20/1d]' out='su_bala.rebin2up' scale=4 mode='simple'
   @endcode

   @par Example:
   Reduces the resolution of HMI images 4096 X 4096 to that of MDI,
   1024 X 1024. Here the input is the HMI images and the output
   is the lower resolution HMI images.
   @code
   rebin2 in='su_bala.rebin2up[2003.10.20/1d]' out='su_bala.rebin2down' scale=0.25 mode='simple'
   @endcode

   @bug
   None known so far.

   @par Code:
   The doxygen code that makes this page is here:
   @verbinclude  ./doxygen_moduletemplate.txt
*/

#include "jsoc.h"
#include "jsoc_main.h"
#include "fstats.h"

char *module_name = "rebin2";

#define DIE(msg) {fflush(stdout);fprintf(stderr,"%s, status=%d\n",msg,status); return(status);}

ModuleArgs_t module_args[] =
{
     {ARG_STRING, "in", "NOT SPECIFIED",  "Input data series."},
     {ARG_STRING, "out", "NOT SPECIFIED",  "Output data series."},
     {ARG_FLOAT, "scale", "NOTSPECIFIED", "Scale factor. negative for shrink factor"},
     {ARG_STRING, "mode", "simple", "conversion type."},
     {ARG_END}
};

int DoIt(void)
  {
  int status = DRMS_SUCCESS;
  DRMS_RecordSet_t *inRS, *outRS;
  int irec, nrecs;
  const char *inQuery = params_get_str(&cmdparams, "in");
  const char *outSeries = params_get_str(&cmdparams, "out");
  const char *mode = params_get_str(&cmdparams, "mode");
  int scale = params_get_int(&cmdparams, "scale");
  int iscale, factor;
  float fscale;
  if (scale < 0.0)
    { // shrinking
    iscale = -scale;
    fscale = 1.0/iscale;
    }
  else
    { // enlarging
    iscale = (int) scale;
    fscale = scale;
    }

  inRS = drms_open_records(drms_env, inQuery, &status);
  if (status || inRS->n == 0)
    DIE("No input data found");
  nrecs = inRS->n;
 if (nrecs == 0)
    DIE("No records found");
  outRS = drms_create_records(drms_env, nrecs, (char *)outSeries, DRMS_PERMANENT, &status);
  if (status)
    DIE("Output recordset not created");

  for (irec=0; irec<nrecs; irec++)
    {
    DRMS_Record_t *outRec, *inRec;
    DRMS_Segment_t *outSeg, *inSeg;
    DRMS_Array_t *inArray, *outArray;
    float *inData, *outData;
    int inx, iny, outx, outy, i, j;
    float val;
    int quality=0;
 
    inRec = inRS->records[irec];
    quality = drms_getkey_int(inRec, "QUALITY", &status);
    if (status || (!status && quality >= 0))
      {
      inSeg = drms_segment_lookupnum(inRec, 0);
      inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
      if (status)
        {
        printf(" No data file found but QUALITY not bad, status=%d\n", status);
        drms_free_array(inArray);
        continue;
        }
  
      int naxis = inArray->naxis;
      int in_nx = inArray->axis[0];
      int in_ny = inArray->axis[1];
      int in_inc = (scale < 0 ? 1 : iscale);
      int out_inc = (scale < 0 ? iscale : 1);
      int out_nx = in_nx * fscale;
      int out_ny = in_ny * fscale;
  
      int outDims[2] = {out_nx, out_ny};
      inData = (float *)inArray->data;
      outArray = drms_array_create(DRMS_TYPE_FLOAT, naxis, outDims, NULL, &status);
      outData = (float *)outArray->data;
  
      if (strcmp(mode, "simple")==0) // setup for better HMI proxy
        {
        if (scale < 0) // need to average input chunk
          {
          for (outy = 0; outy < out_ny; outy += 1)
            for (outx = 0; outx < out_nx; outx += 1)
              {
              double total = 0.0;
              int nn = 0;
              for (j = 0; j < in_inc; j += 1)
                for (i = 0; i < in_inc; i += 1)
                  {
                  val = inData[in_nx*(outy*iscale + j) + outx*iscale + i];
                  if (!drms_ismissing_float(val))
                      {
                       total = total + val; 
                       nn++;
                      }
                  }
              outData[out_nx*outy + outx] = (nn > 0 ? total/nn : DRMS_MISSING_FLOAT); 
              }
          }
        else  // need to replicate input point
          {
          for (iny = 0; iny < in_ny; iny += 1)
            for (inx = 0; inx < in_nx; inx += 1)
              {
              val = inData[in_nx*iny + inx];
              for (j = 0; j < out_inc; j += 1)
                for (i = 0; i < out_inc; i += 1)
                  outData[out_nx*(iny*iscale + j) + inx*iscale + i] = val;
              }
          }
        }
    else
      DIE("invalid conversion mode");
  
    drms_free_array(inArray);
  
    // write data file
    outRec = outRS->records[irec];
    outSeg = drms_segment_lookupnum(outRec, 0);

    // copy all keywords
    drms_copykeys(outRec, inRec, 0, kDRMS_KeyClass_Explicit);
    set_statistics(outSeg, outArray, 1);

    // Now fixup coordinate keywords
    // Only CRPIX1,2 and CDELT1,2 should need repair.
    drms_setkey_double(outRec, "CDELT1", drms_getkey_double(inRec, "CDELT1", NULL)/fscale);
    drms_setkey_double(outRec, "CDELT2", drms_getkey_double(inRec, "CDELT2", NULL)/fscale);
    drms_setkey_double(outRec, "CRPIX1", (drms_getkey_double(inRec, "CRPIX1", NULL)-0.5) * fscale + 0.5);
    drms_setkey_double(outRec, "CRPIX2", (drms_getkey_double(inRec, "CRPIX2", NULL)-0.5) * fscale + 0.5);
    drms_setkey_int(outRec, "SCALE", scale);

  
    // get info for array from segment
    outArray->bzero = outSeg->bzero;
    outArray->bscale = outSeg->bscale;
    outArray->parent_segment = outSeg;
  
    status = drms_segment_write(outSeg, outArray, 0);
    if (status)
      DIE("problem writing file");
    drms_free_array(outArray);
    }
   } //end of "irec" loop

  drms_close_records(inRS, DRMS_FREE_RECORD);
  drms_close_records(outRS, DRMS_INSERT_RECORD);
  return (DRMS_SUCCESS);
  } // end of DoIt
Revision:	1.5
Committed:	Mon May 31 22:50:06 2010 UTC (13 years, 3 months ago) by phil
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Ver_6-0, Ver_6-1, Ver_6-2, Ver_6-3, Ver_6-4, Ver_9-1, Ver_5-14, Ver_5-13, Ver_5-12, Ver_5-11, Ver_5-10, Ver_LATEST, Ver_9-3, Ver_9-41, Ver_9-2, Ver_8-8, Ver_8-2, Ver_8-3, Ver_8-0, Ver_8-1, Ver_8-6, Ver_8-7, Ver_8-4, Ver_8-5, Ver_7-1, Ver_7-0, Ver_9-5, Ver_9-4, Ver_8-10, Ver_8-11, Ver_8-12, Ver_9-0, HEAD
Changes since 1.4:	+94 -159 lines
Log Message:	rewrite to make it work. Now takes int scale >0 for enlarging, <0 to reduce. Fixes WCS keywrods and statistics.
#	Content
1	/**
2	@defgroup rebin2 rebin2 reduce/increase resolution to that of MDI/HMI
3	@ingroup su_util
4
5	@brief Reduce (increase) the resolution of the input data to that of MDI (HMI) by factor of 2.
6
7	@par Synopsis:
8	@code
9	rebin2 in=input data out=output data factor=2*n mode=simple
10	where n is a multiple/fraction of 2.
11	@endcode
12
13	This is a general purpose module that takes a series of input
14	data and modifies its spatial resolution by a factor (multiples
15	or fractions of 2) as required and gives out a set of output
16	data. The method for avaraging (interpolation) can be specified
17	through the input "mode". The current version handles a simple
18	boxcar average. If 'scale' < 0 then the input is reduced in size to 1/\|scale\|.
19
20	Make sure you created the appropriate output series before running
21	the program. For example, su_bala.rebin2up.jsd
22
23	@par Flags:
24	@c none
25	Currently it doesn't have any flags.
26
27	@par GEN_FLAGS:
28	Ubiquitous flags present in every module.
29	@ref jsoc_main
30
31	@param in The input data series.
32	@param out The output series.
33
34	@par Exit_Status:
35	Brief description of abnormal, non-zero, exit values.
36
37	@par Example:
38	Takes a series of 1024 X 1024 MDI full disk magnetogram and
39	produces images with the resolution of HMI, 4096 X 4096.
40
41	@code
42	rebin2 in='mdi.fd_M_96m_lev18[2003.10.20/1d]' out='su_bala.rebin2up' scale=4 mode='simple'
43	@endcode
44
45	@par Example:
46	Reduces the resolution of HMI images 4096 X 4096 to that of MDI,
47	1024 X 1024. Here the input is the HMI images and the output
48	is the lower resolution HMI images.
49	@code
50	rebin2 in='su_bala.rebin2up[2003.10.20/1d]' out='su_bala.rebin2down' scale=0.25 mode='simple'
51	@endcode
52
53	@bug
54	None known so far.
55
56	@par Code:
57	The doxygen code that makes this page is here:
58	@verbinclude ./doxygen_moduletemplate.txt
59	*/
60
61	#include "jsoc.h"
62	#include "jsoc_main.h"
63	#include "fstats.h"
64
65	char *module_name = "rebin2";
66
67	#define DIE(msg) {fflush(stdout);fprintf(stderr,"%s, status=%d\n",msg,status); return(status);}
68
69	ModuleArgs_t module_args[] =
70	{
71	{ARG_STRING, "in", "NOT SPECIFIED", "Input data series."},
72	{ARG_STRING, "out", "NOT SPECIFIED", "Output data series."},
73	{ARG_FLOAT, "scale", "NOTSPECIFIED", "Scale factor. negative for shrink factor"},
74	{ARG_STRING, "mode", "simple", "conversion type."},
75	{ARG_END}
76	};
77
78	int DoIt(void)
79	{
80	int status = DRMS_SUCCESS;
81	DRMS_RecordSet_t inRS, outRS;
82	int irec, nrecs;
83	const char *inQuery = params_get_str(&cmdparams, "in");
84	const char *outSeries = params_get_str(&cmdparams, "out");
85	const char *mode = params_get_str(&cmdparams, "mode");
86	int scale = params_get_int(&cmdparams, "scale");
87	int iscale, factor;
88	float fscale;
89	if (scale < 0.0)
90	{ // shrinking
91	iscale = -scale;
92	fscale = 1.0/iscale;
93	}
94	else
95	{ // enlarging
96	iscale = (int) scale;
97	fscale = scale;
98	}
99
100	inRS = drms_open_records(drms_env, inQuery, &status);
101	if (status \|\| inRS->n == 0)
102	DIE("No input data found");
103	nrecs = inRS->n;
104	if (nrecs == 0)
105	DIE("No records found");
106	outRS = drms_create_records(drms_env, nrecs, (char *)outSeries, DRMS_PERMANENT, &status);
107	if (status)
108	DIE("Output recordset not created");
109
110	for (irec=0; irec<nrecs; irec++)
111	{
112	DRMS_Record_t outRec, inRec;
113	DRMS_Segment_t outSeg, inSeg;
114	DRMS_Array_t inArray, outArray;
115	float inData, outData;
116	int inx, iny, outx, outy, i, j;
117	float val;
118	int quality=0;
119
120	inRec = inRS->records[irec];
121	quality = drms_getkey_int(inRec, "QUALITY", &status);
122	if (status \|\| (!status && quality >= 0))
123	{
124	inSeg = drms_segment_lookupnum(inRec, 0);
125	inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
126	if (status)
127	{
128	printf(" No data file found but QUALITY not bad, status=%d\n", status);
129	drms_free_array(inArray);
130	continue;
131	}
132
133	int naxis = inArray->naxis;
134	int in_nx = inArray->axis[0];
135	int in_ny = inArray->axis[1];
136	int in_inc = (scale < 0 ? 1 : iscale);
137	int out_inc = (scale < 0 ? iscale : 1);
138	int out_nx = in_nx * fscale;
139	int out_ny = in_ny * fscale;
140
141	int outDims[2] = {out_nx, out_ny};
142	inData = (float *)inArray->data;
143	outArray = drms_array_create(DRMS_TYPE_FLOAT, naxis, outDims, NULL, &status);
144	outData = (float *)outArray->data;
145
146	if (strcmp(mode, "simple")==0) // setup for better HMI proxy
147	{
148	if (scale < 0) // need to average input chunk
149	{
150	for (outy = 0; outy < out_ny; outy += 1)
151	for (outx = 0; outx < out_nx; outx += 1)
152	{
153	double total = 0.0;
154	int nn = 0;
155	for (j = 0; j < in_inc; j += 1)
156	for (i = 0; i < in_inc; i += 1)
157	{
158	val = inData[in_nx(outyiscale + j) + outx*iscale + i];
159	if (!drms_ismissing_float(val))
160	{
161	total = total + val;
162	nn++;
163	}
164	}
165	outData[out_nx*outy + outx] = (nn > 0 ? total/nn : DRMS_MISSING_FLOAT);
166	}
167	}
168	else // need to replicate input point
169	{
170	for (iny = 0; iny < in_ny; iny += 1)
171	for (inx = 0; inx < in_nx; inx += 1)
172	{
173	val = inData[in_nx*iny + inx];
174	for (j = 0; j < out_inc; j += 1)
175	for (i = 0; i < out_inc; i += 1)
176	outData[out_nx(inyiscale + j) + inx*iscale + i] = val;
177	}
178	}
179	}
180	else
181	DIE("invalid conversion mode");
182
183	drms_free_array(inArray);
184
185	// write data file
186	outRec = outRS->records[irec];
187	outSeg = drms_segment_lookupnum(outRec, 0);
188
189	// copy all keywords
190	drms_copykeys(outRec, inRec, 0, kDRMS_KeyClass_Explicit);
191	set_statistics(outSeg, outArray, 1);
192
193	// Now fixup coordinate keywords
194	// Only CRPIX1,2 and CDELT1,2 should need repair.
195	drms_setkey_double(outRec, "CDELT1", drms_getkey_double(inRec, "CDELT1", NULL)/fscale);
196	drms_setkey_double(outRec, "CDELT2", drms_getkey_double(inRec, "CDELT2", NULL)/fscale);
197	drms_setkey_double(outRec, "CRPIX1", (drms_getkey_double(inRec, "CRPIX1", NULL)-0.5) * fscale + 0.5);
198	drms_setkey_double(outRec, "CRPIX2", (drms_getkey_double(inRec, "CRPIX2", NULL)-0.5) * fscale + 0.5);
199	drms_setkey_int(outRec, "SCALE", scale);
200
201
202	// get info for array from segment
203	outArray->bzero = outSeg->bzero;
204	outArray->bscale = outSeg->bscale;
205	outArray->parent_segment = outSeg;
206
207	status = drms_segment_write(outSeg, outArray, 0);
208	if (status)
209	DIE("problem writing file");
210	drms_free_array(outArray);
211	}
212	} //end of "irec" loop
213
214	drms_close_records(inRS, DRMS_FREE_RECORD);
215	drms_close_records(outRS, DRMS_INSERT_RECORD);
216	return (DRMS_SUCCESS);
217	} // end of DoIt