37 #define SIZETSUB(X, Y) ((X) > (Y) ? (X-Y) : (Y-X))
40 namespace SourceXtractor {
45 itsVariance = variance_map;
47 itsMaskType = mask_type_flag;
51 if (variance_map!=
nullptr)
63 itsNaxes[0] = (
size_t)image->getWidth();
64 itsNaxes[1] = (
size_t)image->getHeight();
67 SE2BackgroundModeller::~SE2BackgroundModeller() {
69 delete[] itsWhtMeanVals;
77 DetectionAccessor detectionAccessor(itsImage);
78 MaskAccessor maskAccessor(itsMask);
80 size_t gridSize[2] = {0,0};
83 long increment[2]={1,1};
89 size_t subImgNaxes[2] = {0,0};
115 divResult =
std::div(
long(itsNaxes[0]),
long(bckCellSize[0]));
116 gridSize[0] =
size_t(divResult.quot);
121 divResult =
std::div(
long(itsNaxes[1]),
long(bckCellSize[1]));
122 gridSize[1] =
size_t(divResult.quot);
127 nGridPoints = gridSize[0]*gridSize[1];
130 bckMeanVals =
new PIXTYPE[nGridPoints];
131 bckSigVals =
new PIXTYPE[nGridPoints];
133 itsWhtMeanVals =
new PIXTYPE[nGridPoints];
134 whtSigVals =
new PIXTYPE[nGridPoints];
144 for (
size_t yIndex=0; yIndex<gridSize[1]; yIndex++){
147 fpixel[1] = (long)yIndex*bckCellSize[1];
148 lpixel[1] = yIndex < gridSize[1]-1 ? (long)(yIndex+1)*bckCellSize[1] : (long)itsNaxes[1];
151 for (
size_t xIndex=0; xIndex<gridSize[0]; xIndex++){
154 fpixel[0] = (long)xIndex*bckCellSize[0];
155 lpixel[0] = xIndex < gridSize[0]-1 ? (long)(xIndex+1)*bckCellSize[0] : (long)itsNaxes[0];
158 subImgNaxes[0] =(
size_t)(lpixel[0]-fpixel[0]);
159 subImgNaxes[1] =(
size_t)(lpixel[1]-fpixel[1]);
162 bck_model_logger.
debug() <<
"Background cell from fpixel=(" << fpixel[0] <<
"," << fpixel[1] <<
") to lpixel=("<< lpixel[0] <<
"," << lpixel[1] <<
")";
167 getMinIncr(nElements, increment, subImgNaxes);
170 if (nElements!=nData) {
172 gridData =
new PIXTYPE[nElements];
175 delete [] weightData;
176 weightData =
new PIXTYPE[nElements];
183 for (
auto yIndex=fpixel[1]; yIndex<lpixel[1]; yIndex+=increment[1])
184 for (
auto xIndex=fpixel[0]; xIndex<lpixel[0]; xIndex+=increment[0]){
185 gridData[pixIndex++] = (
PIXTYPE)detectionAccessor.getValue(
int(xIndex),
int(yIndex));
189 for (
auto yIndex=fpixel[1]; yIndex<lpixel[1]; yIndex+=increment[1])
190 for (
auto xIndex=fpixel[0]; xIndex<lpixel[0]; xIndex+=increment[0])
191 weightData[pixIndex++] = (
PIXTYPE)detectionAccessor.getValue(
int(xIndex),
int(yIndex));
196 for (
auto yIndex=fpixel[1]; yIndex<lpixel[1]; yIndex+=increment[1])
197 for (
auto xIndex=fpixel[0]; xIndex<lpixel[0]; xIndex+=increment[0], pixIndex++)
199 if (maskAccessor.getValue(
int(xIndex),
int(yIndex)) & itsMaskType){
200 gridData[pixIndex] = -
BIG;
207 oneCell =
new BackgroundCell(gridData, nElements, weightData, weightVarThreshold);
209 oneCell->
getBackgroundValues(bckMeanVals[gridIndex], bckSigVals[gridIndex], itsWhtMeanVals[gridIndex], whtSigVals[gridIndex]);
220 filter(bckMeanVals, bckSigVals, gridSize, filterBoxSize, filterThreshold);
222 filter(itsWhtMeanVals, whtSigVals, gridSize, filterBoxSize, filterThreshold);
228 computeScalingFactor(itsWhtMeanVals, bckSigVals, sigFac, nGridPoints);
235 for (
size_t index=0; index<nGridPoints; index++)
236 bckSigVals[index] *= bckSigVals[index];
243 delete [] whtSigVals;
245 delete [] weightData;
248 void SE2BackgroundModeller::getMinIncr(
size_t &nElements,
long* incr,
const size_t * subImgNaxes)
254 size_t tmpImgNaxes[2] = {subImgNaxes[0], subImgNaxes[1]};
262 nElements = tmpImgNaxes[0]*tmpImgNaxes[1];
267 axisRatio = float(tmpImgNaxes[0]) / float(tmpImgNaxes[1]);
275 if (axisRatio >= 1.0)
281 divResult =
std::div(
long(subImgNaxes[0]),
long(incr[0]));
282 tmpImgNaxes[0] =
size_t(divResult.quot);
287 divResult =
std::div(
long(subImgNaxes[1]),
long(incr[1]));
288 tmpImgNaxes[1] =
size_t(divResult.quot);
294 nElements = tmpImgNaxes[0]*tmpImgNaxes[1];
295 axisRatio = float(tmpImgNaxes[0]) / float(tmpImgNaxes[1]);
299 bck_model_logger.
debug() <<
"New increment=(" << incr[0] <<
"," << incr[1] <<
") sampledPixels=("<< tmpImgNaxes[0] <<
"," << tmpImgNaxes[1] <<
") nElements=" << nElements;
304 void SE2BackgroundModeller::filter(
PIXTYPE* bckVals,
PIXTYPE* sigmaVals,
const size_t* gridSize,
const size_t* filterSize,
const float &filterThreshold)
307 replaceUNDEF(bckVals, sigmaVals, gridSize);
310 filterMedian(bckVals, sigmaVals, gridSize, filterSize, filterThreshold);
315 void SE2BackgroundModeller::replaceUNDEF(
PIXTYPE* bckVals,
PIXTYPE* sigmaVals,
const size_t* gridSize)
324 size_t iAct,jAct, nx,ny, nmin;
331 np = gridSize[0]*gridSize[1];
334 for (
size_t index=0; index< np; index++)
336 if (bckVals[index]<=-
BIG)
363 for (
size_t py=0; py<ny; py++)
365 for (
size_t px=0; px<nx; px++)
371 backMod[iAct]=back[iAct];
374 if (back[iAct]<=-
BIG)
381 for (
size_t y=0;
y<ny;
y++)
383 for (
size_t x=0;
x<nx;
x++)
402 else if (
fabs(dist-distMin)<1
e-05)
414 backMod[iAct] = nmin ? val/nmin: 0.0;
415 sigma[iAct] = nmin ? sval/nmin: 1.0;
421 for (
size_t index=0; index< np; index++)
423 bckVals[index] =backMod[index];
433 void SE2BackgroundModeller::filterMedian(
PIXTYPE* bckVals,
PIXTYPE* sigmaVals,
const size_t* gridSize,
const size_t* filterSize,
const float filterThresh)
435 PIXTYPE *back, *sigma, *sigmat;
441 int i,nx,ny,npx,npx2,npy,npy2,
x,
y;
446 if (filterSize[0]<2 && filterSize[1]<2)
453 nx = (int)gridSize[0];
454 ny = (int)gridSize[1];
458 npx = (int)filterSize[0]/2;
459 npy = (int)filterSize[1]/2;
463 bmask =
new PIXTYPE[(2*npx+1)*(2*npy+1)];
464 smask =
new PIXTYPE[(2*npx+1)*(2*npy+1)];
475 for (
int py=0; py<np; py+=nx)
485 for (
int px=0; px<nx; px++)
497 for (
int dpy = -npy2; dpy<=npy2; dpy+=nx)
500 for (
int dpx = -npx2; dpx <= npx2; dpx++)
503 bmask[i] = back[x+
y];
504 smask[i++] = sigma[x+
y];
510 median = SE2BackgroundUtils::fqMedian(bmask, i);
511 if (
fabs((median-back[px+py]))>=(
PIXTYPE)filterThresh)
514 backFilt[px+py] = median;
515 sigmaFilt[px+py] = SE2BackgroundUtils::fqMedian(smask, i);
520 backFilt[px+py] = back[px+py];
521 sigmaFilt[px+py] = sigma[px+py];
527 for (
int index=0; index<np; index++)
528 back[index] = backFilt[index];
531 for (
int index=0; index<np; index++)
532 sigma[index] = sigmaFilt[index];
537 allSigmaMed = SE2BackgroundUtils::fqMedian(sigmaFilt, np);
541 if (allSigmaMed<=0.0)
543 sigmat = sigmaFilt+np;
544 for (i=np; i-- && *(--sigmat)>0.0;);
545 if (i>=0 && i<(np-1))
546 allSigmaMed = SE2BackgroundUtils::fqMedian(sigmat+1, np-1-i);
566 void SE2BackgroundModeller::computeScalingFactor(
PIXTYPE* whtMeanVals,
PIXTYPE* bckSigVals,
PIXTYPE& sigFac,
const size_t nGridPoints)
577 ratio =
new PIXTYPE[nGridPoints];
582 for (
size_t index=0; index<nGridPoints; index++){
583 if (whtMeanVals[index]>0.0){
585 actRatio = bckSigVals[index] * bckSigVals[index] / whtMeanVals[index];
594 sigFac = SE2BackgroundUtils::fqMedian(ratio, nr);
597 for (lowIndex=0; lowIndex<nr && ratio[lowIndex]<=0.0; lowIndex++);
607 sigFac = SE2BackgroundUtils::fqMedian(ratio+lowIndex, nr-lowIndex);
611 bck_model_logger.
debug() <<
"\tNull or negative global weighting factor: " <<
" | " << lowIndex <<
"defaulted to 1 " << nr;
620 void SE2BackgroundModeller::rescaleThreshold(
PIXTYPE &weightVarThreshold,
const PIXTYPE &weightThreshold)
623 if (weightThreshold<0.0){
624 throw Elements::Exception() <<
"The weight threshold is: " << weightThreshold <<
" but can not be smaller than 0.0";
628 switch (itsWeightTypeFlag){
632 weightVarThreshold = weightThreshold;
637 weightVarThreshold = weightThreshold*weightThreshold;
642 if (weightThreshold>0.0){
643 weightVarThreshold = 1.0/weightThreshold;
646 weightVarThreshold =
BIG;
652 weightVarThreshold =
BIG;
658 PIXTYPE *SE2BackgroundModeller::getWhtMeanVals()
660 return itsWhtMeanVals;
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
void debug(const std::string &logMessage)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y