diff --git a/src/Function2D.cpp b/src/Function2D.cpp index 8bbfd93..5fe5932 100644 --- a/src/Function2D.cpp +++ b/src/Function2D.cpp @@ -1034,3 +1034,87 @@ Matrix Aurora::sub2ind(const Matrix &aVMatrixSize, std::initializer_list 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]; + } + } +} diff --git a/src/Function2D.cu b/src/Function2D.cu index 57f4f46..fa2f253 100644 --- a/src/Function2D.cu +++ b/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]=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]=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<<>>(aData.getData(), rowCount, colCount, aNpeaks, aMinPeakHeight, aMinPeakProminence,peaks,*aOutIndexs); + cudaDeviceSynchronize(); + // cudaFree(peaks); + return Aurora::CudaMatrix::fromRawData(peaks, aNpeaks, colCount); +} \ No newline at end of file diff --git a/src/Function2D.cuh b/src/Function2D.cuh index 2fde164..01b1888 100644 --- a/src/Function2D.cuh +++ b/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__ \ No newline at end of file diff --git a/src/Function2D.h b/src/Function2D.h index c10f62d..36ab9a7 100644 --- a/src/Function2D.h +++ b/src/Function2D.h @@ -178,7 +178,7 @@ namespace Aurora * @return */ Matrix sub2ind(const Matrix &aVMatrixSize, std::initializer_list aSliceIdxs); - + void findPeaksHost(const Matrix & aData, int aNPeaks, float aMinPeakHeight, float MinPeakProminece,int* outIndex); }; #endif // AURORA_FUNCTION2D_H diff --git a/test/Function2D_Cuda_Test.cpp b/test/Function2D_Cuda_Test.cpp index 7af762c..356a5b6 100644 --- a/test/Function2D_Cuda_Test.cpp +++ b/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; +} + diff --git a/test/Function2D_Test.cpp b/test/Function2D_Test.cpp index 822e058..2994997 100644 --- a/test/Function2D_Test.cpp +++ b/test/Function2D_Test.cpp @@ -2,6 +2,8 @@ #include #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; +} +