feat: add findPeaks Function

src/Function2D.cpp

@@ -1034,3 +1034,87 @@ Matrix Aurora::sub2ind(const Matrix &aVMatrixSize, std::initializer_list<Matrix>
     delete [] strides;
     return Matrix::New(output,returnVectorSize,1,1);
 }
+
+void Aurora::findPeaksHost(const Matrix & aData,  int aNPeaks, float aMinPeakHeight, float MinPeakProminece,
+        int* outIndex)
+{
+    int  signalSize = aData.getDimSize(0);
+    int  signalCount = aData.getDimSize(1);
+
+    #pragma omp parallel for
+    for (size_t threadIndex = 0; threadIndex < signalCount; threadIndex++)
+    {
+        float* dataPointer  = aData.getData() + threadIndex*signalSize ;
+        float newPeak = dataPointer[0];
+        float newValley = dataPointer[0];
+        int peakIndex = 0;
+        float higherValley;
+        int indexs[32];
+        float values[32];
+        for (size_t i = 0; i < aNPeaks; i++)
+        {
+            indexs[i] = signalSize;
+            values[i] = 0;
+        }
+        int save_index=0;
+
+        for (int i = 1; i < signalSize - 1; ++i)
+        {
+            // find peaks
+            if (dataPointer[i] > dataPointer[i - 1] && dataPointer[i] > dataPointer[i + 1])
+            {
+                newPeak = dataPointer[i];
+                peakIndex  = i;
+            };
+            // find valley
+            if (dataPointer[i] < dataPointer[i - 1] && dataPointer[i] < dataPointer[i + 1])
+            {
+                higherValley = std::max(newValley, dataPointer[i]);
+                newValley = dataPointer[i];
+                if (newPeak >= aMinPeakHeight)
+                {
+                    float prominece = newPeak - higherValley;
+                    if (prominece >= MinPeakProminece)
+                    {
+                        if (save_index < aNPeaks)
+                        {
+                            values[save_index] = newPeak;
+                            indexs[save_index] = peakIndex;
+                            save_index++;
+                        }
+                        else
+                        {
+                            for (size_t j = 0; j < aNPeaks; j++)
+                            {
+                                if (values[j] < newPeak)
+                                {
+                                    std::swap(values[j], newPeak);
+                                    std::swap(indexs[j], peakIndex);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+        if (save_index>=aNPeaks)
+        {
+            for (size_t i = 0; i < aNPeaks; i++)
+            {
+                for (size_t j = i+1; j < aNPeaks; j++)
+                {
+                    if (indexs[i]>indexs[j])
+                    {
+                        std::swap(values[j], values[i]);
+                        std::swap(indexs[j], indexs[i]);
+                    }
+                }
+            }
+        }
+        for (size_t i = 0; i < aNPeaks; i++)
+        {
+            // aOutPeaks[threadIndex*aNPeaks+i] = values[i];
+            outIndex[threadIndex*aNPeaks+i] = indexs[i];
+        }
+    }
+}

src/Function2D.cu

@@ -1711,7 +1711,7 @@ __global__ void validKernel(const float* aData, const float* aValid, float* aOut
     }
 }
 
-Aurora::CudaMatrix Aurora::valid(const Aurora::CudaMatrix aData, const Aurora::CudaMatrix aValid)
+Aurora::CudaMatrix Aurora::valid(const Aurora::CudaMatrix& aData, const Aurora::CudaMatrix aValid)
 {
     auto validEnable = Aurora::sum(aValid);
     if(validEnable.toHostMatrix()[0] == 0)
@@ -1748,3 +1748,107 @@ Aurora::CudaMatrix Aurora::valid(const Aurora::CudaMatrix aData, const Aurora::C
     delete[] validProcessed;
     return Aurora::CudaMatrix::fromRawData(result, rowCount, validColumnCount);
 }
+
+__global__ void findPeaksKernel(float* aData, int aSignalSize, int aSignalCount, int aNPeaks, float aMinPeakHeight, float MinPeakProminece,
+         float* aOutPeaks, int* outIndex)
+{
+    int threadIndex = blockIdx.x * blockDim.x;
+    float* dataPointer  = aData + threadIndex*aSignalSize ;
+    float newPeak = dataPointer[0];
+    float newValley = dataPointer[0];
+    float higherValley;
+    __shared__ int indexs[32];
+    __shared__ float values[32];
+    for (size_t i = 0; i < aNPeaks; i++)
+    {
+        indexs[i] = aSignalSize;
+        values[i] = 0;
+    }
+    int save_index = 0;
+    int peakIndex = -1;
+    {
+        for(int i=1; i < aSignalSize-1; ++i)
+        {
+            //find peaks
+            if (dataPointer[i]>dataPointer[i-1] && dataPointer[i]>dataPointer[i+1]){
+                newPeak= dataPointer[i];
+                peakIndex = i;
+            };
+            //find valley
+            if (dataPointer[i]<dataPointer[i-1] && dataPointer[i]<dataPointer[i+1]){
+                higherValley = max(newValley,dataPointer[i]);
+                newValley = dataPointer[i];
+                if (newPeak>=aMinPeakHeight)
+                {
+                    float prominece = newPeak-higherValley;
+                    if (prominece>=MinPeakProminece)
+                    {
+                        if (save_index < aNPeaks){
+                            values[save_index] = newPeak;
+                            indexs[save_index] = peakIndex;
+                            save_index ++;    
+                        }
+                        else{
+                            float tempPeak = newPeak;
+                            int tempIndex = peakIndex;
+
+                            for (size_t j = 0; j < aNPeaks; j++)
+                            {
+                                if (values[j]<newPeak)
+                                {
+                                    tempPeak = values[j];
+                                    tempIndex = indexs[j];
+                                    values[j] = newPeak;
+                                    indexs[j] = peakIndex;
+                                    newPeak = tempPeak;
+                                    peakIndex = tempIndex;
+                                }
+                            }
+                        }
+                    }
+                }
+            };
+
+        }
+    }
+    if (save_index>=aNPeaks)
+    {
+        for (size_t i = 0; i < aNPeaks; i++)
+        {
+            for (size_t j = i+1; j < aNPeaks; j++)
+            {
+                if (indexs[i]>indexs[j])
+                {
+                    int temp = indexs[j];
+                    indexs[j] = indexs[i];
+                    indexs[i] = temp;
+                    float tempValue = values[j];
+                    values[j] = values[i];
+                    values[i] = tempValue;
+                }
+            }
+        }
+    }
+    for (size_t i = 0; i < aNPeaks; i++)
+    {
+        aOutPeaks[threadIndex*aNPeaks+i] = values[i];
+        outIndex[threadIndex*aNPeaks+i] = indexs[i];
+    }
+}
+
+Aurora::CudaMatrix Aurora::findPeaks(const Aurora::CudaMatrix& aData, int aNpeaks, float aMinPeakHeight, float aMinPeakProminence,int** aOutIndexs )
+{
+    int rowCount = aData.getDimSize(0);
+    int colCount = aData.getDimSize(1);
+
+    float* peaks = nullptr ;
+    cudaMalloc((void**)&peaks, sizeof(float) * aNpeaks * colCount );
+    cudaMalloc((void**)aOutIndexs, sizeof(int) * aNpeaks * colCount );
+
+    int threadPerBlock = 256;
+    int blockPerGrid = colCount / threadPerBlock + 1;
+    findPeaksKernel<<<colCount, 1>>>(aData.getData(), rowCount, colCount, aNpeaks, aMinPeakHeight, aMinPeakProminence,peaks,*aOutIndexs);
+    cudaDeviceSynchronize();
+    // cudaFree(peaks);
+    return Aurora::CudaMatrix::fromRawData(peaks, aNpeaks, colCount);
+}

src/Function2D.cuh

@@ -87,8 +87,19 @@ namespace Aurora
      */
     CudaMatrix ifft_symmetric(const CudaMatrix &aMatrix,long aLength);
 
-    CudaMatrix valid(const CudaMatrix aData, const CudaMatrix aValid);
+    CudaMatrix valid(const CudaMatrix& aData, const CudaMatrix aValid);
 
+    /**
+     * findPeaks 按列进行峰查找和匹配
+     * @attention 不要给aOutIndexs提前分配内存
+     * @param aData 输入
+     * @param aNpeaks 峰数量
+     * @param aMinPeakHeight 最小高度
+     * @param aMinPeakProminence 最小相对高度
+     * @param aOutIndexs 空指针，会在函数内分配device内存
+     * @return 筛选出的峰高度
+    */
+    CudaMatrix findPeaks(const CudaMatrix& aData, int aNpeaks, float aMinPeakHeight, float aMinPeakProminence, int** aOutIndexs);
 }
 
 #endif // __FUNCTION2D_CUDA_H__

src/Function2D.h

@@ -178,7 +178,7 @@ namespace Aurora
      * @return
      */
     Matrix sub2ind(const Matrix &aVMatrixSize, std::initializer_list<Matrix> aSliceIdxs);
-    
+    void findPeaksHost(const Matrix & aData,  int aNPeaks, float aMinPeakHeight, float MinPeakProminece,int* outIndex);
 };
 
 #endif // AURORA_FUNCTION2D_H

test/Function2D_Cuda_Test.cpp

@@ -5,6 +5,7 @@
 #include "Function.h"
 #include "Matrix.h"
 #include "TestUtility.h"
+#include "MatlabReader.h"
 
 #include "Function2D.h"
 #include "Function2D.cuh"
@@ -18,11 +19,15 @@ protected:
     static void TearDownTestCase(){
     }
     public:
+    Aurora::Matrix mSignal;
+    Aurora::CudaMatrix dmSignal;
     Aurora::Matrix B;
     Aurora::CudaMatrix dB;
 
     void SetUp(){
-
+        MatlabReader m("/home/krad/TestData/peaks.mat");
+        mSignal = m.read("AScan_env_norm");
+        dmSignal = mSignal.toDeviceMatrix();
     }
     void TearDown(){
     }
@@ -997,3 +1002,17 @@ TEST_F(Function2D_Cuda_Test, hilbert) {
         EXPECT_NEAR(ret1[i], ret2.getValue(i), 0.01);
     }
 }
+
+TEST_F(Function2D_Cuda_Test, findPeaks) {
+
+    int* indexs = new int[mSignal.getDimSize(1)*10];
+    auto ret1 = Aurora::findPeaks(dmSignal,10, 0.2, 0.05,indexs);
+
+    auto reH = ret1.toHostMatrix();
+    for(unsigned int i=0; i<10; ++i)
+    {
+        printf("%d,",indexs[i]);
+    }
+    delete [] indexs;
+}
+

test/Function2D_Test.cpp

@@ -2,6 +2,8 @@
 #include <vector>
 #include "TestUtility.h"
 
+#include "MatlabReader.h"
+
 #include "Matrix.h"
 #include "Function.h"
 #include "Function1D.h"
@@ -16,7 +18,11 @@ protected:
     }
     static void TearDownTestCase(){
     }
+    public:
+    Aurora::Matrix mSignal;
     void SetUp(){
+        MatlabReader m("/home/krad/TestData/peaks.mat");
+        mSignal = m.read("AScan_env_norm");
     }
     void TearDown(){
     }
@@ -573,3 +579,15 @@ TEST_F(Function2D_Test, sub2ind) {
 
 }
 
+TEST_F(Function2D_Test, findPeaks) {
+
+    int* indexs = new int[mSignal.getDimSize(1)*10];
+    Aurora::findPeaksHost(mSignal,10, 0.2, 0.05,indexs);
+
+    for(unsigned int i=0; i<10; ++i)
+    {
+        printf("%d,",indexs[i]);
+    }
+    delete [] indexs;
+}
+