SourceXtractorPlusPlus  0.15
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BackgroundHistogram.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 #include <cmath>
26 #include <cstddef>
29 
30 namespace SourceXtractor {
31 
32 BackgroundHistogram::BackgroundHistogram(const double& mean, const double& sigma, const size_t& ndata)
33 {
34  double step;
35 
36  // Astromatic, back.c:back_stat(), falls from the sky!
37  step = sqrt(2.0/M_PI)*QUANTIF_NSIGMA/QUANTIF_AMIN;
38 
39  // set the basic statistic info
40  itsMean = mean;
41  itsSigma = sigma;
42  itsStatNData = ndata;
43 
44  // find the right number of levels
45  itsNLevels = (size_t)(step*itsStatNData+1);
48 
49  // set the scale
50  if (itsSigma>0.0)
51  itsQscale = 2.0*QUANTIF_NSIGMA*(float)itsSigma/(float)itsNLevels;
52  else
53  itsQscale = 1.0;
54 
55  // set the zero level
56  itsQzero = (float)mean - QUANTIF_NSIGMA*(float)sigma;
57 
58  // compute the constant
59  itsCste = 0.499999 - itsQzero/itsQscale;
60 
61  // allocate the histogram
62  itsHisto = new int[itsNLevels];
63  for (size_t index=0; index<itsNLevels; index++)
64  itsHisto[index]=0;
65 }
66 
68  // delete the histogram
69  if (itsHisto)
70  delete[] itsHisto;
71  itsHisto=NULL;
72 }
73 /*
74 void BackgroundHistogram::fillInData(const PIXTYPE* cellData, const size_t ndata)
75 {
76  int histIndex;
77  PIXTYPE pixVal;
78 
79  // go over the data
80  for (size_t index=0; index<ndata; index++)
81  {
82  // take the current value
83  pixVal =cellData[index];
84 
85  // compute the index in the histogram and enhance
86  // the bin if possible
87  histIndex = (int)(pixVal/(PIXTYPE)itsQscale + (PIXTYPE)itsCste);
88  if (histIndex>=0 && histIndex<(int)itsNLevels)
89  itsHisto[histIndex]++;
90  }
91  //logger.info() << "Created the histogram!";
92  return;
93 }
94 */
96 {
97  int histIndex;
98 
99  // compute the index in the histogram and enhance
100  // the bin if possible
101  histIndex = (int)(pixVal/(PIXTYPE)itsQscale + (PIXTYPE)itsCste);
102  if (histIndex>=0 && histIndex<(int)itsNLevels)
103  itsHisto[histIndex]++;
104 
105  return;
106 }
107 
109 {
110  int *histo, *hilow, *hihigh, *histot;
111  unsigned long lowsum, highsum, sum;
112  double ftemp, mea, sig, sig1, med, dpix;
113  int i, n, lcut,hcut, nlevelsm1, pix;
114 
115  int hilowIndex, hihighIndex;
116  //int hcutIndex, lcutIndex;
117  //int hilowVal, hihighVal;
118  unsigned long mySum;
119  double myMea, mySig, myDpix;
120  //int myI, myN, myLcut,myHcut, myNlevelsm1, myPix;
121  int myPix;
122 
123 
124  histo = itsHisto;
125  hcut = itsNLevels-1;
126  lcut = 0;
127  nlevelsm1 = itsNLevels-1;
128  lcut = 0;
129 
130  sig = 10.0*nlevelsm1;
131  sig1 = 1.0;
132  mea = itsMean;
133  med = itsMean;
134  for (n=100; n-- && (sig>=0.1) && (fabs(sig/sig1-1.0)>BACK_EPS);)
135  {
136  sig1 = sig;
137  sum = 0.0;
138  mea = 0.0;
139  sig = 0.0;
140  lowsum = 0;
141  highsum = 0;
142 
143  mySum=0.0;
144  myMea=0.0;
145  mySig=0.0;
146 
147  histot = histo+lcut;
148  hilow = histo+lcut;
149  hihigh = histo+hcut;
150 
151  // save the values and not the pointers
152  //histotIndex=lcut;
153  hilowIndex=lcut;
154  hihighIndex=hcut;
155  //histotVal = histo[histotIndex];
156  //hilowVal = histo[hilowIndex];
157  //hihighVal = histo[hihighIndex];
158 
159  for (i=lcut; i<=hcut; i++)
160  {
161  if (lowsum<highsum)
162  {
163  lowsum += *(hilow++);
164  hilowIndex++;
165  //hilowVal = histo[hilowIndex];
166  }
167  else
168  {
169  highsum += *(hihigh--);
170  hihighIndex--;
171  //hihighVal = histo[hihighIndex];
172  }
173 
174  sum += (pix = *(histot++));
175  mea += (dpix = (double)pix*i);
176  sig += dpix*i;
177 
178  myPix = histo[i];
179  mySum += myPix;
180  myDpix = (double)myPix*i;
181  myMea += myDpix;
182  mySig += myDpix*i;
183  }
184 
185 
186  med = hihigh>=histo?((hihigh-histo)+0.5+((double)highsum-lowsum)/(2.0*(*hilow>*hihigh?*hilow:*hihigh))): 0.0;
187  /*if (hihigh>=histo){
188  printf("%d %d %i %i\n", *hihigh, *hilow, hihighIndex, hihighVal);
189  int iTmp;
190  if (*hilow>*hihigh)
191  iTmp = *hilow;
192  else
193  iTmp = *hihigh;
194  med = (hihigh-histo)+0.5+((double)highsum-lowsum)/(2.0*iTmp);
195  }
196  else {
197  med=0.0;
198  }
199  */
200  if (sum)
201  {
202  mea /= (double)sum;
203  sig = sig/sum - mea*mea;
204  }
205  sig = sig>0.0?sqrt(sig):0.0;
206  lcut = (ftemp=med-3.0*sig)>0.0 ?(int)(ftemp>0.0?ftemp+0.5:ftemp-0.5):0;
207  hcut = (ftemp=med+3.0*sig)<nlevelsm1 ?(int)(ftemp>0.0?ftemp+0.5:ftemp-0.5): nlevelsm1;
208  //printf("cuts: %i <-> %i\n", lcut, hcut);
209  }
210 
211 
212  bckVal = fabs(sig)>0.0? (fabs(itsSigma/(sig*itsQscale)-1) < 0.0 ?
213  itsQzero+mea*itsQscale
214  :(fabs((mea-med)/sig)< 0.3 ?
215  itsQzero+(2.5*med-1.5*mea)*itsQscale
216  :itsQzero+med*itsQscale))
217  :itsQzero+mea*itsQscale;
218  sigmaVal = sig*itsQscale;
219 
220  return;
221 }
222 
224 {
225  // code taken from Astromatic back.c:back_guess();
226  // small adjustments for the class environment;
227  // TODO: re-factor that such that it is understandable!!
228  int *histo, *hilow, *hihigh, *histot;
229  unsigned long lowsum, highsum, sum;
230  double ftemp, mea, sig, sig1, med, dpix;
231  int i, n, lcut,hcut, nlevelsm1, pix;
232 
233 
234  histo = itsHisto;
235  hcut = itsNLevels-1;
236  nlevelsm1 = itsNLevels-1;
237  lcut = 0;
238 
239  sig = 10.0*nlevelsm1;
240  sig1 = 1.0;
241  mea = itsMean;
242  med = itsMean;
243 
244  // iterates the clipping until one of the following conditions is fullfilled:
245  // - sig < 0.1 (meaning that the histogram of the pixels sharply peaks to the mean i.e. all the pixels are in the same histogram bin)
246  // - fabs(sig/sig1-1) <= 1e-4 (meaning that the difference with respect to the previous iteration is small)
247 
248  for (n=100; n-- && (sig>=0.1) && (fabs(sig/sig1-1.0)>BACK_EPS);)
249  {
250  sig1 = sig;
251  sum = mea = sig = 0.0;
252  lowsum = highsum = 0;
253  histot = hilow = histo+lcut;
254  hihigh = histo+hcut;
255  for (i=lcut; i<=hcut; i++)
256  {
257  if (lowsum<highsum)
258  lowsum += *(hilow++);
259  else
260  highsum += *(hihigh--);
261  sum += (pix = *(histot++));
262  mea += (dpix = (double)pix*i);
263  sig += dpix*i;
264  }
265 
266  med = hihigh>=histo?
267  ((hihigh-histo)+0.5+((double)highsum-lowsum)/(2.0*(*hilow>*hihigh?
268  *hilow:*hihigh)))
269  : 0.0;
270  if (sum)
271  {
272  mea /= (double)sum;
273  sig = sig/sum - mea*mea;
274  }
275  sig = sig>0.0?sqrt(sig):0.0;
276  lcut = (ftemp=med-3.0*sig)>0.0 ?(int)(ftemp>0.0?ftemp+0.5:ftemp-0.5):0;
277  hcut = (ftemp=med+3.0*sig)<nlevelsm1 ?(int)(ftemp>0.0?ftemp+0.5:ftemp-0.5)
278  : nlevelsm1;
279  }
280  // The following lines can be translated into:
281  //
282  // if (fabs(sig>0.0)) <=== "If sigma is not zero"
283  // {
284  // if (fabs(itsSigma/(sig*itsQscale)-1) < 0.0) <=== I think this condition is never satisfied
285  // {
286  // bckVal=itsQzero+mea*itsQscale;
287  // }
288  // else
289  // if (fabs((mea-med)/sig) < 0.3) <=== This is a measure of how many outliers are in the distribution.
290  // { If mean and median do not differ too much, we are in the NOT crowded case
291  // bckVal=itsQzero+(2.5*med-1.5*mea)*itsQscale;
292  // }
293  // else <=== This is the crowded case
294  // {
295  // bckVal=itsQzero+med*itsQscale;
296  // }
297  // }
298  // else bckVal=itsQzero+mea*itsQscale; <=== "If sigma is zero"
299  //
300  bckVal = fabs(sig)>0.0? (fabs(itsSigma/(sig*itsQscale)-1) < 0.0 ?
301  itsQzero+mea*itsQscale
302  :(fabs((mea-med)/sig)< 0.3 ?
303  itsQzero+(2.5*med-1.5*mea)*itsQscale
304  :itsQzero+med*itsQscale))
305  :itsQzero+mea*itsQscale;
306 
307  // In the following lines I implement the previous conditions as they are descripted
308  // in the sextractor user's manual (E. Bertin, SExtractor v2.13, pag. 16)
309  //
310  // if (fabs(sig>0.0)) <=== "If sigma is not zero"
311  // {
312  // if (fabs(sig/sig1-1.) < 0.2) <=== This is a measure of how many outliers are in the distribution.
313  // { If mean and median do not differ too much, we are in the NOT crowded case
314  // bckVal=itsQzero+mea*itsQscale;
315  // }
316  // else <=== This is the crowded case
317  // {
318  // bckVal=itsQzero+(2.5*med-1.5*mea)*itsQscale;
319  // }
320  // }
321  // else bckVal=itsQzero+mea*itsQscale; <=== "If sigma is zero"
322  //
323 
324 
325  sigmaVal = sig*itsQscale;
326 
327  return;
328 }
329 } // end of namespace SourceXtractor
330 
#define BACK_EPS
BackgroundHistogram(const double &mean, const double &sigm, const size_t &ndata)
#define QUANTIF_AMIN
#define QUANTIF_NSIGMA
T fabs(T...args)
void getBackGuessMod(PIXTYPE &bckVal, PIXTYPE &sigmaVal)
void getBackGuess(PIXTYPE &bckVal, PIXTYPE &sigmaVal)
void addDatum(const PIXTYPE &pixVal)
T sqrt(T...args)
#define QUANTIF_NMAXLEVELS