function yFiltM = applyRecursGaussFilter(xM,coeffS,dim)
% function yFiltM = applyRecursGaussFilter(xM,coeffS,dim)
%
% This function applies recursive filter of order 4 to xM along the passed dim.
% coeffS is struture with filter coefficients. 
% For example for LoG filter:
%     scan3M: 3D scanArray (for example, planC{indexS.scan}.scanArray - ctOffset).
%     pixelSizeV = [1,1,3]; % voxel size
%     sigmaVal = 3; % sigma
%     dim = 1; % filter along rows
%     derivativeOrder = 'zero';
%     coeffS.sigmad = sigmaVal / pixelSizeV(dim);
%     coeffS = setGaussOrder(coeffS,derivativeOrder);
%     filt3M = applyRecursGaussFilter(scan3M,coeffS,dim);
%
% APA, 6/12/2018

ykPlusM = zeros(size(xM));
ykMinusM = zeros(size(xM));

switch dim
    case 1
        
        % Causal
        ykPlusM(1,:,:) = coeffS.N0 * xM(1,:,:) + coeffS.N1 * xM(1,:,:) + coeffS.N2 * xM(1,:,:)...
            + coeffS.N3 * xM(1,:,:);
        ykPlusM(2,:,:) = coeffS.N0 * xM(2,:,:) + coeffS.N1 * xM(1,:,:) + coeffS.N2 * xM(1,:,:)...
            + coeffS.N3 * xM(1,:,:);
        ykPlusM(3,:,:) = coeffS.N0 * xM(3,:,:) + coeffS.N1 * xM(2,:,:) + coeffS.N2 * xM(1,:,:)...
            + coeffS.N3 * xM(1,:,:);
        ykPlusM(4,:,:) = coeffS.N0 * xM(4,:,:) + coeffS.N1 * xM(3,:,:) + coeffS.N2 * xM(2,:,:)...
            + coeffS.N3 * xM(1,:,:);
        
        ykPlusM(1,:,:) = ykPlusM(1,:,:) - coeffS.BN1 * xM(1,:,:) - coeffS.BN2 * xM(1,:,:) - coeffS.BN3 * xM(1,:,:)...
            - coeffS.BN4 * xM(1,:,:);
        ykPlusM(2,:,:) = ykPlusM(2,:,:) + -coeffS.D1 * ykPlusM(1,:,:) - coeffS.BN2 * xM(1,:,:) - coeffS.BN3 * xM(1,:,:)...
            - coeffS.BN4 * xM(1,:,:);
        ykPlusM(3,:,:) = ykPlusM(3,:,:)  + -coeffS.D1 * ykPlusM(1,:,:) -+ coeffS.D2 * ykPlusM(2,:,:) - coeffS.BN3 * xM(1,:,:)...
            - coeffS.BN4 * xM(1,:,:);
        ykPlusM(4,:,:) = ykPlusM(4,:,:) + -coeffS.D1 * ykPlusM(1,:,:) -+ coeffS.D2 * ykPlusM(2,:,:) -+ coeffS.D3 * ykPlusM(3,:,:)...
            - coeffS.BN4 * xM(1,:,:);
        
        for ind = 5:size(xM,dim)
            ykPlusM(ind,:,:) = coeffS.N0 * xM(ind,:,:) + coeffS.N1 * xM(ind-1,:,:) + coeffS.N2 * xM(ind-2,:,:)...
                + coeffS.N3 * xM(ind-3,:,:) -+ coeffS.D1 * ykPlusM(ind-1,:,:) -+ coeffS.D2 * ykPlusM(ind-2,:,:) ...
                -+ coeffS.D3 * ykPlusM(ind-3,:,:) -+ coeffS.D4 * ykPlusM(ind-4,:,:);
        end
        
        % Anticausal
        ykMinusM(end,:,:) = coeffS.M1 * xM(end,:,:) + coeffS.M2 * xM(end,:,:) + coeffS.M3 * xM(end,:,:)...
            + coeffS.M4 * xM(end,:,:);
        ykMinusM(end-1,:,:) = coeffS.M1 * xM(end,:,:) + coeffS.M2 * xM(end,:,:) + coeffS.M3 * xM(end,:,:)...
            + coeffS.M4 * xM(end,:,:);
        ykMinusM(end-2,:,:) = coeffS.M1 * xM(end-1,:,:) + coeffS.M2 * xM(end-1,:,:) + coeffS.M3 * xM(end,:,:)...
            + coeffS.M4 * xM(end,:,:);
        ykMinusM(end-3,:,:) = coeffS.M1 * xM(end-2,:,:) + coeffS.M2 * xM(end-2,:,:) + coeffS.M3 * xM(end-1,:,:)...
            + coeffS.M4 * xM(end,:,:);
        
        ykMinusM(end,:,:) = ykMinusM(end,:,:) - coeffS.BM1 * xM(end,:,:) - coeffS.BM2 * xM(end,:,:) - coeffS.BM3 * xM(end,:,:)...
            - coeffS.BM4 * xM(end,:,:);
        ykMinusM(end-1,:,:) = ykMinusM(end-1,:,:) + -coeffS.D1 * ykMinusM(end,:,:) - coeffS.BM2 * xM(end,:,:) - coeffS.BM3 * xM(end,:,:)...
            - coeffS.BM4 * xM(end,:,:);
        ykMinusM(end-2,:,:) = ykMinusM(end-2,:,:) + -coeffS.D1 * ykMinusM(end-1,:,:) +- coeffS.D2 * ykMinusM(end,:,:) - coeffS.BM3 * xM(end,:,:)...
            - coeffS.BM4 * xM(end,:,:);
        ykMinusM(end-3,:,:) = ykMinusM(end-3,:,:) + -coeffS.D1 * ykMinusM(end-2,:,:) +- coeffS.D2 * ykMinusM(end-1,:,:) +- coeffS.D3 * ykMinusM(end,:,:)...
            - coeffS.BM4 * xM(end,:,:);
        
        for ind = size(xM,1)-4:-1:1
            ykMinusM(ind,:,:) = coeffS.M1 * xM(ind+1,:,:) + coeffS.M2 * xM(ind+2,:,:) + coeffS.M3 * xM(ind+3,:,:)...
                + coeffS.M4 * xM(ind+4,:,:) - coeffS.D1 * ykMinusM(ind+1,:,:) - coeffS.D2 * ykMinusM(ind+2,:,:) ...
                - coeffS.D3 * ykMinusM(ind+3,:,:) - coeffS.D4 * ykMinusM(ind+4,:,:);
        end
        
    case 2

        % Causal
        ykPlusM(:,1,:) = coeffS.N0 * xM(:,1,:) + coeffS.N1 * xM(:,1,:) + coeffS.N2 * xM(:,1,:)...
            + coeffS.N3 * xM(:,1,:);
        ykPlusM(:,2,:) = coeffS.N0 * xM(:,2,:) + coeffS.N1 * xM(:,1,:) + coeffS.N2 * xM(:,1,:)...
            + coeffS.N3 * xM(:,1,:);
        ykPlusM(:,3,:) = coeffS.N0 * xM(:,3,:) + coeffS.N1 * xM(:,2,:) + coeffS.N2 * xM(:,1,:)...
            + coeffS.N3 * xM(:,1,:);
        ykPlusM(:,4,:) = coeffS.N0 * xM(:,4,:) + coeffS.N1 * xM(:,3,:) + coeffS.N2 * xM(:,2,:)...
            + coeffS.N3 * xM(:,1,:);
        
        ykPlusM(:,1,:) = ykPlusM(:,1,:) - coeffS.BN1 * xM(:,1,:) - coeffS.BN2 * xM(:,1,:) - coeffS.BN3 * xM(:,1,:)...
            - coeffS.BN4 * xM(:,1,:);
        ykPlusM(:,2,:) = ykPlusM(:,2,:) + -coeffS.D1 * ykPlusM(:,1,:) - coeffS.BN2 * xM(:,1,:) - coeffS.BN3 * xM(:,1,:)...
            - coeffS.BN4 * xM(:,1,:);
        ykPlusM(:,3,:) = ykPlusM(:,3,:) + -coeffS.D1 * ykPlusM(:,2,:) -+ coeffS.D2 * ykPlusM(:,1,:) - coeffS.BN3 * xM(:,1,:)...
            - coeffS.BN4 * xM(:,1,:);
        ykPlusM(:,4,:) = ykPlusM(:,4,:) + -coeffS.D1 * ykPlusM(:,3,:) -+ coeffS.D2 * ykPlusM(:,2,:) -+ coeffS.D3 * ykPlusM(:,1,:)...
            - coeffS.BN4 * xM(:,1,:);
        
        for ind = 5:size(xM,dim)
            ykPlusM(:,ind,:) = coeffS.N0 * xM(:,ind,:) + coeffS.N1 * xM(:,ind-1,:) + coeffS.N2 * xM(:,ind-2,:)...
                + coeffS.N3 * xM(:,ind-3,:) -+ coeffS.D1 * ykPlusM(:,ind-1,:) -+ coeffS.D2 * ykPlusM(:,ind-2,:) ...
                -+ coeffS.D3 * ykPlusM(:,ind-3,:) -+ coeffS.D4 * ykPlusM(:,ind-4,:);
        end
        
        % Anticausal
        ykMinusM(:,end,:) = coeffS.M1 * xM(:,end,:) + coeffS.M2 * xM(:,end,:) + coeffS.M3 * xM(:,end,:)...
            + coeffS.M4 * xM(:,end,:);
        ykMinusM(:,end-1,:) = coeffS.M1 * xM(:,end,:) + coeffS.M2 * xM(:,end,:) + coeffS.M3 * xM(:,end,:)...
            + coeffS.M4 * xM(:,end,:);
        ykMinusM(:,end-2,:) = coeffS.M1 * xM(:,end-1,:) + coeffS.M2 * xM(:,end-1,:) + coeffS.M3 * xM(:,end,:)...
            + coeffS.M4 * xM(:,end,:);
        ykMinusM(:,end-3,:) = coeffS.M1 * xM(:,end-2,:) + coeffS.M2 * xM(:,end-2,:) + coeffS.M3 * xM(:,end-1,:)...
            + coeffS.M4 * xM(:,end,:);
        
        ykMinusM(:,end,:) = ykMinusM(:,end,:) - coeffS.BM1 * xM(:,end,:) - coeffS.BM2 * xM(:,end,:) - coeffS.BM3 * xM(:,end,:)...
            - coeffS.BM4 * xM(:,end,:);
        ykMinusM(:,end-1,:) = ykMinusM(:,end-1,:) + -coeffS.D1 * ykMinusM(:,end,:) - coeffS.BM2 * xM(:,end,:) - coeffS.BM3 * xM(:,end,:)...
            - coeffS.BM4 * xM(:,end,:);
        ykMinusM(:,end-2,:) = ykMinusM(:,end-2,:) + -coeffS.D1 * ykMinusM(:,end-1,:) +- coeffS.D2 * ykMinusM(:,end,:) - coeffS.BM3 * xM(:,end,:)...
            - coeffS.BM4 * xM(:,end,:);
        ykMinusM(:,end-3,:) = ykMinusM(:,end-3,:) + -coeffS.D1 * ykMinusM(:,end-2,:) +- coeffS.D2 * ykMinusM(:,end-1,:) +- coeffS.D3 * ykMinusM(:,end,:)...
            - coeffS.BM4 * xM(:,end,:);
        
        for ind = size(xM,dim)-4:-1:1
            ykMinusM(:,ind,:) = coeffS.M1 * xM(:,ind+1,:) + coeffS.M2 * xM(:,ind+2,:) + coeffS.M3 * xM(:,ind+3,:)...
                + coeffS.M4 * xM(:,ind+4,:) - coeffS.D1 * ykMinusM(:,ind+1,:) - coeffS.D2 * ykMinusM(:,ind+2,:) ...
                - coeffS.D3 * ykMinusM(:,ind+3,:) - coeffS.D4 * ykMinusM(:,ind+4,:);
        end

        
    case 3
        
        % Causal
        ykPlusM(:,:,1) = coeffS.N0 * xM(:,:,1) + coeffS.N1 * xM(:,:,1) + coeffS.N2 * xM(:,:,1)...
            + coeffS.N3 * xM(:,:,1);
        ykPlusM(:,:,2) = coeffS.N0 * xM(:,:,2) + coeffS.N1 * xM(:,:,1) + coeffS.N2 * xM(:,:,1)...
            + coeffS.N3 * xM(:,:,1);
        ykPlusM(:,:,3) = coeffS.N0 * xM(:,:,3) + coeffS.N1 * xM(:,:,2) + coeffS.N2 * xM(:,:,1)...
            + coeffS.N3 * xM(:,:,1);
        ykPlusM(:,:,4) = coeffS.N0 * xM(:,:,4) + coeffS.N1 * xM(:,:,3) + coeffS.N2 * xM(:,:,2)...
            + coeffS.N3 * xM(:,:,1);
        
        ykPlusM(:,:,1) = ykPlusM(:,:,1) - coeffS.BN1 * xM(:,:,1) - coeffS.BN2 * xM(:,:,1) - coeffS.BN3 * xM(:,:,1)...
            - coeffS.BN4 * xM(:,:,1);
        ykPlusM(:,:,2) = ykPlusM(:,:,2) + -coeffS.D1 * ykPlusM(:,:,1) - coeffS.BN2 * xM(:,:,1) - coeffS.BN3 * xM(:,:,1)...
            - coeffS.BN4 * xM(:,:,1);
        ykPlusM(:,:,3) = ykPlusM(:,:,3) + -coeffS.D1 * ykPlusM(:,:,1) -+ coeffS.D2 * ykPlusM(:,:,2) - coeffS.BN3 * xM(:,:,1)...
            - coeffS.BN4 * xM(:,:,1);
        ykPlusM(:,:,4) = ykPlusM(:,:,4) + -coeffS.D1 * ykPlusM(:,:,1) -+ coeffS.D2 * ykPlusM(:,:,2) -+ coeffS.D3 * ykPlusM(:,:,3)...
            - coeffS.BN4 * xM(:,:,1);
        
        for ind = 5:size(xM,dim)
            ykPlusM(:,:,ind) = coeffS.N0 * xM(:,:,ind) + coeffS.N1 * xM(:,:,ind-1) + coeffS.N2 * xM(:,:,ind-2)...
                + coeffS.N3 * xM(:,:,ind-3) -+ coeffS.D1 * ykPlusM(:,:,ind-1) -+ coeffS.D2 * ykPlusM(:,:,ind-2) ...
                -+ coeffS.D3 * ykPlusM(:,:,ind-3) -+ coeffS.D4 * ykPlusM(:,:,ind-4);
        end
        
        % Anticausal
        ykMinusM(:,:,end) = coeffS.M1 * xM(:,:,end) + coeffS.M2 * xM(:,:,end) + coeffS.M3 * xM(:,:,end)...
            + coeffS.M4 * xM(:,:,end);
        ykMinusM(:,:,end-1) = coeffS.M1 * xM(:,:,end) + coeffS.M2 * xM(:,:,end) + coeffS.M3 * xM(:,:,end)...
            + coeffS.M4 * xM(:,:,end);
        ykMinusM(:,:,end-2) = coeffS.M1 * xM(:,:,end-1) + coeffS.M2 * xM(:,:,end-1) + coeffS.M3 * xM(:,:,end)...
            + coeffS.M4 * xM(:,:,end);
        ykMinusM(:,:,end-3) = coeffS.M1 * xM(:,:,end-2) + coeffS.M2 * xM(:,:,end-2) + coeffS.M3 * xM(:,:,end-1)...
            + coeffS.M4 * xM(:,:,end);
        
        ykMinusM(:,:,end) = ykMinusM(:,:,end) - coeffS.BM1 * xM(:,:,end) - coeffS.BM2 * xM(:,:,end) - coeffS.BM3 * xM(:,:,end)...
            - coeffS.BM4 * xM(:,:,end);
        ykMinusM(:,:,end-1) = ykMinusM(:,:,end-1) + -coeffS.D1 * ykMinusM(:,:,end) - coeffS.BM2 * xM(:,:,end) - coeffS.BM3 * xM(:,:,end)...
            - coeffS.BM4 * xM(:,:,end);
        ykMinusM(:,:,end-2) = ykMinusM(:,:,end-2) + -coeffS.D1 * ykMinusM(:,:,end-1) +- coeffS.D2 * ykMinusM(:,:,end) - coeffS.BM3 * xM(:,:,end)...
            - coeffS.BM4 * xM(:,:,end);
        ykMinusM(:,:,end-3) = ykMinusM(:,:,end-3) + -coeffS.D1 * ykMinusM(:,:,end-2) +- coeffS.D2 * ykMinusM(:,:,end-1) +- coeffS.D3 * ykMinusM(:,:,end)...
            - coeffS.BM4 * xM(:,:,end);
        
        for ind = size(xM,dim)-4:-1:1
            ykMinusM(:,:,ind) = coeffS.M1 * xM(:,:,ind+1) + coeffS.M2 * xM(:,:,ind+2) + coeffS.M3 * xM(:,:,ind+3)...
                + coeffS.M4 * xM(:,:,ind+4) - coeffS.D1 * ykMinusM(:,:,ind+1) - coeffS.D2 * ykMinusM(:,:,ind+2) ...
                - coeffS.D3 * ykMinusM(:,:,ind+3) - coeffS.D4 * ykMinusM(:,:,ind+4);
        end
        


end


yFiltM = ykPlusM + ykMinusM;

