From 0f158107e9e48e3466f9546d6ec7d75d3e72a915 Mon Sep 17 00:00:00 2001 From: Catherine Olschanowky Date: Wed, 6 Nov 2019 23:03:38 -0700 Subject: [PATCH] Added support for flat sections to be identified as peaks --- src/GaussianFitter.cpp | 80 ++++++++++++++++++++++++++++---- src/GaussianFitter_unittests.cpp | 35 +++++++++----- 2 files changed, 94 insertions(+), 21 deletions(-) diff --git a/src/GaussianFitter.cpp b/src/GaussianFitter.cpp index e881596..f591c15 100644 --- a/src/GaussianFitter.cpp +++ b/src/GaussianFitter.cpp @@ -614,19 +614,81 @@ int GaussianFitter::guess_peaks(std::vector* results, //are pointing to space used in LidarVolume results->clear(); - //Sign of gradient: - // = 1 for increasing - // = 0 for level AND PREVIOUSLY INCREASING (so potential wide peak) - // = -1 for decreasing OR level, but previously decreasing - //A sharp peak is identified by grad=1 -> grad=-1 - //A wide peak is identified by grad=0 -> grad=-1 std::vector peak_guesses_loc; //Store peaks x-values here std::vector peak_guesses_amp; - int wideStart = -1; //The start of any current wide peak int prev_grad = -1; int grad = -1; - for(int i = 0; i<(int)ampData.size()-1; i++){ + int firstDiffs[(int)ampData.size()]; + int secondDiffs[(int)ampData.size()]; + firstDiffs[1] = ampData[1] - ampData[0]; + for(int i = 2; i<(int)ampData.size(); i++){ + firstDiffs[i] = ampData[i] - ampData[i-1]; + secondDiffs[i] = firstDiffs[i] - firstDiffs[i-1]; + } + for(int i = 2; i<(int)ampData.size(); i++){ + if(grad == 0){ + if(firstDiffs[i] > 0){ + // flat to positive (record peak if not a trough) + if(ampData[i]>noise_level && prev_grad == 1){ + //handle flat peaks + int j = i-2; + for( ; ampData[j] == ampData[i-1] ; j--); + // there is no body to the above for loop -- just looking + // for the stopping point + float lenOfFlat = i-1 -j -1.; + //this is truncating, we need + //to fix this + float loc = i-1 - lenOfFlat/2.; + //record the peak + peak_guesses_loc.push_back(loc); + peak_guesses_amp.push_back(ampData[i-1]); + } + grad = 1; + }else if(firstDiffs[i] < 0 ){ + // flat to negative (record peak) + if(ampData[i]>noise_level){ + //handle flat peaks + int j = i-2; + for( ; ampData[j] == ampData[i-1] ; j--); + // there is no body to the above for loop -- just looking + // for the stopping point + float lenOfFlat = i-1 -j -1.; + //this is truncating, we need + //to fix this + float loc = i-1 - lenOfFlat/2.; + //record the peak + peak_guesses_loc.push_back(loc); + peak_guesses_amp.push_back(ampData[i-1]); + } + grad = -1; + } + + }else if(grad == 1){ + prev_grad = 1; + if(firstDiffs[i] == 0 ){ + grad = 0; + }else if(firstDiffs[i] < 0){ + // positive to negative (record peak) + if(ampData[i]>noise_level){ + //record the peak + peak_guesses_loc.push_back(float(i-1)); + peak_guesses_amp.push_back(ampData[i-1]); + } + grad = -1; + } + }else{ // gradient is -1 + prev_grad = -1; + if(firstDiffs[i] == 0 ){ + // down to flat + grad = 0; + }else if(firstDiffs[i] > 0){ + // down to up + grad = 1; + } + } + } + /*for(int i = 0; i<(int)ampData.size()-1; i++){ // were we sloping up before? if (grad == 1){ // sloping down @@ -659,7 +721,7 @@ int GaussianFitter::guess_peaks(std::vector* results, grad = 1; } } - } + }*/ //Figure out the size of ampData size_t n = ampData.size(); diff --git a/src/GaussianFitter_unittests.cpp b/src/GaussianFitter_unittests.cpp index c39b723..d2b3195 100644 --- a/src/GaussianFitter_unittests.cpp +++ b/src/GaussianFitter_unittests.cpp @@ -72,15 +72,20 @@ TEST_F(GaussianFitterTest, NayaniClipped1_guess){ std::vector peaks; int count = fitter.guess_peaks(&peaks, ampData, idxData); - EXPECT_EQ(2,peaks.size()); + ASSERT_EQ(3,peaks.size()); EXPECT_EQ(200,peaks.at(0)->amp); - EXPECT_EQ(34 ,peaks.at(1)->amp); EXPECT_NEAR(18.5, peaks.at(0)->location, .25); - EXPECT_EQ(38, peaks.at(1)->location); EXPECT_NEAR(10.5, peaks.at(0)->fwhm, 1); - EXPECT_NEAR(6, peaks.at(1)->fwhm, 1); - EXPECT_EQ(2, count); + EXPECT_EQ(43 ,peaks.at(1)->amp); + EXPECT_NEAR(30.5, peaks.at(1)->location,.25); + EXPECT_NEAR(4, peaks.at(1)->fwhm, 1); + + EXPECT_EQ(34 ,peaks.at(2)->amp); + EXPECT_EQ(38, peaks.at(2)->location); + EXPECT_NEAR(6, peaks.at(2)->fwhm, 1); + + EXPECT_EQ(3, count); } @@ -102,7 +107,7 @@ TEST_F(GaussianFitterTest, NayaniClipped2_guess){ ASSERT_EQ(2,peaks.size()); EXPECT_EQ(235,peaks.at(0)->amp); - EXPECT_EQ(18, peaks.at(0)->location); + EXPECT_NEAR(17.5, peaks.at(0)->location,1); EXPECT_NEAR(7.8, peaks.at(0)->fwhm, 1); EXPECT_EQ(29, peaks.at(1)->location); @@ -186,12 +191,15 @@ TEST_F(GaussianFitterTest, NayaniClipped4_guess){ std::vector peaks; int count = fitter.guess_peaks(&peaks, ampData, idxData); - EXPECT_EQ(1,peaks.size()); + ASSERT_EQ(2,peaks.size()); EXPECT_EQ(240,peaks.at(0)->amp); EXPECT_EQ(20, peaks.at(0)->location); EXPECT_NEAR(7.4, peaks.at(0)->fwhm, 1); - EXPECT_EQ(1, count); + EXPECT_EQ(16,peaks.at(1)->amp); + EXPECT_NEAR(29.5, peaks.at(1)->location,.25); + EXPECT_NEAR(8, peaks.at(1)->fwhm, 1); + EXPECT_EQ(2, count); } @@ -276,15 +284,18 @@ TEST_F(GaussianFitterTest, NayaniClipped7_guess){ fitter.noise_level = 9; int count = fitter.guess_peaks(&peaks, ampData, idxData); - ASSERT_EQ(2,peaks.size()); + ASSERT_EQ(3,peaks.size()); EXPECT_EQ(98,peaks.at(0)->amp); EXPECT_EQ(168,peaks.at(1)->amp); EXPECT_EQ(18, peaks.at(0)->location); EXPECT_EQ(31, peaks.at(1)->location); EXPECT_NEAR(5.2, peaks.at(0)->fwhm, 1); EXPECT_NEAR(10, peaks.at(1)->fwhm, 1); + EXPECT_EQ(13,peaks.at(2)->amp); + EXPECT_NEAR(40.5, peaks.at(2)->location,.25); + EXPECT_NEAR(11, peaks.at(2)->fwhm, 1); - EXPECT_EQ(2, count); + EXPECT_EQ(3, count); } @@ -1681,10 +1692,10 @@ TEST_F(GaussianFitterTest, problem_waveform_5_find){ fitter.smoothing_expt(&Data); int count = fitter.find_peaks(&peaks,ampData,idxData, 200); - EXPECT_EQ(2,peaks.size()); + ASSERT_EQ(2,peaks.size()); + EXPECT_NEAR(19.5, peaks.at(0)->location,.25); EXPECT_EQ(186, peaks.at(0)->amp); EXPECT_EQ(14, peaks.at(1)->amp); - EXPECT_NEAR(19.5, peaks.at(0)->location,.25); EXPECT_EQ(30, peaks.at(1)->location); EXPECT_NEAR(5.7, peaks.at(0)->fwhm, 1); EXPECT_NEAR(8, peaks.at(1)->fwhm, 1);