diff --git a/ChangeLog.md b/ChangeLog.md index f0ef8f0..65c2de0 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -26,3 +26,4 @@ ### 6.2 * 曲线平滑 * 加载与保存数据 +* 高斯寻峰 diff --git a/MCA/ControlView.cpp b/MCA/ControlView.cpp index 7a832d5..136ce68 100644 --- a/MCA/ControlView.cpp +++ b/MCA/ControlView.cpp @@ -1,12 +1,5 @@ // ControlView.cpp: 实现文件 // - -#include -#include -#include -#include -#include - #include "pch.h" #include "MCA.h" diff --git a/MCA/ControlView.h b/MCA/ControlView.h index 5469944..f49b269 100644 --- a/MCA/ControlView.h +++ b/MCA/ControlView.h @@ -1,10 +1,5 @@ #pragma once -#include -#include -#include -#include - #include "MCADoc.h" // CControlView 窗体视图 diff --git a/MCA/DetailView.cpp b/MCA/DetailView.cpp index b60b566..b70c22d 100644 --- a/MCA/DetailView.cpp +++ b/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); diff --git a/MCA/MCA.rc b/MCA/MCA.rc index 9d8e7d3..bc63674 100644 Binary files a/MCA/MCA.rc and b/MCA/MCA.rc differ diff --git a/MCA/MCA.vcxproj b/MCA/MCA.vcxproj index 6275421..02d1e7f 100644 --- a/MCA/MCA.vcxproj +++ b/MCA/MCA.vcxproj @@ -187,6 +187,7 @@ + diff --git a/MCA/MCA.vcxproj.filters b/MCA/MCA.vcxproj.filters index 9585ff3..433bb0f 100644 --- a/MCA/MCA.vcxproj.filters +++ b/MCA/MCA.vcxproj.filters @@ -51,6 +51,9 @@ 头文件 + + 头文件 + @@ -137,5 +140,14 @@ 资源文件 + + 资源文件 + + + 资源文件 + + + 资源文件 + \ No newline at end of file diff --git a/MCA/MainFrm.cpp b/MCA/MainFrm.cpp index ec6e05f..338c3fc 100644 --- a/MCA/MainFrm.cpp +++ b/MCA/MainFrm.cpp @@ -2,10 +2,17 @@ // MainFrm.cpp: CMainFrame 类的实现 // +#include +#include +#include +#include +#include + #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); +} \ No newline at end of file diff --git a/MCA/MainFrm.h b/MCA/MainFrm.h index 46fb03e..17a07e9 100644 --- a/MCA/MainFrm.h +++ b/MCA/MainFrm.h @@ -4,6 +4,11 @@ #pragma once +#include +#include +#include +#include + #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]{}; }; diff --git a/MCA/Matrix.h b/MCA/Matrix.h new file mode 100644 index 0000000..e3bb818 --- /dev/null +++ b/MCA/Matrix.h @@ -0,0 +1,234 @@ +#ifndef DESCSS_Matrix_h +#define DESCSS_Matrix_h + +#include +#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 \ No newline at end of file diff --git a/MCA/Resource.h b/MCA/Resource.h index b77b58f..61a304f 100644 --- a/MCA/Resource.h +++ b/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 diff --git a/README.md b/README.md index bb9194e..08d9612 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ - [x] 对数坐标轴 - [x] 读取与保存 - [x] 模拟信号 -- [ ] 高斯拟合寻峰 +- [x] 高斯拟合寻峰 - [x] 曲线平滑 ## Bugs @@ -23,3 +23,4 @@ - [x] 保存文件乱码 - [x] 退出询问保存 - [x] 读取文件后刷新视图无效 +- [ ] 禁止更改窗口大小 \ No newline at end of file