function Out = RTtrilateral(image, numIterations, sigmaS, sigmaR, sigmaI, sigmaP, NN, returnAllIts,...
                            dispEvery, toPause, handles)
%RTTRILATERAL Returns the image denoised by trilateral filtering,
%  it is integrated with the RT denoising demo
%Images are given and returned in uint8 array format.
%Syntax: Out = RTtrilateral(image, sigmaS, sigmaR, sigmaI, sigmaP, N, returnAllIts,...
%                           dispEvery, toPause, handles)
%  IMAGE is the input image
%  SIGMAS is a positive double with default value = 5 and range [0.3 10]
%  SIGMAR is a positive double with default value = 20 and range [10 100]
%  SIGMAI is a positive double with default value = 30 and range [20 70]
%  SIGMAP is a positive double with default value = 30 and range [20 70]
%  N is a positive integer which gives the size of the window used for filtering
%    its typical range is [1 10]
%    --n.b. N = 1 means 3x3, N=2 means 5x5, etc.
%  NUMITERATIONS is a positive integer with default = 1, range [1 5]
%  RETURNALLITS is a switch: when set to 1, the result after each iteration is saved
%    and a cell array of length NUMITERATIONS containing each result is returned.;
%    when set to 0, only the final result is returned.
%  DISPEVERY is a positive integer specifying how many iterations the
%    algorithm is allowed to run before the displayed result is updated.
%    e.g. set it to 1 to see the result of each iteration
%         set it to NUMITERATIONS to see only the final result
%  TOPAUSE is a switch: if set to 1, the program waits for the user to
%    press a key after each time a new iteration result is displayed
%  HANDLES is the handles structure of the main RT GUI window


%initialization
[numRows, numCols] = size(image);
avgitNum = 0;
avgitDen = 0;

%create combined impulsivity lookup table to speed calculation
indices = cumsum(ones(255*8+1,1));
w_PLookup = 1-exp(((indices-1).^2)./(-8*sigmaP^2));

for k = 1:numIterations
    %initialize iteration
    beginTime = clock;
    numerator = zeros(numRows, numCols);
    denominator = zeros(numRows, numCols);
    imageK = double(image);    
    
    %calculate matrices of neighbors in the given direction  ("reflect" at borders of image)
    N = [imageK(2,:); imageK(1:numRows-1, :)];
    S = [imageK(2:numRows, :); imageK(numRows-1, :)];
    E = [imageK(:, 2:numCols), imageK(:, numCols-1)];
    W = [imageK(:,2), imageK(:, 1:numCols-1)];
    NW = [[imageK(2,2), imageK(2,1:numCols-1)];...
          [imageK(1:numRows-1,2), imageK(1:numRows-1, 1:numCols-1)]];
    SE = [[imageK(2:numRows, 2:numCols), imageK(2:numRows, numCols-1)];...
          [imageK(numRows-1, 2:numCols), imageK(numRows-1, numCols-1)]];
    NE = [[imageK(2, 2:numCols), imageK(2, numCols-1)];...
          [imageK(1:numRows-1, 2:numCols), imageK(1:numRows-1, numCols-1)]];
    SW = [[imageK(2:numRows, 2), imageK(2:numRows, 1:numCols-1)]; ...
          [imageK(numRows-1, 2), imageK(numRows-1, 1:numCols-1)]];
    
    %calculate differences ordered by rank
    loopDiff = cat(3, abs(imageK-N), abs(imageK-NE), abs(imageK-E), abs(imageK-SE), ...
                      abs(imageK-S), abs(imageK-SW), abs(imageK-W), abs(imageK-NW)); 
    sorted = sort(loopDiff, 3);
    
    %calculate ROAD and exponent of the impulsive weight
    ROAD = sum(sorted(:, :, 1:4), 3); 
    powW_I = (ROAD.^2)./(-2*sigmaI^2);
    
    %calculate weighted mean-- the numerator and denominator are calculated by accumulation
    for m = -NN:NN
        %enable interruption of method by closing the result window
        drawnow
        handles = guidata(handles.figure1);
        if (handles.interrupt)
            Out = 0;
            return
        end
        for n = -NN:NN           
            if (m == 0 & n == 0)  
                %calculate the contribution of the point in question  
                tempW_P = w_PLookup(ROAD+ROAD+1);
                tempW = exp(powW_I.*tempW_P);
                numerator = numerator + tempW.*imageK;
                denominator = denominator + tempW;
                
            else         
                %calculate the contributions of neighboring points
                
                %the 'point in question'
                XpartRows = 1-m*(m<0):numRows - m*(m>0);
                XpartCols = 1-n*(n<0):numCols - n*(n>0);
                
                %the neighboring point m down and n to the right
                YpartRows = 1 + m*(m>0): numRows + m*(m<0);
                YpartCols = 1 + n*(n>0): numCols + n*(n<0);
                Ypart = imageK(YpartRows, YpartCols);
                                
                %calculate the exponent of the spatial weight              
                spatialDiffSquared = (m^2 + n^2);              
                tempPowW_S = -spatialDiffSquared/(2*sigmaS^2);
                
                %calculate the exponent of the radiometric weight
                radioDiff = imageK(XpartRows, XpartCols) - Ypart;     
                tempPowW_R = (radioDiff.^2)./(-2*sigmaR^2);
                
                %calculate the exponent of the impulsive weight
                tempPowW_I = powW_I(YpartRows, YpartCols);              
                              
                %calculate combined impulsivity
                tempW_P = w_PLookup(ROAD(YpartRows, YpartCols)+...
                                    ROAD(XpartRows, XpartCols) +1);
                               
                %calculate total weight
                tempW = zeros(numRows, numCols);
                tempW(XpartRows, XpartCols) = ...
                    exp(tempPowW_S + tempPowW_R + ((tempPowW_I-tempPowW_R).*tempW_P));
                
                %calculate contribution of neighboring point
                Y = zeros(numRows, numCols);
                Y(XpartRows, XpartCols) = Ypart;
                denominator = denominator + tempW;
                numerator = numerator + tempW .* Y;
            end           
        end
    end
    %calculate new image
    image = numerator./denominator;

    if (returnAllIts)
        images{k} = uint8(round(image));  %should be round
    end
    
    %display results of iteration
    endTime = clock;
    avgitNum = avgitNum + etime(endTime, beginTime);
    avgitDen = avgitDen + 1;
    avgit = avgitNum./avgitDen;
    drawnow
    handles = guidata(handles.figure1);
    if (handles.interrupt)
        Out = 0;
        return
    end
    RTiterDisp(image, k, numIterations, avgit, dispEvery, toPause, handles);
end
if (returnAllIts)
    Out = images;
else
    Out = uint8(round(image));
end