OpalFilter.cpp

Go to the documentation of this file.

//
// Class OpalFilter
//   The class for the Filter Object.
//   A FILTER definition is used to define a filter which can be applied
//   to a 1D histogram in order to get rid of noise.
//
// Copyright (c) 2008 - 2022, Christof Metzger-Kraus
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "Utilities/OpalFilter.h"
 
#include "AbsBeamline/ElementBase.h"
#include "AbstractObjects/OpalData.h"
#include "Attributes/Attributes.h"
#include "Filters/Filters.h"
#include "Physics/Physics.h"
#include "Utilities/OpalException.h"
 
#include "Utility/IpplInfo.h"
 
#include <cmath>
#include <unordered_map>
 
#define NPOINTS_DEFAULT 129
#define NLEFT_DEFAULT 64
#define NRIGHT_DEFAULT 64
#define POLYORDER_DEFAULT 1
 
extern Inform* gmsg;
 
// The attributes of class OpalFilter.
namespace {
    enum {
        TYPE,       // The type of filter
        NFREQ,      // Number of frequencies in fixedFFTLowPass filter
        THRESHOLD,  // Relative threshold for amplitude of frequency in relativeFFTLowPass
        NPOINTS,    // Number of points in Savitzky-Golay filter
        NLEFT,      // Number of points to the left in S-G filter
        NRIGHT,     // Number of points to the right in S-G filter
        POLYORDER,  // Polynomial order in S-G Filter
        SIZE
    };
}
 
OpalFilter::OpalFilter()
    : Definition(
        SIZE, "FILTER",
        "The \"FILTER\" statement defines a 1 dimensional filter to be "
        "applied on histogram."),
      filter_m(0) {
    itsAttr[TYPE] = Attributes::makePredefinedString(
        "TYPE", "Specifies the type of filter.",
        {"SAVITZKY-GOLAY", "FIXEDFFTLOWPASS", "RELATIVEFFTLOWPASS", "STENCIL"});
 
    itsAttr[NFREQ] =
        Attributes::makeReal("NFREQ", "Number of frequencies to use in fixedFFTLowPass filter", 9.);
 
    itsAttr[THRESHOLD] = Attributes::makeReal(
        "THRESHOLD", "Relative threshold for amplitude of frequencies in relativeFFTLowPass filter",
        1.e-6);
 
    itsAttr[NPOINTS] = Attributes::makeReal(
        "NPOINTS", "Number of points in Savitzky-Golay filter", NPOINTS_DEFAULT);
 
    itsAttr[NLEFT] = Attributes::makeReal(
        "NLEFT", "Number of points to the left in Savitzky-Golay filter", NLEFT_DEFAULT);
 
    itsAttr[NRIGHT] = Attributes::makeReal(
        "NRIGHT", "Number of points to the right in Savitzky-Golay filter", NRIGHT_DEFAULT);
 
    itsAttr[POLYORDER] = Attributes::makeReal(
        "POLYORDER", "Polynomial order for local fit-function in Savitzky-Golay filter",
        POLYORDER_DEFAULT);
 
    registerOwnership(AttributeHandler::STATEMENT);
 
    OpalFilter* defFilter = clone("UNNAMED_FILTER");
    defFilter->builtin    = true;
 
    try {
        defFilter->update();
        OpalData::getInstance()->define(defFilter);
    } catch (...) {
        delete defFilter;
    }
}
 
OpalFilter::OpalFilter(const std::string& name, OpalFilter* parent)
    : Definition(name, parent), filter_m(0) {
}
 
OpalFilter::~OpalFilter() {
    if (filter_m)
        delete filter_m;
}
 
bool OpalFilter::canReplaceBy(Object* object) {
    // Can replace only by another WAKE.
    return dynamic_cast<OpalFilter*>(object) != 0;
}
 
OpalFilter* OpalFilter::clone(const std::string& name) {
    return new OpalFilter(name, this);
}
 
void OpalFilter::execute() {
    update();
}
 
OpalFilter* OpalFilter::find(const std::string& name) {
    OpalFilter* filter = dynamic_cast<OpalFilter*>(OpalData::getInstance()->find(name));
 
    if (filter == 0) {
        throw OpalException("OpalFilter::find()", "OpalFilter \"" + name + "\" not found.");
    }
    return filter;
}
 
void OpalFilter::update() {
    // Set default name.
    if (getOpalName().empty())
        setOpalName("UNNAMED_FILTER");
}
 
void OpalFilter::initOpalFilter() {
    if (filter_m == 0) {
        *gmsg << "* ************* F I L T E R "
                 "************************************************************"
              << endl;
        *gmsg << "OpalFilter::initOpalFilterfunction " << endl;
        *gmsg
            << "* "
               "**********************************************************************************"
            << endl;
 
        static const std::unordered_map<std::string, FilterType> stringFilterType_s = {
            {"SAVITZKY-GOLAY", FilterType::SAVITZKYGOLAY},
            {"FIXEDFFTLOWPASS", FilterType::FIXEDFFTLOWPASS},
            {"RELATIVEFFTLOWPASS", FilterType::RELATIVEFFTLOWPASS},
            {"STENCIL", FilterType::STENCIL}};
 
        FilterType type = stringFilterType_s.at(Attributes::getString(itsAttr[TYPE]));
        switch (type) {
            case FilterType::SAVITZKYGOLAY: {
                int num_points       = (int)(Attributes::getReal(itsAttr[NPOINTS]));
                int num_points_left  = (int)(Attributes::getReal(itsAttr[NLEFT]));
                int num_points_right = (int)(Attributes::getReal(itsAttr[NRIGHT]));
                int polynomial_order = std::abs((int)(Attributes::getReal(itsAttr[POLYORDER])));
 
                Inform svg("Savitzky-Golay: ");
                if (num_points_left < 0) {
                    svg << "Number of points to the left negative; using default (" << NLEFT_DEFAULT
                        << ");" << endl;
                    num_points_left = NLEFT_DEFAULT;
                }
                if (num_points_right < 0) {
                    svg << "Number of points to the right negative; using default ("
                        << NRIGHT_DEFAULT << ");" << endl;
                    num_points_right = NRIGHT_DEFAULT;
                }
                if (num_points < num_points_left + num_points_right) {
                    svg << "Total number of points small than sum of the ones to the left and to "
                           "the right plus 1; using default (NLEFT + NRIGHT + 1);"
                        << endl;
                    num_points = num_points_left + num_points_right + 1;
                }
                if (polynomial_order > num_points_left + num_points_right) {
                    svg << "Polynomial order bigger than sum of points to the left and to the "
                           "right; using default (NLEFT + NRIGHT);"
                        << endl;
                    polynomial_order = num_points_left + num_points_right;
                }
 
                filter_m = new SavitzkyGolayFilter(
                    num_points, num_points_left, num_points_right, polynomial_order);
                break;
            }
            case FilterType::FIXEDFFTLOWPASS: {
                filter_m =
                    new FixedFFTLowPassFilter(std::abs((int)(Attributes::getReal(itsAttr[NFREQ]))));
                break;
            }
            case FilterType::RELATIVEFFTLOWPASS: {
                filter_m =
                    new RelativeFFTLowPassFilter(std::abs(Attributes::getReal(itsAttr[THRESHOLD])));
                break;
            }
            case FilterType::STENCIL: {
                filter_m = new StencilFilter();
                break;
            }
            default: {
                filter_m = 0;
                *ippl::Info << "No filter attached" << endl;
                break;
            }
        }
    }
}
 
void OpalFilter::print(std::ostream& os) const {
    os << "* ************* F I L T E R ********************************************************\n"
       << "* FILTER         " << getOpalName() << '\n'
       << "* TYPE           " << Attributes::getString(itsAttr[TYPE]) << '\n'
       << "* NFREQ          " << Attributes::getReal(itsAttr[NFREQ]) << '\n'
       << "* THRESHOLD      " << Attributes::getReal(itsAttr[THRESHOLD]) << '\n'
       << "* NPOINTS        " << Attributes::getReal(itsAttr[NPOINTS]) << '\n'
       << "* NLEFT          " << Attributes::getReal(itsAttr[NLEFT]) << '\n'
       << "* NRIGHT         " << Attributes::getReal(itsAttr[NRIGHT]) << '\n'
       << "* POLYORDER      " << Attributes::getReal(itsAttr[POLYORDER]) << '\n'
       << "* ********************************************************************************** "
       << std::endl;
}