Gaussian fitting c ++ implementation, Gaussian fitting c

Gaussian fitting c ++ implementation, Gaussian fitting c

Gaussian function is a very important mathematical function. The density function of the normal distribution we are familiar with is Gaussian function, also known as Gaussian distribution. Normal Distribution is undoubtedly the most important distribution in probability theory and mathematical statistics.

The problem now is how to find a Gaussian function to fit some points set if some points are given!

Of course, the fitting method is the least square method, and the fitting function form is:

Y = a * exp (-(x-B)/c) ^ 2 );

There are three parameters in total, a, B, and c. However, this exponential function is difficult to fit, so the logarithm transformation is used to convert it into a quadratic function, as shown below:

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

In this way, the other z = ln (y), A =-(1/c ^ 2), B = (2b/c ^ 2), C = ln () -(B/c) ^ 2; then we can use the least square method to fit the data,

Here I still use the Singular Value Decomposition Method to Solve the least square solution. However, before that, I first processed the point set. I did not fit all the points set, because the effect was very poor, after all, a logarithm transformation is performed,

As a result, the flat point has a great influence on the curve, So I first determined the peak point and then fitted the peak point. It turns out that my idea is correct, and the fitting effect has improved a lot!

The VC implementation result is as follows:

The core program implementation is as follows:

// GaussFit. h

/*************************************** * ********************************* Version: function Description: The Least Squares Gaussian fitting is provided for some column points on the plane. The least squares solution is obtained using the Singular Value Decomposition method as the Gaussian parameter. Call method: gaussfit (arrayx, arrayy, int n, float box, miny); parameter description: arrayx [n]. Each value is an arrayx point on the X axis: arrayy [n], where each value is a point on the Y axis n: number of points box: box [3]. The three parameters of the Gaussian function are a, B, c; miny: translation in the y direction. 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) {ret Urn p1.x> p2.x ;}}; class GaussFit {public: GaussFit (void );~ GaussFit (void); void gaussfit (int * arrayx, int * arrayy, int n, float * box, int & miny); private: 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 for point. x sort priority_queue <int> py; // used to calculate the point smaller than m. yint m = n/10; m = m> 0? M: 1; for (I = 0; I <n; I ++) {GPOINT tmp (arrayx [I], arrayy [I]); gp. push (tmp); if (py. size () <m) {py. push (arrayy [I]);} else if (py. top ()> arrayy [I]) {py. pop (); py. push (arrayy [I]);} miny = py. top (); // replace y with y smaller than m to prevent 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 mean Y = 0; for (int I = 0; I <n; I ++) {meany + = Pointy [I];} meany/= n; // calculate the peak values of 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 <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 <s2) {VG. clear (); for (int j = 0; j <Vtmp. size (); j ++) {GPOINT tmp1 (Vtmp [j]. x, Vtmp [j]. y); VG. push_back (tmp1) ;}// perform Gaussian fitting on the peak value int size = VG. size (); 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] = e Xp (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 = gm Iv (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 <= 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); I F (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> 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 = (f Loat) 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 ;}