From 1dc0d2ea551cc32e00c03dc916d4a3e0bd2de13d Mon Sep 17 00:00:00 2001 From: Robert Butora <robert.butora@inaf.it> Date: Fri, 9 May 2025 17:08:58 +0200 Subject: [PATCH] cutout: modifies SodaImpl to use vlkb-cutout instead on two calls to vlkb-overlap + vlkb-imcopy --- data-access/engine/src/vlkb/src/imcopy.cpp | 878 +++++++++--------- data-access/engine/src/vlkb/src/main.cpp | 18 +- .../servlet/src/main/java/cutout/Soda.java | 2 +- .../src/main/java/cutout/SodaImpl.java | 276 ++---- .../main/java/cutout/engine/cli/ExecCmd.java | 16 +- .../java/cutout/webapi/ServletCutout.java | 7 +- .../src/main/java/mcutout/VlkbCli.java | 6 +- 7 files changed, 523 insertions(+), 680 deletions(-) diff --git a/data-access/engine/src/vlkb/src/imcopy.cpp b/data-access/engine/src/vlkb/src/imcopy.cpp index ac45274..c4df34c 100644 --- a/data-access/engine/src/vlkb/src/imcopy.cpp +++ b/data-access/engine/src/vlkb/src/imcopy.cpp @@ -23,23 +23,22 @@ using namespace std; int stream_cutout(string pathname, int extnum, string region) { - vector<uint_bounds> bnds; + vector<uint_bounds> bnds; - int ov_code = vlkb_overlap(pathname, region, bnds); + int ov_code = vlkb_overlap(pathname, region, bnds); if((ov_code >= 2) && (ov_code <= 5)) - { - string pixfilter{ to_cfitsio_format(bnds) }; - - imcopy(pathname, extnum, pixfilter); - return EXIT_SUCCESS; - - } - else if((ov_code == 1) || (ov_code == 6)) - { - // no overlap - return EXIT_SUCCESS; - } + { + string pixfilter{ to_cfitsio_format(bnds) }; + + imcopy(pathname, extnum, pixfilter); + // pixfilter must be within NAXIS[] otherwise cfitsio-error + return 0; // EXIT_SUCCESS + } + else if((ov_code == 1) || (ov_code == 6)) + { + return 1; // no overlap EXIT_SUCCESS + } else { throw runtime_error("overlap code invalid: " + to_string(ov_code)); @@ -47,107 +46,106 @@ int stream_cutout(string pathname, int extnum, string region) } - // ---------------------------------------------------------------------------------- int fits_copy_image_section2( - fitsfile *fptr, /* I - pointer to input image */ - fitsfile *newptr, /* I - pointer to output image */ - char *expr, /* I - Image section expression */ - int *status); + fitsfile *fptr, /* I - pointer to input image */ + fitsfile *newptr, /* I - pointer to output image */ + char *expr, /* I - Image section expression */ + int *status); void imcopy(std::string filename, int extnum, std::string pixfilter) { - LOG_trace(__func__); - - LOG_STREAM << filename << " EXT: " << extnum << endl; - LOG_STREAM << "pixfilter: " << pixfilter << endl; - - int status = 0; - - pixfilter.erase( remove(pixfilter.begin(), pixfilter.end(), '['), pixfilter.end() ); - pixfilter.erase( remove(pixfilter.begin(), pixfilter.end(), ']'), pixfilter.end() ); - - LOG_STREAM << "filter expr: " << pixfilter << endl; - - char *expr = (char*)pixfilter.c_str(); /* I - Image section expression */ - - fitsfile *fptr; /* I - pointer to input image */ - fitsfile *newfptr; /* I - pointer to output image */ - - const char * cfilename = filename.c_str(); - - fits_open_diskfile(&fptr, cfilename, READONLY, &status); - if (status) - { - string errmsg{ fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, status) }; - throw runtime_error("fits_open_file " + filename + " failed: " + errmsg); - } - - /* cfitsio(v4.2.0) manual: 10.2.3 Notes about the stream filetype driver - * "-" writes to stdout: but first loads all file into mem - * "stream://" keeps NIOBUF=40 buffers, and header must fit into it */ - fits_create_file(&newfptr, "stream://", &status); - if (status) - { - string errmsg{ fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, status) }; - - int tstatus = 0; - fits_close_file(fptr, &tstatus); - if(tstatus) LOG_STREAM << fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, tstatus) << tstatus<< endl; - - throw runtime_error("fits_open_file to stream failed: " + errmsg); - } - - /* err codes 1..1000 reserved by cfitsio */ - #define HEADER_TOO_LONG (10000) - #define MAX_CARD_COUNT (1400) - /* see notes about stream driver in cfitsio manual */ - fits_copy_image_section2(fptr, newfptr, expr, &status); - if (status) - { - int tstatus = 0; - fits_close_file(fptr, &tstatus); - if(tstatus) LOG_STREAM << fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, tstatus) << tstatus<< endl; - tstatus = 0; - fits_close_file(newfptr, &tstatus); - if(tstatus) LOG_STREAM << fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, tstatus) << tstatus<< endl; - - if(status == HEADER_TOO_LONG) - { - throw runtime_error("Cutout from " + filename - + " failed: header is too long. Direct streaming in current build support max " - + to_string(MAX_CARD_COUNT) + " cards in header."); - } - else - { - string errmsg{fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, status)}; - throw runtime_error("fits_copy_image_section " + filename - + " to cut-file with " + string{expr} + " failed: " + errmsg); - } - } - - fits_close_file(fptr, &status); - if (status) - { - string errmsg{fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, status)}; - - int tstatus = 0; - fits_close_file(newfptr, &tstatus); - if(tstatus) LOG_STREAM << fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, tstatus) << tstatus<< endl; - - throw runtime_error("fits_close_file " + filename + " failed: " + errmsg); - } - - fits_close_file(newfptr, &status); - if (status) - { - string errmsg{fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, status)}; - throw runtime_error("fits_close_file cut failed: " + errmsg); - } + LOG_trace(__func__); + + LOG_STREAM << filename << " EXT: " << extnum << endl; + LOG_STREAM << "pixfilter: " << pixfilter << endl; + + int status = 0; + + pixfilter.erase( remove(pixfilter.begin(), pixfilter.end(), '['), pixfilter.end() ); + pixfilter.erase( remove(pixfilter.begin(), pixfilter.end(), ']'), pixfilter.end() ); + + LOG_STREAM << "filter expr: " << pixfilter << endl; + + char *expr = (char*)pixfilter.c_str(); /* I - Image section expression */ + + fitsfile *fptr; /* I - pointer to input image */ + fitsfile *newfptr; /* I - pointer to output image */ + + const char * cfilename = filename.c_str(); + + fits_open_diskfile(&fptr, cfilename, READONLY, &status); + if (status) + { + string errmsg{ fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, status) }; + throw runtime_error("fits_open_file " + filename + " failed: " + errmsg); + } + + /* cfitsio(v4.2.0) manual: 10.2.3 Notes about the stream filetype driver + * "-" writes to stdout: but first loads all file into mem + * "stream://" keeps NIOBUF=40 buffers, and header must fit into it */ + fits_create_file(&newfptr, "stream://", &status); + if (status) + { + string errmsg{ fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, status) }; + + int tstatus = 0; + fits_close_file(fptr, &tstatus); + if(tstatus) LOG_STREAM << fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, tstatus) << tstatus<< endl; + + throw runtime_error("fits_open_file to stream failed: " + errmsg); + } + + /* err codes 1..1000 reserved by cfitsio */ +#define HEADER_TOO_LONG (10000) +#define MAX_CARD_COUNT (1400) + /* see notes about stream driver in cfitsio manual */ + fits_copy_image_section2(fptr, newfptr, expr, &status); + if (status) + { + int tstatus = 0; + fits_close_file(fptr, &tstatus); + if(tstatus) LOG_STREAM << fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, tstatus) << tstatus<< endl; + tstatus = 0; + fits_close_file(newfptr, &tstatus); + if(tstatus) LOG_STREAM << fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, tstatus) << tstatus<< endl; + + if(status == HEADER_TOO_LONG) + { + throw runtime_error("Cutout from " + filename + + " failed: header is too long. Direct streaming in current build support max " + + to_string(MAX_CARD_COUNT) + " cards in header."); + } + else + { + string errmsg{fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, status)}; + throw runtime_error("fits_copy_image_section " + filename + + " to cut-file with " + string{expr} + " failed: " + errmsg); + } + } + + fits_close_file(fptr, &status); + if (status) + { + string errmsg{fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, status)}; + + int tstatus = 0; + fits_close_file(newfptr, &tstatus); + if(tstatus) LOG_STREAM << fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, tstatus) << tstatus<< endl; + + throw runtime_error("fits_close_file " + filename + " failed: " + errmsg); + } + + fits_close_file(newfptr, &status); + if (status) + { + string errmsg{fitsfiles::cfitsio_errmsg(__FILE__, __LINE__, status)}; + throw runtime_error("fits_close_file cut failed: " + errmsg); + } } @@ -157,340 +155,340 @@ void imcopy(std::string filename, int extnum, std::string pixfilter) * Also consider to duplicate the buffer and run reader * on parallel-thread to avoid reader-write wait for each other. */ int fits_copy_image_section2( - fitsfile *fptr, /* I - pointer to input image */ - fitsfile *newptr, /* I - pointer to output image */ - char *expr, /* I - Image section expression */ - int *status) + fitsfile *fptr, /* I - pointer to input image */ + fitsfile *newptr, /* I - pointer to output image */ + char *expr, /* I - Image section expression */ + int *status) { - /* copies an image section from the input file to a new output HDU */ - - int bitpix, naxis, numkeys, nkey; - long naxes[] = {1,1,1,1,1,1,1,1,1}, smin, smax, sinc; - long fpixels[] = {1,1,1,1,1,1,1,1,1}; - long lpixels[] = {1,1,1,1,1,1,1,1,1}; - long incs[] = {1,1,1,1,1,1,1,1,1}; - char *cptr, keyname[FLEN_KEYWORD], card[FLEN_CARD]; - int ii, tstatus, anynull; - long minrow, maxrow, minslice, maxslice, mincube, maxcube; - long firstpix; - long ncubeiter, nsliceiter, nrowiter, kiter, jiter, iiter; - int klen, kk, jj; - long outnaxes[9], outsize, buffsize; - double *buffer, crpix, cdelt; - - if (*status > 0) - return(*status); - - /* get the size of the input image */ - fits_get_img_type(fptr, &bitpix, status); - fits_get_img_dim(fptr, &naxis, status); - if (fits_get_img_size(fptr, naxis, naxes, status) > 0) - return(*status); - - if (naxis < 1 || naxis > 4) - { - ffpmsg( - "Input image either had NAXIS = 0 (NULL image) or has > 4 dimensions"); - return(*status = BAD_NAXIS); - } - - /* create output image with same size and type as the input image */ - /* Will update the size later */ - fits_create_img(newptr, bitpix, naxis, naxes, status); - - /* copy all other non-structural keywords from the input to output file */ - fits_get_hdrspace(fptr, &numkeys, NULL, status); - - int cards_written = 4; - for (nkey = 4; nkey <= numkeys; nkey++) /* skip the first few keywords */ - { - fits_read_record(fptr, nkey, card, status); - - /* cfitsio(v4.2.0) manual: 10.2.3 Notes about the stream filetype driver - * stream:// driver has internal buffer NIOBUF=40 blocks e.g. 1400 cards - * TYP_COMM_KEY: skip HISTORY COMMENT and empty keys */ - if(numkeys > MAX_CARD_COUNT) - if (fits_get_keyclass(card) == TYP_COMM_KEY) continue; - - if (fits_get_keyclass(card) > TYP_CMPRS_KEY) - { - if(++cards_written > MAX_CARD_COUNT) return(HEADER_TOO_LONG); - - /* write the record to the output file */ - fits_write_record(newptr, card, status); - } - } - - if (*status > 0) - { - ffpmsg("error copying header from input image to output image"); - return(*status); - } - - /* parse the section specifier to get min, max, and inc for each axis */ - /* and the size of each output image axis */ - - cptr = expr; - for (ii=0; ii < naxis; ii++) - { - if (fits_get_section_range(&cptr, &smin, &smax, &sinc, status) > 0) - { - ffpmsg("error parsing the following image section specifier:"); - ffpmsg(expr); - return(*status); - } - - if (smax == 0) - smax = naxes[ii]; /* use whole axis by default */ - else if (smin == 0) - smin = naxes[ii]; /* use inverted whole axis */ - - if (smin > naxes[ii] || smax > naxes[ii]) - { - ffpmsg("image section exceeds dimensions of input image:"); - ffpmsg(expr); - return(*status = BAD_NAXIS); - } - - fpixels[ii] = smin; - lpixels[ii] = smax; - incs[ii] = sinc; - - if (smin <= smax) - outnaxes[ii] = (smax - smin + sinc) / sinc; - else - outnaxes[ii] = (smin - smax + sinc) / sinc; - - /* modify the NAXISn keyword */ - fits_make_keyn("NAXIS", ii + 1, keyname, status); - fits_modify_key_lng(newptr, keyname, outnaxes[ii], NULL, status); - - /* modify the WCS keywords if necessary */ - - if (fpixels[ii] != 1 || incs[ii] != 1) - { - for (kk=-1;kk<26; kk++) /* modify any alternate WCS keywords */ - { - /* read the CRPIXn keyword if it exists in the input file */ - fits_make_keyn("CRPIX", ii + 1, keyname, status); - - if (kk != -1) { - klen = (int)strlen(keyname); - keyname[klen]= (char)((int)'A' + kk); - keyname[klen + 1] = '\0'; - } - - tstatus = 0; - if (fits_read_key(fptr, TDOUBLE, keyname, - &crpix, NULL, &tstatus) == 0) - { - /* calculate the new CRPIXn value */ - if (fpixels[ii] <= lpixels[ii]) { - crpix = (crpix - double(fpixels[ii])) / double(incs[ii]) + 1.0; - /* crpix = (crpix - (fpixels[ii] - 1.0) - .5) / incs[ii] + 0.5; */ - } else { - crpix = (double(fpixels[ii]) - crpix) / double(incs[ii]) + 1.0; - /* crpix = (fpixels[ii] - (crpix - 1.0) - .5) / incs[ii] + 0.5; */ - } - - /* modify the value in the output file */ - fits_modify_key_dbl(newptr, keyname, crpix, 15, NULL, status); - - if (incs[ii] != 1 || fpixels[ii] > lpixels[ii]) - { - /* read the CDELTn keyword if it exists in the input file */ - fits_make_keyn("CDELT", ii + 1, keyname, status); - - if (kk != -1) { - klen = (int)strlen(keyname); - keyname[klen]=(char)((int)'A' + kk); - keyname[klen + 1] = '\0'; - } - - tstatus = 0; - if (fits_read_key(fptr, TDOUBLE, keyname, - &cdelt, NULL, &tstatus) == 0) - { - /* calculate the new CDELTn value */ - if (fpixels[ii] <= lpixels[ii]) - cdelt = cdelt * double(incs[ii]); - else - cdelt = cdelt * double(-incs[ii]); - - /* modify the value in the output file */ - fits_modify_key_dbl(newptr, keyname, cdelt, 15, NULL, status); - } - - /* modify the CDi_j keywords if they exist in the input file */ - - fits_make_keyn("CD1_", ii + 1, keyname, status); - - if (kk != -1) { - klen = (int)strlen(keyname); - keyname[klen]=(char)((int)'A' + kk); - keyname[klen + 1] = '\0'; - } - - for (jj=0; jj < 9; jj++) /* look for up to 9 dimensions */ - { - keyname[2] = (char)((int)'1' + jj); - - tstatus = 0; - if (fits_read_key(fptr, TDOUBLE, keyname, - &cdelt, NULL, &tstatus) == 0) - { - /* calculate the new CDi_j value */ - if (fpixels[ii] <= lpixels[ii]) - cdelt = cdelt * double(incs[ii]); - else - cdelt = cdelt * double(-incs[ii]); - - /* modify the value in the output file */ - fits_modify_key_dbl(newptr, keyname, cdelt, 15, NULL, status); - } - } - - } /* end of if (incs[ii]... loop */ - } /* end of fits_read_key loop */ - } /* end of for (kk loop */ - } - } /* end of main NAXIS loop */ - - if (ffrdef(newptr, status) > 0) /* force the header to be scanned */ - { - return(*status); - } - - /* turn off any scaling of the pixel values */ - fits_set_bscale(fptr, 1.0, 0.0, status); - fits_set_bscale(newptr, 1.0, 0.0, status); - - /* to reduce memory foot print, just read/write image 1 row at a time */ - - outsize = outnaxes[0]; - buffsize = (abs(bitpix) / 8) * outsize; - - buffer = (double *) malloc(buffsize); /* allocate memory for the image row */ - if (!buffer) - { - ffpmsg("fits_copy_image_section: no memory for image section"); - return(*status = MEMORY_ALLOCATION); - } - /* read the image section then write it to the output file */ - - minrow = fpixels[1]; - maxrow = lpixels[1]; - if (minrow > maxrow) { - nrowiter = (minrow - maxrow + incs[1]) / incs[1]; - } else { - nrowiter = (maxrow - minrow + incs[1]) / incs[1]; - } - - minslice = fpixels[2]; - maxslice = lpixels[2]; - if (minslice > maxslice) { - nsliceiter = (minslice - maxslice + incs[2]) / incs[2]; - } else { - nsliceiter = (maxslice - minslice + incs[2]) / incs[2]; - } - - mincube = fpixels[3]; - maxcube = lpixels[3]; - if (mincube > maxcube) { - ncubeiter = (mincube - maxcube + incs[3]) / incs[3]; - } else { - ncubeiter = (maxcube - mincube + incs[3]) / incs[3]; - } - - firstpix = 1; - for (kiter = 0; kiter < ncubeiter; kiter++) - { - if (mincube > maxcube) { - fpixels[3] = mincube - (kiter * incs[3]); - } else { - fpixels[3] = mincube + (kiter * incs[3]); - } - - lpixels[3] = fpixels[3]; - - for (jiter = 0; jiter < nsliceiter; jiter++) - { - if (minslice > maxslice) { - fpixels[2] = minslice - (jiter * incs[2]); - } else { - fpixels[2] = minslice + (jiter * incs[2]); - } - - lpixels[2] = fpixels[2]; - - for (iiter = 0; iiter < nrowiter; iiter++) - { - if (minrow > maxrow) { - fpixels[1] = minrow - (iiter * incs[1]); - } else { - fpixels[1] = minrow + (iiter * incs[1]); - } - - lpixels[1] = fpixels[1]; - - if (bitpix == 8) - { - ffgsvb(fptr, 1, naxis, naxes, fpixels, lpixels, incs, 0, - (unsigned char *) buffer, &anynull, status); - - ffpprb(newptr, 1, firstpix, outsize, (unsigned char *) buffer, status); - } - else if (bitpix == 16) - { - ffgsvi(fptr, 1, naxis, naxes, fpixels, lpixels, incs, 0, - (short *) buffer, &anynull, status); - - ffppri(newptr, 1, firstpix, outsize, (short *) buffer, status); - } - else if (bitpix == 32) - { - ffgsvk(fptr, 1, naxis, naxes, fpixels, lpixels, incs, 0, - (int *) buffer, &anynull, status); - - ffpprk(newptr, 1, firstpix, outsize, (int *) buffer, status); - } - else if (bitpix == -32) - { - ffgsve(fptr, 1, naxis, naxes, fpixels, lpixels, incs, FLOATNULLVALUE, - (float *) buffer, &anynull, status); - - ffppne(newptr, 1, firstpix, outsize, (float *) buffer, FLOATNULLVALUE, status); - } - else if (bitpix == -64) - { - ffgsvd(fptr, 1, naxis, naxes, fpixels, lpixels, incs, DOUBLENULLVALUE, - buffer, &anynull, status); - - ffppnd(newptr, 1, firstpix, outsize, buffer, DOUBLENULLVALUE, - status); - } - else if (bitpix == 64) - { - ffgsvjj(fptr, 1, naxis, naxes, fpixels, lpixels, incs, 0, - (LONGLONG *) buffer, &anynull, status); - - ffpprjj(newptr, 1, firstpix, outsize, (LONGLONG *) buffer, status); - } - - - firstpix += outsize; - } - } - } - - free(buffer); /* finished with the memory */ - - if (*status > 0) - { - ffpmsg("fits_copy_image_section: error copying image section"); - return(*status); - } - - return(*status); + /* copies an image section from the input file to a new output HDU */ + + int bitpix, naxis, numkeys, nkey; + long naxes[] = {1,1,1,1,1,1,1,1,1}, smin, smax, sinc; + long fpixels[] = {1,1,1,1,1,1,1,1,1}; + long lpixels[] = {1,1,1,1,1,1,1,1,1}; + long incs[] = {1,1,1,1,1,1,1,1,1}; + char *cptr, keyname[FLEN_KEYWORD], card[FLEN_CARD]; + int ii, tstatus, anynull; + long minrow, maxrow, minslice, maxslice, mincube, maxcube; + long firstpix; + long ncubeiter, nsliceiter, nrowiter, kiter, jiter, iiter; + int klen, kk, jj; + long outnaxes[9], outsize, buffsize; + double *buffer, crpix, cdelt; + + if (*status > 0) + return(*status); + + /* get the size of the input image */ + fits_get_img_type(fptr, &bitpix, status); + fits_get_img_dim(fptr, &naxis, status); + if (fits_get_img_size(fptr, naxis, naxes, status) > 0) + return(*status); + + if (naxis < 1 || naxis > 4) + { + ffpmsg( + "Input image either had NAXIS = 0 (NULL image) or has > 4 dimensions"); + return(*status = BAD_NAXIS); + } + + /* create output image with same size and type as the input image */ + /* Will update the size later */ + fits_create_img(newptr, bitpix, naxis, naxes, status); + + /* copy all other non-structural keywords from the input to output file */ + fits_get_hdrspace(fptr, &numkeys, NULL, status); + + int cards_written = 4; + for (nkey = 4; nkey <= numkeys; nkey++) /* skip the first few keywords */ + { + fits_read_record(fptr, nkey, card, status); + + /* cfitsio(v4.2.0) manual: 10.2.3 Notes about the stream filetype driver + * stream:// driver has internal buffer NIOBUF=40 blocks e.g. 1400 cards + * TYP_COMM_KEY: skip HISTORY COMMENT and empty keys */ + if(numkeys > MAX_CARD_COUNT) + if (fits_get_keyclass(card) == TYP_COMM_KEY) continue; + + if (fits_get_keyclass(card) > TYP_CMPRS_KEY) + { + if(++cards_written > MAX_CARD_COUNT) return(HEADER_TOO_LONG); + + /* write the record to the output file */ + fits_write_record(newptr, card, status); + } + } + + if (*status > 0) + { + ffpmsg("error copying header from input image to output image"); + return(*status); + } + + /* parse the section specifier to get min, max, and inc for each axis */ + /* and the size of each output image axis */ + + cptr = expr; + for (ii=0; ii < naxis; ii++) + { + if (fits_get_section_range(&cptr, &smin, &smax, &sinc, status) > 0) + { + ffpmsg("error parsing the following image section specifier:"); + ffpmsg(expr); + return(*status); + } + + if (smax == 0) + smax = naxes[ii]; /* use whole axis by default */ + else if (smin == 0) + smin = naxes[ii]; /* use inverted whole axis */ + + if (smin > naxes[ii] || smax > naxes[ii]) + { + ffpmsg("image section exceeds dimensions of input image:"); + ffpmsg(expr); + return(*status = BAD_NAXIS); + } + + fpixels[ii] = smin; + lpixels[ii] = smax; + incs[ii] = sinc; + + if (smin <= smax) + outnaxes[ii] = (smax - smin + sinc) / sinc; + else + outnaxes[ii] = (smin - smax + sinc) / sinc; + + /* modify the NAXISn keyword */ + fits_make_keyn("NAXIS", ii + 1, keyname, status); + fits_modify_key_lng(newptr, keyname, outnaxes[ii], NULL, status); + + /* modify the WCS keywords if necessary */ + + if (fpixels[ii] != 1 || incs[ii] != 1) + { + for (kk=-1;kk<26; kk++) /* modify any alternate WCS keywords */ + { + /* read the CRPIXn keyword if it exists in the input file */ + fits_make_keyn("CRPIX", ii + 1, keyname, status); + + if (kk != -1) { + klen = (int)strlen(keyname); + keyname[klen]= (char)((int)'A' + kk); + keyname[klen + 1] = '\0'; + } + + tstatus = 0; + if (fits_read_key(fptr, TDOUBLE, keyname, + &crpix, NULL, &tstatus) == 0) + { + /* calculate the new CRPIXn value */ + if (fpixels[ii] <= lpixels[ii]) { + crpix = (crpix - double(fpixels[ii])) / double(incs[ii]) + 1.0; + /* crpix = (crpix - (fpixels[ii] - 1.0) - .5) / incs[ii] + 0.5; */ + } else { + crpix = (double(fpixels[ii]) - crpix) / double(incs[ii]) + 1.0; + /* crpix = (fpixels[ii] - (crpix - 1.0) - .5) / incs[ii] + 0.5; */ + } + + /* modify the value in the output file */ + fits_modify_key_dbl(newptr, keyname, crpix, 15, NULL, status); + + if (incs[ii] != 1 || fpixels[ii] > lpixels[ii]) + { + /* read the CDELTn keyword if it exists in the input file */ + fits_make_keyn("CDELT", ii + 1, keyname, status); + + if (kk != -1) { + klen = (int)strlen(keyname); + keyname[klen]=(char)((int)'A' + kk); + keyname[klen + 1] = '\0'; + } + + tstatus = 0; + if (fits_read_key(fptr, TDOUBLE, keyname, + &cdelt, NULL, &tstatus) == 0) + { + /* calculate the new CDELTn value */ + if (fpixels[ii] <= lpixels[ii]) + cdelt = cdelt * double(incs[ii]); + else + cdelt = cdelt * double(-incs[ii]); + + /* modify the value in the output file */ + fits_modify_key_dbl(newptr, keyname, cdelt, 15, NULL, status); + } + + /* modify the CDi_j keywords if they exist in the input file */ + + fits_make_keyn("CD1_", ii + 1, keyname, status); + + if (kk != -1) { + klen = (int)strlen(keyname); + keyname[klen]=(char)((int)'A' + kk); + keyname[klen + 1] = '\0'; + } + + for (jj=0; jj < 9; jj++) /* look for up to 9 dimensions */ + { + keyname[2] = (char)((int)'1' + jj); + + tstatus = 0; + if (fits_read_key(fptr, TDOUBLE, keyname, + &cdelt, NULL, &tstatus) == 0) + { + /* calculate the new CDi_j value */ + if (fpixels[ii] <= lpixels[ii]) + cdelt = cdelt * double(incs[ii]); + else + cdelt = cdelt * double(-incs[ii]); + + /* modify the value in the output file */ + fits_modify_key_dbl(newptr, keyname, cdelt, 15, NULL, status); + } + } + + } /* end of if (incs[ii]... loop */ + } /* end of fits_read_key loop */ + } /* end of for (kk loop */ + } + } /* end of main NAXIS loop */ + + if (ffrdef(newptr, status) > 0) /* force the header to be scanned */ + { + return(*status); + } + + /* turn off any scaling of the pixel values */ + fits_set_bscale(fptr, 1.0, 0.0, status); + fits_set_bscale(newptr, 1.0, 0.0, status); + + /* to reduce memory foot print, just read/write image 1 row at a time */ + + outsize = outnaxes[0]; + buffsize = (abs(bitpix) / 8) * outsize; + + buffer = (double *) malloc(buffsize); /* allocate memory for the image row */ + if (!buffer) + { + ffpmsg("fits_copy_image_section: no memory for image section"); + return(*status = MEMORY_ALLOCATION); + } + /* read the image section then write it to the output file */ + + minrow = fpixels[1]; + maxrow = lpixels[1]; + if (minrow > maxrow) { + nrowiter = (minrow - maxrow + incs[1]) / incs[1]; + } else { + nrowiter = (maxrow - minrow + incs[1]) / incs[1]; + } + + minslice = fpixels[2]; + maxslice = lpixels[2]; + if (minslice > maxslice) { + nsliceiter = (minslice - maxslice + incs[2]) / incs[2]; + } else { + nsliceiter = (maxslice - minslice + incs[2]) / incs[2]; + } + + mincube = fpixels[3]; + maxcube = lpixels[3]; + if (mincube > maxcube) { + ncubeiter = (mincube - maxcube + incs[3]) / incs[3]; + } else { + ncubeiter = (maxcube - mincube + incs[3]) / incs[3]; + } + + firstpix = 1; + for (kiter = 0; kiter < ncubeiter; kiter++) + { + if (mincube > maxcube) { + fpixels[3] = mincube - (kiter * incs[3]); + } else { + fpixels[3] = mincube + (kiter * incs[3]); + } + + lpixels[3] = fpixels[3]; + + for (jiter = 0; jiter < nsliceiter; jiter++) + { + if (minslice > maxslice) { + fpixels[2] = minslice - (jiter * incs[2]); + } else { + fpixels[2] = minslice + (jiter * incs[2]); + } + + lpixels[2] = fpixels[2]; + + for (iiter = 0; iiter < nrowiter; iiter++) + { + if (minrow > maxrow) { + fpixels[1] = minrow - (iiter * incs[1]); + } else { + fpixels[1] = minrow + (iiter * incs[1]); + } + + lpixels[1] = fpixels[1]; + + if (bitpix == 8) + { + ffgsvb(fptr, 1, naxis, naxes, fpixels, lpixels, incs, 0, + (unsigned char *) buffer, &anynull, status); + + ffpprb(newptr, 1, firstpix, outsize, (unsigned char *) buffer, status); + } + else if (bitpix == 16) + { + ffgsvi(fptr, 1, naxis, naxes, fpixels, lpixels, incs, 0, + (short *) buffer, &anynull, status); + + ffppri(newptr, 1, firstpix, outsize, (short *) buffer, status); + } + else if (bitpix == 32) + { + ffgsvk(fptr, 1, naxis, naxes, fpixels, lpixels, incs, 0, + (int *) buffer, &anynull, status); + + ffpprk(newptr, 1, firstpix, outsize, (int *) buffer, status); + } + else if (bitpix == -32) + { + ffgsve(fptr, 1, naxis, naxes, fpixels, lpixels, incs, FLOATNULLVALUE, + (float *) buffer, &anynull, status); + + ffppne(newptr, 1, firstpix, outsize, (float *) buffer, FLOATNULLVALUE, status); + } + else if (bitpix == -64) + { + ffgsvd(fptr, 1, naxis, naxes, fpixels, lpixels, incs, DOUBLENULLVALUE, + buffer, &anynull, status); + + ffppnd(newptr, 1, firstpix, outsize, buffer, DOUBLENULLVALUE, + status); + } + else if (bitpix == 64) + { + ffgsvjj(fptr, 1, naxis, naxes, fpixels, lpixels, incs, 0, + (LONGLONG *) buffer, &anynull, status); + + ffpprjj(newptr, 1, firstpix, outsize, (LONGLONG *) buffer, status); + } + + + firstpix += outsize; + } + } + } + + free(buffer); /* finished with the memory */ + + if (*status > 0) + { + ffpmsg("fits_copy_image_section: error copying image section"); + return(*status); + } + + return(*status); } diff --git a/data-access/engine/src/vlkb/src/main.cpp b/data-access/engine/src/vlkb/src/main.cpp index 5ed5a64..7da51fb 100644 --- a/data-access/engine/src/vlkb/src/main.cpp +++ b/data-access/engine/src/vlkb/src/main.cpp @@ -259,10 +259,16 @@ int cmd_cutout(int argc, char * argv[]) if (argc != 4) { std::cerr - << "Usage: cutout <filename.fits> <extnum> <region>\n" - << "\n" - << "Calculate overlap between HDU in file and region.\n\nregion in JSON form of VO-SODA params. For example 'POS=CIRCLE 21.4458 -1.373 0.1' :\n" - " \'{\"pos\":{\"circle\":{\"lat\":-1.373,\"lon\":21.4458,\"radius\":0.1},\"system\":\"ICRS\"},\"service\":\"SUBIMG\"}\'\n"; + << "Usage: cutout <filename.fits> <extnum> <region>" + << endl + << endl + << "Calculate overlap between HDU in file and region. Returns 0 (ok), 1 (no overlap) or EXIT_FAILURE" + << endl + << endl + << "region in JSON form of VO-SODA params. For example 'POS=CIRCLE 21.4458 -1.373 0.1' :" + << endl + << " \'{\"pos\":{\"circle\":{\"lat\":-1.373,\"lon\":21.4458,\"radius\":0.1},\"system\":\"ICRS\"},\"service\":\"SUBIMG\"}\'" + << endl; return EXIT_FAILURE; } else @@ -272,8 +278,8 @@ int cmd_cutout(int argc, char * argv[]) string region{argv[3]}; int rc = stream_cutout(pathname, extnum, region); - if(!rc) - return EXIT_SUCCESS; + if( (rc == 0) || (rc == 1)) + return rc; // OK or no overlap } return EXIT_FAILURE; } diff --git a/data-access/servlet/src/main/java/cutout/Soda.java b/data-access/servlet/src/main/java/cutout/Soda.java index fa10472..7d510f8 100644 --- a/data-access/servlet/src/main/java/cutout/Soda.java +++ b/data-access/servlet/src/main/java/cutout/Soda.java @@ -14,7 +14,7 @@ import vo.parameter.*; public interface Soda { - public void doStream(String relPathname, int hdunum, + public int doStream(String relPathname, int hdunum, Pos pos, Band band, Time time, Pol pol, String pixels, OutputStream outputStream) throws IOException, InterruptedException; diff --git a/data-access/servlet/src/main/java/cutout/SodaImpl.java b/data-access/servlet/src/main/java/cutout/SodaImpl.java index 2bc94f1..b0374bd 100644 --- a/data-access/servlet/src/main/java/cutout/SodaImpl.java +++ b/data-access/servlet/src/main/java/cutout/SodaImpl.java @@ -43,237 +43,67 @@ class SodaImpl implements Soda } - public void doStream(String relPathname, int hdunum, + public int doStream(String relPathname, int hdunum, Pos pos, Band band, Time time, Pol pol, String pixels, OutputStream outputStream) throws IOException, InterruptedException { - LOGGER.fine("trace"); - Instant start = Instant.now(); - - boolean has_overlap = false; - boolean pixels_valid = (pixels != null); - - String boundsString = ""; - - final boolean isDavCall = relPathname.startsWith("http://") || relPathname.startsWith("https://"); - final boolean isAbsPath = relPathname.startsWith("/"); // Resolver removes schema from file:// - // file:///some_abs_path -> /soma_abs_path - // file://some_rel_path -> some_rel_path - - String absPathname = (isDavCall || isAbsPath) ? relPathname : (fitsPaths.surveys() +"/"+ relPathname); - - if( !pixels_valid ) - { - ByteArrayOutputStream bos = new ByteArrayOutputStream(); - if(bos == null) - throw new AssertionError("byte output stream for bounds was not created, is null"); - - JsonEncoder jReq = new JsonEncoder(); - jReq.add(pos); - jReq.add(band); - jReq.add(time); - jReq.add(pol); - String coordString = jReq.toString(); - LOGGER.finest("coordString: " + coordString); - - /* calc bounds */ - - String[] cmdBounds = new String[4]; - cmdBounds[0] = "/usr/local/bin/vlkb"; - cmdBounds[1] = "overlap"; - cmdBounds[2] = absPathname; - cmdBounds[3] = coordString; - - LOGGER.finest(String.join(" ", cmdBounds)); - - ExecCmd execBounds = new ExecCmd(); - execBounds.doRun(bos, cmdBounds); - LOGGER.finest("execBounds exitValue: " + execBounds.exitValue); - - boolean has_result = (execBounds.exitValue == 0); - - if(has_result) - { - boundsString = new String(bos.toByteArray()); - - /* FIXME disable GRID not supported - boundsString = replaceWithGrid(boundsString, pos, band, time, pol); - LOGGER.finest("boundsString(with GRID): " + boundsString); - */ - LOGGER.warning("GRID support was disabled"); - has_overlap = !((boundsString != null) && boundsString.trim().isEmpty()); - - if( !has_overlap ) - { - throw new IllegalArgumentException( - "region in file does not overlap with region defined by SODA parameters"); - } - } - bos.close(); - - Instant boundsDone = Instant.now(); - LOGGER.finer("EXECTIME boundsDone: " + Duration.between(start, boundsDone)); - } - - if(has_overlap || pixels_valid) - { - /* cutout -> outputStream */ - - String pixFilterString = pixels_valid ? pixels : boundsString; - - String[] cmdCut = new String[5]; - cmdCut[0] = "/usr/local/bin/vlkb"; - cmdCut[1] = isDavCall ? "imcopydav" : "imcopy"; - cmdCut[2] = absPathname; - cmdCut[3] = String.valueOf(hdunum-1); - cmdCut[4] = pixFilterString; - - LOGGER.finest(String.join(" ", cmdCut)); - - if(outputStream == null) - LOGGER.finest("supplied outputStream for cut-file is null"); - - ExecCmd execCut = new ExecCmd(); - execCut.doRun(outputStream, cmdCut); - - LOGGER.finest("execCut exitValue: " + execCut.exitValue); - - boolean cut_successful = (execCut.exitValue == 0); - - if(!cut_successful) - { - throw new IllegalArgumentException("cut by pixels not completed for pixels : " + pixFilterString); - } - - Instant cutDone = Instant.now(); - LOGGER.finer("EXECTIME cutDone: " + Duration.between(start, cutDone)); - } - else - { - throw new IllegalArgumentException( - "overlap computation could not be completed with the given arguments"); - } - } - - - private String genRegionForVlkbOverlapCmd(Pos pos, Band band) - { LOGGER.fine("trace"); - String region = ""; - - if(pos != null) - { - String skySystem = pos.system.name(); - - if(pos.shape.equals("CIRCLE")) - { - double l = pos.circle.lon; - double b = pos.circle.lat; - double r = pos.circle.radius; - region = region + "skysystem=" + skySystem + "&l=" + String.valueOf(l) + "&b=" + String.valueOf(b) - + "&r=" + String.valueOf(r); - } - else if(pos.shape.equals("RANGE")) - { - double l = (pos.range.lon1 + pos.range.lon2)/2.0; - double b = (pos.range.lat1 + pos.range.lat2)/2.0; - double dl = (pos.range.lon2 - pos.range.lon1); - double db = (pos.range.lat2 - pos.range.lat1); - region = region + "skysystem=" + skySystem + "&l=" + String.valueOf(l) + "&b=" + String.valueOf(b) - + "&dl=" + String.valueOf(dl) + "&db=" + String.valueOf(db); - } - else - { - LOGGER.finest("FIXME here Exception: POLYGON not supported or pos.shape invalid: " + pos.shape); - } - - } - - if(band != null) - { - String specSystem = band.system.name(); - double vl = band.min; - double vu = band.max; - - region =region + "specsystem=" + specSystem + "&vl=" + String.valueOf(vl) + "&vu=" + String.valueOf(vu); - } - - return region; - } - -/* - private String replaceWithGrid(String wcsBounds, Pos pos, Band band, Time time, Pol pol) - { - // remove end-of-line (was added by vlkb_ast.cpp: cout << ... << endl) - String lineSeparator = System.lineSeparator(); - wcsBounds = wcsBounds.replace(lineSeparator, ""); - LOGGER.finest("BOUNDS: " + wcsBounds); + if(outputStream == null) + LOGGER.finest("supplied outputStream for cut-file is null"); + // FIXME throw excpetion here - // replace in wcsBounds those bounds where pos,band,time or pol has system=GRID + //boolean pixels_valid = (pixels != null); + //final boolean isDavCall = relPathname.startsWith("http://") + // || relPathname.startsWith("https://"); + //final boolean isAbsPath = relPathname.startsWith("/"); // Resolver removes schema from file:// + // file:///some_abs_path -> /soma_abs_path + // file://some_rel_path -> some_rel_path + //String absPathname = (isDavCall || isAbsPath) ? relPathname + // : (fitsPaths.surveys() +"/"+ relPathname); - String[] substr = wcsBounds.split("(?=AXES)", 2); - for(String ss : substr) LOGGER.finest("boundsResult: " + ss); - - String boundsString = substr[0]; - boolean noOverlap = ((boundsString != null) && boundsString.trim().isEmpty()); - if(noOverlap) - { - boundsString = ""; // no overlap - } - else - { - String axesTypes = ""; - if(substr.length > 1) - { - axesTypes = substr[1].replace("AXES"," ").trim(); - LOGGER.finest("AXES TYPES: " + axesTypes); - - String[] bnds = normalize(boundsString); - String[] types = normalize(axesTypes); - // assert: bnds.length == types.length - LOGGER.finest("boundsCount: " + bnds.length + " typesCount: " + types.length); - - if(bnds.length == types.length) - boundsString = replaceBounds(bnds, types, pos, band); - } - } - return boundsString; - } - - private String replaceBounds(String[] bnds, String[] types, Pos pos, Band band) - { - int ix; - for(ix=0; ix<bnds.length; ix++) - { - if( types[ix].equals("LON") && ((pos != null) && (pos.system == Pos.System.GRID)) ) - { - bnds[ix] = pos.lonBoundsString(); - } - else if(types[ix].equals("LAT") && ((pos != null) && (pos.system == Pos.System.GRID))) - { - bnds[ix] = pos.latBoundsString(); - } - else if(types[ix].equals("BAND") && ((band != null) && (band.system == Band.System.GRID))) - { - bnds[ix] = band.boundsString(); - } - } - - LOGGER.finest("replaced: " + String.join(" ", bnds)) ; - - return "[" + String.join(" ", bnds) + "]"; - } + boolean has_overlap = false; + String boundsString = ""; + JsonEncoder jReq = new JsonEncoder(); + jReq.add(pos); + jReq.add(band); + jReq.add(time); + jReq.add(pol); + String coordString = jReq.toString(); + LOGGER.finest("coordString: " + coordString); + + String[] cmd = new String[5]; + cmd[0] = "/usr/local/bin/vlkb"; + cmd[1] = "cutout"; + cmd[2] = fitsPaths.surveys() + "/" + relPathname; + cmd[3] = String.valueOf(hdunum); + cmd[4] = coordString; + LOGGER.finest(String.join(" ", cmd)); + + StringBuilder errorStrBuilder = new StringBuilder(); + + ExecCmd exec = new ExecCmd(); + exec.doRun(outputStream, cmd, errorStrBuilder); + int rc = exec.exitValue; + LOGGER.finest("exec cutout exitValue: " + rc); + LOGGER.finer("EXECTIME cutoutDone: " + Duration.between(start, Instant.now())); + + if(rc == 0) + { + return rc; // OK + } + else if(rc == 1) + { + return rc; // OK, but no overlap + } + else + { + throw new IllegalArgumentException( + "overlap computation could not be completed with the given arguments. " + + errorStrBuilder.toString()); + } + } - // MAKE SURE vlkb overlap returns space delimited bounds: [ a:b c:d ] - // normalize: strip [,] if any, and split into array by space - private String[] normalize(String spaceStr) - { - String other = spaceStr.replace("[","").replace("]",""); - LOGGER.finest("normalize: " + other); - return other.split("\\s+"); - } -*/ } diff --git a/data-access/servlet/src/main/java/cutout/engine/cli/ExecCmd.java b/data-access/servlet/src/main/java/cutout/engine/cli/ExecCmd.java index cee6e77..61a87d3 100644 --- a/data-access/servlet/src/main/java/cutout/engine/cli/ExecCmd.java +++ b/data-access/servlet/src/main/java/cutout/engine/cli/ExecCmd.java @@ -11,10 +11,14 @@ class StreamGobbler extends Thread InputStream is; String type; OutputStream os; + StringBuilder sb; - StreamGobbler(InputStream is, String type) + StreamGobbler(InputStream is, String type, StringBuilder sb) { - this(is, type, null); + this.is = is; + this.type = type; + this.os = null; + this.sb = sb; } StreamGobbler(InputStream is, String type, OutputStream redirect) @@ -48,7 +52,9 @@ class StreamGobbler extends Thread } else { - LOGGER.finest(type + ">" + new String(buffer, 0, nread, "utf-8")); + String str = new String(buffer, 0, nread, "utf-8"); + LOGGER.finest(type + ">" + str);//new String(buffer, 0, nread, "utf-8")); + sb.append(str); } } @@ -69,7 +75,7 @@ class ExecCmd public int exitValue; - public void doRun(OutputStream outputStream, String[] cmd) + public void doRun(OutputStream outputStream, String[] cmd, StringBuilder error) throws IOException, InterruptedException { LOGGER.fine("trace CMD: " + Arrays.toString(cmd)); @@ -82,7 +88,7 @@ class ExecCmd Process proc = rt.exec(cmd); // any error message? - StreamGobbler errorGobbler = new StreamGobbler(proc.getErrorStream(), "ERROR"); + StreamGobbler errorGobbler = new StreamGobbler(proc.getErrorStream(), "ERROR", error); // any output? StreamGobbler outputGobbler = new StreamGobbler(proc.getInputStream(), "OUTPUT", outputStream); diff --git a/data-access/servlet/src/main/java/cutout/webapi/ServletCutout.java b/data-access/servlet/src/main/java/cutout/webapi/ServletCutout.java index 2d441d6..630cb6c 100644 --- a/data-access/servlet/src/main/java/cutout/webapi/ServletCutout.java +++ b/data-access/servlet/src/main/java/cutout/webapi/ServletCutout.java @@ -115,7 +115,7 @@ public class ServletCutout extends HttpServlet } - protected void doCutoutStream(String id, Pos pos, Band band, Time time, Pol pol, String pixels, + protected int doCutoutStream(String id, Pos pos, Band band, Time time, Pol pol, String pixels, OutputStream respOutputStream) throws IOException, InterruptedException { LOGGER.fine("trace " + pos); @@ -128,7 +128,7 @@ public class ServletCutout extends HttpServlet resolver.resolve(id); - soda.doStream(resolver.relPathname(), resolver.hdunum(), pos, band, time, pol, pixels, respOutputStream); + return soda.doStream(resolver.relPathname(), resolver.hdunum(), pos, band, time, pol, pixels, respOutputStream); } @@ -270,7 +270,8 @@ public class ServletCutout extends HttpServlet if(respFormat.startsWith("application/fits")) { response.setContentType(respFormat); - doCutoutStream(id, pos, band, time, pol, pixels, respOutputStream); + int rc = doCutoutStream(id, pos, band, time, pol, pixels, respOutputStream); + if(rc == 1) response.setStatus(HttpServletResponse.SC_NO_CONTENT); } else if(respFormat.startsWith("application/x-vlkb+xml")) { diff --git a/data-access/servlet/src/main/java/mcutout/VlkbCli.java b/data-access/servlet/src/main/java/mcutout/VlkbCli.java index 95cb73e..5cb4af6 100644 --- a/data-access/servlet/src/main/java/mcutout/VlkbCli.java +++ b/data-access/servlet/src/main/java/mcutout/VlkbCli.java @@ -347,8 +347,9 @@ class VlkbCli implements Vlkb cmdBounds[1] = "nullvals"; cmdBounds[2] = absPathname; + StringBuilder errorStrBuilder = new StringBuilder(); ExecCmd exec = new ExecCmd(); - exec.doRun(bos, cmdBounds); + exec.doRun(bos, cmdBounds, errorStrBuilder); LOGGER.finest("exec NullVals exitValue: " + exec.exitValue); bos.close(); @@ -378,7 +379,8 @@ class VlkbCli implements Vlkb } else { - throw new AssertionError("'vlkb nullvals' exited without results for: " + absPathname); + throw new AssertionError("'vlkb nullvals' exited with error: " + + errorStrBuilder.toString() + " for: " + absPathname); } } -- GitLab