Gaussian fit C + + implementation

The Gaussian function is a mathematically important function, and the density function of our familiar normal distribution is the Gaussian function, also called the Gaussian distribution. The normal distribution is undoubtedly one of the most important distributions in probability theory and mathematical statistics.

The question now is how to find a Gaussian function to fit the point set if we give some point sets!

Of course, the fitting method is also the least squares, and the fitted function is:

Y=a*exp (-((x-b)/C) ^2);

A total of three parameters, A, B, C. However, this exponential function fitting is difficult to achieve, so the logarithmic transformation is used to two functions, as follows:

ln (y) =-(1/c^2) *x^2+ (2b/c^2) *x +ln (a)-(B/C) ^2;

In this way, another z=ln (y), a=-(1/c^2), b= (2b/c^2), C=ln (A)-(B/C) ^2; can be fitted with the least squares method,

Here I still use the singular value decomposition method to solve the least-squares solution, but before I do the point set to deal with, I did not fit all the points set, because this effect is very poor, after all, do a logarithmic transformation,

The impact of the flat point on the curve is too large, so I first judged the peak point, and then to the peak point to fit, it turns out my idea is correct, the fitting effect improved a lot!

VC implementation results are as follows:

The core program is implemented as follows:

GaussFit.h

/************************************************************************* Version: 2014-1-06 function Description: The Gaussian fitting of the least squares is given for some column points on the plane, and the least squares solution is solved by the singular value decomposition method as the Gaussian parameter. Call form: Gaussfit (Arrayx, Arrayy,int n,float box,miny); parameter description: Arrayx:arrayx[n], each value is x-axis a point arrayx:arrayy[n], each   The values are y-axis a point n: the number of points box:box[3], 3 parameters of the Gaussian function, respectively, is a,b,c; Miny:y in the direction of the translation, the actual fitting function is Y=a*exp (-((x-b)/C) ^2) +miny*************************************************************** /#pragma oncestruct gpoint{int x;int y; Gpoint (int x,int y): x (x), Y (y) {};friend bool operator < (gpoint p1,gpoint p2) {return p1.x>p2.x;}}; Class Gaussfit{public:gaussfit (void), ~gaussfit (void), void gaussfit (int *arrayx, int *arrayy,int n,float *box,int & Miny);p rivate:int SVD (float *a,int m,int n,float b[],float x[],float ESP); int gmiv (float a[],int m,int n,float b[],float    X[],float aa[],float eps,float u[],float v[],int ka); int Ginv (float a[],int m,int n,float aa[],float eps,float u[],float v[],int ka); int muav (float a[],inT m,int n,float u[],float v[],float eps,int ka);}; 

GaussFit.cpp

#include "StdAfx.h" #include "GaussFit.h" #include <cmath> #include <queue> #include <vector>using namespace Std; Gaussfit::gaussfit (void) {}gaussfit::~gaussfit (void) {}void gaussfit::gaussfit (int *arrayx, int *arrayy,int n,float * Box,int &miny) {float *a1=new float[n*3];float *b1=new float[n];float *pointx=new float[n];float *Pointy=new float[    N];int Maxy=0,midx,minx;miny=int_max;    Const double min_eps = 1e-10;            int i;priority_queue<gpoint> GP;                Used to sort Point.x priority_queue<int>py; Used to calculate the Point.yint m=n/10;m=m>0?m:1;for (i=0;i<n;i++) {gpoint tmp (arrayx[i],arrayy[i)), Gp.push (TMP); Py.size () <m) {Py.push (arrayy[i]);} else if (Py.top () >arrayy[i]) {py.pop ();p Y.push (Arrayy[i]);}}         Miny=py.top ();        Replace the smallest y with the small Y, to prevent the anomaly from Minx=gp.top (). x; for (i=0;i<n;i++) {gpoint tmp=gp.top (); Gp.pop (); pointx[i]= (Tmp.x-minx) *1.0; pointy[i]= (Tmp.y-miny) *1.0;if (pointy[i]>maxy) {maxy=pointy[i];midx=i;} else if (pointy[i]<0) {pointy[i]=0;}} float meany=0;for (int i=0;i<n;i++) {meany+=pointy[i];} meany/=n;//Statistics Peak vector<gpoint>vg,vtmp;for (int i=0;i<n;i++) {if (pointy[i]>meany) {gpoint tmp (Pointx[i], Pointy[i]); Vtmp.push_back (TMP);} Else{int s1=vg.size (), s2=vtmp.size (), if (S1&LT;S2) {vg.clear (); for (int j=0;j<vtmp.size (); j + +) {Gpoint TMP1 (vtmp[j ].X,VTMP[J].Y);    Vg.push_back (TMP1);}}   Vtmp.clear ();}} int s1=vg.size (), s2=vtmp.size (), if (S1&LT;S2) {vg.clear (); for (int j=0;j<vtmp.size (); j + +) {Gpoint tmp1 (vtmp[j].x, VTMP[J].Y); Vg.push_back (TMP1);}}  Gaussian fit int size=vg.size () on peak, if (size>0) {for (i = 0; i < n; i++) {int step= (i) *3;float px,py;px = Vg[i%size].x;py = VG[I%SIZE].Y; B1[i] = log (py); A1[step] = 1.0; A1[step + 1] = px; A1[step + 2] = px * PX;} float *x1=new float[3]; SVD (A1,n,3,b1,x1,min_eps), if (x1[2]<0) {box[2]=sqrt ( -1.0/x1[2]); Box[1]=x1[1]*box[2]*box[2]*0.5;box[0]=exp (x1[ 0]+box[1]*box[1]/(box[2]*box[2])); Box[1]+=minx;} Else{box[0]=box[1]=box[2]=-1;} delete []x1;} Else{box[0]=box[1]=box[2]=-1;} delete []a1;   delete []b1;delete []pointx;delete []pointy;}    int GAUSSFIT::SVD (float *a,int m,int n,float b[],float x[],float ESP) {float *aa;float *u;float *v;    Aa=new Float[n*m];u=new Float[m*m];      V=new Float[n*n];   int ka;   int flag;   if (m>n) {ka=m+1;   }else {ka=n+1;           } flag=gmiv (A,m,n,b,x,aa,esp,u,v,ka);    delete []aa;delete []u;delete []v; return (flag);}     int Gaussfit::gmiv (float a[],int m,int n,float b[],float x[],float aa[],float eps,float u[],float v[],int ka) {int i,j;    I=GINV (A,m,n,aa,eps,u,v,ka);    if (i<0) return (-1);        for (i=0; i<=n-1; i++) {x[i]=0.0;      for (j=0; j<=m-1; j + +) X[i]=x[i]+aa[i*m+j]*b[j];  } return (1); }int Gaussfit::ginv (float a[],int m,int n,float aa[],float eps,float u[],float v[],int ka) {//int muav (float a[],int    M,int n,float u[],float v[],float eps,int ka); int i,j,k,l,t,p,q,f;    I=MUAV (A,m,n,u,v,eps,ka);    if (i<0) return (-1);    J=n;    if (m<n) j=m;    J=j-1; k=0;   while ((K&LT;=J) && (a[k*n+k]!=0.0)) k=k+1;    K=k-1;        For (i=0, i<=n-1; i++) for (j=0; j<=m-1; J + +) {t=i*m+j; aa[t]=0.0;            for (l=0; l<=k; l++) {f=l*n+i; p=j*m+l; q=l*n+l;          AA[T]=AA[T]+V[F]*U[P]/A[Q];  }} return (1); }int Gaussfit::muav (float a[],int m,int n,float u[],float v[],float eps,int ka) {int i,j,k,l,it,ll,kk,ix,iy,mm,nn,iz,m1,    ks    float d,dd,t,sm,sm1,em1,sk,ek,b,c,shh,fg[2],cs[2];    float *s,*e,*w;   void PPP ();     void SSs ();     void PPP (float a[],float e[],float s[],float v[],int m,int N);    void sss (float fg[],float cs[]), s= (float *) malloc (ka*sizeof (float));    E= (float *) malloc (ka*sizeof (float));    w= (float *) malloc (ka*sizeof (float)); it=60;    K=n;    if (m-1<n) k=m-1;    L=m;    if (n-2<m) l=n-2;    if (l<0) l=0;    Ll=k;    if (l>k) ll=l;                if (ll>=1) {for (kk=1; kk<=ll; kk++) {if (kk<=k) {d=0.0; for (I=KK; i<=m;i++) {ix= (i-1) *n+kk-1; D=d+a[ix]*a[ix];}                s[kk-1]= (float) sqrt (d);                    if (s[kk-1]!=0.0) {ix= (kk-1) *n+kk-1;                        if (a[ix]!=0.0) {s[kk-1]= (float) fabs (s[kk-1]);                      if (a[ix]<0.0) s[kk-1]=-s[kk-1];                        } for (I=KK; i<=m; i++) {iy= (i-1) *n+kk-1;                      A[IY]=A[IY]/S[KK-1];                  } A[ix]=1.0f+a[ix];              } S[kk-1]=-s[kk-1]; } if (n>=kk+1) {for (j=kk+1; j<=n; J + +) {if (kk<=k) && (S[kk-1]                        !=0.0)) {d=0.0;                            for (I=KK; i<=m; i++) {ix= (i-1) *n+kk-1;                            Iy= (i-1) *n+j-1;                          D=d+a[ix]*a[iy]; } d=-d/a[(Kk-1) *n+kk-1];                            for (I=KK; i<=m; i++) {ix= (i-1) *n+j-1;                            Iy= (i-1) *n+kk-1;                          A[ix]=a[ix]+d*a[iy];                  }} e[j-1]=a[(Kk-1) *n+j-1]; }} if (Kk<=k) {for (I=KK; i<=m; i++) {ix= (i-1) *m+kk-1; iy= (                    i-1) *n+kk-1;                  U[ix]=a[iy];                }} if (kk<=l) {d=0.0;                for (i=kk+1; i<=n; i++) d=d+e[i-1]*e[i-1];                e[kk-1]= (float) sqrt (d);                        if (e[kk-1]!=0.0) {if (e[kk]!=0.0) {e[kk-1]= (float) fabs (e[kk-1]);                      if (e[kk]<0.0) e[kk-1]=-e[kk-1];                    } for (i=kk+1; i<=n; i++) e[i-1]=e[i-1]/e[kk-1];          E[KK]=1.0F+E[KK];        } E[kk-1]=-e[kk-1];                    if ((kk+1<=m) && (e[kk-1]!=0.0)) {for (i=kk+1; i<=m; i++) w[i-1]=0.0; for (j=kk+1; j<=n; j + +) for (i=kk+1; i<=m; i++) w[i-1]=w[i-1]+e[j-1]*a[(                    i-1) *n+j-1];                          for (j=kk+1; j<=n; j + +) for (i=kk+1; i<=m; i++) {ix= (i-1) *n+j-1;                        A[IX]=A[IX]-W[I-1]*E[J-1]/E[KK];              }} for (i=kk+1; i<=n; i++) v[(i-1) *n+kk-1]=e[i-1];    }}} Mm=n;    if (m+1<n) mm=m+1;    if (k<n) s[k]=a[k*n+k];    if (m<mm) s[mm-1]=0.0;    if (l+1<mm) e[l]=a[l*n+mm-1];    e[mm-1]=0.0;    Nn=m;    if (m>n) nn=n;            if (nn>=k+1) {for (j=k+1; j<=nn; J + +) {for (i=1; i<=m; i++) u[(i-1) *m+j-1]=0.0;         u[(j-1) *m+j-1]=1.0; }} if (K>=1) {for (ll=1; ll<=k; ll++) {kk=k-ll+1; iz= (kk-1) *m+kk-1;                       if (s[kk-1]!=0.0) {if (nn>=kk+1) for (j=kk+1; j<=nn; J + +) {d=0.0;                          for (I=KK; i<=m; i++) {ix= (i-1) *m+kk-1;                          Iy= (i-1) *m+j-1;                        D=d+u[ix]*u[iy]/u[iz];                      } d=-d;                          for (I=KK; i<=m; i++) {ix= (i-1) *m+j-1;                          Iy= (i-1) *m+kk-1;                        U[ix]=u[ix]+d*u[iy];                  }} for (I=KK; i<=m; i++) {ix= (i-1) *m+kk-1; U[ix]=-u[ix];}                  U[iz]=1.0f+u[iz];              if (kk-1>=1) for (i=1; i<=kk-1; i++) u[(i-1) *m+kk-1]=0.0; } else {for (i=1; i<=m;i++) u[(i-1) *m+kk-1]=0.0;              u[(kk-1) *m+kk-1]=1.0;        }}} for (ll=1; ll<=n; ll++) {kk=n-ll+1; iz=kk*n+kk-1;                if ((kk<=l) && (e[kk-1]!=0.0)) {for (j=kk+1; j<=n; J + +) {d=0.0;                    for (i=kk+1; i<=n; i++) {ix= (i-1) *n+kk-1; iy= (i-1) *n+j-1;                  D=d+v[ix]*v[iy]/v[iz];                } d=-d;                    for (i=kk+1; i<=n; i++) {ix= (i-1) *n+j-1; iy= (i-1) *n+kk-1;                  V[ix]=v[ix]+d*v[iy];        }}} for (I=1; i<=n; i++) v[(i-1) *n+kk-1]=0.0;      v[iz-n]=1.0;    } for (I=1; i<=m; i++) for (j=1; j<=n; j + +) a[(i-1) *n+j-1]=0.0; m1=mm;    it=60;            while (1==1) {if (mm==0) {PPP (A,E,S,V,M,N); Free (s); Free (e); Free (w);          return (1);          } if (it==0) {PPP (A,E,S,V,M,N);  Free (s); Free (e); Free (w);          Return (-1);            } kk=mm-1;while ((kk!=0) && (Fabs (e[kk-1])!=0.0) {d= (float) (Fabs (s[kk-1]) +fabs (S[KK]));            dd= (float) fabs (e[kk-1]);            if (dd>eps*d) kk=kk-1;          else e[kk-1]=0.0;            } if (kk==mm-1) {kk=kk+1;                if (s[kk-1]<0.0) {s[kk-1]=-s[kk-1];              for (I=1; i<=n; i++) {ix= (i-1) *n+kk-1; V[ix]=-v[ix];}                } while (Kk!=m1 && (S[kk-1]<s[kk])) {d=s[kk-1]; s[kk-1]=s[kk]; s[kk]=d;                      if (kk<n) for (i=1; i<=n; i++) {ix= (i-1) *n+kk-1; iy= (i-1) *n+kk; D=v[ix]; V[ix]=v[iy];                    V[iy]=d; } if (kk<m) for (i=1; i<=m; i++) {ix= (i-1) *m+kk-1; iy= (i-1) *m+kk                      ; D=u[ix]; U[ix]=u[iy];                    U[iy]=d; } kk=kk+1;            } it=60;          Mm=mm-1;            } else {ks=mm;                while ((KS&GT;KK) && (Fabs (s[ks-1])!=0.0) {d=0.0;                if (ks!=mm) d=d+ (float) fabs (e[ks-1]);                if (ks!=kk+1) d=d+ (float) fabs (e[ks-2]);                dd= (float) fabs (s[ks-1]);                if (dd>eps*d) ks=ks-1;              else s[ks-1]=0.0;                } if (KS==KK) {kk=kk+1;                d= (float) fabs (s[mm-1]);                t= (float) fabs (s[mm-2]);                if (t>d) d=t;                t= (float) fabs (e[mm-2]);                if (t>d) d=t;                t= (float) fabs (s[kk-1]);                if (t>d) d=t;                t= (float) fabs (e[kk-1]);                if (t>d) d=t; sm=s[mm-1]/d;                sm1=s[mm-2]/d;                em1=e[mm-2]/d; sk=s[kk-1]/d;                ek=e[kk-1]/d;                B= ((SM1+SM) * (SM1-SM) +em1*em1)/2.0f; C=sm*em1; C=c*c;    shh=0.0;            if ((b!=0.0) | | |                  (c!=0.0))                    {shh= (float) sqrt (b*b+c);                    if (b<0.0) shh=-shh;                  shh=c/(B+SHH);                } fg[0]= (SK+SM) * (SK-SM)-shh;                Fg[1]=sk*ek;                    for (I=KK; i<=mm-1; i++) {sss (FG,CS);                    if (I!=KK) e[i-2]=fg[0];                    FG[0]=CS[0]*S[I-1]+CS[1]*E[I-1];                    E[I-1]=CS[0]*E[I-1]-CS[1]*S[I-1];                    Fg[1]=cs[1]*s[i];                    S[i]=cs[0]*s[i]; if ((cs[0]!=1.0) | | |                      (cs[1]!=0.0))                          for (j=1; j<=n; J + +) {ix= (j-1) *n+i-1;                          Iy= (j-1) *n+i;                          D=cs[0]*v[ix]+cs[1]*v[iy];                          V[iy]=-cs[1]*v[ix]+cs[0]*v[iy];                        V[ix]=d;                    } sss (FG,CS);                    S[I-1]=FG[0];  Fg[0]=cs[0]*e[i-1]+cs[1]*s[i];                  S[i]=-cs[1]*e[i-1]+cs[0]*s[i];                    Fg[1]=cs[1]*e[i];                    E[i]=cs[0]*e[i]; if (i<m) if ((cs[0]!=1.0) | |                        (cs[1]!=0.0))                            for (j=1; j<=m; J + +) {ix= (j-1) *m+i-1;                            Iy= (j-1) *m+i;                            D=cs[0]*u[ix]+cs[1]*u[iy];                            U[iy]=-cs[1]*u[ix]+cs[0]*u[iy];                          U[ix]=d;                }} E[mm-2]=fg[0];              it=it-1;                    } else {if (ks==mm) {kk=kk+1; FG[1]=E[MM-2];                    e[mm-2]=0.0;                        for (LL=KK; ll<=mm-1; ll++) {i=mm+kk-ll-1;                        FG[0]=S[I-1];                        SSS (FG,CS);                        S[I-1]=FG[0];    if (I!=KK) {fg[1]=-cs[1]*e[i-2];                        E[I-2]=CS[0]*E[I-2]; } if ((cs[0]!=1.0) | | |                          (cs[1]!=0.0))                              for (j=1; j<=n; J + +) {ix= (j-1) *n+i-1;                              Iy= (j-1) *n+mm-1;                              D=cs[0]*v[ix]+cs[1]*v[iy];                              V[iy]=-cs[1]*v[ix]+cs[0]*v[iy];                            V[ix]=d;                    }}}} else {kk=ks+1;                    FG[1]=E[KK-2];                    e[kk-2]=0.0;                        for (I=KK; i<=mm; i++) {fg[0]=s[i-1];                        SSS (FG,CS);                        S[I-1]=FG[0];                        FG[1]=-CS[1]*E[I-1];                        E[I-1]=CS[0]*E[I-1]; if ((cs[0]!=1.0) | | |                          (cs[1]!=0.0))   for (j=1; j<=m; J + +) {ix= (j-1) *m+i-1;                           Iy= (j-1) *m+kk-2;                              D=cs[0]*u[ix]+cs[1]*u[iy];                              U[iy]=-cs[1]*u[ix]+cs[0]*u[iy];                            U[ix]=d;   }}}}}} free (s), free (e), free (w);  return (1);    } void PPP (float a[],float e[],float s[],float v[],int m,int n) {int i,j,p,q;    float D;    if (m>=n) i=n;    else i=m;        for (j=1; j<=i-1; J + +) {a[(j-1) *n+j-1]=s[j-1];      a[(j-1) *n+j]=e[j-1];    } a[(I-1) *n+i-1]=s[i-1];    if (m<n) a[(i-1) *n+i]=e[i-1];        For (I=1, i<=n-1; i++) for (j=i+1; j<=n; J + +) {p= (i-1) *n+j-1; q= (j-1) *n+i-1; D=V[P]; V[P]=V[Q];      V[q]=d;  } return;    } void SSS (float fg[],float cs[]) {float r,d;    if ((Fabs (fg[0]) +fabs (fg[1]) ==0.0) {cs[0]=1.0; cs[1]=0.0; d=0.0;}        else {d= (float) sqrt (fg[0]*fg[0]+fg[1]*fg[1]);     if (Fabs (fg[0]) >fabs (fg[1]))     {d= (float) fabs (d);          if (fg[0]<0.0) d=-d;            } if (Fabs (fg[1]) >=fabs (fg[0]) {d= (float) fabs (d);          if (fg[1]<0.0) d=-d; } cs[0]=fg[0]/d;      cs[1]=fg[1]/d;    } r=1.0;    if (Fabs (fg[0]) >fabs (fg[1])) r=cs[1];    else if (cs[0]!=0.0) r=1.0f/cs[0]; Fg[0]=d;    Fg[1]=r;  Return }

Gaussian fit C + + implementation