한동안 만들지 않았던 장비를 이번에 새로 2대나 만들게 되었는데 예전에 사용했던 Gaussian Fitting 함수를 찾지 못해서 애를 먹다가 회사서버 하드디스크 구석에 쳐박혀 숨어있던 코드를 찾아 내었습니다.
필요하신 분들이 계실것 같아서 올려 봅니다.
사용법은
double DataX[DATA_CNT], DataY[DATA_CNT], gaussY[DATA_CNT] ;
GaussDeviation(DataY, gaussY, DATA_CNT, peak_value, peak_pos, pdiam, pdev, 0, 0.5) ;
이런식으로 호출하시면 됩니다.
물론 DataX[DATA_CNT], DataY[DATA_CNT] 에 값이 들어 있어야 겠죠..
함수를 호출하면 gaussY 버퍼로 fitting 된 값들이 들어갑니다.
이걸 그래프로 그려 보면 아무리 노이즈가 심한 데이타도 깨끗하게 그려집니다.
#include <algorithm>
using namespace std;
int __fastcall GaussDeviation(double *buf, double *gbuf, int len,
double &peak, double &pkpos, double &pdiam, double &pdev, double base, double cut_level)
{
int i, p1, p2, gbuf_len ;
double pix=0, maxi=0 ;
double sumix = 0, sumi = 0 ; //long type에서 double로 변경, sumi가 0이 되어 gaussianfit하지 않음.
double thresh,thresh1, maxinten ;
double radius=0, center, r2=0., oldr2=0 ;
double *temp ;
temp = max_element(&buf[0], &buf[len]) ;
maxi = *temp ;
//------------ 추가 코드 --------------------
temp = min_element(&buf[0], &buf[len]) ;
base = *temp ;
//-------------------------------------------
maxinten = maxi - base ;
thresh = maxinten * cut_level ;
peak = 0 ;
pkpos = 0 ;
pdiam = 0 ;
pdev = 0 ;
if(len == 0) return 0 ;
for(i=1 ; i<len-1 ; i++)
{
if(buf[i] < thresh) continue ;
pix = buf[i] - thresh ;
sumi += pix ;
sumix += i * pix ;
}
for(p1=1 ; buf[p1]<=thresh && p1<len-1 ; p1++) ; // find first edge
for(p2=len-1 ; buf[p2]<=thresh && p2>p1 ; p2--) ; // find second edge
if(fabs(buf[p1] - thresh) > fabs(buf[p1+1] - thresh)) p1 = p1 + 1 ;
if(fabs(buf[p2] - thresh) > fabs(buf[p2-1] - thresh)) p2 = p2 - 1 ;
center = p1 + (p2 - p1) / 2.0 ;
radius = (p2 - p1 + 1) / 2.0 ;
memset(gbuf, 0, sizeof(double)*len) ;
if(sumi < 1) // If no points, just give center
{
sumix = len / 2 ;
sumi = 1 ;
radius = 0 ;
}
if(sumi >= 1 && radius > 1)
{
double a=0, b=0 ;
if(cut_level > 0)
{
for(i=0 ; i<50 ; i++)
{
GaussianFit(buf,gbuf,len,base,maxinten,radius,center,&a,&b,&r2) ;
maxinten += a ;
radius += b ;
if(r2 < 0.) r2=0. ; // On horrible curves, r2 can go negative!
if(r2<0.005|| (i>1 && a<0.001 && b<0.001)) break ;
if(i>1 && fabs(oldr2-r2)/r2 < .0002) break ; // r2 didn't change
oldr2 = r2 ;
}
}
else
{
for(i=0 ; i<5 ; i++)
{
GaussianFit(buf,gbuf,len,base,maxinten,radius,center,&a,&b,&r2) ;
maxinten += a ;
radius += b ;
if(r2<0.) r2=0. ; // On horrible curves, r2 can go negative!
if(r2<0.005 || (i>1 && a<0.001 && b<0.001)) break ;
if(i>1 && fabs(oldr2-r2)/r2 < .0002) break ; // r2 didn't change
oldr2 = r2 ;
}
}
}
// radius = (p2-p1) / 2. ; // Use the image's radius, not the gaussian's
pdiam = p1 ; //radius*2; 0808..
pdev = r2 ;
peak = maxinten ;
for(i=0 ; i<len ; i++)
{
if(gbuf[i] >= maxinten) pkpos = i ;
}
gbuf_len = p2 - p1 + 1 ;
return gbuf_len ;
}
//---------------------------------------------------------------------------
void __fastcall GaussianFit(double *pix, double *gbuf,int npix, double mininten,
double A, double B, double T, double *a, double *b, double *r2)
{
/* pix is the row of pixel data
gbuf is output of fit gaussian
npix is the number of pixels in the row
A is the current amplitude estimate (0 < A < 256)
B is the current width estimate (in pixels)
T is the centroid (assumed to be accurate/constant; 0 < T < npix)
a is the calculated correction for A
b is the calculated correction for B
r2 is the r?goodness-of-fit for the revised estimators
*/
/* note that it would be intuitive for this function to return r2 */
/* SUMMATION VARIABLES */
/* for each pixel : */
double x; /* the distance of the pixel from the centroid = i-T */
double u; /* the partial derivative of the Gaussian estimate w.r.t. A */
double v; /* the partial derivative of the Gaussian estimate w.r.t. B */
double f; /* the current Gaussian estimate */
double F; /* the deviation of the pixel value from the estimate */
/* overall: */
double C=0; /* the sum of u?*/
double D=0; /* the sum of uv */
double E=0; /* the sum of Fu */
double G=0; /* the sum of v?*/
double H=0; /* the sum of Fv */
double J=0; /* the sum of F?*/
/* additional variables from Coherents formulas (without names) */
double sumf=0; /* the sum of revised estimates */
double sumF=0; /* the sum of errors in the revised estimates */
double avgf; /* the average new estimate */
double var; /* the variance between a pixel value and the average */
double sumvar2=0; /* the sum of the squared variances */
int i ;
for(i=0 ; i<npix ; i++) /* this is the slowest part of this function */
{
x = (double)(i-(int)T) ;
f = A * exp(-2.*((x*x)/(B*B))) ; /* (1) */
F = (double)(pix[i]-mininten)-f ; /* (5) */
if(A==0) A = 1 ;
u = f / A ; /* (3) */
if(B==0) B = 1 ;
v = 4.* f *(x*x)/(B*B*B) ; /* (4) */
C += u*u ;
D += u*v ;
E += F*u ;
G += v*v ;
H += F*v ;
J += F*F ;
}
/* new coefficient estimates */
double temp ;
temp = (D*D-C*G) ;
if(temp==0) temp = 1 ;
*b = (E*D-C*H)/ temp ; /* (8) */
*a = -(E*G-D*H)/temp ; /* (9) */
A += *a;
B += *b;
/* analyze estimates using new coefficients */
J = 0 ;
for(i=0 ; i<npix ; i++)
{
x = (double)(i-(int)T) ;
if(B==0) B = 1 ;
f = A*exp(-2.*((x*x)/(B*B))) ;
F = (double)(pix[i]-mininten) - f ;
gbuf[i] = mininten + f ; /* PLOT IDEAL GAUSSIAN */
sumf += f ;
sumF += F ;
J += F*F ;
}
if(npix == 0) npix = 1 ;
avgf = sumf / npix ;
for(i=0 ; i<npix ; i++)
{
var = (double)(pix[i]-mininten) - avgf ;
sumvar2 += var * var ;
}
if(sumvar2==0) sumvar2 = 1 ;
*r2 = 1. - sqrt(J / sumvar2) ;
}
//---------------------------------------------------------------------------
