35 #define SIZETSUB(X, Y) ((X) > (Y) ? (X-Y) : (Y-X))
43 itsVariance = variance_map;
45 itsMaskType = mask_type_flag;
49 if (variance_map!=
nullptr)
61 itsNaxes[0] = (
size_t)image->getWidth();
62 itsNaxes[1] = (
size_t)image->getHeight();
65 SE2BackgroundModeller::~SE2BackgroundModeller() {
67 delete[] itsWhtMeanVals;
72 size_t gridSize[2] = {0,0};
75 long increment[2]={1,1};
81 size_t subImgNaxes[2] = {0,0};
107 divResult =
std::div(
long(itsNaxes[0]),
long(bckCellSize[0]));
108 gridSize[0] =
size_t(divResult.quot);
113 divResult =
std::div(
long(itsNaxes[1]),
long(bckCellSize[1]));
114 gridSize[1] =
size_t(divResult.quot);
119 nGridPoints = gridSize[0]*gridSize[1];
122 bckMeanVals =
new PIXTYPE[nGridPoints];
123 bckSigVals =
new PIXTYPE[nGridPoints];
125 itsWhtMeanVals =
new PIXTYPE[nGridPoints];
126 whtSigVals =
new PIXTYPE[nGridPoints];
136 for (
size_t yIndex=0; yIndex<gridSize[1]; yIndex++){
139 fpixel[1] = (long)yIndex*bckCellSize[1];
140 lpixel[1] = yIndex < gridSize[1]-1 ? (long)(yIndex+1)*bckCellSize[1] : (long)itsNaxes[1];
143 for (
size_t xIndex=0; xIndex<gridSize[0]; xIndex++){
146 fpixel[0] = (long)xIndex*bckCellSize[0];
147 lpixel[0] = xIndex < gridSize[0]-1 ? (long)(xIndex+1)*bckCellSize[0] : (long)itsNaxes[0];
150 subImgNaxes[0] =(
size_t)(lpixel[0]-fpixel[0]);
151 subImgNaxes[1] =(
size_t)(lpixel[1]-fpixel[1]);
156 getMinIncr(nElements, increment, subImgNaxes);
159 if (nElements!=nData) {
161 gridData =
new PIXTYPE[nElements];
164 delete [] weightData;
165 weightData =
new PIXTYPE[nElements];
172 for (
auto yIndex=fpixel[1]; yIndex<lpixel[1]; yIndex+=increment[1])
173 for (
auto xIndex=fpixel[0]; xIndex<lpixel[0]; xIndex+=increment[0])
174 gridData[pixIndex++] = (
PIXTYPE)itsImage->getValue(
int(xIndex),
int(yIndex));
178 for (
auto yIndex=fpixel[1]; yIndex<lpixel[1]; yIndex+=increment[1])
179 for (
auto xIndex=fpixel[0]; xIndex<lpixel[0]; xIndex+=increment[0])
180 weightData[pixIndex++] = (
PIXTYPE)itsVariance->getValue(
int(xIndex),
int(yIndex));
185 for (
auto yIndex=fpixel[1]; yIndex<lpixel[1]; yIndex+=increment[1])
186 for (
auto xIndex=fpixel[0]; xIndex<lpixel[0]; xIndex+=increment[0], pixIndex++)
188 if (itsMask->getValue(
int(xIndex),
int(yIndex)) & itsMaskType){
189 gridData[pixIndex] = -
BIG;
196 oneCell =
new BackgroundCell(gridData, nElements, weightData, weightVarThreshold);
198 oneCell->
getBackgroundValues(bckMeanVals[gridIndex], bckSigVals[gridIndex], itsWhtMeanVals[gridIndex], whtSigVals[gridIndex]);
209 filter(bckMeanVals, bckSigVals, gridSize, filterBoxSize, filterThreshold);
211 filter(itsWhtMeanVals, whtSigVals, gridSize, filterBoxSize, filterThreshold);
217 computeScalingFactor(itsWhtMeanVals, bckSigVals, sigFac, nGridPoints);
224 for (
size_t index=0; index<nGridPoints; index++)
225 bckSigVals[index] *= bckSigVals[index];
232 delete [] whtSigVals;
234 delete [] weightData;
237 void SE2BackgroundModeller::getMinIncr(
size_t &nElements,
long* incr,
size_t * subImgNaxes)
243 nElements = subImgNaxes[0]*subImgNaxes[1];
255 axisRatio = float(subImgNaxes[0]) / float(subImgNaxes[1]);
264 if (axisRatio >= 1.0)
270 divResult =
std::div(
long(subImgNaxes[0]),
long(incr[0]));
271 subImgNaxes[0] =
size_t(divResult.quot);
276 divResult =
std::div(
long(subImgNaxes[1]),
long(incr[1]));
277 subImgNaxes[1] =
size_t(divResult.quot);
283 nElements = subImgNaxes[0]*subImgNaxes[1];
284 axisRatio = float(subImgNaxes[0]) / float(subImgNaxes[1]);
339 void SE2BackgroundModeller::filter(
PIXTYPE* bckVals,
PIXTYPE* sigmaVals,
const size_t* gridSize,
const size_t* filterSize,
const float &filterThreshold)
342 replaceUNDEF(bckVals, sigmaVals, gridSize);
345 filterMedian(bckVals, sigmaVals, gridSize, filterSize, filterThreshold);
350 void SE2BackgroundModeller::replaceUNDEF(
PIXTYPE* bckVals,
PIXTYPE* sigmaVals,
const size_t* gridSize)
359 size_t iAct,jAct, nx,ny, nmin;
366 np = gridSize[0]*gridSize[1];
369 for (
size_t index=0; index< np; index++)
371 if (bckVals[index]<=-
BIG)
398 for (
size_t py=0; py<ny; py++)
400 for (
size_t px=0; px<nx; px++)
406 backMod[iAct]=back[iAct];
409 if (back[iAct]<=-
BIG)
416 for (
size_t y=0;
y<ny;
y++)
418 for (
size_t x=0;
x<nx;
x++)
437 else if (
fabs(dist-distMin)<1
e-05)
449 backMod[iAct] = nmin ? val/nmin: 0.0;
450 sigma[iAct] = nmin ? sval/nmin: 1.0;
456 for (
size_t index=0; index< np; index++)
458 bckVals[index] =backMod[index];
468 void SE2BackgroundModeller::filterMedian(
PIXTYPE* bckVals,
PIXTYPE* sigmaVals,
const size_t* gridSize,
const size_t* filterSize,
const float filterThresh)
470 PIXTYPE *back, *sigma, *sigmat;
476 int i,nx,ny,npx,npx2,npy,npy2,
x,
y;
481 if (filterSize[0]<2 && filterSize[1]<2)
488 nx = (int)gridSize[0];
489 ny = (int)gridSize[1];
493 npx = (int)filterSize[0]/2;
494 npy = (int)filterSize[1]/2;
498 bmask =
new PIXTYPE[(2*npx+1)*(2*npy+1)];
499 smask =
new PIXTYPE[(2*npx+1)*(2*npy+1)];
510 for (
int py=0; py<np; py+=nx)
520 for (
int px=0; px<nx; px++)
532 for (
int dpy = -npy2; dpy<=npy2; dpy+=nx)
535 for (
int dpx = -npx2; dpx <= npx2; dpx++)
538 bmask[i] = back[
x+
y];
539 smask[i++] = sigma[
x+
y];
545 median = SE2BackgroundUtils::fqMedian(bmask, i);
546 if (
fabs((median-back[px+py]))>=(
PIXTYPE)filterThresh)
549 backFilt[px+py] = median;
550 sigmaFilt[px+py] = SE2BackgroundUtils::fqMedian(smask, i);
555 backFilt[px+py] = back[px+py];
556 sigmaFilt[px+py] = sigma[px+py];
562 for (
int index=0; index<np; index++)
563 back[index] = backFilt[index];
566 for (
int index=0; index<np; index++)
567 sigma[index] = sigmaFilt[index];
572 allSigmaMed = SE2BackgroundUtils::fqMedian(sigmaFilt, np);
576 if (allSigmaMed<=0.0)
578 sigmat = sigmaFilt+np;
579 for (i=np; i-- && *(--sigmat)>0.0;);
580 if (i>=0 && i<(np-1))
581 allSigmaMed = SE2BackgroundUtils::fqMedian(sigmat+1, np-1-i);
601 void SE2BackgroundModeller::computeScalingFactor(
PIXTYPE* whtMeanVals,
PIXTYPE* bckSigVals,
PIXTYPE& sigFac,
const size_t nGridPoints)
617 for (
size_t index=0; index<nGridPoints; index++){
618 if (whtMeanVals[index]>0.0){
620 actRatio = bckSigVals[index] * bckSigVals[index] / whtMeanVals[index];
629 sigFac = SE2BackgroundUtils::fqMedian(
ratio, nr);
632 for (lowIndex=0; lowIndex<nr &&
ratio[lowIndex]<=0.0; lowIndex++);
642 sigFac = SE2BackgroundUtils::fqMedian(
ratio+lowIndex, nr-lowIndex);
646 bck_model_logger.
debug() <<
"\tNull or negative global weighting factor: " <<
" | " << lowIndex <<
"defaulted to 1 " << nr;
655 void SE2BackgroundModeller::rescaleThreshold(
PIXTYPE &weightVarThreshold,
const PIXTYPE &weightThreshold)
658 if (weightThreshold<0.0){
659 throw Elements::Exception() <<
"The weight threshold is: " << weightThreshold <<
" but can not be smaller than 0.0";
663 switch (itsWeightTypeFlag){
667 weightVarThreshold = weightThreshold;
672 weightVarThreshold = weightThreshold*weightThreshold;
677 if (weightThreshold>0.0){
678 weightVarThreshold = 1.0/weightThreshold;
681 weightVarThreshold =
BIG;
687 weightVarThreshold =
BIG;
693 PIXTYPE *SE2BackgroundModeller::getWhtMeanVals()
695 return itsWhtMeanVals;