27 #include <boost/filesystem.hpp>
37 namespace SourceXtractor {
43 itsNaxes[0] = naxes[0];
44 itsNaxes[1] = naxes[1];
45 itsGridCellSize[0] = gridCellSize[0];
46 itsGridCellSize[1] = gridCellSize[1];
47 itsNGrid[0] = nGrid[0];
48 itsNGrid[1] = nGrid[1];
49 itsNGridPoints = nGrid[0] * nGrid[1];
51 itsGridData = gridData;
52 itsDerivData = makeSplineDeriv(itsNGrid, gridData);
54 itsBackLine =
new PIXTYPE[itsNaxes[0]];
56 itsMedianValue = computeMedian(itsGridData, itsNGridPoints);
60 itsGridData = loadModelFromFits(modelFile);
61 itsDerivData = makeSplineDeriv(itsNGrid, itsGridData);
65 return itsGridCellSize;
73 return itsNGridPoints;
81 return itsMedianValue;
96 delete[] itsDerivData;
101 delete[] itsBackLine;
107 long naxes[2] = { 0, 0 };
108 long fpixel[2] = { 1, 1 };
109 fitsfile *outFits = NULL;
114 if (fitsName.string().size() < 1) {
118 if (boost::filesystem::exists(fitsName)) {
120 boost::filesystem::remove(fitsName);
129 fits_create_file(&outFits, fitsName.string().c_str(), &status);
137 naxes[0] = itsNGrid[0];
138 naxes[1] = itsNGrid[1];
139 fits_create_img(outFits, -32, 2, naxes, &status);
141 throw Elements::Exception() <<
"Problems creating an image extension for: " << fitsName <<
"!";
146 fits_write_pixnull(outFits, TFLOAT, fpixel, (
long) itsNGrid[0] * itsNGrid[1], itsGridData, &undefNumber, &status);
148 fits_report_error(stderr, status);
149 throw Elements::Exception() <<
"Problems writing data to the FITS file: " << fitsName <<
"!";
158 fits_write_key(outFits, TLONG,
"NXIMG1", &itsNaxes[0],
"x-dimension of the original image", &status);
160 fits_report_error(stderr, status);
161 throw Elements::Exception() <<
"Problems writing keyword 'NXIMG1' to the FITS file: " << fitsName <<
"!";
164 fits_write_key(outFits, TLONG,
"NXIMG2", &itsNaxes[1],
"y-dimension of the original image", &status);
166 fits_report_error(stderr, status);
167 throw Elements::Exception() <<
"Problems writing keyword 'NXIMG2' to the FITS file: " << fitsName <<
"!";
170 fits_write_key(outFits, TLONG,
"BCKGRID1", &itsGridCellSize[0],
"background cell size in x", &status);
172 fits_report_error(stderr, status);
173 throw Elements::Exception() <<
"Problems writing keyword 'BCKGRID1' to the FITS file: " << fitsName <<
"!";
176 fits_write_key(outFits, TLONG,
"BCKGRID2", &itsGridCellSize[1],
"background cell size in y", &status);
178 fits_report_error(stderr, status);
179 throw Elements::Exception() <<
"Problems writing keyword 'BCKGRID2' to the FITS file: " << fitsName <<
"!";
182 fits_write_key(outFits, TFLOAT,
"MEDIAN", &itsMedianValue,
"median value in grid", &status);
184 fits_report_error(stderr, status);
185 throw Elements::Exception() <<
"Problems writing keyword 'MEDIAN' to the FITS file: " << fitsName <<
"!";
190 fits_close_file(outFits, &status);
199 long naxes[2] = { 0, 0 };
200 long fpixel[2] = { 1, 1 };
201 fitsfile *outFits = NULL;
207 if (fitsName.string().size() < 1) {
211 if (boost::filesystem::exists(fitsName)) {
213 boost::filesystem::remove(fitsName);
223 fits_create_file(&outFits, fitsName.string().c_str(), &status);
232 naxes[0] = itsNaxes[0];
233 naxes[1] = itsNaxes[1];
234 fits_create_img(outFits, -32, 2, naxes, &status);
236 throw Elements::Exception() <<
"Problems creating an image extension for: " << fitsName <<
"!";
242 pixBuffer =
new PIXTYPE[itsNaxes[0]];
245 for (
size_t yIndex = 0; yIndex < itsNaxes[1]; yIndex++) {
247 splineLine(pixBuffer, yIndex, 0, itsNaxes[0]);
251 fpixel[1] = (long) yIndex + 1;
254 fits_write_pixnull(outFits, TFLOAT, fpixel, (
long) itsNaxes[0], pixBuffer, &undefNumber, &status);
256 fits_report_error(stderr, status);
257 throw Elements::Exception() <<
"Problems writing data to the FITS file: " << fitsName <<
"!";
264 fits_close_file(outFits, &status);
278 if (itsBackLineY!=y){
279 splineLine (itsBackLine, y, 0, itsNaxes[0]);
282 rValue = itsBackLine[
x];
287 int i, j,
x, yl, nbx, nbxm1, nby, nx, ystep, changepoint;
288 float dx, dx0,
dy, dy3, cdx, cdy, cdy3, temp, xstep, *node, *nodep, *dnode, *blo, *bhi, *dblo, *dbhi, *u;
298 dy = (float) y / itsGridCellSize[1] - 0.5;
299 dy -= (yl = (int) dy);
303 }
else if (yl >= nby - 1) {
310 dy3 = (dy * dy * dy -
dy);
311 cdy3 = (cdy * cdy * cdy - cdy);
313 blo = itsGridData + ystep;
315 dblo = itsDerivData + ystep;
318 node =
new float[nbx];
321 *(nodep++) = cdy * *(blo++) + dy * *(bhi++) + cdy3 * *(dblo++) + dy3 * *(dbhi++);
325 dnode =
new float[nbx];
328 u =
new float[nbxm1];
331 for (x = nbxm1; --
x; nodep++) {
332 temp = -1 / (*(dnode++) + 4);
334 temp *= *(u++) - 6 * (*(nodep + 1) + *(nodep - 1) - 2 * *nodep);
338 for (x = nbx - 2; x--;) {
340 *dnode = (*dnode * temp + *(u--)) / 6.0;
348 dnode = itsDerivData;
353 nx = itsGridCellSize[0];
355 changepoint = nx / 2;
356 dx = (xstep - 1) / 2;
357 dx0 = ((nx + 1) % 2) * xstep / 2;
362 for (x = i = 0, j = xStart; j--; i++, dx += xstep) {
363 if (i == changepoint && x > 0 && x < nbxm1) {
377 for (j = width; j--; i++, dx += xstep) {
378 if (i == changepoint && x > 0 && x < nbxm1) {
388 *(backline++) = (
PIXTYPE) (cdx * (*blo + (cdx * cdx - 1) * *dblo) + dx * (*bhi + (dx * dx - 1) * *dbhi));
395 for (j = width; j--;)
397 *(backline++) = (
PIXTYPE) *node;
408 int x,
y, nbx, nby, nbym1;
409 PIXTYPE *dmap, *dmapt, *mapt, *u, temp;
416 dmap =
new PIXTYPE[nGrid[0] * nGrid[1]];
418 for (x = 0; x < nbx; x++) {
429 for (y = 1; y < nbym1; y++, mapt += nbx) {
430 temp = -1 / (*dmapt + 4);
432 *(dmapt += nbx) = temp;
433 temp *= *(u++) - 6 * (*(mapt + nbx) + *(mapt - nbx) - 2 * *mapt);
437 *(dmapt += nbx) = 0.0;
438 for (y = nby - 2; y--;) {
441 *dmapt = (*dmapt * temp + *(u--)) / 6.0;
457 long fpixel[2] = { 1, 1 };
458 char someComment[80];
460 size_t gridSize[2] = { 0, 0 };
466 fitsfile* inputFits = NULL;
469 fits_open_image(&inputFits, modelFile.string().c_str(), READONLY, &status);
478 fits_get_img_param(inputFits, 2, &bitpix, &naxis, naxes, &status);
480 throw Elements::Exception() <<
"Problem reading parameters from image: " << modelFile <<
"!";
485 throw Elements::Exception() <<
"The image: " << modelFile <<
" has " << naxis <<
"!=2 dimensions!";
490 itsNGrid[0] = (
size_t) naxes[0];
491 itsNGrid[1] = (
size_t) naxes[1];
494 fits_read_key(inputFits, TLONG,
"NXIMG1", &itsNaxes[0], someComment, &status);
496 if (status == KEY_NO_EXIST || status == VALUE_UNDEFINED) {
497 throw Elements::Exception() <<
"Keyword 'NXIMG1' does not exist in image: " << modelFile <<
"!";
501 throw Elements::Exception() <<
"Problem reading keyword 'NXIMG1' from image: " << modelFile <<
"!";
505 fits_read_key(inputFits, TLONG,
"NXIMG2", &itsNaxes[1], someComment, &status);
507 if (status == KEY_NO_EXIST || status == VALUE_UNDEFINED) {
508 throw Elements::Exception() <<
"Keyword 'NXIMG2' does not exist in image: " << modelFile <<
"!";
512 throw Elements::Exception() <<
"Problem reading keyword 'NXIMG2' from image: " << modelFile <<
"!";
518 fits_read_key(inputFits, TLONG,
"BCKGRID1", &itsGridCellSize[0], someComment, &status);
520 if (status == KEY_NO_EXIST || status == VALUE_UNDEFINED) {
521 throw Elements::Exception() <<
"Keyword 'BCKGRID1' does not exist in image: " << modelFile <<
"!";
525 throw Elements::Exception() <<
"Problem reading keyword 'BCKGRID1' from image: " << modelFile <<
"!";
529 fits_read_key(inputFits, TLONG,
"BCKGRID2", &itsGridCellSize[1], someComment, &status);
531 if (status == KEY_NO_EXIST || status == VALUE_UNDEFINED) {
532 throw Elements::Exception() <<
"Keyword 'BCKGRID2' does not exist in image: " << modelFile <<
"!";
536 throw Elements::Exception() <<
"Problem reading keyword 'BCKGRID2' from image: " << modelFile <<
"!";
542 fits_read_key(inputFits, TFLOAT,
"MEDIAN", &itsMedianValue, someComment, &status);
544 if (status == KEY_NO_EXIST || status == VALUE_UNDEFINED) {
545 throw Elements::Exception() <<
"Keyword 'MEDIAN' does not exist in image: " << modelFile <<
"!";
549 throw Elements::Exception() <<
"Problem reading keyword 'MEDIAN' from image: " << modelFile <<
"!";
555 divResult =
std::div(
long(itsNaxes[0]),
long(itsGridCellSize[0]));
556 gridSize[0] =
size_t(divResult.quot);
559 if (gridSize[0] != itsNGrid[0]) {
560 throw Elements::Exception() <<
"Inferred and actual x-size differ in : " << modelFile <<
"!";
565 divResult =
std::div(
long(itsNaxes[1]),
long(itsGridCellSize[1]));
566 gridSize[1] =
size_t(divResult.quot);
569 if (gridSize[1] != itsNGrid[1]) {
570 throw Elements::Exception() <<
"Inferred and actual y-size differ in " << modelFile <<
"!";
574 itsNGridPoints = itsNGrid[0] * itsNGrid[1];
578 throw Elements::Exception() <<
"The model in: " << modelFile <<
" is too large with " << itsNGrid[0]*itsNGrid[1] <<
" pixels!";
585 gridData =
new PIXTYPE[itsNGridPoints];
589 fits_read_pix(inputFits, TFLOAT, fpixel, itsNGrid[0] * itsNGrid[1], &undefNumber, gridData, &anynul, &status);
591 throw Elements::Exception() <<
"Problem reading pixel data from image: " << modelFile <<
"!";
597 throw Elements::Exception() <<
"There are undefined values in the model: " << modelFile <<
"!";
603 fits_close_file(inputFits, &status);
619 tmpArray =
new PIXTYPE[nGridPoints];
620 for (
size_t index = 0; index < nGridPoints; index++) {
621 tmpArray[index] = itsGridData[index];
628 if (::
isnan((
double) median)) {
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
std::shared_ptr< EngineParameter > dx
std::shared_ptr< EngineParameter > dy