add: Gauss Fit

ChangeLog.md

@@ -26,3 +26,4 @@
 ### 6.2
 * 曲线平滑
 * 加载与保存数据
+* 高斯寻峰

MCA/ControlView.cpp

@@ -1,12 +1,5 @@
 // ControlView.cpp: 实现文件
 //
-
-#include <math.h>
-#include <chrono>
-#include <fstream>
-#include <sstream>
-#include <iostream>
-
 #include "pch.h"
 
 #include "MCA.h"

MCA/ControlView.h

@@ -1,10 +1,5 @@
 #pragma once
 
-#include <chrono>
-#include <fstream>
-#include <sstream>
-#include <iostream>
-
 #include "MCADoc.h"
 
 // CControlView 窗体视图

MCA/DetailView.cpp

@@ -123,6 +123,7 @@ void CDetailView::OnLButtonDown(UINT nFlags, CPoint point)
 	if (pView == NULL)   pView = (CControlView*)pDoc->GetView(RUNTIME_CLASS(CControlView));
 	int nX = point.x / ZF + pView->m_nCursor1;
 	pView->m_nCursorROI = nX;
+	pView->m_nCursorROICount = pDoc->m_nChannelCount[nX];
 	pView->UpdateData(FALSE);
 	pDoc->UpdateAllViews(NULL);
 	CView::OnLButtonDown(nFlags, point);

MCA/MCA.vcxproj

@@ -187,6 +187,7 @@
     <ClInclude Include="DetailView.h" />
     <ClInclude Include="framework.h" />
     <ClInclude Include="MainFrm.h" />
+    <ClInclude Include="Matrix.h" />
     <ClInclude Include="MCADoc.h" />
     <ClInclude Include="MCAView.h" />
     <ClInclude Include="MCA.h" />

MCA/MCA.vcxproj.filters

@@ -51,6 +51,9 @@
     <ClInclude Include="DetailView.h">
       <Filter>头文件</Filter>
     </ClInclude>
+    <ClInclude Include="Matrix.h">
+      <Filter>头文件</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="MCA.cpp">
@@ -137,5 +140,14 @@
     <Image Include="res\auto.ico">
       <Filter>资源文件</Filter>
     </Image>
+    <Image Include="res\left.ico">
+      <Filter>资源文件</Filter>
+    </Image>
+    <Image Include="res\pause.ico">
+      <Filter>资源文件</Filter>
+    </Image>
+    <Image Include="res\right.ico">
+      <Filter>资源文件</Filter>
+    </Image>
   </ItemGroup>
 </Project>

MCA/MainFrm.cpp

@@ -2,10 +2,17 @@
 // MainFrm.cpp: CMainFrame 类的实现
 //
 
+#include <math.h>
+#include <chrono>
+#include <fstream>
+#include <sstream>
+#include <iostream>
+
 #include "pch.h"
 #include "framework.h"
 
 #include "MCA.h"
+#include "Matrix.h"
 #include "MainFrm.h"
 #include "TotalView.h"
 #include "DetailView.h"
@@ -42,6 +49,7 @@ BEGIN_MESSAGE_MAP(CMainFrame, CFrameWnd)
 	ON_COMMAND(ID_DATA_ORIGIN, &CMainFrame::OnDataOrigin)
 	ON_COMMAND(ID_DATA_3, &CMainFrame::OnData3)
 	ON_COMMAND(ID_DATA_5, &CMainFrame::OnData5)
+	ON_COMMAND(ID_DATA_PEAK, &CMainFrame::OnDataPeak)
 
 	ON_UPDATE_COMMAND_UI(ID_STA_START, &CMainFrame::OnUpdateStaStart)
 	ON_UPDATE_COMMAND_UI(ID_STA_STOP, &CMainFrame::OnUpdateStaStop)
@@ -649,3 +657,77 @@ void CMainFrame::OnUpdateData5(CCmdUI* pCmdUI)
 {
 	pCmdUI->Enable(m_bFSmoothFlag);
 }
+
+void CMainFrame::OnDataPeak()
+{
+	if (pDoc == nullptr)   pDoc = (CMCADoc*)GetActiveDocument();
+	if (pView == nullptr)   pView = (CControlView*)pDoc->GetView(RUNTIME_CLASS(CControlView));
+
+	int L = pView->m_nCursor1;
+	int R = pView->m_nCursor2;
+	int n = R - L + 1;
+
+	if (n > 256) {
+		AfxMessageBox((CString)"ROI区域太大（需≤256）");
+		return ;
+	}
+
+	int cnt = 0;
+	double amplitude, mean, sigma, b0, b1, b2;
+
+	for (int i = L; i <= R; i++)
+		if (pDoc->m_nChannelSmooth[i] != 0)   cnt++;
+
+	matrix mZ(cnt, 1);
+	matrix mX(cnt, 3);
+	matrix mB(3, 1);
+
+	cnt = 0;
+	for (int i = L; i <= R; i++) {
+		if (pDoc->m_nChannelSmooth[i] == 0)   continue;
+		vX[cnt] = i, vY[cnt] = pDoc->m_nChannelSmooth[i];
+		mZ(cnt, 0) = log(pDoc->m_nChannelSmooth[i]);
+		mX(cnt, 0) = 1, mX(cnt, 1) = i, mX(cnt, 2) = i * i;
+		cnt++;
+	}
+
+	mB = (mX.T() * mX).I() * mX.T() * mZ;
+	b0 = mB(0, 0), b1 = mB(1, 0), b2 = mB(2, 0);
+
+	sigma = -1 / b2;
+	mean = b1 * sigma / 2;
+	amplitude = exp(b0 + mean * mean / sigma);
+	sigma = sqrt(sigma);
+
+	matrix mJ(cnt, 3);
+	matrix mH(3, 3);
+	matrix mF(cnt, 1);
+	matrix mG(3, 1);
+	matrix mS(3, 1);
+	matrix mL(3, 3);
+	matrix mU(3, 3);
+
+	double maxV;
+	int* P = new int[3]();
+	double* b = new double[3]();
+
+	for (int i = 0; i < 10; i++) {
+		mJ = Jacobian(cnt, &vX[0], amplitude, mean, sigma);
+		mH = mJ.T() * mJ;
+		for (int j = 0; j < cnt; j++)   mF(j, 0) = vY[j] - Gaussian(vX[j], amplitude, mean, sigma);
+		mG = mJ.T() * mF;
+		for (int j = 0; j < 3; j++)   b[j] = -mG(j, 0);
+
+		if (!LUPDescomposition(mH, mL, mU, P))   return;
+		mS = LUPSolve(mL, mU, P, b);
+		amplitude += mS(0, 0), mean += mS(1, 0), sigma += mS(2, 0);
+
+		maxV = 0;
+		for (int j = 0; j < 3; j++)   maxV = max(maxV, abs(mS(j, 0)));
+		if (maxV < 1e-5)   break;
+	}
+
+	CString strText;
+	strText.Format(_T("幅值%.3f\n均值%.3f\n标准差%.3f"), amplitude, mean, sigma);
+	AfxMessageBox(strText);
+}

MCA/MainFrm.h

@@ -4,6 +4,11 @@
 
 #pragma once
 
+#include <chrono>
+#include <fstream>
+#include <sstream>
+#include <iostream>
+
 #include "MCADoc.h"
 #include "ControlView.h"
 
@@ -83,6 +88,7 @@ public:
     afx_msg void OnDataOrigin();
     afx_msg void OnData3();
     afx_msg void OnData5();
+    afx_msg void OnDataPeak();
 
     // 点击状态
     BOOL m_bStartFlag = TRUE;
@@ -110,6 +116,10 @@ public:
     afx_msg void OnUpdateDataOrigin(CCmdUI* pCmdUI);
     afx_msg void OnUpdateData3(CCmdUI* pCmdUI);
     afx_msg void OnUpdateData5(CCmdUI* pCmdUI);
+
+    // 寻峰
+    double vX[256]{};
+    double vY[256]{};
 };
 
 

MCA/Matrix.h

@@ -0,0 +1,234 @@
+#ifndef DESCSS_Matrix_h
+#define DESCSS_Matrix_h
+
+#include <iostream>
+#define swap2(x, y) ((x)=(x)^(y),(y)=(x)^(y),(x)=(x)^(y))
+
+class matrix {
+public:
+    matrix(int n, int m);
+    matrix(const matrix&);
+    ~matrix();
+
+    void Input(double* p);
+    void Output();
+    int rows() const { return n; };
+    int cols() const { return m; };
+
+    double& operator()(int i, int j) const;
+    matrix& operator=(const matrix&);
+    matrix operator+(const matrix&) const;
+    matrix operator*(const matrix&) const;
+    matrix operator+(const double&) const;
+    matrix operator*(const double&) const;
+    matrix operator/(const double&) const;
+
+    matrix I();
+    matrix T();
+
+protected:
+    int n, m;
+    double** element;
+};
+
+matrix::matrix(int n, int m) {
+    if (m < 0 || n < 0) std::cout << "Rows and Cols must be >= 0 " << std::endl;
+    if ((m == 0 || n == 0) && (m != 0 || n != 0))
+        std::cout << "Either both or neither rows and columns should be zero " << std::endl;
+    this->m = m, this->n = n;
+    element = new double*[n];
+    for (int i = 0; i < n; i++) element[i] = new double[m];
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) element[i][j] = 0;
+}
+
+matrix::matrix(const matrix& A) {
+    n = A.n, m = A.m;
+    element = new double*[n];
+    for (int i = 0; i < n; i++) element[i] = new double[m];
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) element[i][j] = A.element[i][j];
+}
+
+matrix::~matrix() {}
+
+matrix& matrix::operator=(const matrix& A) {
+    if (this != &A) {
+        delete[] element;
+        n = A.n, m = A.m;
+        element = new double*[n];
+        for (int i = 0; i < n; i++) element[i] = new double[m];
+        for (int i = 0; i < n; i++)
+            for (int j = 0; j < m; j++) element[i][j] = A.element[i][j];
+    }
+    return *this;
+}
+
+void matrix::Input(double* p) {
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) this->element[i][j] = *(p++);
+}
+
+void matrix::Output() {
+    for (int i = 0; i < n; i++) {
+        for (int j = 0; j < m; j++) std::cout << element[i][j] << " ";
+        std::cout << std::endl;
+    }
+}
+
+double& matrix::operator()(int i, int j) const {
+    if (i < 0 || i > n - 1 || j < 0 || j > m - 1) std::cout << "Matrix Index Out Of Bounds " << std::endl;
+    return element[i][j];
+}
+
+matrix matrix::operator+(const matrix& A) const {
+    if (n != A.n || m != A.m) std::cout << "Matrix Size is Out of batch " << std::endl;
+    matrix w(n, m);
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) w.element[i][j] = element[i][j] + A.element[i][j];
+    return w;
+}
+
+matrix matrix::operator+(const double& C) const {
+    matrix w(n, m);
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) w.element[i][j] = element[i][j] + C;
+    return w;
+}
+
+matrix matrix::operator*(const matrix& A) const {
+    if (m != A.n) std::cout << "Matrix Style is Out of batch " << std::endl;
+    matrix w(n, A.m);
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < A.m; j++) {
+            w.element[i][j] = 0;
+            for (int k = 0; k < m; k++) w.element[i][j] += element[i][k] * A.element[k][j];
+        }
+    return w;
+}
+
+matrix matrix::operator*(const double& C) const {
+    matrix w(n, m);
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) w.element[i][j] = element[i][j] * C;
+    return w;
+}
+
+matrix matrix::operator/(const double& C) const {
+    matrix w(n, m);
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) w.element[i][j] = element[i][j] / C;
+    return w;
+}
+
+BOOL LUPDescomposition(matrix mA, matrix& mL, matrix& mU, int* P)
+{
+    int n = mA.rows();
+    int row = 0, tmp;
+    double L, U, tmp2;
+    for (int i = 0; i < n; i++)   P[i] = i;
+    for (int i = 0; i < n - 1; i++) {
+        tmp2 = 0.;
+        for(int j = i; j < n; j++)
+            if (std::abs(mA(j, i)) > tmp2) {
+                tmp2 = std::abs(mA(j, i));
+                row = j;
+            }
+
+        if (tmp2 == 0) {
+            AfxMessageBox((CString)"矩阵奇异，无法计算");
+            return FALSE;
+        }
+
+        tmp = P[i], P[i] = P[row], P[row] = tmp;
+
+        for (int j = 0; j < n; j++)
+            tmp2 = mA(i, j), mA(i, j) = mA(row, j), mA(row, j) = tmp2;
+
+        U = mA(i, i), L = 0;
+        for (int j = i + 1; j < n; j++) {
+            L = mA(j, i) / U;
+            mA(j, i) = L;
+            for (int k = i + 1; k < n; k++)
+                mA(j, k) = mA(j, k) - mA(i, k) * L;
+        }
+    }
+    for (int i = 0; i < n; i++) {
+        mL(i, i) = 1;
+        for (int j = 0; j < i; j++)   mL(i, j) = mA(i, j);
+        for (int j = i; j < n; j++)   mU(i, j) = mA(i, j);
+    }
+    return TRUE;
+}
+
+matrix LUPSolve(matrix& mL, matrix& mU, int* P, double* b) {
+    int n = mL.rows();
+    matrix mX(n, 1);
+    double* y = new double[n]();
+
+    for (int i = 0; i < n; i++) {
+        y[i] = b[P[i]];
+        for (int j = 0; j < i; j++)   y[i] -= mL(i, j) * y[j];
+    }
+
+    for (int i = n - 1; i >= 0; i--) {
+        mX(i, 0) = y[i];
+        for (int j = n - 1; j > i; j--)   mX(i, 0) -= mU(i, j) * mX(j, 0);
+        mX(i, 0) /= mU(i, i);
+    }
+
+    return mX;
+}
+
+// Inverse
+matrix matrix::I() {
+    matrix mX(n, 1);
+    matrix mL(n, m);
+    matrix mU(n, m);
+    matrix mMir(n, m);
+    matrix mInv(n, m);
+    int* P = new int[n]();
+    double* b = new double[n]();
+
+    if (n != m) {
+        AfxMessageBox((CString)"无法计算非方阵");
+        return mInv;
+    }
+
+    for (int j = 0; j < n; j++)
+            for (int k = 0; k < m; k++) mMir.element[j][k] = element[j][k];
+    LUPDescomposition(mMir, mL, mU, P);
+
+    for (int i = 0; i < n; i++) {
+        for (int j = 0; j < n; j++)   b[j] = 0;
+        b[i] = 1;
+
+        mX = LUPSolve(mL, mU, P, b);
+        for (int j = 0; j < n; j++)   mInv(j, i) = mX(j, 0);
+    }
+    return mInv;
+}
+
+// Transposition
+matrix matrix::T() {
+    matrix w(m, n);
+    for (int i = 0; i < m; i++)
+        for (int j = 0; j < n; j++) w.element[i][j] = element[j][i];
+    return w;
+}
+
+double Gaussian(double x, double A, double mean, double sigma) {
+    return A * exp(-(x - mean) * (x - mean) / (2 * sigma * sigma));
+}
+
+matrix Jacobian(int n, double* X, double A, double mean, double sigma) {
+    matrix J(n, 3);
+    for (int i = 0; i < n; i++) {
+        J(i, 0) = -Gaussian(X[i], A, mean, sigma) / A;
+        J(i, 1) = -(X[i] - mean) * Gaussian(X[i], A, mean, sigma) / (sigma * sigma);
+        J(i, 2) = -(X[i] - mean) * (X[i] - mean) * Gaussian(X[i], A, mean, sigma) / (sigma * sigma * sigma);
+    }
+    return J;
+}
+
+#endif

MCA/Resource.h

@@ -55,13 +55,16 @@
 #define ID_STA_CLEAR                    32837
 #define ID_OPT_SAVE                     32840
 #define ID_OPT_OPEN                     32841
+#define ID_32842                        32842
+#define ID_                             32843
+#define ID_DATA_PEAK                    32844
 
 // Next default values for new objects
 // 
 #ifdef APSTUDIO_INVOKED
 #ifndef APSTUDIO_READONLY_SYMBOLS
 #define _APS_NEXT_RESOURCE_VALUE        325
-#define _APS_NEXT_COMMAND_VALUE         32842
+#define _APS_NEXT_COMMAND_VALUE         32847
 #define _APS_NEXT_CONTROL_VALUE         1019
 #define _APS_NEXT_SYMED_VALUE           310
 #endif

README.md

@@ -11,7 +11,7 @@
 - [x] 对数坐标轴
 - [x] 读取与保存
 - [x] 模拟信号
-- [ ] 高斯拟合寻峰
+- [x] 高斯拟合寻峰
 - [x] 曲线平滑
 
 ## Bugs
@@ -23,3 +23,4 @@
 - [x] 保存文件乱码
 - [x] 退出询问保存
 - [x] 读取文件后刷新视图无效
+- [ ] 禁止更改窗口大小