diff --git a/Source/Data/CMakeLists.txt b/Source/Data/CMakeLists.txt index 441d34e5d..dabd31bf4 100644 --- a/Source/Data/CMakeLists.txt +++ b/Source/Data/CMakeLists.txt @@ -4,6 +4,7 @@ set (DATA_HEADERFILES SpectrumAnalysis/KTAmplitudeDistribution.hh SpectrumAnalysis/KTAnalyticAssociateData.hh + SpectrumAnalysis/KTBilateralFilteredData.hh SpectrumAnalysis/KTCluster1DData.hh SpectrumAnalysis/KTCorrelationData.hh SpectrumAnalysis/KTCorrelationTSData.hh @@ -67,6 +68,7 @@ set (DATA_HEADERFILES set (DATA_SOURCEFILES SpectrumAnalysis/KTAmplitudeDistribution.cc SpectrumAnalysis/KTAnalyticAssociateData.cc + SpectrumAnalysis/KTBilateralFilteredData.cc SpectrumAnalysis/KTCluster1DData.cc SpectrumAnalysis/KTCorrelationData.cc SpectrumAnalysis/KTCorrelationTSData.cc diff --git a/Source/Data/SpectrumAnalysis/KTBilateralFilteredData.cc b/Source/Data/SpectrumAnalysis/KTBilateralFilteredData.cc new file mode 100644 index 000000000..33abd2310 --- /dev/null +++ b/Source/Data/SpectrumAnalysis/KTBilateralFilteredData.cc @@ -0,0 +1,41 @@ +/* + * KTBilateralFilteredData.cc + * + * Created on: Mar 7, 2018 + * Author: buzinsky + */ + +#include "KTBilateralFilteredData.hh" + +namespace Katydid +{ + const std::string KTBilateralFilteredFSDataFFTW::sName("blfed-fs-fftw"); + + KTBilateralFilteredFSDataFFTW::KTBilateralFilteredFSDataFFTW() + {} + KTBilateralFilteredFSDataFFTW::~KTBilateralFilteredFSDataFFTW() + {} + + KTBilateralFilteredFSDataFFTW& KTBilateralFilteredFSDataFFTW::SetNComponents(unsigned components) + { + unsigned oldSize = fSpectra.size(); + if (components < oldSize) + { + for (unsigned iComponent = components; iComponent < oldSize; ++iComponent) + { + delete fSpectra[iComponent]; + } + } + fSpectra.resize(components); + if (components > oldSize) + { + for (unsigned iComponent = oldSize; iComponent < components; ++iComponent) + { + fSpectra[iComponent] = NULL; + } + } + return *this; + } + + +} /* namespace Katydid */ diff --git a/Source/Data/SpectrumAnalysis/KTBilateralFilteredData.hh b/Source/Data/SpectrumAnalysis/KTBilateralFilteredData.hh new file mode 100644 index 000000000..44e095074 --- /dev/null +++ b/Source/Data/SpectrumAnalysis/KTBilateralFilteredData.hh @@ -0,0 +1,30 @@ +/** + @file KTBilateralFilteredData.hh + @brief Contains KTBilateralFilteredFSDataPolar and KTBilateralFilteredFSDataFFTW + @details + @author: N. Buzinsky + @date: Mar 7, 2018 + */ + +#ifndef KTBILATERALFILTEREDDATA_HH_ +#define KTBILATERALFILTEREDDATA_HH_ + +#include "KTFrequencySpectrumDataFFTW.hh" + +namespace Katydid +{ + class KTBilateralFilteredFSDataFFTW : public KTFrequencySpectrumDataFFTWCore, public Nymph::KTExtensibleData< KTBilateralFilteredFSDataFFTW > + { + public: + KTBilateralFilteredFSDataFFTW(); + virtual ~KTBilateralFilteredFSDataFFTW(); + + KTBilateralFilteredFSDataFFTW& SetNComponents(unsigned components); + + public: + static const std::string sName; + + }; + +} /* namespace Katydid */ +#endif /* KTBILATERALFILTEREDDATA_HH_ */ diff --git a/Source/IO/BasicROOTFileWriter/KTBasicROOTTypeWriterSpectrumAnalysis.cc b/Source/IO/BasicROOTFileWriter/KTBasicROOTTypeWriterSpectrumAnalysis.cc index 89052e536..767f8eb6e 100644 --- a/Source/IO/BasicROOTFileWriter/KTBasicROOTTypeWriterSpectrumAnalysis.cc +++ b/Source/IO/BasicROOTFileWriter/KTBasicROOTTypeWriterSpectrumAnalysis.cc @@ -9,6 +9,7 @@ #include "KT2ROOT.hh" #include "KTAnalyticAssociateData.hh" +#include "KTBilateralFilteredData.hh" #include "KTCorrelationData.hh" #include "KTCorrelationTSData.hh" #include "KTFrequencySpectrumPolar.hh" @@ -77,6 +78,7 @@ namespace Katydid fWriter->RegisterSlot("wv-dist", this, &KTBasicROOTTypeWriterSpectrumAnalysis::WriteWignerVilleDataDistribution); fWriter->RegisterSlot("wv-2d", this, &KTBasicROOTTypeWriterSpectrumAnalysis::WriteWV2DData); fWriter->RegisterSlot("kd-tree-ss", this, &KTBasicROOTTypeWriterSpectrumAnalysis::WriteKDTreeSparseSpectrogram); + fWriter->RegisterSlot("blf-fs-fftw", this, &KTBasicROOTTypeWriterSpectrumAnalysis::WriteBilateralFilteredFSDataFFTW); #ifdef ENABLE_TUTORIAL fWriter->RegisterSlot("lpf-fs-polar", this, &KTBasicROOTTypeWriterSpectrumAnalysis::WriteLowPassFilteredFSDataPolar); fWriter->RegisterSlot("lpf-fs-fftw", this, &KTBasicROOTTypeWriterSpectrumAnalysis::WriteLowPassFilteredFSDataFFTW); @@ -756,6 +758,36 @@ namespace Katydid return; } + void KTBasicROOTTypeWriterSpectrumAnalysis::WriteBilateralFilteredFSDataFFTW(Nymph::KTDataPtr data) + { + if (! data) return; + + uint64_t sliceNumber = data->Of().GetSliceNumber(); + + KTBilateralFilteredFSDataFFTW& fsData = data->Of(); + unsigned nComponents = fsData.GetNComponents(); + + if (! fWriter->OpenAndVerifyFile()) return; + + for (unsigned iChannel=0; iChannel> histName; + TH1D* powerSpectrum = KT2ROOT::CreateMagnitudeHistogram(spectrum, histName); + powerSpectrum->SetDirectory(fWriter->GetFile()); + powerSpectrum->Write(); + KTDEBUG(publog, "Histogram <" << histName << "> written to ROOT file"); + } + } + return; + } + + #ifdef ENABLE_TUTORIAL void KTBasicROOTTypeWriterSpectrumAnalysis::WriteLowPassFilteredFSDataPolar(Nymph::KTDataPtr data) { @@ -786,6 +818,7 @@ namespace Katydid return; } + void KTBasicROOTTypeWriterSpectrumAnalysis::WriteLowPassFilteredFSDataFFTW(Nymph::KTDataPtr data) { if (! data) return; diff --git a/Source/IO/BasicROOTFileWriter/KTBasicROOTTypeWriterSpectrumAnalysis.hh b/Source/IO/BasicROOTFileWriter/KTBasicROOTTypeWriterSpectrumAnalysis.hh index 451c03963..d0bad2683 100644 --- a/Source/IO/BasicROOTFileWriter/KTBasicROOTTypeWriterSpectrumAnalysis.hh +++ b/Source/IO/BasicROOTFileWriter/KTBasicROOTTypeWriterSpectrumAnalysis.hh @@ -89,6 +89,13 @@ namespace Katydid public: void WriteKDTreeSparseSpectrogram(Nymph::KTDataPtr data); + //************************ + // Bilateral Filter Data + //************************ + + public: + void WriteBilateralFilteredFSDataFFTW(Nymph::KTDataPtr data); + #ifdef ENABLE_TUTORIAL //************************ // LPF Tutorial Data diff --git a/Source/SpectrumAnalysis/CMakeLists.txt b/Source/SpectrumAnalysis/CMakeLists.txt index 53542ad99..0083eb113 100644 --- a/Source/SpectrumAnalysis/CMakeLists.txt +++ b/Source/SpectrumAnalysis/CMakeLists.txt @@ -49,12 +49,14 @@ if (FFTW_FOUND) set (SPECTRUMANALYSIS_HEADERFILES ${SPECTRUMANALYSIS_HEADERFILES} KTAnalyticAssociator.hh + KTBilateralFilter.hh KTConvolution.hh KTWignerVille.hh ) set (SPECTRUMANALYSIS_SOURCEFILES ${SPECTRUMANALYSIS_SOURCEFILES} KTAnalyticAssociator.cc + KTBilateralFilter.cc KTConvolution.cc KTWignerVille.cc ) diff --git a/Source/SpectrumAnalysis/KTBilateralFilter.cc b/Source/SpectrumAnalysis/KTBilateralFilter.cc new file mode 100644 index 000000000..494d57c38 --- /dev/null +++ b/Source/SpectrumAnalysis/KTBilateralFilter.cc @@ -0,0 +1,149 @@ +/* + * KTBilateralFilter.cc + * + * Created on: Mar 7, 2018 + * Author: N. Buzinsky + */ + +#include "KTBilateralFilter.hh" + +#include "KTFrequencySpectrumDataFFTW.hh" +#include "KTFrequencySpectrumFFTW.hh" + +#include "KTBilateralFilteredData.hh" + +#include "KTLogger.hh" +#include "KTMath.hh" + +#include "KTMultiFSDataFFTW.hh" +#include "param.hh" + +#include +#include + +namespace Katydid +{ + KTLOGGER(gclog, "KTBilateralFilter"); + + // Register the processor + KT_REGISTER_PROCESSOR(KTBilateralFilter, "bilateral-filter"); + + KTBilateralFilter::KTBilateralFilter(const std::string& name) : + KTProcessor(name), + fSigmaPixels(2.), + fSigmaRange(0.05), + fFSFFTWSignal("fs-fftw", this), + fFSFFTWSlot("str-fs-fftw", this, &KTBilateralFilter::Filter, &fFSFFTWSignal) + { + } + + KTBilateralFilter::~KTBilateralFilter() + { + } + + bool KTBilateralFilter::Configure(const scarab::param_node* node) + { + if (node == NULL) return false; + + SetSigmaPixels(node->get_value("sigma-pixels", GetSigmaPixels())); + SetSigmaRange(node->get_value("sigma-range", GetSigmaRange())); + + return true; + } + + bool KTBilateralFilter::Filter(KTMultiFSDataFFTW& fsData) + { + const int nComponents = fsData.GetNComponents(); + + KTFrequencySpectrumDataFFTW& newData = fsData.Of< KTFrequencySpectrumDataFFTW >().SetNComponents(nComponents); + //KTBilateralFilteredFSDataFFTW& newData = fsData.Of< KTBilateralFilteredFSDataFFTW>().SetNComponents(nComponents); + + if (fSigmaRange == 0 || fSigmaPixels == 0) + { + KTERROR(gclog, "Bilateral filter has zero for SigmaRange"); + } + + for (unsigned iComponent=0; iComponentsize(); + const int nFrequencyBins = (*frequencySpectrum)(0)->size(); + + const int centerTimeBin = (nTimeBins%2) ? ((nTimeBins-1)/2): (nTimeBins - 2)/2; + + const int nSigmaPixels = 3; + + int timeRange[2] = { std::max(0, centerTimeBin - 3* int(fSigmaPixels)), std::min( nTimeBins , centerTimeBin + 3 * int(fSigmaRange)) }; + int frequencyRange[2]; + + KTFrequencySpectrumFFTW* newSpectrum = new KTFrequencySpectrumFFTW(nFrequencyBins, (*frequencySpectrum)(0)->GetRangeMin(), (*frequencySpectrum)(0)->GetRangeMax()); + + fftw_complex filterWeightNumerator; + double filterWeightDenominator; + + fftw_complex pixelPower[2]; + double pixelWeight; + + for (int iCenterFrequencyBin = 0; iCenterFrequencyBin < nFrequencyBins; ++iCenterFrequencyBin) + { + frequencyRange[0] = std::max(0, iCenterFrequencyBin - 3 * int(fSigmaPixels)); + frequencyRange[1] = std::min( nFrequencyBins, iCenterFrequencyBin + 3 * int(fSigmaPixels)); + + filterWeightNumerator[0] = 0; + filterWeightNumerator[1] = 0; + filterWeightDenominator = 0; + + for(int iTimeBin = timeRange[0] ; iTimeBin < timeRange[1]; ++iTimeBin) + { + for(int iFrequencyBin = frequencyRange[0]; iFrequencyBin < frequencyRange[1]; ++iFrequencyBin) + { + pixelPower[0][0] = (*(*frequencySpectrum)(centerTimeBin))(iCenterFrequencyBin)[0]; + pixelPower[0][1] = (*(*frequencySpectrum)(centerTimeBin))(iCenterFrequencyBin)[1]; + pixelPower[1][0] = (*(*frequencySpectrum)(iTimeBin))(iFrequencyBin)[0]; + pixelPower[1][1] = (*(*frequencySpectrum)(iTimeBin))(iFrequencyBin)[1]; + + pixelWeight = GaussianWeightRange( pixelPower[0] , pixelPower[1] ) * GaussianWeightPixels(centerTimeBin, iCenterFrequencyBin, iTimeBin, iFrequencyBin ); + + filterWeightNumerator[0] += pixelWeight * pixelPower[1][0]; + filterWeightNumerator[1] += pixelWeight * pixelPower[1][1]; + filterWeightDenominator += pixelWeight; + } + } + + (*newSpectrum)(iCenterFrequencyBin)[0] = filterWeightNumerator[0] / filterWeightDenominator; + (*newSpectrum)(iCenterFrequencyBin)[1] = filterWeightNumerator[1] / filterWeightDenominator; + } + + return newSpectrum; + } + + double KTBilateralFilter::GaussianWeightRange(const fftw_complex &I1, const fftw_complex &I2) const + { + return exp( - (pow(I1[0] - I2[0], 2.) + pow(I1[1] - I2[1], 2.)) / (2. * pow(fSigmaRange, 2.) ) ); + } + + double KTBilateralFilter::GaussianWeightPixels(const double &i, const double &j,const double &k,const double &l) const + { + return exp( - ( pow(i-k, 2.) + pow(j - l, 2. )) / (2. * pow(fSigmaPixels, 2.) ) ); + } + +} /* namespace Katydid */ diff --git a/Source/SpectrumAnalysis/KTBilateralFilter.hh b/Source/SpectrumAnalysis/KTBilateralFilter.hh new file mode 100644 index 000000000..0905826c0 --- /dev/null +++ b/Source/SpectrumAnalysis/KTBilateralFilter.hh @@ -0,0 +1,88 @@ +/** + @file KTBilateralFilter.hh + @brief Contains KTBilateralFilter + @details Applies a direct bilateral filter to time-frequency plane + @author: N. Buzinsky + @date: Mar 7, 2018 + */ + +#ifndef KTBILATERALFILTER_HH_ +#define KTBILATERALFILTER_HH_ + +#include "KTProcessor.hh" + +#include "KTMemberVariable.hh" +#include "KTSlot.hh" +#include "KTMultiFSDataFFTW.hh" + +#include + +namespace scarab +{ + class param_node; +} + +namespace Katydid +{ + class KTFrequencySpectrumDataFFTW; + class KTFrequencySpectrumFFTW; + + /*! + @class KTBilateralFilter + @author N. Buzinsky + + @brief Applies a bilateral filter + + @details + Approximate implementation of a bilateral filter processor: https://en.wikipedia.org/wiki/Bilateral_filter + + Configuration name: "bilateral-filter" + + Available configuration values: + - "sigma-pixels": double -- width of gaussian for averaging, in pixels + - "sigma-range": double -- width of gaussian for averaging, in terms of differences in power + + Slots: + - "str-fs-fftw": void (KTDataPtr) -- Applies a bilateral filter; Requires KTMultiFSDataFFTW; Adds KTBilateralFilteredFSDataFFTW; Emits signal "fs-fftw" + + Signals: + - "fs-fftw": void (KTDataPtr) -- Emitted upon bilateral filtering; Guarantees KTBilateralFilteredFSDataFFTW. + */ + class KTBilateralFilter : public Nymph::KTProcessor + { + public: + KTBilateralFilter(const std::string& name = "bilateral-filter"); + virtual ~KTBilateralFilter(); + + bool Configure(const scarab::param_node* node); + + MEMBERVARIABLE(double, SigmaPixels); + MEMBERVARIABLE(double, SigmaRange); + + public: + bool Filter(KTMultiFSDataFFTW& fsData); + + KTFrequencySpectrumFFTW* Filter(const KTMultiFSFFTW * frequencySpectrum) const; + + //*************** + // Signals + //*************** + + private: + Nymph::KTSignalData fFSFFTWSignal; + + //*************** + // Slots + //*************** + + private: + Nymph::KTSlotDataOneType< KTMultiFSDataFFTW > fFSFFTWSlot; + double GaussianWeightRange(const fftw_complex &I1, const fftw_complex &I2) const; //Weight due to difference in powers between pixels + double GaussianWeightPixels(const double &i, const double &j,const double &k,const double &l) const; //Weight between pixels at (i,j) and (k,l) due to spacing + + + }; + +} + /* namespace Katydid */ +#endif /* KTBILATERALFILTER_HH_ */