SourceXtractorPlusPlus  0.8
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SplineModel.cpp
Go to the documentation of this file.
1 
17 /*
18  * Created on Jan 05, 2015
19  * @author: mkuemmel@usm.lmu.de
20  *
21  * Date: $Date$
22  * Revision: $Revision$
23  * Author: $Author$
24  */
25 
26 #include <math.h>
27 #include <boost/filesystem.hpp> // for boost path type
28 #include "fitsio.h"
29 #include "ElementsKernel/Exception.h" // for Elements Exception
33 #include <iostream>
35 //#define QMALLOC(ptr, typ, nel) ptr = (typ *)malloc((size_t)(nel)*sizeof(typ))
36 
37 namespace SourceXtractor {
38 
39 using namespace std;
40 
41 SplineModel::SplineModel (const size_t* naxes, const size_t* gridCellSize, const size_t* nGrid, PIXTYPE* gridData) {
42 
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];
50 
51  itsGridData = gridData;
52  itsDerivData = makeSplineDeriv(itsNGrid, gridData);
53 
54  itsBackLine = new PIXTYPE[itsNaxes[0]];
55 
56  itsMedianValue = computeMedian(itsGridData, itsNGridPoints);
57 }
58 
59 SplineModel::SplineModel (const boost::filesystem::path modelFile) {
60  itsGridData = loadModelFromFits(modelFile);
61  itsDerivData = makeSplineDeriv(itsNGrid, itsGridData);
62 }
63 
65  return itsGridCellSize;
66 }
67 
68 size_t * SplineModel::getNGrid () {
69  return itsNGrid;
70 }
71 
73  return itsNGridPoints;
74 }
75 
77  return itsNaxes;
78 }
79 
81  return itsMedianValue;
82 }
83 
85  return itsGridData;
86 }
87 
89  // delete the grid data
90  if (itsGridData)
91  delete[] itsGridData;
92  itsGridData = NULL;
93 
94  // delete the derivative data
95  if (itsDerivData)
96  delete[] itsDerivData;
97  itsDerivData = NULL;
98 
99  // delete the line data
100  if (itsBackLine)
101  delete[] itsBackLine;
102  itsBackLine=NULL;
103 }
104 
105 void SplineModel::gridToFits (boost::filesystem::path& fitsName, const bool overwrite) {
106  int status = 0;
107  long naxes[2] = { 0, 0 };
108  long fpixel[2] = { 1, 1 };
109  fitsfile *outFits = NULL;
110  PIXTYPE undefNumber = -BIG;
111 
112  // delete any previous versions
113  // of this file if allowed
114  if (fitsName.string().size() < 1) {
115  throw Elements::Exception() << "No output filename provided!";
116  //Utils::throwElementsException(std::string("No output filename provided!"));
117  } else {
118  if (boost::filesystem::exists(fitsName)) {
119  if (overwrite)
120  boost::filesystem::remove(fitsName);
121  else {
122  throw Elements::Exception() << "The file:" << fitsName << " does already exist!";
123  //Utils::throwElementsException(std::string("The file:") + fitsName.generic_string() + std::string(" does already exist!"));
124  }
125  }
126  }
127 
128  // create a new FITS image
129  fits_create_file(&outFits, fitsName.string().c_str(), &status);
130  if (status) {
131  throw Elements::Exception() << "Problems opening the new FITS file:" << fitsName << "!";
132  //Utils::throwElementsException(std::string("Problems opening the new FITS file: ") + fitsName.generic_string() + std::string("!"));
133  }
134 
135  // create an image HDU as the
136  // primary extension
137  naxes[0] = itsNGrid[0];
138  naxes[1] = itsNGrid[1];
139  fits_create_img(outFits, -32, 2, naxes, &status);
140  if (status) {
141  throw Elements::Exception() << "Problems creating an image extension for: " << fitsName << "!";
142  //Utils::throwElementsException(std::string("Problems creating an image extension for: ") + fitsName.generic_string() + std::string("!"));
143  }
144 
145  // fill in the pixels into the (new) image
146  fits_write_pixnull(outFits, TFLOAT, fpixel, (long) itsNGrid[0] * itsNGrid[1], itsGridData, &undefNumber, &status);
147  if (status) {
148  fits_report_error(stderr, status);
149  throw Elements::Exception() << "Problems writing data to the FITS file: " << fitsName << "!";
150  //Utils::throwElementsException(std::string("Problems writing data to the FITS file: ") + fitsName.generic_string() + std::string("!"));
151  }
152 
153  // store the metadata for the spline model in keywords:
154  // NXIMG1 = itsNaxes[0] the x-dimension of the base image
155  // NXIMG2 = itsNaxes[1] the y-dimension of the base image
156  // BCKGRID1 = itsGridCellSize[0] the grid cell size in x
157  // BCKGRID2 = itsGridCellSize[1] the grid cell size in y
158  fits_write_key(outFits, TLONG, "NXIMG1", &itsNaxes[0], "x-dimension of the original image", &status);
159  if (status) {
160  fits_report_error(stderr, status);
161  throw Elements::Exception() << "Problems writing keyword 'NXIMG1' to the FITS file: " << fitsName << "!";
162  //Utils::throwElementsException(std::string("Problems writing keyword 'NXIMG1' to the FITS file: ") + fitsName.generic_string()+ std::string("!"));
163  }
164  fits_write_key(outFits, TLONG, "NXIMG2", &itsNaxes[1], "y-dimension of the original image", &status);
165  if (status) {
166  fits_report_error(stderr, status);
167  throw Elements::Exception() << "Problems writing keyword 'NXIMG2' to the FITS file: " << fitsName << "!";
168  //Utils::throwElementsException(std::string("Problems writing keyword 'NXIMG2' to the FITS file: ") + fitsName.generic_string()+ std::string("!"));
169  }
170  fits_write_key(outFits, TLONG, "BCKGRID1", &itsGridCellSize[0], "background cell size in x", &status);
171  if (status) {
172  fits_report_error(stderr, status);
173  throw Elements::Exception() << "Problems writing keyword 'BCKGRID1' to the FITS file: " << fitsName << "!";
174  //Utils::throwElementsException(std::string("Problems writing keyword 'BCKGRID1' to the FITS file: ") + fitsName.generic_string()+ std::string("!"));
175  }
176  fits_write_key(outFits, TLONG, "BCKGRID2", &itsGridCellSize[1], "background cell size in y", &status);
177  if (status) {
178  fits_report_error(stderr, status);
179  throw Elements::Exception() << "Problems writing keyword 'BCKGRID2' to the FITS file: " << fitsName << "!";
180  //Utils::throwElementsException(std::string("Problems writing keyword 'BCKGRID2' to the FITS file: ") + fitsName.generic_string()+ std::string("!"));
181  }
182  fits_write_key(outFits, TFLOAT, "MEDIAN", &itsMedianValue, "median value in grid", &status);
183  if (status) {
184  fits_report_error(stderr, status);
185  throw Elements::Exception() << "Problems writing keyword 'MEDIAN' to the FITS file: " << fitsName << "!";
186  //Utils::throwElementsException(std::string("Problems writing keyword 'MEDIAN' to the FITS file: ") + fitsName.generic_string()+ std::string("!"));
187  }
188 
189  // close the FITS file
190  fits_close_file(outFits, &status);
191  if (status) {
192  throw Elements::Exception() << "Problems closing the FITS file: " << fitsName << "!";
193  //Utils::throwElementsException(std::string("Problems closing the FITS file: ") + fitsName.generic_string() + std::string("!"));
194  }
195 }
196 
197 void SplineModel::toFits (boost::filesystem::path& fitsName, const bool overwrite) {
198  int status = 0;
199  long naxes[2] = { 0, 0 };
200  long fpixel[2] = { 1, 1 };
201  fitsfile *outFits = NULL;
202  PIXTYPE undefNumber = -BIG;
203  PIXTYPE* pixBuffer = NULL;
204 
205  // delete any previous versions
206  // of this file if allowed
207  if (fitsName.string().size() < 1) {
208  throw Elements::Exception() << "No output filename provided!";
209  //Utils::throwElementsException(std::string("No output filename provided!"));
210  } else {
211  if (boost::filesystem::exists(fitsName)) {
212  if (overwrite)
213  boost::filesystem::remove(fitsName);
214  else {
215  throw Elements::Exception() << "The file:" << fitsName << " does already exist!";
216  //Utils::throwElementsException(
217  // std::string("The file:") + fitsName.generic_string() + std::string(" does already exist!"));
218  }
219  }
220  }
221 
222  // create a new FITS image
223  fits_create_file(&outFits, fitsName.string().c_str(), &status);
224  if (status) {
225  throw Elements::Exception() << "Problems opening the new FITS file:" << fitsName << "!";
226  //Utils::throwElementsException(
227  // std::string("Problems opening the new FITS file:") + fitsName.generic_string() + std::string("!"));
228  }
229 
230  // create an image HDU as the
231  // primary extension
232  naxes[0] = itsNaxes[0];
233  naxes[1] = itsNaxes[1];
234  fits_create_img(outFits, -32, 2, naxes, &status);
235  if (status) {
236  throw Elements::Exception() << "Problems creating an image extension for: " << fitsName << "!";
237  //Utils::throwElementsException(
238  // std::string("Problems creating an image extension for: ") + fitsName.generic_string() + std::string("!"));
239  }
240 
241  // allocate memory for the buffer
242  pixBuffer = new PIXTYPE[itsNaxes[0]];
243 
244  // go over all rows
245  for (size_t yIndex = 0; yIndex < itsNaxes[1]; yIndex++) {
246  // evaluate the spline for one row
247  splineLine(pixBuffer, yIndex, 0, itsNaxes[0]);
248 
249  // set the new starting position
250  // at the beginning of the row
251  fpixel[1] = (long) yIndex + 1;
252 
253  // fill in the pixels of the row into the (new) image
254  fits_write_pixnull(outFits, TFLOAT, fpixel, (long) itsNaxes[0], pixBuffer, &undefNumber, &status);
255  if (status) {
256  fits_report_error(stderr, status);
257  throw Elements::Exception() << "Problems writing data to the FITS file: " << fitsName << "!";
258  //Utils::throwElementsException(
259  // std::string("Problems writing data to the FITS file: ") + fitsName.generic_string() + std::string("!"));
260  }
261  }
262 
263  // close the FITS file
264  fits_close_file(outFits, &status);
265  if (status) {
266  throw Elements::Exception() << "Problems closing the FITS file: " << fitsName << "!";
267  //Utils::throwElementsException(
268  // std::string("Problems writing data to the FITS file: ") + fitsName.generic_string() + std::string("!"));
269  }
270 
271  // release the memory
272  if (pixBuffer)
273  delete[] pixBuffer;
274 }
275 
277  PIXTYPE rValue;
278  if (itsBackLineY!=y){
279  splineLine (itsBackLine, y, 0, itsNaxes[0]);
280  itsBackLineY=y;
281  }
282  rValue = itsBackLine[x];
283  return rValue;
284 }
285 
286 void SplineModel::splineLine (PIXTYPE *line, const size_t y, const size_t xStart, const size_t width) {
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;
289  PIXTYPE *backline;
290 
291  backline = line;
292 
293  nbx = itsNGrid[0];
294  nbxm1 = nbx - 1;
295  nby = itsNGrid[1];
296 
297  if (nby > 1) {
298  dy = (float) y / itsGridCellSize[1] - 0.5;
299  dy -= (yl = (int) dy);
300  if (yl < 0) {
301  yl = 0;
302  dy -= 1.0;
303  } else if (yl >= nby - 1) {
304  yl = nby - 2;
305  dy += 1.0;
306  }
307 
308  /*-- Interpolation along y for each node */
309  cdy = 1 - dy;
310  dy3 = (dy * dy * dy - dy);
311  cdy3 = (cdy * cdy * cdy - cdy);
312  ystep = nbx * yl;
313  blo = itsGridData + ystep;
314  bhi = blo + nbx;
315  dblo = itsDerivData + ystep;
316  dbhi = dblo + nbx;
317  //QMALLOC(node, float, nbx); /* Interpolated background */
318  node = new float[nbx]; /* Interpolated background */
319  nodep = node;
320  for (x = nbx; x--;)
321  *(nodep++) = cdy * *(blo++) + dy * *(bhi++) + cdy3 * *(dblo++) + dy3 * *(dbhi++);
322 
323  /*-- Computation of 2nd derivatives along x */
324  //QMALLOC(dnode, float, nbx); /* 2nd derivative along x */
325  dnode = new float[nbx]; /* 2nd derivative along x */
326  if (nbx > 1) {
327  //QMALLOC(u, float, nbxm1); /* temporary array */
328  u = new float[nbxm1]; /* temporary array */
329  *dnode = *u = 0.0; /* "natural" lower boundary condition */
330  nodep = node + 1;
331  for (x = nbxm1; --x; nodep++) {
332  temp = -1 / (*(dnode++) + 4);
333  *dnode = temp;
334  temp *= *(u++) - 6 * (*(nodep + 1) + *(nodep - 1) - 2 * *nodep);
335  *u = temp;
336  }
337  *(++dnode) = 0.0; /* "natural" upper boundary condition */
338  for (x = nbx - 2; x--;) {
339  temp = *(dnode--);
340  *dnode = (*dnode * temp + *(u--)) / 6.0;
341  }
342  delete[] u;
343  dnode--;
344  }
345  } else {
346  /*-- No interpolation and no new 2nd derivatives needed along y */
347  node = itsGridData;
348  dnode = itsDerivData;
349  }
350 
351  /*-- Interpolation along x */
352  if (nbx > 1) {
353  nx = itsGridCellSize[0];
354  xstep = 1.0 / nx;
355  changepoint = nx / 2;
356  dx = (xstep - 1) / 2; /* dx of the first pixel in the row */
357  dx0 = ((nx + 1) % 2) * xstep / 2; /* dx of the 1st pixel right to a bkgnd node */
358  blo = node;
359  bhi = node + 1;
360  dblo = dnode;
361  dbhi = dnode + 1;
362  for (x = i = 0, j = xStart; j--; i++, dx += xstep) {
363  if (i == changepoint && x > 0 && x < nbxm1) {
364  blo++;
365  bhi++;
366  dblo++;
367  dbhi++;
368  dx = dx0;
369  }
370  cdx = 1 - dx;
371  if (i == nx) {
372  x++;
373  i = 0;
374  }
375  }
376 
377  for (j = width; j--; i++, dx += xstep) {
378  if (i == changepoint && x > 0 && x < nbxm1) {
379  blo++;
380  bhi++;
381  dblo++;
382  dbhi++;
383  dx = dx0;
384  }
385  cdx = 1 - dx;
386 
387  //*(line++) -= (*(backline++) = (PIXTYPE)(cdx*(*blo+(cdx*cdx-1)**dblo) + dx*(*bhi+(dx*dx-1)**dbhi)));
388  *(backline++) = (PIXTYPE) (cdx * (*blo + (cdx * cdx - 1) * *dblo) + dx * (*bhi + (dx * dx - 1) * *dbhi));
389  if (i == nx) {
390  x++;
391  i = 0;
392  }
393  }
394  } else
395  for (j = width; j--;)
396  //*(line++) -= (*(backline++) = (PIXTYPE)*node);
397  *(backline++) = (PIXTYPE) *node;
398 
399  if (nby > 1) {
400  delete[] node;
401  delete[] dnode;
402  }
403 
404  return;
405 }
406 
407 PIXTYPE* SplineModel::makeSplineDeriv (const size_t* nGrid, PIXTYPE* gridData) {
408  int x, y, nbx, nby, nbym1;
409  PIXTYPE *dmap, *dmapt, *mapt, *u, temp;
410  //float *dmap,*dmapt,*mapt, *u, temp;
411  nbx = nGrid[0];
412  nby = nGrid[1];
413  nbym1 = nby - 1;
414  //QMALLOC(dmap, float, nGrid[0]*nGrid[1]);
415 
416  dmap = new PIXTYPE[nGrid[0] * nGrid[1]];
417 
418  for (x = 0; x < nbx; x++) {
419  mapt = gridData + x;
420  dmapt = dmap + x;
421 
422  if (nby > 1) {
423  //QMALLOC(u, float, nbym1); // temporary array
424  u = new PIXTYPE[nbym1]; // temporary array
425  *dmapt = *u = 0.0; // "natural" lower boundary condition
426 
427  mapt += nbx;
428 
429  for (y = 1; y < nbym1; y++, mapt += nbx) {
430  temp = -1 / (*dmapt + 4);
431 
432  *(dmapt += nbx) = temp;
433  temp *= *(u++) - 6 * (*(mapt + nbx) + *(mapt - nbx) - 2 * *mapt);
434 
435  *u = temp;
436  }
437  *(dmapt += nbx) = 0.0; /* "natural" upper boundary condition */
438  for (y = nby - 2; y--;) {
439  temp = *dmapt;
440  dmapt -= nbx;
441  *dmapt = (*dmapt * temp + *(u--)) / 6.0;
442  }
443  //free(u);
444  delete[] u;
445  } else
446  *dmapt = 0.0;
447  }
448  return dmap;
449 }
450 
451 PIXTYPE* SplineModel::loadModelFromFits (const boost::filesystem::path modelFile) {
452  int status = 0;
453  int bitpix;
454  int naxis;
455  int anynul = 0;
456  long naxes[2];
457  long fpixel[2] = { 1, 1 };
458  char someComment[80];
459 
460  size_t gridSize[2] = { 0, 0 };
461  ldiv_t divResult;
462 
463  PIXTYPE* gridData;
464  PIXTYPE undefNumber = -BIG;
465 
466  fitsfile* inputFits = NULL;
467 
468  // open the FITS file and store the pointer; give some logger info
469  fits_open_image(&inputFits, modelFile.string().c_str(), READONLY, &status);
470  if (status) {
471  throw Elements::Exception() << "Problem opening the image: " << modelFile << "!";
472  //Utils::throwElementsException(
473  // std::string("Problem opening the image: ") + modelFile.generic_string() + std::string("!"));
474  }
475 
476  // get the basic image information;
477  // make sure the image has two dimensions
478  fits_get_img_param(inputFits, 2, &bitpix, &naxis, naxes, &status);
479  if (status) {
480  throw Elements::Exception() << "Problem reading parameters from image: " << modelFile << "!";
481  //Utils::throwElementsException(
482  // std::string("Problem reading parameters from image: ") + modelFile.generic_string() + std::string("!"));
483  }
484  if (naxis != 2) {
485  throw Elements::Exception() << "The image: " << modelFile << " has " << naxis <<"!=2 dimensions!";
486  //Utils::throwElementsException(
487  // std::string("The image: ") + modelFile.generic_string() + std::string(" has ") + tostr(naxis)
488  // + std::string("!=2 dimensions!"));
489  }
490  itsNGrid[0] = (size_t) naxes[0];
491  itsNGrid[1] = (size_t) naxes[1];
492 
493  // read in the keywords NXIMG1 and NXIMG2
494  fits_read_key(inputFits, TLONG, "NXIMG1", &itsNaxes[0], someComment, &status);
495  if (status) {
496  if (status == KEY_NO_EXIST || status == VALUE_UNDEFINED) {
497  throw Elements::Exception() << "Keyword 'NXIMG1' does not exist in image: " << modelFile << "!";
498  //Utils::throwElementsException(
499  // std::string("Keyword 'NXIMG1' does not exist in image: ") + modelFile.generic_string() + std::string("!"));
500  }
501  throw Elements::Exception() << "Problem reading keyword 'NXIMG1' from image: " << modelFile << "!";
502  //Utils::throwElementsException(
503  // std::string("Problem reading keyword 'NXIMG1' from image: ") + modelFile.generic_string() + std::string("!"));
504  }
505  fits_read_key(inputFits, TLONG, "NXIMG2", &itsNaxes[1], someComment, &status);
506  if (status) {
507  if (status == KEY_NO_EXIST || status == VALUE_UNDEFINED) {
508  throw Elements::Exception() << "Keyword 'NXIMG2' does not exist in image: " << modelFile << "!";
509  //Utils::throwElementsException(
510  // std::string("Keyword 'NXIMG2' does not exist in image: ") + modelFile.generic_string() + std::string("!"));
511  }
512  throw Elements::Exception() << "Problem reading keyword 'NXIMG2' from image: " << modelFile << "!";
513  //Utils::throwElementsException(
514  // std::string("Problem reading keyword 'NXIMG2' from image: ") + modelFile.generic_string() + std::string("!"));
515  }
516 
517  // read in the keywords BCKGRID1 and BCKGRID2
518  fits_read_key(inputFits, TLONG, "BCKGRID1", &itsGridCellSize[0], someComment, &status);
519  if (status) {
520  if (status == KEY_NO_EXIST || status == VALUE_UNDEFINED) {
521  throw Elements::Exception() << "Keyword 'BCKGRID1' does not exist in image: " << modelFile << "!";
522  //Utils::throwElementsException(
523  // std::string("Keyword 'BCKGRID1' does not exist in image: ") + modelFile.generic_string() + std::string("!"));
524  }
525  throw Elements::Exception() << "Problem reading keyword 'BCKGRID1' from image: " << modelFile << "!";
526  //Utils::throwElementsException(
527  // std::string("Problem reading keyword 'BCKGRID1' from image: ") + modelFile.generic_string() + std::string("!"));
528  }
529  fits_read_key(inputFits, TLONG, "BCKGRID2", &itsGridCellSize[1], someComment, &status);
530  if (status) {
531  if (status == KEY_NO_EXIST || status == VALUE_UNDEFINED) {
532  throw Elements::Exception() << "Keyword 'BCKGRID2' does not exist in image: " << modelFile << "!";
533  //Utils::throwElementsException(
534  // std::string("Keyword 'BCKGRID2' does not exist in image: ") + modelFile.generic_string() + std::string("!"));
535  }
536  throw Elements::Exception() << "Problem reading keyword 'BCKGRID2' from image: " << modelFile << "!";
537  //Utils::throwElementsException(
538  // std::string("Problem reading keyword 'BCKGRID2' from image: ") + modelFile.generic_string() + std::string("!"));
539  }
540 
541  // read in the keyword 'MEDIAN'
542  fits_read_key(inputFits, TFLOAT, "MEDIAN", &itsMedianValue, someComment, &status);
543  if (status) {
544  if (status == KEY_NO_EXIST || status == VALUE_UNDEFINED) {
545  throw Elements::Exception() << "Keyword 'MEDIAN' does not exist in image: " << modelFile << "!";
546  //Utils::throwElementsException(
547  // std::string("Keyword 'MEDIAN' does not exist in image: ") + modelFile.generic_string() + std::string("!"));
548  }
549  throw Elements::Exception() << "Problem reading keyword 'MEDIAN' from image: " << modelFile << "!";
550  //Utils::throwElementsException(
551  // std::string("Problem reading keyword 'MEDIAN' from image: ") + modelFile.generic_string() + std::string("!"));
552  }
553 
554  // make sure that the model x-dimension and the image dimension are identical
555  divResult = std::div(long(itsNaxes[0]), long(itsGridCellSize[0]));
556  gridSize[0] = size_t(divResult.quot);
557  if (divResult.rem)
558  gridSize[0] += 1;
559  if (gridSize[0] != itsNGrid[0]) {
560  throw Elements::Exception() << "Inferred and actual x-size differ in : " << modelFile << "!";
561  //Utils::throwElementsException(
562  // std::string("Inferred and actual x-size differ in : ") + modelFile.generic_string() + std::string("!"));
563  }
564  // make sure that the model y-dimension and the image dimension are identical
565  divResult = std::div(long(itsNaxes[1]), long(itsGridCellSize[1]));
566  gridSize[1] = size_t(divResult.quot);
567  if (divResult.rem)
568  gridSize[1] += 1;
569  if (gridSize[1] != itsNGrid[1]) {
570  throw Elements::Exception() << "Inferred and actual y-size differ in " << modelFile << "!";
571  //Utils::throwElementsException(
572  // std::string("Inferred and actual y-size differ in : ") + modelFile.generic_string() + std::string("!"));
573  }
574  itsNGridPoints = itsNGrid[0] * itsNGrid[1];
575 
576  // make sure not too much R=memory is deployed
577  if (itsNGrid[0] * itsNGrid[1] > BACK_BUFSIZE) {
578  throw Elements::Exception() << "The model in: " << modelFile << " is too large with " << itsNGrid[0]*itsNGrid[1] << " pixels!";
579  //Utils::throwElementsException(
580  // std::string("The model in: ") + modelFile.generic_string() + std::string(" is too large with ")
581  // + tostr(itsNGridPoints) + std::string("pixels!"));
582  }
583 
584  // allocate memory for the data
585  gridData = new PIXTYPE[itsNGridPoints];
586 
587  // read all the pixels in;
588  // make sure there are no undefined values
589  fits_read_pix(inputFits, TFLOAT, fpixel, itsNGrid[0] * itsNGrid[1], &undefNumber, gridData, &anynul, &status);
590  if (status) {
591  throw Elements::Exception() << "Problem reading pixel data from image: " << modelFile << "!";
592  //Utils::throwElementsException(
593  // std::string("Problem reading pixel data from image: ") + modelFile.generic_string() + std::string("!"));
594  }
595  if (anynul) {
596 
597  throw Elements::Exception() << "There are undefined values in the model: " << modelFile << "!";
598  //Utils::throwElementsException(
599  // std::string("There are undefined values in the model: ") + modelFile.generic_string() + std::string("!"));
600  }
601 
602  // close the FITS file
603  fits_close_file(inputFits, &status);
604  if (status) {
605  throw Elements::Exception() << "Problem closing the image: " << modelFile << "!";
606  //Utils::throwElementsException(
607  // std::string("Problem closing the image: ") + modelFile.generic_string() + std::string("!"));
608  }
609  inputFits = NULL;
610 
611  return gridData;
612 }
613 
614 PIXTYPE SplineModel::computeMedian (PIXTYPE* itsGridData, const size_t nGridPoints) {
615  PIXTYPE median;
616  PIXTYPE* tmpArray;
617 
618  // create a tmp array and transfer the data
619  tmpArray = new PIXTYPE[nGridPoints];
620  for (size_t index = 0; index < nGridPoints; index++) {
621  tmpArray[index] = itsGridData[index];
622  }
623 
624  // compute the median
625  median = SE2BackgroundUtils::fqMedian(tmpArray, nGridPoints);
626 
627  // the median can not be NaN, throw an Exception
628  if (::isnan((double) median)) {
629  throw Elements::Exception() << "The median of the background is NaN. This can not be!";
630  //Utils::throwElementsException(std::string("The median of the background is NaN. This can not be!"));
631  }
632 
633  // mop up and return
634  delete[] tmpArray;
635  return median;
636 }
637 
638 } // end of namespace SourceXtractor
#define BACK_BUFSIZE
PIXTYPE computeMedian(PIXTYPE *gridData, const size_t nGridPoints)
PIXTYPE * makeSplineDeriv(const size_t *nGrid, PIXTYPE *gridData)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
void toFits(boost::filesystem::path &fitsName, const bool overwrite=true)
PIXTYPE * loadModelFromFits(const boost::filesystem::path)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
void splineLine(PIXTYPE *line, const size_t y, const size_t xStart, const size_t width)
PIXTYPE getValue(size_t x, size_t y)
T div(T...args)
void gridToFits(boost::filesystem::path &fitsName, const bool overwrite=true)
static PIXTYPE fqMedian(PIXTYPE *ra, size_t n)
T isnan(T...args)
SplineModel(const size_t *naxes, const size_t *gridCellSize, const size_t *nGrid, PIXTYPE *gridData)
Definition: SplineModel.cpp:41
std::shared_ptr< EngineParameter > dx
#define BIG
std::shared_ptr< EngineParameter > dy