SourceXtractorPlusPlus  0.11
Please provide a description of the project.
SE2BackgroundModeller.cpp
Go to the documentation of this file.
1 
17 /*
18  * Created on Jan 05, 2015
19  * @author: mkuemmel@usm.lmu.de
20  * Date: $Date$
21  * Revision: $Revision$
22  * Author: $Author$
23  */
24 #include <iostream>
25 #include <cstdlib>
26 
27 #include "fitsio.h"
28 
36 #define SIZETSUB(X, Y) ((X) > (Y) ? (X-Y) : (Y-X))
37 using namespace std;
38 
39 namespace SourceXtractor {
40 
41 SE2BackgroundModeller::SE2BackgroundModeller(std::shared_ptr<DetectionImage> image, std::shared_ptr<WeightImage> variance_map, std::shared_ptr<Image<unsigned char>> mask, const unsigned char mask_type_flag)
42 {
43  itsImage = image;
44  itsVariance = variance_map;
45  itsMask = mask;
46  itsMaskType = mask_type_flag;
47  //itsWeightTypeFlag = weight_type_flag;
48 
49  //check for variance
50  if (variance_map!=nullptr)
51  itsHasVariance=true;
52  else
53  itsHasVariance=false;
54 
55  //check for mask
56  if (mask!=nullptr)
57  itsHasMask=true;
58  else
59  itsHasMask=false;
60 
61  // store the image size
62  itsNaxes[0] = (size_t)image->getWidth();
63  itsNaxes[1] = (size_t)image->getHeight();
64 }
65 
66 SE2BackgroundModeller::~SE2BackgroundModeller() {
67  if (itsWhtMeanVals)
68  delete[] itsWhtMeanVals;
69 }
70 
71 void SE2BackgroundModeller::createSE2Models(std::shared_ptr<TypedSplineModelWrapper<SeFloat>> &bckPtr, std::shared_ptr<TypedSplineModelWrapper<SeFloat>> &varPtr, PIXTYPE &sigFac, const size_t *bckCellSize, const WeightImage::PixelType varianceThreshold, const size_t *filterBoxSize, const float &filterThreshold)
72 {
73  size_t gridSize[2] = {0,0};
74  size_t nGridPoints=0;
75 
76  long increment[2]={1,1};
77  long fpixel[2];
78  long lpixel[2];
79 
80  size_t nElements=0;
81  size_t nData=0;
82  size_t subImgNaxes[2] = {0,0};
83 
84  PIXTYPE* gridData=NULL;
85  PIXTYPE* weightData=NULL;
86 
87  //PIXTYPE undefNumber=-BIG;
88 
89  BackgroundCell* oneCell=NULL;
90 
91  PIXTYPE* bckMeanVals=NULL;
92  PIXTYPE* bckSigVals=NULL;
93  PIXTYPE* whtSigVals=NULL;
94 
95  size_t gridIndex=0;
96 
97  ldiv_t divResult;
98 
99  PIXTYPE weightVarThreshold=(PIXTYPE)varianceThreshold;
100 
101  // not necessary, since all is in variance
102  // re-scale the weight threshold
103  //if (itsHasVariance){
104  // rescaleThreshold(weightVarThreshold, weightThreshold);
105  //}
106 
107  // get the number of grid cells in x
108  divResult = std::div(long(itsNaxes[0]), long(bckCellSize[0]));
109  gridSize[0] = size_t(divResult.quot);
110  if (divResult.rem)
111  gridSize[0] += 1;
112 
113  // get the number of grid cells in y
114  divResult = std::div(long(itsNaxes[1]), long(bckCellSize[1]));
115  gridSize[1] = size_t(divResult.quot);
116  if (divResult.rem)
117  gridSize[1] += 1;
118 
119  // compute the number of grid points
120  nGridPoints = gridSize[0]*gridSize[1];
121 
122  // create the arrays
123  bckMeanVals = new PIXTYPE[nGridPoints];
124  bckSigVals = new PIXTYPE[nGridPoints];
125  if (itsHasVariance){
126  itsWhtMeanVals = new PIXTYPE[nGridPoints];
127  whtSigVals = new PIXTYPE[nGridPoints];
128  }
129 
130  // give some feedback on the parameters
131  bck_model_logger.debug() << "\tBackground cell size=("<<bckCellSize[0]<<"," << bckCellSize[1]<< ")";
132  bck_model_logger.debug() << "\tFilter box size=("<<filterBoxSize[0]<<"," << filterBoxSize[1]<< ")";
133  bck_model_logger.debug() << "\tThe bad pixel threshold is: "<< weightVarThreshold;
134 
135  // iterate over cells in y
136  gridIndex=0;
137  for (size_t yIndex=0; yIndex<gridSize[1]; yIndex++){
138 
139  // set the boundaries in y
140  fpixel[1] = (long)yIndex*bckCellSize[1];
141  lpixel[1] = yIndex < gridSize[1]-1 ? (long)(yIndex+1)*bckCellSize[1] : (long)itsNaxes[1];
142 
143  // iterate over cells in x
144  for (size_t xIndex=0; xIndex<gridSize[0]; xIndex++){
145 
146  // set the boundaries in x
147  fpixel[0] = (long)xIndex*bckCellSize[0];
148  lpixel[0] = xIndex < gridSize[0]-1 ? (long)(xIndex+1)*bckCellSize[0] : (long)itsNaxes[0];
149 
150  // compute the length of the cell sub-image
151  subImgNaxes[0] =(size_t)(lpixel[0]-fpixel[0]);
152  subImgNaxes[1] =(size_t)(lpixel[1]-fpixel[1]);
153 
154  // some feedback on the corners of the image to be treated
155  bck_model_logger.debug() << "Background cell from fpixel=(" << fpixel[0] << "," << fpixel[1] << ") to lpixel=("<< lpixel[0] << "," << lpixel[1] << ")";
156 
157  // compute the increments to perhaps limit the number
158  // of pixels read in, the total number of elements
159  // and the numbers read in x and y
160  getMinIncr(nElements, increment, subImgNaxes);
161 
162  // define or re-define the buffers
163  if (nElements!=nData) {
164  delete [] gridData;
165  gridData = new PIXTYPE[nElements];
166  //if (itsInputWeight){
167  if (itsHasVariance){
168  delete [] weightData;
169  weightData = new PIXTYPE[nElements];
170  }
171  nData=nElements;
172  }
173 
174  // load in the image data
175  long pixIndex=0;
176  for (auto yIndex=fpixel[1]; yIndex<lpixel[1]; yIndex+=increment[1])
177  for (auto xIndex=fpixel[0]; xIndex<lpixel[0]; xIndex+=increment[0]){
178  gridData[pixIndex++] = (PIXTYPE)itsImage->getValue(int(xIndex), int(yIndex));
179  }
180  if (itsHasVariance){
181  long pixIndex=0;
182  for (auto yIndex=fpixel[1]; yIndex<lpixel[1]; yIndex+=increment[1])
183  for (auto xIndex=fpixel[0]; xIndex<lpixel[0]; xIndex+=increment[0])
184  weightData[pixIndex++] = (PIXTYPE)itsVariance->getValue(int(xIndex), int(yIndex));
185  }
186  if (itsHasMask){
187  // load in the image data
188  long pixIndex=0;
189  for (auto yIndex=fpixel[1]; yIndex<lpixel[1]; yIndex+=increment[1])
190  for (auto xIndex=fpixel[0]; xIndex<lpixel[0]; xIndex+=increment[0], pixIndex++)
191  //if (!itsMask->getValue(int(xIndex), int(yIndex)))
192  if (itsMask->getValue(int(xIndex), int(yIndex)) & itsMaskType){
193  gridData[pixIndex] = -BIG;
194  bck_model_logger.debug() << "\tReplacing data value";
195  }
196  }
197 
198  // create a background cell compute and store the values
199  // and then delete the cell again
200  oneCell = new BackgroundCell(gridData, nElements, weightData, weightVarThreshold);
201  if (itsHasVariance)
202  oneCell->getBackgroundValues(bckMeanVals[gridIndex], bckSigVals[gridIndex], itsWhtMeanVals[gridIndex], whtSigVals[gridIndex]);
203  else
204  oneCell->getBackgroundValues(bckMeanVals[gridIndex], bckSigVals[gridIndex]);
205  delete oneCell;
206 
207  // enhance the grid index
208  gridIndex++;
209  }
210  }
211 
212  // do some filtering on the data
213  filter(bckMeanVals, bckSigVals, gridSize, filterBoxSize, filterThreshold);
214  if (itsHasVariance){
215  filter(itsWhtMeanVals, whtSigVals, gridSize, filterBoxSize, filterThreshold);
216  }
217 
218  if (itsHasVariance){
219  //if (itsHasVariance && (itsWeightTypeFlag & (VAR_FIELD|WEIGHT_FIELD))){
220  // compute the scaling factor
221  computeScalingFactor(itsWhtMeanVals, bckSigVals, sigFac, nGridPoints);
222  }
223  else{
224  sigFac=0.0;
225  }
226 
227  // convert the grid of rms values to variance
228  for (size_t index=0; index<nGridPoints; index++)
229  bckSigVals[index] *= bckSigVals[index];
230 
231  // create the splined interpolation images for background and variance
232  bckPtr = TypedSplineModelWrapper<SeFloat>::create(itsNaxes, bckCellSize, gridSize, bckMeanVals);
233  varPtr = TypedSplineModelWrapper<SeFloat>::create(itsNaxes, bckCellSize, gridSize, bckSigVals);
234 
235  // release memory
236  delete [] whtSigVals;
237  delete [] gridData;
238  delete [] weightData;
239 }
240 
241 void SE2BackgroundModeller::getMinIncr(size_t &nElements, long* incr, const size_t * subImgNaxes)
242 {
243  float axisRatio;
244  ldiv_t divResult;
245 
246  // copy the original dimensions
247  size_t tmpImgNaxes[2] = {subImgNaxes[0], subImgNaxes[1]};
248 
249  // set the default increment
250  incr[0] = 1;
251  incr[1] = 1;
252 
253  // compute the number of elements;
254  // check if something needs to be done at all
255  nElements = tmpImgNaxes[0]*tmpImgNaxes[1];
256  if (nElements <= BACK_BUFSIZE)
257  return;
258 
259  // compute the axis ratio
260  axisRatio = float(tmpImgNaxes[0]) / float(tmpImgNaxes[1]);
261 
262  // iterate until the number is small enough
263  while (nElements > BACK_BUFSIZE)
264  {
265  // change the increments such
266  // that pixels sampled in x and y
267  // are about equal
268  if (axisRatio >= 1.0)
269  incr[0] += 1;
270  else
271  incr[1] += 1;
272 
273  // get the number of pixels sampled in x
274  divResult = std::div(long(subImgNaxes[0]), long(incr[0]));
275  tmpImgNaxes[0] = size_t(divResult.quot);
276  if (divResult.rem)
277  tmpImgNaxes[0] += 1;
278 
279  // get the number of pixels sampled in y
280  divResult = std::div(long(subImgNaxes[1]), long(incr[1]));
281  tmpImgNaxes[1] = size_t(divResult.quot);
282  if (divResult.rem)
283  tmpImgNaxes[1] += 1;
284 
285  // re-compute the number of elements and the
286  // axis ratio
287  nElements = tmpImgNaxes[0]*tmpImgNaxes[1];
288  axisRatio = float(tmpImgNaxes[0]) / float(tmpImgNaxes[1]);
289  }
290 
291  // some feedback for the new increments and sizes
292  bck_model_logger.debug() << "New increment=(" << incr[0] << "," << incr[1] << ") sampledPixels=("<< tmpImgNaxes[0] << "," << tmpImgNaxes[1] << ") nElements=" << nElements;
293 
294  return;
295 }
296 
297 void SE2BackgroundModeller::filter(PIXTYPE* bckVals, PIXTYPE* sigmaVals,const size_t* gridSize, const size_t* filterSize, const float &filterThreshold)
298 {
299  // replace undefined values
300  replaceUNDEF(bckVals, sigmaVals, gridSize);
301 
302  // do a median filtering of the values
303  filterMedian(bckVals, sigmaVals, gridSize, filterSize, filterThreshold);
304 
305  return;
306 }
307 
308 void SE2BackgroundModeller::replaceUNDEF(PIXTYPE* bckVals, PIXTYPE* sigmaVals,const size_t* gridSize)
309 {
310  PIXTYPE *back=NULL;
311  PIXTYPE *sigma=NULL;
312  PIXTYPE *backMod=NULL;
313  PIXTYPE val;
314  PIXTYPE sval;
315  float dist=0.;
316  float distMin=0;
317  size_t iAct,jAct, nx,ny, nmin;
318  size_t np;
319  bool hasHoles=false;
320 
321  // take the sizing information
322  nx = gridSize[0];
323  ny = gridSize[1];
324  np = gridSize[0]*gridSize[1];
325 
326  // check whether something needs to be done
327  for (size_t index=0; index< np; index++)
328  {
329  if (bckVals[index]<=-BIG)
330  {
331  hasHoles=true;
332  break;
333  }
334  }
335 
336  // nothing to do,
337  // just go back
338  if (!hasHoles)
339  return;
340 
341  // give some feedback that undefined data is replaced
342  bck_model_logger.debug() << "\tReplacing undefined data";
343 
344  // vector for the final
345  // background values
346  backMod = new PIXTYPE[np];
347 
348  // go to the array start
349  back = bckVals;
350  sigma = sigmaVals;
351 
352  // look for `bad' meshes and interpolate them if necessary
353  val = 0.0;
354  sval = 0.0;
355  iAct=0;
356  for (size_t py=0; py<ny; py++)
357  {
358  for (size_t px=0; px<nx; px++)
359  {
360  // compute the current index
361  iAct = py*nx+px;
362 
363  // transfer the value
364  backMod[iAct]=back[iAct];
365 
366  // check for undefined data
367  if (back[iAct]<=-BIG)
368  {
369 
370  // seek the closest data points,
371  // search over the entire array
372  distMin = BIG;
373  nmin = 0;
374  for (size_t y=0; y<ny; y++)
375  {
376  for (size_t x=0; x<nx; x++)
377  {
378  // compute the current index
379  jAct = y*nx+x;
380 
381  if (back[jAct]>-BIG)
382  {
383  // compute the pixel distance
384  dist = (float)SIZETSUB(x,px)*SIZETSUB(x,px)+SIZETSUB(y,py)*SIZETSUB(y,py);
385 
386  // check for a new minimum distance
387  if (dist<distMin)
388  {
389  // start new average values
390  val = back[jAct];
391  sval = sigma[jAct];
392  nmin = 1;
393  distMin = dist;
394  }
395  else if (fabs(dist-distMin)<1e-05)
396  {
397  // add equal distance pixel for averaging
398  val += back[jAct];
399  sval += sigma[jAct];
400  nmin++;
401  }
402  }
403  }
404  }
405  // take the mean of the closest
406  // defined data points
407  backMod[iAct] = nmin ? val/nmin: 0.0;
408  sigma[iAct] = nmin ? sval/nmin: 1.0;
409  }
410  }
411  }
412 
413  // push the modified values back
414  for (size_t index=0; index< np; index++)
415  {
416  bckVals[index] =backMod[index];
417  }
418 
419  // release the memory
420  if (backMod)
421  delete [] backMod;
422 
423  return;
424 }
425 
426 void SE2BackgroundModeller::filterMedian(PIXTYPE* bckVals, PIXTYPE* sigmaVals, const size_t* gridSize, const size_t* filterSize, const float filterThresh)
427 {
428  PIXTYPE *back, *sigma, *sigmat;
429  PIXTYPE* backFilt=NULL;
430  PIXTYPE* sigmaFilt=NULL;
431  PIXTYPE* bmask=NULL;
432  PIXTYPE* smask=NULL;
433  PIXTYPE allSigmaMed, median; //allBckMed
434  int i,nx,ny,npx,npx2,npy,npy2,x,y;
435  int np;
436 
437  // check whether something needs to be done at all
438  // note that the code would run nevertheless
439  if (filterSize[0]<2 && filterSize[1]<2)
440  return;
441 
442  // Note: I am converting the "size_t" to int's since
443  // there are computations done down.
444 
445  // take the sizing information
446  nx = (int)gridSize[0];
447  ny = (int)gridSize[1];
448  np = nx*ny;
449 
450  // store the filter size
451  npx = (int)filterSize[0]/2;
452  npy = (int)filterSize[1]/2;
453  npy *= nx;
454 
455  // allocate space for the work area
456  bmask = new PIXTYPE[(2*npx+1)*(2*npy+1)];
457  smask = new PIXTYPE[(2*npx+1)*(2*npy+1)];
458 
459  // allocate space for filtered arrays
460  backFilt = new PIXTYPE[np];
461  sigmaFilt = new PIXTYPE[np];
462 
463  // store the arrays locally
464  back = bckVals;
465  sigma = sigmaVals;
466 
467  // go over all y
468  for (int py=0; py<np; py+=nx)
469  {
470  // limit the filter box in y
471  npy2 = np - py - nx;
472  if (npy2>npy)
473  npy2 = npy;
474  if (npy2>py)
475  npy2 = py;
476 
477  // go over all x
478  for (int px=0; px<nx; px++)
479  {
480  // limit the filter box in x
481  npx2 = nx - px - 1;
482  if (npx2>npx)
483  npx2 = npx;
484  if (npx2>px)
485  npx2 = px;
486 
487  // store all values in the box
488  // in an array
489  i=0;
490  for (int dpy = -npy2; dpy<=npy2; dpy+=nx)
491  {
492  y = py+dpy;
493  for (int dpx = -npx2; dpx <= npx2; dpx++)
494  {
495  x = px+dpx;
496  bmask[i] = back[x+y];
497  smask[i++] = sigma[x+y];
498  }
499  }
500 
501  // compute the median, check
502  // whether the median is above the threshold
503  median = SE2BackgroundUtils::fqMedian(bmask, i);
504  if (fabs((median-back[px+py]))>=(PIXTYPE)filterThresh)
505  {
506  // use the median values
507  backFilt[px+py] = median;
508  sigmaFilt[px+py] = SE2BackgroundUtils::fqMedian(smask, i);
509  }
510  else
511  {
512  // use the original value
513  backFilt[px+py] = back[px+py];
514  sigmaFilt[px+py] = sigma[px+py];
515  }
516  }
517  }
518 
519  // transfer the filtered background values back
520  for (int index=0; index<np; index++)
521  back[index] = backFilt[index];
522 
523  // transfer the filtered sigma values back
524  for (int index=0; index<np; index++)
525  sigma[index] = sigmaFilt[index];
526 
527  // compute the median values for the background
528  // and the sigma
529  //allBckMed = SE2BackgroundUtils::fqMedian(backFilt, np);
530  allSigmaMed = SE2BackgroundUtils::fqMedian(sigmaFilt, np);
531 
532  // NOTE: I don't understand what that does.
533  // Why should the median sigma be <.0?
534  if (allSigmaMed<=0.0)
535  {
536  sigmat = sigmaFilt+np;
537  for (i=np; i-- && *(--sigmat)>0.0;);
538  if (i>=0 && i<(np-1))
539  allSigmaMed = SE2BackgroundUtils::fqMedian(sigmat+1, np-1-i);
540  else
541  {
542  //if (field->flags&(DETECT_FIELD|MEASURE_FIELD))
543  // warning("Image contains mainly constant data; ",
544  // "I'll try to cope with that...");
545  //field->backsig = 1.0;
546  allSigmaMed = 1.0;
547  }
548  }
549 
550  // release memory
551  delete [] sigmaFilt;
552  delete [] backFilt;
553  delete [] bmask;
554  delete [] smask;
555 
556  return;
557 }
558 
559 void SE2BackgroundModeller::computeScalingFactor(PIXTYPE* whtMeanVals, PIXTYPE* bckSigVals, PIXTYPE& sigFac, const size_t nGridPoints)
560 {
561 
562  size_t nr = 0;
563  size_t lowIndex=0;
564  PIXTYPE* ratio=NULL;
565 
566  PIXTYPE actRatio;
567  //PIXTYPE tmp;
568 
569  // allocate memory
570  ratio = new PIXTYPE[nGridPoints];
571 
572  // form the list of ratios between measured sigma values
573  // and the sigma's properly derived from the weights
574 
575  for (size_t index=0; index<nGridPoints; index++){
576  if (whtMeanVals[index]>0.0){
577  //actRatio = bckSigVals[index] / sqrt(whtMeanVals[index]);
578  actRatio = bckSigVals[index] * bckSigVals[index] / whtMeanVals[index]; // scaling factor for the variance image
579  if (actRatio>0.0){
580  ratio[nr]=actRatio;
581  nr++;
582  }
583  }
584  }
585 
586  // use the median ratio as scaling factor
587  sigFac = SE2BackgroundUtils::fqMedian(ratio, nr);
588 
589  // count leading 0.0 values in the list
590  for (lowIndex=0; lowIndex<nr && ratio[lowIndex]<=0.0; lowIndex++);
591 
592  // re-calculated the median, omitting
593  // the leading 0.0 values; I can't see
594  // how there should be leading zeros
595  // TODO: check whether this make sense
596  if (lowIndex>0){
597  if (lowIndex<nr){
598  // make a log message
599  bck_model_logger.debug() << "\tRe-calculating the scaling due to leading zeros";
600  sigFac = SE2BackgroundUtils::fqMedian(ratio+lowIndex, nr-lowIndex);
601  }
602  else {
603  //warning("Null or negative global weighting factor:","defaulted to 1");
604  bck_model_logger.debug() << "\tNull or negative global weighting factor: " << " | " << lowIndex << "defaulted to 1 " << nr;
605  sigFac = 1.0;
606  }
607  }
608 
609  delete [] ratio;
610  return;
611 }
612 
613 void SE2BackgroundModeller::rescaleThreshold(PIXTYPE &weightVarThreshold, const PIXTYPE &weightThreshold)
614 {
615  // the threshold needs to be larger than zero
616  if (weightThreshold<0.0){
617  throw Elements::Exception() << "The weight threshold is: " << weightThreshold << " but can not be smaller than 0.0";
618  }
619 
620  // check the type flag
621  switch (itsWeightTypeFlag){
622 
623  case VAR_FIELD:
624  // just copy for variance
625  weightVarThreshold = weightThreshold;
626  break;
627 
628  case RMS_FIELD:
629  // square for rms
630  weightVarThreshold = weightThreshold*weightThreshold;
631  break;
632 
633  case WEIGHT_FIELD:
634  // give default or invert for weight
635  if (weightThreshold>0.0){
636  weightVarThreshold = 1.0/weightThreshold;
637  }
638  else {
639  weightVarThreshold = BIG;
640  }
641  break;
642 
643  default:
644  // this is the default
645  weightVarThreshold = BIG;
646  break;
647  }
648  return;
649 }
650 
651 PIXTYPE *SE2BackgroundModeller::getWhtMeanVals()
652 {
653  return itsWhtMeanVals;
654 }
655 } // end of namespace SourceXtractor
WEIGHT_FIELD
#define WEIGHT_FIELD
Definition: BackgroundDefine.h:41
SourceXtractor::PIXTYPE
float PIXTYPE
Definition: BackgroundDefine.h:30
std::shared_ptr
STL class.
std::fabs
T fabs(T... args)
RMS_FIELD
#define RMS_FIELD
Definition: BackgroundDefine.h:39
SIZETSUB
#define SIZETSUB(X, Y)
Definition: SE2BackgroundModeller.cpp:36
std::div
T div(T... args)
SourceXtractor::Image< unsigned char >
BackgroundDefine.h
BACK_BUFSIZE
#define BACK_BUFSIZE
Definition: BackgroundDefine.h:50
Utils.h
SourceXtractor
Definition: Aperture.h:30
std::ratio
Exception.h
Elements::Exception
VAR_FIELD
#define VAR_FIELD
Definition: BackgroundDefine.h:40
Elements::Logging::debug
void debug(const std::string &logMessage)
SourceXtractor::BackgroundCell
Definition: BackgroundCell.h:32
e
constexpr double e
SourceXtractor::BackgroundCell::getBackgroundValues
void getBackgroundValues(PIXTYPE &meanVal, PIXTYPE &sigmaVal)
Definition: BackgroundCell.cpp:64
BIG
#define BIG
Definition: BackgroundDefine.h:32
SourceXtractor::bck_model_logger
static Elements::Logging bck_model_logger
Definition: Utils.h:25
BackgroundCell.h
x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
Definition: MoffatModelFittingTask.cpp:93
std
STL namespace.
SE2BackgroundUtils.h
SE2BackgroundModeller.h
std::size_t
TypedSplineModelWrapper.h
y
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
Definition: MoffatModelFittingTask.cpp:93
SourceXtractor::TypedSplineModelWrapper
Definition: TypedSplineModelWrapper.h:36