Feature/bilateral filter

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

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 */

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_ */

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<KTSliceHeader>().GetSliceNumber();
+
+        KTBilateralFilteredFSDataFFTW& fsData = data->Of<KTBilateralFilteredFSDataFFTW>();
+        unsigned nComponents = fsData.GetNComponents();
+
+        if (! fWriter->OpenAndVerifyFile()) return;
+
+        for (unsigned iChannel=0; iChannel<nComponents; iChannel++)
+        {
+            const KTFrequencySpectrumFFTW* spectrum = fsData.GetSpectrumFFTW(iChannel);
+            if (spectrum != NULL)
+            {
+                stringstream conv;
+                conv << "histBLFFSfftw_" << sliceNumber << "_" << iChannel;
+                string histName;
+                conv >> 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;

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

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
     )

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 <iostream>
+#include <algorithm>
+
+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; iComponent<nComponents; ++iComponent)
+        {
+            KTFrequencySpectrumFFTW* newSpectrum = Filter(fsData.GetSpectra(iComponent));
+
+            if (newSpectrum == NULL)
+            {
+                KTERROR(gclog, "Bilateral filter of spectrum " << iComponent << " failed for some reason. Continuing processing.");
+                continue;
+            }
+
+            KTDEBUG(gclog, "Computed bilateral filter");
+            newData.SetSpectrum(newSpectrum, iComponent);
+        }
+        KTINFO(gclog, "Completed bilateral filter of " << nComponents << " frequency spectra (FFTW)");
+
+        return true;
+    }
+
+    KTFrequencySpectrumFFTW* KTBilateralFilter::Filter(const KTMultiFSFFTW * frequencySpectrum) const
+    {
+        KTDEBUG(gclog, "Creating new FS for filtered data");
+
+        const int nTimeBins = frequencySpectrum->size();
+        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 */