1 tplarson 1.1 /*
2 this module (optionally) rebins its input data and then (optionally) convolves with a gaussian and subsamples
3 */
4
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <time.h>
8 #include <math.h>
9 #include <sys/time.h>
10 #include <sys/resource.h>
11 #include "jsoc_main.h"
12 #include "drms_dsdsapi.h"
13 #include "fresize.h"
|
14 tplarson 1.6 #include "fstats.h"
|
15 tplarson 1.1
16 #define ARRLENGTH(ARR) (sizeof(ARR)/sizeof(ARR[0]))
17 #define QUAL_NODATA (0x80000000)
18 #define kNOTSPECIFIED "not specified"
19
20 char *module_name = "jrebinsmooth";
|
21 tplarson 1.6 char *cvsinfo_jrebinsmooth = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jrebinsmooth.c,v 1.5 2013/10/17 23:35:13 tplarson Exp $";
|
22 tplarson 1.1
23 ModuleArgs_t module_args[] =
24 {
25 {ARG_STRING, "in", "", "input data records"},
26 {ARG_STRING, "out", "", "output data series"},
27 {ARG_STRING, "interout", kNOTSPECIFIED, "output data series for intermediate data (binned only)"},
28 {ARG_STRING, "segin", kNOTSPECIFIED, "name of input segment if not using segment 0"},
29 {ARG_STRING, "segout", kNOTSPECIFIED, "name of output segment if not using segment 0"},
30 {ARG_STRING, "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"},
31 {ARG_STRING, "srclink", "SRCDATA", "name of link to source data"},
32 {ARG_INT, "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
33 {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 */
34
35 {ARG_INT, "IBIN", "0", "flag for binning"},
36 {ARG_INT, "NBIN", "-1", "factor for binning"},
37 {ARG_INT, "BIN_XOFF", "0", "offset in pixels in x direction to apply before binning"},
38 {ARG_INT, "BIN_YOFF", "0", "offset in pixels in y direction to apply before binning"},
39 {ARG_FLOAT, "BIN_FILL", "0.0", "value to use outside valid area"},
40
41 {ARG_INT, "IGAUSS", "0", "flag for convolving with gaussian"},
42 {ARG_FLOAT, "GAUSS_SIG", "-1.0", "width of gaussian"},
43 tplarson 1.1 {ARG_INT, "GAUSS_WID", "-1", "half width of kernel"},
44 {ARG_INT, "GAUSS_NSUB", "-1", "distance between sampled points in pixels"},
45 {ARG_INT, "GAUSS_XOFF", "0", "offset in pixels to add to x index of input image"},
46 {ARG_INT, "GAUSS_YOFF", "0", "offset in pixels to add to y index of input image"},
47 {ARG_FLOAT, "GAUSS_FILL", "0.0", "value to use outside valid area"},
48
49 {ARG_END}
50 };
51
52 #include "saveparm.c"
53 #include "timing.c"
54 #include "set_history.c"
55
56 int DoIt(void)
57 {
58
59 int newstat = 0;
60 int status = DRMS_SUCCESS;
61 int quality;
62 TIME trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
63 char trecstr[100];
64 tplarson 1.1 struct fresize_struct fress;
65 int xpixin, ypixin, xpix1, ypix1, xpix2, ypix2;
|
66 tplarson 1.6 float x0in, y0in, xscalein, yscalein, x01, y01, xscale1, yscale1, x02, y02, xscale2, yscale2;
|
67 tplarson 1.1 float *inptr, *outptr;
68
69 DRMS_RecordSet_t *inrecset = NULL;
70 DRMS_Segment_t *segin = NULL;
71 DRMS_Segment_t *segout = NULL;
72 DRMS_Array_t *inarr = NULL;
73 DRMS_Array_t *outarr = NULL;
74 DRMS_Record_t *inrec = NULL;
75 DRMS_Record_t *outrec = NULL;
76
77 int length[2];
78 DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
79 DRMS_RecLifetime_t lifetime;
80
81 int nrecs, irec, error;
82 long long histrecnum=-1;
83
84 int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
85 int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);
86
87 char *inrecquery = (char *)cmdparams_save_str(&cmdparams, "in", &newstat);
88 tplarson 1.1 char *outseries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat);
89 char *segnamein = (char *)cmdparams_save_str(&cmdparams, "segin", &newstat);
90 char *segnameout = (char *)cmdparams_save_str(&cmdparams, "segout", &newstat);
91 int seginflag = strcmp(kNOTSPECIFIED, segnamein);
92 int segoutflag = strcmp(kNOTSPECIFIED, segnameout);
93 char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);
94 char *srclinkname = (char *)cmdparams_save_str(&cmdparams, "srclink", &newstat);
95 char *interout = (char *)cmdparams_save_str(&cmdparams, "interout", &newstat);
96 int interflag = strcmp(kNOTSPECIFIED, interout);
97 float *interptr=NULL;
98
99 int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
100 int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
101 if (permflag)
102 lifetime = DRMS_PERMANENT;
103 else
104 lifetime = DRMS_TRANSIENT;
105
106 int ibin = cmdparams_save_int(&cmdparams, "IBIN", &newstat);
107 int nbin = cmdparams_save_int(&cmdparams, "NBIN", &newstat);
108 int bxoff = cmdparams_save_int(&cmdparams, "BIN_XOFF", &newstat);
109 tplarson 1.1 int byoff = cmdparams_save_int(&cmdparams, "BIN_YOFF", &newstat);
110 float bfill = cmdparams_save_float(&cmdparams, "BIN_FILL", &newstat);
111
112 int igauss = cmdparams_save_int(&cmdparams, "IGAUSS", &newstat);
113 float sigma = cmdparams_save_float(&cmdparams, "GAUSS_SIG", &newstat);
114 int hwidth = cmdparams_save_int(&cmdparams, "GAUSS_WID", &newstat);
115 int nsub = cmdparams_save_int(&cmdparams, "GAUSS_NSUB", &newstat);
116 int gxoff = cmdparams_save_int(&cmdparams, "GAUSS_XOFF", &newstat);
117 int gyoff = cmdparams_save_int(&cmdparams, "GAUSS_YOFF", &newstat);
118 float gfill = cmdparams_save_float(&cmdparams, "GAUSS_FILL", &newstat);
119
120 if (newstat)
121 {
122 fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
123 cpsave_decode_error(newstat);
124 return 1;
125 }
126 else if (savestrlen != strlen(savestr))
127 {
128 fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr));
129 return 1;
130 tplarson 1.1 }
131
132 if (ibin == 0 && igauss == 0)
133 {
134 fprintf(stderr,"ERROR: must set at least one of IBIN and IGAUSS\n");
135 return 1;
136 }
137
138 if (interflag != 0 && (ibin == 0 || igauss == 0))
139 {
140 fprintf(stderr,"ERROR: cannot specify interout unless both IBIN and IGAUSS are set\n");
141 return 1;
142 }
143
144 if (ibin != 0 && nbin < 0)
145 {
146 fprintf(stderr,"ERROR: must specify NBIN > 0 when IBIN is set\n");
147 return 1;
148 }
149
150 if (igauss != 0 && (sigma < 0.0 || hwidth < 0 || nsub < 0))
151 tplarson 1.1 {
152 fprintf(stderr,"ERROR: must specify nonnegative value for each of GAUSS_SIG, GAUSS_WID, and GAUSS_NSUB when IGAUSS is set\n");
153 return 1;
154 }
155
156 if (ibin == 0)
157 nbin=1;
158 if (igauss == 0)
159 nsub=1;
160
161
162 DRMS_Record_t *tempoutrec = drms_create_record(drms_env, outseries, DRMS_TRANSIENT, &status);
163
164 if (status != DRMS_SUCCESS || tempoutrec == NULL)
165 {
166 fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outseries, status);
167 return 1;
168 }
169
170 DRMS_Link_t *histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);
171 DRMS_Link_t *srclink = hcon_lookup_lower(&tempoutrec->links, srclinkname);
172 tplarson 1.1
|
173 tplarson 1.6 // char *cvsinfo = strdup("$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jrebinsmooth.c,v 1.5 2013/10/17 23:35:13 tplarson Exp $");
|
174 tplarson 1.1 if (histlink != NULL)
175 {
|
176 tplarson 1.4 histrecnum=set_history(histlink);
|
177 tplarson 1.1 if (histrecnum < 0)
178 {
179 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
180 return 1;
181 }
182 }
183 else
184 {
185 fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
186 }
187
188 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
189
190 inrecset = drms_open_records(drms_env, inrecquery, &status);
191 nrecs = inrecset->n;
192 if (status != DRMS_SUCCESS || inrecset == NULL)
193 {
194 fprintf(stderr, "ERROR: problem opening input recordset: status = %d\n", status);
195 return 1;
196 }
197
198 tplarson 1.1 if (nrecs == 0)
199 {
200 printf("ERROR: input recordset contains no records\n");
201 drms_close_records(inrecset, DRMS_FREE_RECORD);
202 return 1;
203 }
204
205 if (verbflag)
206 printf("input recordset opened, nrecs = %d\n",nrecs);
207
208 char *inchecklist[] = {"T_REC", "QUALITY"};
209 inrec=inrecset->records[0];
210 int itest;
211 for (itest=0; itest < ARRLENGTH(inchecklist); itest++)
212 {
213 DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]);
214 if (!inkeytest)
215 {
216 fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
217 drms_close_records(inrecset, DRMS_FREE_RECORD);
218 return 1;
219 tplarson 1.1 }
220 }
221
222
223 init_fresize_gaussian(&fress, sigma, hwidth, nsub);
224 // init_fresize_gaussian_fft(&fress, sigma, hwidth, nsub, nxin, nyin);
225
|
226 tplarson 1.6 DRMS_RecordSet_t *outrecset = drms_create_records(drms_env, nrecs, outseries, lifetime, &status);
227 if (status != DRMS_SUCCESS || outrecset == NULL)
228 {
229 fprintf(stderr,"ERROR: unable to create record in output dataseries, status = %d\n", status);
230 return 1;
231 }
232
233 DRMS_RecordSet_t *interoutrecset;
234 if (interflag)
235 {
236 interoutrecset = drms_create_records(drms_env, nrecs, interout, lifetime, &status);
237 if (status != DRMS_SUCCESS || interoutrecset == NULL)
238 {
239 fprintf(stderr,"ERROR: unable to create record in intermediate output dataseries, status = %d\n", status);
240 return 1;
241 }
242 }
243
|
244 tplarson 1.1 error=0;
|
245 tplarson 1.6 int nodata=0;
|
246 tplarson 1.1 for (irec=0; irec<nrecs; irec++)
247 {
248
249 inrec = inrecset->records[irec];
|
250 tplarson 1.6 outrec = outrecset->records[irec];
|
251 tplarson 1.1
252 if (histlink)
253 drms_setlink_static(outrec, histlinkname, histrecnum);
254 if (srclink)
255 drms_setlink_static(outrec, srclinkname, inrec->recnum);
256
257 trec = drms_getkey_time(inrec, "T_REC", &status);
258 if (status != DRMS_SUCCESS)
259 {
260 fprintf(stderr, "SKIP: error getting primekey T_REC: irec = %d, status = %d, recnum = %lld\n", irec, status, inrec->recnum);
261 continue;
262 }
263 sprint_time(trecstr, trec, "TAI", 0);
264
265 quality = drms_getkey_int(inrec, "QUALITY", &status);
266 if (status != DRMS_SUCCESS)
267 {
268 fprintf(stderr, "SKIP: error getting keyword QUALITY: irec = %d, status = %d, T_REC = %s, recnum = %lld\n", irec, status, trecstr, inrec->recnum);
269 continue;
270 }
271
272 tplarson 1.1 if (quality & QUAL_NODATA)
273 {
|
274 tplarson 1.6 nodata=1;
275 if (verbflag > 1)
276 fprintf(stderr,"SKIP: record rejected based on quality = %d: irec = %d, T_REC = %s, recnum = %lld\n", quality, irec, trecstr, inrec->recnum);
277 // continue;
|
278 tplarson 1.1 }
279 else
|
280 tplarson 1.6 {
281 nodata=0;
|
282 tplarson 1.1
|
283 tplarson 1.6 if (seginflag)
284 segin = drms_segment_lookup(inrec, segnamein);
285 else
286 segin = drms_segment_lookupnum(inrec, 0);
|
287 tplarson 1.1
|
288 tplarson 1.6 if (segin)
289 inarr = drms_segment_read(segin, usetype, &status);
|
290 tplarson 1.1
|
291 tplarson 1.6 if (segin == NULL || inarr == NULL || status != DRMS_SUCCESS)
292 {
293 fprintf(stderr, "ERROR: problem with input segment or array: irec = %d, status = %d, T_REC = %s, recnum = %lld \n", irec, status, trecstr, inrec->recnum);
294 if (inarr)
295 drms_free_array(inarr);
296 error=1;
297 break;
298 }
299
300 xpixin=inarr->axis[0];
301 ypixin=inarr->axis[1];
302 xpix1=xpixin/nbin;
303 ypix1=ypixin/nbin;
304 xpix2=xpix1/nsub;
305 ypix2=ypix1/nsub;
306
307 x0in = drms_getkey_float(inrec, "CRPIX1", &status) - 1;
308 y0in = drms_getkey_float(inrec, "CRPIX2", &status) - 1;
309 xscalein = drms_getkey_float(inrec, "CDELT1", &status);
310 yscalein = drms_getkey_float(inrec, "CDELT2", &status);
311
312 tplarson 1.6 x01 = (x0in - bxoff + 0.5)/nbin - 0.5;
313 y01 = (y0in - byoff + 0.5)/nbin - 0.5;
314 xscale1=xscalein*nbin;
315 yscale1=yscalein*nbin;
316
317 if (ibin)
318 {
319 interptr=(float *)malloc(xpix1*ypix1*sizeof(float));
320 fbin(
|
321 tplarson 1.1 (float *)inarr->data,
322 interptr,
323 xpixin,
324 ypixin,
325 xpixin,
326 xpix1,
327 ypix1,
328 xpix1,
329 nbin,
330 bxoff,
331 byoff,
|
332 tplarson 1.6 bfill);
333 }
|
334 tplarson 1.1 }
335
|
336 tplarson 1.6 if (interflag && nodata)
337 {
338 DRMS_Record_t *outrec = interoutrecset->records[irec];
339 DRMS_Link_t *histlink = hcon_lookup_lower(&outrec->links, histlinkname);
340 DRMS_Link_t *srclink = hcon_lookup_lower(&outrec->links, srclinkname);
341 if (histlink)
342 drms_setlink_static(outrec, histlinkname, histrecnum);
343 if (srclink)
344 drms_setlink_static(outrec, srclinkname, inrec->recnum);
345 drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
346 tnow = (double)time(NULL);
347 tnow += UNIX_epoch;
348 drms_setkey_time(outrec, "DATE", tnow);
349
350 }
351 else if (interflag)
|
352 tplarson 1.1 {
|
353 tplarson 1.6 // DRMS_Record_t *outrec = drms_create_record(drms_env, interout, lifetime, &status);
354 DRMS_Record_t *outrec = interoutrecset->records[irec];
|
355 tplarson 1.1 DRMS_Link_t *histlink = hcon_lookup_lower(&outrec->links, histlinkname);
356 DRMS_Link_t *srclink = hcon_lookup_lower(&outrec->links, srclinkname);
357 if (histlink)
358 drms_setlink_static(outrec, histlinkname, histrecnum);
359 if (srclink)
360 drms_setlink_static(outrec, srclinkname, inrec->recnum);
361 if (segoutflag)
362 segout = drms_segment_lookup(outrec, segnameout);
363 else
364 segout = drms_segment_lookupnum(outrec, 0);
365 length[0]=xpix1;
366 length[1]=ypix1;
367 outarr = drms_array_create(usetype, 2, length, interptr, &status);
368 outarr->bzero=segout->bzero;
369 outarr->bscale=segout->bscale;
|
370 tplarson 1.6 status=drms_segment_write(segout, outarr, 0);
|
371 tplarson 1.1 drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
|
372 tplarson 1.5 drms_setkey_int(outrec, "TOTVALS", xpix1*ypix1);
373 set_statistics(segout, outarr, 1);
|
374 tplarson 1.1 drms_setkey_time(outrec, "T_REC", trec);
375 drms_setkey_int(outrec, "QUALITY", quality);
376 drms_setkey_float(outrec, "CRPIX1", x01+1);
377 drms_setkey_float(outrec, "CRPIX2", y01+1);
378 drms_setkey_float(outrec, "CDELT1", xscale1);
379 drms_setkey_float(outrec, "CDELT2", yscale1);
380 tnow = (double)time(NULL);
381 tnow += UNIX_epoch;
382 drms_setkey_time(outrec, "DATE", tnow);
|
383 tplarson 1.6 // drms_close_record(outrec, DRMS_INSERT_RECORD);
384 }
385
386 if (nodata)
387 {
388 drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
389 tnow = (double)time(NULL);
390 tnow += UNIX_epoch;
391 drms_setkey_time(outrec, "DATE", tnow);
392 continue;
|
393 tplarson 1.1 }
394
395 length[0]=xpix2;
396 length[1]=ypix2;
397
398 if (segoutflag)
399 segout = drms_segment_lookup(outrec, segnameout);
400 else
401 segout = drms_segment_lookupnum(outrec, 0);
402
403 if (!igauss)
404 outarr = drms_array_create(usetype, 2, length, interptr, &status);
405 else
406 outarr = drms_array_create(usetype, 2, length, NULL, &status);
407
408 if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
409 {
410 fprintf(stderr, "ERROR: problem with output segment or array: irec = %d, status = %d\n", irec, status);
411 drms_free_array(inarr);
412 if (outarr)
413 drms_free_array(outarr);
414 tplarson 1.1 error=1;
415 break;
416 }
417
418 if (!ibin)
419 inptr=(float *)inarr->data;
420 else
421 inptr=interptr;
422
423 if (igauss)
424 {
425 fresize(
426 &fress, // Must have been initialized by init_fresize_XXX
427 inptr,
428 (float *)(outarr->data),
429 xpix1, // Size of input image
430 ypix1, // Size of input image
431 xpix1, // Leading dimension of input image. nleadin>=nxin
432 xpix2, // Size of xin, yin and image_out
433 ypix2, // Size of xin, yin and image_out
434 xpix2, // Leading dimension. nlead>=nx
435 tplarson 1.1 gxoff, // Offset in x direction
436 gyoff, // Offset in y direction
437 gfill // Value to use if outside area
438 );
439 }
440
|
441 tplarson 1.6 x02 = (x01 - gxoff)/nsub;
442 y02 = (y01 - gyoff)/nsub;
443 xscale2=xscale1*nsub;
444 yscale2=yscale1*nsub;
|
445 tplarson 1.1
446 outarr->bzero=segout->bzero;
447 outarr->bscale=segout->bscale;
|
448 tplarson 1.6 status=drms_segment_write(segout, outarr, 0);
|
449 tplarson 1.1 drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
|
450 tplarson 1.5 drms_setkey_int(outrec, "TOTVALS", xpix2*ypix2);
451 set_statistics(segout, outarr, 1);
|
452 tplarson 1.1 drms_setkey_time(outrec, "T_REC", trec);
453 drms_setkey_int(outrec, "QUALITY", quality);
454 drms_setkey_float(outrec, "CRPIX1", x02+1);
455 drms_setkey_float(outrec, "CRPIX2", y02+1);
456 drms_setkey_float(outrec, "CDELT1", xscale2);
457 drms_setkey_float(outrec, "CDELT2", yscale2);
458 tnow = (double)time(NULL);
459 tnow += UNIX_epoch;
460 drms_setkey_time(outrec, "DATE", tnow);
461
|
462 tplarson 1.6 // drms_close_record(outrec, DRMS_INSERT_RECORD);
|
463 tplarson 1.1 drms_free_array(inarr);
464 drms_free_array(outarr);
465 if (interptr != NULL)
466 free(interptr);
467
468 }
469
470 drms_close_records(inrecset, DRMS_FREE_RECORD);
|
471 tplarson 1.6 if (error)
472 {
473 drms_close_records(outrecset, DRMS_FREE_RECORD);
474 if (interflag)
475 drms_close_records(interoutrecset, DRMS_FREE_RECORD);
476 }
477 else
478 {
479 drms_close_records(outrecset, DRMS_INSERT_RECORD);
480 if (interflag)
481 drms_close_records(interoutrecset, DRMS_INSERT_RECORD);
482 }
|
483 tplarson 1.1
484 free_fresize(&fress);
485
486 if (verbflag)
487 {
488 if (!error)
489 printf("module %s successful completion\n", cmdparams.argv[0]);
490 }
491
492 return 0;
493 }
|
Return to
jrebinsmooth.c
CVS log
Up to [Development] / JSOC / proj / globalhs / apps