Stencil.cpp

Go to the documentation of this file.

#include "Filters/Stencil.h"

using namespace std;

void StencilFilter::apply(vector<double> &LineDensity) {
    int size = LineDensity.size();
    vector<double> temp1(size + 8, 0.0);
    vector<double> temp2(5, 0.0);

    for(int i = 0; i < size; ++i)
        temp1[i+4] = LineDensity[i];

    for(int i = 0; i < 4; ++ i) {
        temp1[3-i] = temp1[5+i];
        temp1[size+5+i] = temp1[size+3-i];
    }

    for(int i = 0; i < size; ++i) {
        temp2[i%5] = (7. * (temp1[i] + temp1[i+8]) + 24. * (temp1[i+2] + temp1[i+6]) + 34.*temp1[i+4]) / 96.;
        temp1[i] = temp2[(i+1)%5];
    }

    temp1[2] = temp1[6];
    temp1[3] = temp1[5];
    temp1[size+5] = temp1[size+1];
    temp1[size+4] = temp1[size+2];

    for(int i = 0; i < size; ++i)
        LineDensity[i] = (7. * (temp1[i+2] + temp1[i+6]) + 24. * (temp1[i+3] + temp1[i+5]) + 34.*temp1[i+4]) / 96.;
}

void StencilFilter::calc_derivative(vector<double> &LineDensity, const double &h) {
    vector<double> temp(LineDensity.begin(), LineDensity.end());
    const double N = LineDensity.size();
    LineDensity[0] = (-25.*temp[0] + 48.*temp[1] - 36.*temp[2] + 16.*temp[3] + temp[4]) / (12.*h);
    LineDensity[1] = (-3.*temp[0] - 10.*temp[1] + 18.*temp[2] - 6.*temp[3] + temp[4]) / (12.*h);
    for(int i = 2; i < N - 2; ++ i) {
        LineDensity[i] = (temp[i-2] - 8.*temp[i-1] + 8.*temp[i+1] - temp[i+2]) / (12.*h);
    }
    LineDensity[N-2] = (-temp[N-5] + 6.*temp[N-4] - 18.*temp[N-3] + 10.*temp[N-2] + 3.*temp[N-1]) / (12.*h);
    LineDensity[N-1] = (3.*temp[N-5] - 16.*temp[N-4] + 36.*temp[N-3] - 48.*temp[N-2] + 25.*temp[N-1]) / (12.*h);
    apply(LineDensity);
}