| 1 |
/* adapted from proj/lev1.5_aia/aia_lev1p5.c */ |
| 2 |
/* restrictions: should not resize by "too much" since does simple bicubic mapping |
| 3 |
and large changes will not preserve all information. Also, the rotation is done |
| 4 |
about CRPIX1 and CRPIX2 which must be contained in the image, or at least |
| 5 |
the final location must be contained in the image after rotation. A patch from |
| 6 |
HMI will result in all missing data. |
| 7 |
*/ |
| 8 |
|
| 9 |
#include <string.h> |
| 10 |
#include "jsoc_main.h" |
| 11 |
#include "drms.h" |
| 12 |
#include "drms_names.h" |
| 13 |
|
| 14 |
#define NOT_SPECIFIED "***Not Specified***" |
| 15 |
#define DIE(msg) {fprintf(stderr,"$$$$ %s: %s\n", module_name, msg); return 1;} |
| 16 |
|
| 17 |
#define NBINS 1048576 |
| 18 |
static int hist[NBINS]; |
| 19 |
|
| 20 |
int image_magrotate(void *, int nin, int min, int data_type_input, float theta, float mag, |
| 21 |
float dx, float dy, void **outarray, int *nx, int *ny, int regridtype_input, int stretch_lines); |
| 22 |
|
| 23 |
int drms_appendcomment(DRMS_Record_t *rec, const char *comment, int nl); |
| 24 |
|
| 25 |
ModuleArgs_t module_args[] = |
| 26 |
{ |
| 27 |
{ARG_STRING, "in", NOT_SPECIFIED, "Input series query"}, |
| 28 |
{ARG_STRING, "out", NOT_SPECIFIED, "Output series"}, |
| 29 |
// {ARG_STRING, "despike", "0", "remove cosmic ray hits"}, |
| 30 |
{ARG_INT, "rescale", "0", "rescale to fixed plate scale"}, |
| 31 |
{ARG_INT, "regrid", "1", "regrid type 0: nearest neighbor, 1: bicubic"}, |
| 32 |
{ARG_INT, "center_to", "0", "center to 0: Sun center, 1: First image, 2: No change"}, |
| 33 |
{ARG_DOUBLE, "scale_to", "0.6", "rescale to new CDELT scale"}, |
| 34 |
{ARG_STRING, "do_stretchmarks", "0", "replicate pixels created else use missing value"}, |
| 35 |
{ARG_STRING, "requestid", NOT_SPECIFIED, "RequestID if used for export"}, |
| 36 |
{ARG_FLAG, "c", "0", "crop image to RSUN_OBS-1/4 pixel"}, |
| 37 |
{ARG_FLAG, "h", "0", "Print usage message and quit"}, |
| 38 |
{ARG_FLAG, "v", "0", "verbose flag"}, |
| 39 |
{ARG_END} |
| 40 |
}; |
| 41 |
|
| 42 |
char *module_name = "jsoc_resize"; |
| 43 |
|
| 44 |
|
| 45 |
// Definitions at top of file |
| 46 |
// --------------------------------------------------------------------- |
| 47 |
|
| 48 |
#define Deg2Rad (M_PI/180.0) |
| 49 |
#define Rad2arcsec (3600.0/Deg2Rad) |
| 50 |
#define arcsec2Rad (Deg2Rad/3600.0) |
| 51 |
#define Rad2Deg (180.0/M_PI) |
| 52 |
|
| 53 |
struct ObsInfo_struct |
| 54 |
{ |
| 55 |
// from observation info |
| 56 |
TIME t_obs; |
| 57 |
double rsun_obs, obs_vr, obs_vw, obs_vn; |
| 58 |
double crpix1, crpix2, cdelt1, cdelt2, crota2; |
| 59 |
double crval1, crval2; |
| 60 |
double cosa, sina; |
| 61 |
double obs_b0; |
| 62 |
// observed point |
| 63 |
int i,j; |
| 64 |
// parameters for observed point |
| 65 |
double x,y,r; |
| 66 |
double rho; |
| 67 |
double lon; |
| 68 |
double lat; |
| 69 |
double sinlat, coslat; |
| 70 |
double sig; |
| 71 |
double mu; |
| 72 |
double chi; |
| 73 |
double obs_v; |
| 74 |
}; |
| 75 |
|
| 76 |
typedef struct ObsInfo_struct ObsInfo_t; |
| 77 |
|
| 78 |
void rebinArraySF(DRMS_Array_t *out, DRMS_Array_t *in); |
| 79 |
int upNcenter(DRMS_Array_t *arr, ObsInfo_t *ObsLoc); |
| 80 |
int crop_image(DRMS_Array_t *arr, ObsInfo_t *ObsLoc); |
| 81 |
const char *get_input_recset(DRMS_Env_t *drms_env, const char *in, DRMS_Record_t **template); |
| 82 |
|
| 83 |
ObsInfo_t *GetObsInfo(DRMS_Segment_t *seg, ObsInfo_t *pObsLoc, int *rstatus); |
| 84 |
|
| 85 |
// --------------------------------------------------------------------- |
| 86 |
|
| 87 |
int verbose; |
| 88 |
|
| 89 |
int nice_intro(int help) |
| 90 |
{ |
| 91 |
int usage = cmdparams_get_int(&cmdparams, "h", NULL) != 0; |
| 92 |
verbose = cmdparams_get_int(&cmdparams, "v", NULL) != 0; |
| 93 |
if (usage || help) { |
| 94 |
printf("jsoc_resize {-h} {-v} dsinp=series_record_spec dsout=output_series\n" |
| 95 |
" -h: print this message\n" |
| 96 |
" -v: verbose\n" |
| 97 |
// "despike=0, do not despike, default=0, despike\n" |
| 98 |
"rescale=0, do not rescale to fixed plate scale, default=0\n" |
| 99 |
"regrid=0, nearest neighbor, default=1, bicubic\n" |
| 100 |
"center_to=0, centering target, 0=Sun, 1=First_image, 2=no change; default=0\n" |
| 101 |
"scale_to=0.6, rescale to fixed plate scale, default=0.6 arcsec/pixel\n" |
| 102 |
"do_stretchmarks=0, fill in empty pixels created, default=0\n" |
| 103 |
"in=<recordset query> as <series>{[record specifier]} - required\n" |
| 104 |
"out=<series> - required\n"); |
| 105 |
return(1); |
| 106 |
} |
| 107 |
return(0); |
| 108 |
} |
| 109 |
|
| 110 |
int set_statistics(DRMS_Segment_t *seg, DRMS_Array_t *arr, int mode); |
| 111 |
|
| 112 |
int DoIt () |
| 113 |
{ |
| 114 |
if (nice_intro(0)) return(0); |
| 115 |
|
| 116 |
int irec, iseg, nrecs, nsegs, status, is_aia=0; |
| 117 |
char now_str[100]; |
| 118 |
float crpix1, crpix2, cdelt1, cdelt2, crota2, x0, y0; |
| 119 |
float mag = 1.0, dx, dy; |
| 120 |
DRMS_Record_t *inprec, *outrec; |
| 121 |
DRMS_Record_t *template = NULL; |
| 122 |
DRMS_RecordSet_t *inprs; |
| 123 |
DRMS_Keyword_t *inpkey = NULL, *outkey = NULL; |
| 124 |
DRMS_Array_t *inparr=NULL, *outarr=NULL; |
| 125 |
DRMS_Segment_t *inpseg, *outseg; |
| 126 |
const char *inQuery; |
| 127 |
|
| 128 |
const char *requestid = cmdparams_get_str(&cmdparams, "requestid", NULL); |
| 129 |
char *dsinp = strdup(cmdparams_get_str(&cmdparams, "in", NULL)); |
| 130 |
const char *dsout = cmdparams_get_str(&cmdparams, "out", NULL); |
| 131 |
if (strcmp(dsinp, NOT_SPECIFIED)==0) DIE("\"in\" argument is required"); |
| 132 |
if (strcmp(dsout, NOT_SPECIFIED)==0) |
| 133 |
{ |
| 134 |
char newout[DRMS_MAXNAMELEN]; |
| 135 |
char *inp, inseries[DRMS_MAXNAMELEN]; |
| 136 |
strncpy(inseries, dsinp, DRMS_MAXNAMELEN); |
| 137 |
inp = index(inseries, '['); |
| 138 |
if (inp) |
| 139 |
*inp = '\0'; |
| 140 |
strncpy(newout, inseries, DRMS_MAXNAMELEN); |
| 141 |
strncat(newout, "_resized", DRMS_MAXNAMELEN); |
| 142 |
dsout = strdup(newout); |
| 143 |
} |
| 144 |
else |
| 145 |
dsout = strdup(dsout); |
| 146 |
inQuery = get_input_recset(drms_env, dsinp, &template); |
| 147 |
// int despike = cmdparams_get_int(&cmdparams, "despike", NULL); |
| 148 |
int rescale = cmdparams_get_int(&cmdparams, "rescale", NULL); |
| 149 |
double scale_to = cmdparams_get_double(&cmdparams, "scale_to", NULL); |
| 150 |
int regridtype = cmdparams_get_int(&cmdparams, "regrid", NULL); |
| 151 |
int center_to = cmdparams_get_int(&cmdparams, "center_to", NULL); |
| 152 |
int do_stretchmarks = cmdparams_get_int(&cmdparams, "do_stretchmarks", NULL); |
| 153 |
int do_crop = cmdparams_get_int(&cmdparams, "c", NULL) != 0; |
| 154 |
|
| 155 |
double scale_by; |
| 156 |
char comment[4096]; |
| 157 |
|
| 158 |
if (strstr(dsinp, "aia")) is_aia = 1; |
| 159 |
inprs = drms_open_records(drms_env, inQuery, &status); |
| 160 |
if (status) DIE("cant open recordset query"); |
| 161 |
if (dsinp != inQuery && *inQuery == '@') |
| 162 |
unlink(inQuery+1); |
| 163 |
drms_stage_records(inprs, 1, 0); |
| 164 |
nrecs = inprs->n; |
| 165 |
printf("%d records\n", nrecs); |
| 166 |
float usedx = 0; |
| 167 |
float usedy = 0; |
| 168 |
float usemag = 1.0; |
| 169 |
|
| 170 |
int iSet; |
| 171 |
int iSubRec; |
| 172 |
int nSubRecs = 0; |
| 173 |
|
| 174 |
|
| 175 |
/* We need the record-set specfication for each record, so loop by recordset subset. */ |
| 176 |
for (iSet = 0, irec = 0; iSet < inprs->ss_n; iSet++) |
| 177 |
{ |
| 178 |
int hasSeglist = 0; |
| 179 |
|
| 180 |
/* For each record, we need to know if the segments it contains came from a record-set specification that had a segment list. |
| 181 |
* To do that, we have to search for '{' in the specification that resolved to this record. So we need to spec for this record, |
| 182 |
* which is in inprs->ss_queries[iSet]. */ |
| 183 |
if (strchr(inprs->ss_queries[iSet], '{')) |
| 184 |
{ |
| 185 |
hasSeglist = 1; |
| 186 |
} |
| 187 |
|
| 188 |
nSubRecs = drms_recordset_getssnrecs(inprs, iSet, &status); |
| 189 |
|
| 190 |
for (iSubRec = 0; iSubRec < nSubRecs; iSubRec++, irec++) |
| 191 |
{ |
| 192 |
int quality; |
| 193 |
TIME t_rec, t_obs; |
| 194 |
ObsInfo_t *ObsLoc; |
| 195 |
char inQuery[DRMS_MAXQUERYLEN]; |
| 196 |
|
| 197 |
inprec = inprs->records[(inprs->ss_starts)[iSet] + iSubRec]; |
| 198 |
|
| 199 |
t_obs = drms_getkey_time(inprec, "T_OBS", &status); if (status) DIE("T_OBS not found!"); |
| 200 |
if (t_obs < 0) |
| 201 |
continue; |
| 202 |
if (is_aia) |
| 203 |
{ |
| 204 |
int qualmask = 0x800100f0; |
| 205 |
int quallev0; |
| 206 |
quallev0 = drms_getkey_int(inprec, "QUALLEV0", &status); if (status) DIE("QUALLEV0 not found!"); |
| 207 |
if (quallev0 & qualmask) |
| 208 |
continue; |
| 209 |
} |
| 210 |
quality = drms_getkey_int(inprec, "QUALITY", &status); |
| 211 |
if (verbose) |
| 212 |
fprintf(stderr,"rec %d of %d, quality=%X\n",irec,nrecs,quality); |
| 213 |
if (quality < 0) |
| 214 |
continue; |
| 215 |
|
| 216 |
outrec = drms_create_record(drms_env, (char*)dsout, DRMS_PERMANENT, &status); |
| 217 |
if (status) DIE("cant create recordset"); |
| 218 |
status = drms_copykeys(outrec, inprec, 0, kDRMS_KeyClass_Explicit); |
| 219 |
if (status) DIE("Error in drms_copykeys()"); |
| 220 |
if (strcmp(requestid, NOT_SPECIFIED) != 0) |
| 221 |
drms_setkey_string(outrec, "requestid", requestid); |
| 222 |
drms_setkey_time(outrec, "DATE", CURRENT_SYSTEM_TIME); |
| 223 |
drms_sprint_rec_query(inQuery, inprec); |
| 224 |
drms_setkey_string(outrec, "SOURCE", inQuery); |
| 225 |
|
| 226 |
if (is_aia) |
| 227 |
{ |
| 228 |
if (!drms_keyword_lookup(inprec, "T_REC", 1) && drms_keyword_lookup(outrec, "T_REC", 1)) |
| 229 |
{ |
| 230 |
double tr_step; |
| 231 |
long long tr_index; |
| 232 |
TIME tr_epoch; |
| 233 |
if ( 0 == drms_setkey_time(outrec, "T_REC", t_obs)) |
| 234 |
{ |
| 235 |
tr_index = drms_getkey_longlong(outrec, "T_REC_index", &status); if (status) DIE("T_REC_index not found!"); |
| 236 |
tr_step = drms_getkey_double(outrec, "T_REC_step", &status); if (status) DIE("T_REC_step not found!"); |
| 237 |
tr_epoch = drms_getkey_time(outrec, "T_REC_epoch", &status); if (status) DIE("T_REC_epoch not found!"); |
| 238 |
t_rec = tr_epoch + tr_index*tr_step; |
| 239 |
drms_setkey_time(outrec, "T_REC", t_rec); |
| 240 |
} |
| 241 |
} |
| 242 |
if (!drms_keyword_lookup(inprec, "GAEX_OBS", 1)) |
| 243 |
drms_setkey_double(outrec, "GAEX_OBS", drms_getkey_double(inprec, "GCIEC_X", NULL)); |
| 244 |
if (!drms_keyword_lookup(inprec, "GAEY_OBS", 1)) |
| 245 |
drms_setkey_double(outrec, "GAEY_OBS", drms_getkey_double(inprec, "GCIEC_Y", NULL)); |
| 246 |
if (!drms_keyword_lookup(inprec, "GAEZ_OBS", 1)) |
| 247 |
drms_setkey_double(outrec, "GAEZ_OBS", drms_getkey_double(inprec, "GCIEC_Z", NULL)); |
| 248 |
if (!drms_keyword_lookup(inprec, "HAEX_OBS", 1)) |
| 249 |
drms_setkey_double(outrec, "HAEX_OBS", drms_getkey_double(inprec, "HCIEC_X", NULL)); |
| 250 |
if (!drms_keyword_lookup(inprec, "HAEY_OBS", 1)) |
| 251 |
drms_setkey_double(outrec, "HAEY_OBS", drms_getkey_double(inprec, "HCIEC_Y", NULL)); |
| 252 |
if (!drms_keyword_lookup(inprec, "HAEZ_OBS", 1)) |
| 253 |
drms_setkey_double(outrec, "HAEZ_OBS", drms_getkey_double(inprec, "HCIEC_Z", NULL)); |
| 254 |
} |
| 255 |
|
| 256 |
crpix1 = drms_getkey_double(inprec, "CRPIX1", &status); if (status) DIE("CRPIX1 not found!"); |
| 257 |
crpix2 = drms_getkey_double(inprec, "CRPIX2", &status); if (status) DIE("CRPIX2 not found!"); |
| 258 |
cdelt1 = drms_getkey_double(inprec, "CDELT1", &status); if (status) DIE("CDELT1 not found!"); |
| 259 |
if (rescale) |
| 260 |
{ |
| 261 |
mag = fabs(cdelt1 / scale_to); |
| 262 |
cdelt1 /= mag; |
| 263 |
} |
| 264 |
crota2 = drms_getkey_double(inprec, "CROTA2", &status); if (status) crota2 = 0.0; // WCS default |
| 265 |
|
| 266 |
/* Segment loop. */ |
| 267 |
HIterator_t *segIter = NULL; |
| 268 |
DRMS_Segment_t *orig = NULL; |
| 269 |
|
| 270 |
/* inpseg could be a linked segment, but we need to obtain the segment struct for the original series. |
| 271 |
* There is no pointer from the linked segment to the original segment. And we do not have the name or |
| 272 |
* number of the original segment. Or we didn't, so I made drms_record_nextseg2(). */ |
| 273 |
while ((inpseg = drms_record_nextseg2(inprec, &segIter, 1, &orig)) != NULL) |
| 274 |
{ |
| 275 |
/* CAVEAT: There is no check to see if the segment on which this code is operating is a suitable |
| 276 |
* segment for this module to process. For example, there is the assumption that the segment |
| 277 |
* being processed has two dimensions. Rejecting incompatible segments is something to implement |
| 278 |
* in the future. */ |
| 279 |
void *output_array = NULL; |
| 280 |
int i, ix, iy, n, m, dtyp, nx, ny, npix; |
| 281 |
|
| 282 |
{ |
| 283 |
inparr = drms_segment_read(inpseg, DRMS_TYPE_FLOAT, &status); if (status) DIE("drms_segment_read failed!"); |
| 284 |
|
| 285 |
/* The original intent was for the output series' segments to parallel the input series' segments. |
| 286 |
* However, the names were allowed to vary. So the segment descriptions match, and they are in the |
| 287 |
* same order. However, the names given to these descriptions might vary. Therefore, we need to |
| 288 |
* look-up output segments by segment number, not name. |
| 289 |
* |
| 290 |
* Use the segment number of the input segment (the original input segment, in case a link was followed). |
| 291 |
*/ |
| 292 |
outseg = drms_segment_lookupnum(outrec, orig->info->segnum); |
| 293 |
|
| 294 |
if (!outseg) |
| 295 |
{ |
| 296 |
DIE("Cant get output segment"); |
| 297 |
} |
| 298 |
|
| 299 |
int outdims[2]; |
| 300 |
dtyp = 3; |
| 301 |
n = inparr->axis[0]; |
| 302 |
m = inparr->axis[1]; |
| 303 |
|
| 304 |
if (is_aia) // set outer rows and cols to 0 |
| 305 |
{ |
| 306 |
float *fval; |
| 307 |
fval = inparr->data; |
| 308 |
for (ix = 0; ix<n; ix++) |
| 309 |
{ |
| 310 |
fval[ix] = 0.0; |
| 311 |
fval[(m - 1)*n + ix] = 0.0; |
| 312 |
} |
| 313 |
for (iy = 0; iy<m; iy++) |
| 314 |
{ |
| 315 |
fval[iy*n] = 0.0; |
| 316 |
fval[iy*n + n - 1] = 0.0; |
| 317 |
} |
| 318 |
} |
| 319 |
else |
| 320 |
{ |
| 321 |
ObsLoc = GetObsInfo(inpseg, NULL, &status); |
| 322 |
upNcenter(inparr, ObsLoc); |
| 323 |
if (do_crop) |
| 324 |
crop_image(inparr, ObsLoc); |
| 325 |
crota2 = ObsLoc->crota2; |
| 326 |
crpix1 = ObsLoc->crpix1; |
| 327 |
crpix2 = ObsLoc->crpix2; |
| 328 |
} |
| 329 |
|
| 330 |
if (center_to == 0) |
| 331 |
{ |
| 332 |
usedx = (n + 1.0)/2.0; // center for corner is (1,1) as is crpix |
| 333 |
usedy = (m + 1.0)/2.0; |
| 334 |
} |
| 335 |
else // center_to > 0 |
| 336 |
{ |
| 337 |
if (irec == 0 || center_to == 2) |
| 338 |
{ |
| 339 |
usedx = crpix1; |
| 340 |
usedy = crpix2; |
| 341 |
usemag = cdelt1; |
| 342 |
} |
| 343 |
mag = cdelt1/usemag; |
| 344 |
cdelt1 = usemag; |
| 345 |
} |
| 346 |
dx = crpix1 - usedx; |
| 347 |
dy = crpix2 - usedy; |
| 348 |
status = image_magrotate( (float *)inparr->data, n, m, dtyp, crota2, |
| 349 |
mag, dx, dy, &output_array, &nx, &ny, regridtype, do_stretchmarks); |
| 350 |
if (verbose) |
| 351 |
{ |
| 352 |
fprintf(stderr,"image_magrotate called with: " |
| 353 |
"data=%p, n=%d, m=%d, dtyp=%d, crota2=%f, mag=%f, dx=%f, dy=%f, regridtype=%d, do_stretch=%d\n", |
| 354 |
(float *)inparr->data, n, m, dtyp, crota2, mag, dx, dy,regridtype, do_stretchmarks); |
| 355 |
fprintf(stderr,"image_magrotate returned: " |
| 356 |
"into out=%p, nx=%d, ny=%d, status=%d\n", output_array, nx, ny, status); |
| 357 |
} |
| 358 |
sprintf(comment,"resize: scale_to=%f, mag=%f, dx=%f, dy=%f, regridtype=%d, do_stretch=%d, center_to=%s", |
| 359 |
scale_to, mag, dx, dy,regridtype, do_stretchmarks, |
| 360 |
(center_to==0 ? "solar disk" : (center_to==1 ? "First frame" : "No change"))); |
| 361 |
drms_appendcomment(outrec, comment, 0); |
| 362 |
if (status) DIE("image_magrotate failed!"); |
| 363 |
outdims[0] = nx; |
| 364 |
outdims[1] = ny; |
| 365 |
outarr = drms_array_create(inparr->type, 2, outdims, output_array, &status); |
| 366 |
if (status) DIE("drms_array_create failed!"); |
| 367 |
|
| 368 |
drms_setkey_float(outrec, "CROTA2", 0.0); |
| 369 |
drms_setkey_float(outrec, "CRPIX1", (nx + 1.0)*0.5); |
| 370 |
drms_setkey_float(outrec, "CRPIX2", (ny + 1.0)*0.5); |
| 371 |
drms_setkey_float(outrec, "CDELT1", cdelt1); |
| 372 |
drms_setkey_float(outrec, "CDELT2", cdelt1); |
| 373 |
set_statistics(outseg, outarr, 1); |
| 374 |
// get info for array from segment |
| 375 |
outarr->bzero = outseg->bzero; |
| 376 |
outarr->bscale = outseg->bscale; |
| 377 |
if (requestid == NOT_SPECIFIED) |
| 378 |
drms_segment_write(outseg, outarr, 0); |
| 379 |
else |
| 380 |
drms_segment_writewithkeys(outseg, outarr, 0); |
| 381 |
if (inparr) drms_free_array(inparr); |
| 382 |
if (outarr) drms_free_array(outarr); |
| 383 |
} |
| 384 |
|
| 385 |
/* If the user did not provide a segment specifier (which means that all segment structs are in the rec->segments container, |
| 386 |
* and there was no seglist), then process only the first segment. */ |
| 387 |
if (drms_record_numsegments(inprec) == drms_record_numsegments(template) && !hasSeglist) |
| 388 |
{ |
| 389 |
break; |
| 390 |
} |
| 391 |
} /* end segment loop */ |
| 392 |
|
| 393 |
if (segIter) |
| 394 |
{ |
| 395 |
hiter_destroy(&segIter); |
| 396 |
} |
| 397 |
|
| 398 |
drms_close_record(outrec, DRMS_INSERT_RECORD); |
| 399 |
} |
| 400 |
} /* end record loop */ |
| 401 |
return 0; |
| 402 |
} |
| 403 |
|
| 404 |
|
| 405 |
|
| 406 |
// definitions at bottom of file |
| 407 |
|
| 408 |
// ---------------------------------------------------------------------- |
| 409 |
|
| 410 |
/* center whith whole pixel shifts and rotate by 180 if needed */ |
| 411 |
/* Only apply center if it will not result in an image crop. I.e. not ever |
| 412 |
for AIA, and not for HMI or MDI or other if a shift of more than 20 arcsec |
| 413 |
is implied */ |
| 414 |
int upNcenter(DRMS_Array_t *arr, ObsInfo_t *ObsLoc) |
| 415 |
{ |
| 416 |
int nx, ny, ix, iy, i, j, xoff, yoff, max_off; |
| 417 |
double rot, x0, y0, midx, midy; |
| 418 |
float *data; |
| 419 |
float *data2; |
| 420 |
if (!arr || !ObsLoc) |
| 421 |
return(1); |
| 422 |
data = arr->data; |
| 423 |
nx = arr->axis[0]; |
| 424 |
ny = arr->axis[1]; |
| 425 |
x0 = ObsLoc->crpix1 - 1; |
| 426 |
y0 = ObsLoc->crpix2 - 1; |
| 427 |
midx = (nx-1.0)/2.0; |
| 428 |
midy = (ny-1.0)/2.0; |
| 429 |
if ((rot = fabs(ObsLoc->crota2)) > 179 && rot < 181) |
| 430 |
{ |
| 431 |
// rotate image by 180 degrees by a flip flip |
| 432 |
float val; |
| 433 |
int half = ny / 2; |
| 434 |
int odd = ny & 1; |
| 435 |
if (odd) half++; |
| 436 |
for (iy=0; iy<half; iy++) |
| 437 |
{ |
| 438 |
for (ix=0; ix<nx; ix++) |
| 439 |
{ |
| 440 |
i = iy*nx + ix; |
| 441 |
j = (ny - 1 - iy)*nx + (nx - 1 - ix); |
| 442 |
val = data[i]; |
| 443 |
data[i] = data[j]; |
| 444 |
data[j] = val; |
| 445 |
} |
| 446 |
} |
| 447 |
x0 = nx - 1 - x0; |
| 448 |
y0 = ny - 1 - y0; |
| 449 |
rot = ObsLoc->crota2 - 180.0; |
| 450 |
if (rot < -90.0) rot += 360.0; |
| 451 |
ObsLoc->crota2 = rot; |
| 452 |
} |
| 453 |
// Center to nearest pixel - if OK to do so |
| 454 |
xoff = round(x0 - midx); |
| 455 |
yoff = round(y0 - midy); |
| 456 |
max_off = 20.0 / ObsLoc->cdelt1; |
| 457 |
if (arr->parent_segment && |
| 458 |
arr->parent_segment->record && |
| 459 |
arr->parent_segment->record->seriesinfo && |
| 460 |
arr->parent_segment->record->seriesinfo->seriesname && |
| 461 |
strncasecmp(arr->parent_segment->record->seriesinfo->seriesname, "aia", 3) && |
| 462 |
abs(xoff) < max_off && abs(yoff) < max_off) |
| 463 |
{ |
| 464 |
if (abs(xoff) >= 1 || abs(yoff) >= 1) |
| 465 |
{ |
| 466 |
data2 = malloc(4*nx*ny); |
| 467 |
for (iy=0; iy<ny; iy++) |
| 468 |
{ |
| 469 |
int jy = iy + yoff; |
| 470 |
for (ix=0; ix<nx; ix++) |
| 471 |
{ |
| 472 |
int jx = ix + xoff; |
| 473 |
int idx = jy*nx + jx; |
| 474 |
int idx2 = iy*nx + ix; |
| 475 |
if (jx<0 || jx>=nx || jy<0 || jy>=ny) |
| 476 |
data2[idx2] = DRMS_MISSING_FLOAT; |
| 477 |
else |
| 478 |
data2[idx2] = data[idx]; |
| 479 |
} |
| 480 |
} |
| 481 |
x0 -= xoff; |
| 482 |
y0 -= yoff; |
| 483 |
free(data); |
| 484 |
arr->data = data2; |
| 485 |
} |
| 486 |
} |
| 487 |
// update center location |
| 488 |
ObsLoc->crpix1 = x0 + 1; |
| 489 |
ObsLoc->crpix2 = y0 + 1; |
| 490 |
return(0); |
| 491 |
} |
| 492 |
|
| 493 |
// ---------------------------------------------------------------------- |
| 494 |
|
| 495 |
int crop_image(DRMS_Array_t *arr, ObsInfo_t *ObsLoc) |
| 496 |
{ |
| 497 |
int nx, ny, ix, iy, i, j, xoff, yoff; |
| 498 |
double x0, y0; |
| 499 |
double rsun = ObsLoc->rsun_obs/ObsLoc->cdelt1; |
| 500 |
double scale, crop_limit2; |
| 501 |
float *data; |
| 502 |
if (!arr || !ObsLoc) |
| 503 |
return(1); |
| 504 |
data = arr->data; |
| 505 |
nx = arr->axis[0]; |
| 506 |
ny = arr->axis[1]; |
| 507 |
x0 = ObsLoc->crpix1 - 1; |
| 508 |
y0 = ObsLoc->crpix2 - 1; |
| 509 |
scale = 1.0/rsun; |
| 510 |
// crop_limit = 0.99975; // 1 - 1/4000, 1/2 HMI pixel. |
| 511 |
crop_limit2 = 0.99950; // square of 1 - 1/4000, 1/2 HMI pixel. |
| 512 |
for (iy=0; iy<ny; iy++) |
| 513 |
for (ix=0; ix<nx; ix++) |
| 514 |
{ |
| 515 |
double x, y, R2; |
| 516 |
float *Ip = data + iy*nx + ix; |
| 517 |
if (drms_ismissing_float(*Ip)) |
| 518 |
continue; |
| 519 |
x = ((double)ix - x0) * scale; /* x,y in pixel coords */ |
| 520 |
y = ((double)iy - y0) * scale; |
| 521 |
R2 = x*x + y*y; |
| 522 |
if (R2 > crop_limit2) |
| 523 |
*Ip = DRMS_MISSING_FLOAT; |
| 524 |
} |
| 525 |
return(0); |
| 526 |
} |
| 527 |
|
| 528 |
// ---------------------------------------------------------------------- |
| 529 |
|
| 530 |
#define CHECK(keyname) {if (status) {fprintf(stderr,"Keyword failure to find: %s, status=%d\n",keyname,status); *rstatus=status; return(NULL);}} |
| 531 |
|
| 532 |
ObsInfo_t *GetObsInfo(DRMS_Segment_t *seg, ObsInfo_t *pObsLoc, int *rstatus) |
| 533 |
{ |
| 534 |
TIME t_prev; |
| 535 |
DRMS_Record_t *rec; |
| 536 |
TIME t_obs; |
| 537 |
double dv; |
| 538 |
ObsInfo_t *ObsLoc; |
| 539 |
int status; |
| 540 |
|
| 541 |
if (!seg || !(rec = seg->record)) |
| 542 |
{ *rstatus = 1; return(NULL); } |
| 543 |
|
| 544 |
ObsLoc = (pObsLoc ? pObsLoc : (ObsInfo_t *)malloc(sizeof(ObsInfo_t))); |
| 545 |
if (!pObsLoc) |
| 546 |
memset(ObsLoc, 0, sizeof(ObsInfo_t)); |
| 547 |
|
| 548 |
t_prev = ObsLoc->t_obs; |
| 549 |
t_obs = drms_getkey_time(rec, "T_OBS", &status); CHECK("T_OBS"); |
| 550 |
|
| 551 |
if (t_obs <= 0.0) |
| 552 |
{ *rstatus = 2; return(NULL); } |
| 553 |
|
| 554 |
if (t_obs != t_prev) |
| 555 |
{ |
| 556 |
ObsLoc->crpix1 = drms_getkey_double(rec, "CRPIX1", &status); CHECK("CRPIX1"); |
| 557 |
ObsLoc->crpix2 = drms_getkey_double(rec, "CRPIX2", &status); CHECK("CRPIX2"); |
| 558 |
ObsLoc->crval1 = drms_getkey_double(rec, "CRVAL1", &status); CHECK("CRVAL1"); |
| 559 |
ObsLoc->crval2 = drms_getkey_double(rec, "CRVAL2", &status); CHECK("CRVAL2"); |
| 560 |
ObsLoc->cdelt1 = drms_getkey_double(rec, "CDELT1", &status); CHECK("CDELT1"); |
| 561 |
ObsLoc->cdelt2 = drms_getkey_double(rec, "CDELT2", &status); CHECK("CDELT1"); |
| 562 |
ObsLoc->crota2 = drms_getkey_double(rec, "CROTA2", &status); if (status) ObsLoc->crota2 = 0.0; // WCS default |
| 563 |
ObsLoc->sina = sin(ObsLoc->crota2*Deg2Rad); |
| 564 |
ObsLoc->cosa = sqrt (1.0 - ObsLoc->sina*ObsLoc->sina); |
| 565 |
ObsLoc->rsun_obs = drms_getkey_double(rec, "RSUN_OBS", &status); |
| 566 |
if (status) |
| 567 |
{ |
| 568 |
double dsun_obs = drms_getkey_double(rec, "DSUN_OBS", &status); CHECK("DSUN_OBS"); |
| 569 |
ObsLoc->rsun_obs = asin(696000000.0/dsun_obs)/arcsec2Rad; |
| 570 |
} |
| 571 |
ObsLoc->obs_vr = drms_getkey_double(rec, "OBS_VR", &status); CHECK("OBS_VR"); |
| 572 |
ObsLoc->obs_vw = drms_getkey_double(rec, "OBS_VW", &status); CHECK("OBS_VW"); |
| 573 |
ObsLoc->obs_vn = drms_getkey_double(rec, "OBS_VN", &status); CHECK("OBS_VN"); |
| 574 |
ObsLoc->obs_b0 = drms_getkey_double(rec, "CRLT_OBS", &status); CHECK("CRLT_OBS"); |
| 575 |
ObsLoc->t_obs = t_obs; |
| 576 |
} |
| 577 |
*rstatus = 0; |
| 578 |
return(ObsLoc); |
| 579 |
} |
| 580 |
|
| 581 |
// ---------------------------------------------------------------------- |
| 582 |
|
| 583 |
// In cases known to not have compact slotted series and cadence is specified |
| 584 |
// generate explicit recordset list of closest good record to desired grid |
| 585 |
// First get vector of times and quality |
| 586 |
// Then if vector is not OK, quit. |
| 587 |
// then: make temp file to hold recordset list |
| 588 |
// start with first time to define desired grid, |
| 589 |
// make array of desired times. |
| 590 |
// make empty array of recnums |
| 591 |
// search vector for good images nearest desired times |
| 592 |
// for each found time, write record query |
| 593 |
|
| 594 |
|
| 595 |
#define DIE_get_recset(msg) {fprintf(stderr,"$$$$ %s: %s\n", module_name, msg); return NULL;} |
| 596 |
const char *get_input_recset(DRMS_Env_t *drms_env, const char *inQuery, DRMS_Record_t **template) |
| 597 |
{ |
| 598 |
static char newInQuery[DRMS_MAXSERIESNAMELEN+2]; |
| 599 |
int epoch_given = cmdparams_exists(&cmdparams, "epoch"); |
| 600 |
TIME epoch, t_epoch; |
| 601 |
DRMS_Array_t *data; |
| 602 |
DRMS_Record_t *inTemplate; |
| 603 |
TIME t_start, t_stop, t_now, t_want, t_diff, this_t_diff; |
| 604 |
int status = 1; |
| 605 |
int nrecs, irec; |
| 606 |
int nslots, islot; |
| 607 |
long long *recnums; |
| 608 |
TIME *t_this, half; |
| 609 |
TIME cadence; |
| 610 |
double *drecnum, *dquality; |
| 611 |
int quality; |
| 612 |
long long recnum; |
| 613 |
char keylist[DRMS_MAXQUERYLEN]; |
| 614 |
char filename[DRMS_MAXSERIESNAMELEN]; |
| 615 |
char *tmpdir; |
| 616 |
FILE *tmpfile; |
| 617 |
char newIn[DRMS_MAXQUERYLEN]; |
| 618 |
char seriesname[DRMS_MAXQUERYLEN]; |
| 619 |
char *lbracket; |
| 620 |
char *at = index(inQuery, '@'); |
| 621 |
int npkeys; |
| 622 |
char *timekeyname; |
| 623 |
double t_step; |
| 624 |
|
| 625 |
strcpy(seriesname, inQuery); |
| 626 |
lbracket = index(seriesname,'['); |
| 627 |
if (lbracket) *lbracket = '\0'; |
| 628 |
inTemplate = drms_template_record(drms_env, seriesname, &status); |
| 629 |
if (template) |
| 630 |
{ |
| 631 |
*template = inTemplate; |
| 632 |
} |
| 633 |
if (status || !inTemplate) DIE_get_recset("Input series can not be found"); |
| 634 |
|
| 635 |
// Now find the prime time keyword name |
| 636 |
npkeys = inTemplate->seriesinfo->pidx_num; |
| 637 |
timekeyname = NULL; |
| 638 |
if (npkeys > 0) |
| 639 |
{ |
| 640 |
int i; |
| 641 |
for (i=0; i<npkeys; i++) |
| 642 |
{ |
| 643 |
DRMS_Keyword_t *pkey, *skey; |
| 644 |
pkey = inTemplate->seriesinfo->pidx_keywords[i]; |
| 645 |
if (pkey->info->recscope > 1) |
| 646 |
pkey = drms_keyword_slotfromindex(pkey); |
| 647 |
if (pkey->info->type != DRMS_TYPE_TIME) |
| 648 |
continue; |
| 649 |
if(i > 0) DIE_get_recset("Input series must have TIME keyword first, for now..."); |
| 650 |
timekeyname = pkey->info->name; |
| 651 |
t_step = drms_keyword_getdouble(drms_keyword_stepfromslot(pkey), &status); |
| 652 |
if (status) DIE_get_recset("problem getting t_step"); |
| 653 |
t_epoch = drms_keyword_getdouble(drms_keyword_epochfromslot(pkey), &status); |
| 654 |
if (status) DIE_get_recset("problem getting t_epoch"); |
| 655 |
} |
| 656 |
} |
| 657 |
else |
| 658 |
DIE_get_recset("Must have time prime key"); |
| 659 |
epoch = epoch_given ? params_get_time(&cmdparams, "epoch") : t_epoch; |
| 660 |
|
| 661 |
if (at && *at && ((strncmp(inQuery,"aia.lev1[", 9)==0 || |
| 662 |
strncmp(inQuery,"hmi.lev1[", 9)==0 || |
| 663 |
strncmp(inQuery,"aia.lev1_nrt2[",14)==0 || |
| 664 |
strncmp(inQuery,"hmi.lev1_nrt[", 13)==0 ) || |
| 665 |
epoch_given)) |
| 666 |
{ |
| 667 |
char *ip=(char *)inQuery, *op=newIn, *p; |
| 668 |
long n, mul; |
| 669 |
while ( *ip && ip<at ) |
| 670 |
*op++ = *ip++; |
| 671 |
ip++; // skip the '@' |
| 672 |
n = strtol(ip, &p, 10); // get digits only |
| 673 |
if (*p == 's') mul = 1; |
| 674 |
else if (*p == 'm') mul = 60; |
| 675 |
else if (*p == 'h') mul = 3600; |
| 676 |
else if (*p == 'd') mul = 86400; |
| 677 |
else |
| 678 |
DIE_get_recset("cant make sense of @xx cadence spec"); |
| 679 |
cadence = n * mul; |
| 680 |
ip = ++p; // skip cadence multiplier |
| 681 |
while ( *ip ) |
| 682 |
*op++ = *ip++; |
| 683 |
*op = '\0'; |
| 684 |
half = cadence/2.0; |
| 685 |
sprintf(keylist, "%s,QUALITY,recnum", timekeyname); |
| 686 |
data = drms_record_getvector(drms_env, newIn, keylist, DRMS_TYPE_DOUBLE, 0, &status); |
| 687 |
if (!data || status) |
| 688 |
{ |
| 689 |
fprintf(stderr, "status=%d\n", status); |
| 690 |
DIE_get_recset("getkey_vector failed\n"); |
| 691 |
} |
| 692 |
nrecs = data->axis[1]; |
| 693 |
irec = 0; |
| 694 |
t_this = (TIME *)data->data; |
| 695 |
dquality = (double *)data->data + 1*nrecs; |
| 696 |
drecnum = (double *)data->data + 2*nrecs; |
| 697 |
if (epoch_given) |
| 698 |
{ |
| 699 |
int s0 = (t_this[0] - epoch)/cadence; |
| 700 |
TIME t0 = s0*cadence + epoch; |
| 701 |
t_start = (t0 < t_this[0] ? t0 + cadence : t0); |
| 702 |
} |
| 703 |
else |
| 704 |
t_start = t_this[0]; |
| 705 |
t_stop = t_this[nrecs-1]; |
| 706 |
nslots = (t_stop - t_start + cadence/2)/cadence; |
| 707 |
recnums = (long long *)malloc(nslots*sizeof(long long)); |
| 708 |
for (islot=0; islot<nslots; islot++) |
| 709 |
recnums[islot] = 0; |
| 710 |
islot = 0; |
| 711 |
t_want = t_start; |
| 712 |
t_diff = 1.0e9; |
| 713 |
for (irec = 0; irec<nrecs; irec++) |
| 714 |
{ |
| 715 |
t_now = t_this[irec]; |
| 716 |
quality = (int)dquality[irec] & 0xFFFFFFFF; |
| 717 |
recnum = (long long)drecnum[irec]; |
| 718 |
this_t_diff = fabs(t_now - t_want); |
| 719 |
if (quality < 0) |
| 720 |
continue; |
| 721 |
if (t_now <= (t_want-half)) |
| 722 |
continue; |
| 723 |
while (t_now > (t_want+half)) |
| 724 |
{ |
| 725 |
islot++; |
| 726 |
if (islot >= nslots) |
| 727 |
break; |
| 728 |
t_want = t_start + cadence * islot; |
| 729 |
this_t_diff = fabs(t_now - t_want); |
| 730 |
t_diff = 1.0e8; |
| 731 |
} |
| 732 |
if (this_t_diff <= t_diff) |
| 733 |
recnums[islot] = recnum; |
| 734 |
t_diff = fabs(t_now - t_want); |
| 735 |
} |
| 736 |
if (islot+1 < nslots) |
| 737 |
nslots = islot+1; // take what we got. |
| 738 |
tmpdir = getenv("TMPDIR"); |
| 739 |
if (!tmpdir) tmpdir = "/tmp"; |
| 740 |
sprintf(filename, "%s/%sXXXXXX", tmpdir, module_name); |
| 741 |
mkstemp(filename); |
| 742 |
tmpfile = fopen(filename,"w"); |
| 743 |
for (islot=0; islot<nslots; islot++) |
| 744 |
if (recnums[islot]) |
| 745 |
fprintf(tmpfile, "%s[:#%lld]\n", seriesname, recnums[islot]); |
| 746 |
fclose(tmpfile); |
| 747 |
free(recnums); |
| 748 |
drms_free_array(data); |
| 749 |
sprintf(newInQuery,"@%s", filename); |
| 750 |
return(newInQuery); |
| 751 |
} |
| 752 |
else |
| 753 |
return(inQuery); |
| 754 |
} |
