function profile = clutter_bresenham2_BJTU(a, b, mat, mat2, res, br_height)
% This function calculates the profile by using the dual Bresenham
% algorithm, which can be found in Internet or literature.
% INPUT:
% a - start value with [column, row]
% b - final value with [column, row]
% map - raster database map of height [m]
% map2 - raster database of various clutters
% res - resolution of the raster database [?m,?m]
% br_height - crossing bridge height relative to the ground (track) [m]
%
% OUTPUT:
% profile - Matrix with the profile bodies
%  Column 1: x-values
%  Column 2: y-values
%  Column 3: Values of mat
%  Column 4: Distance to the 1st profile point
% 
% Author: Dr. Ke Guan
% Date: 20.07.2016

dx = b(1) - a(1);
dy = b(2) - a(2);

adx = abs(dx);  ady = abs(dy);      
sdx = sign(dx); sdy = sign(dy);     

if adx > ady                
    pdx = sdx; pdy = 0;     
    ddx = sdx; ddy = sdy;   
    es  = ady; el  = adx;   
else                        
    pdx = 0;   pdy = sdy;   
    ddx = sdx; ddy = sdy;   
    es  = adx; el  = ady;   
end

profile = [el + 1; 4];

x = a(1);
y = a(2);

profile(1, 1) = x;
profile(1, 2) = y;
profile(1, 3) = mat(y, x);
profile(1, 4) = 0;

error = el / 2;

for i = 1:el
    error = error - es;      
    
    if error < 0   
        error = error + el;  
        x = x + ddx;          
        y = y + ddy;
    else         
        x = x + pdx;          
        y = y + pdy;
    end
      
    profile(i + 1, 1) = x;  % Coordinate point (x, y)
    profile(i + 1, 2) = y;
       
    if ismember(mat2(y,x),[71])   %1st bridge
    profile(i+1,3) = mat(y,x) + br_height(1);  
        elseif ismember(mat2(y,x),[72])   %2nd bridge
    profile(i+1,3) = mat(y,x) + br_height(2);  
        elseif ismember(mat2(y,x),[73])   %3rd bridge
    profile(i+1,3) = mat(y,x) + br_height(3);  
        elseif ismember(mat2(y,x),[74])   %4th bridge
    profile(i+1,3) = mat(y,x) + br_height(4);  
        elseif ismember(mat2(y,x),[75])   %5th bridge
    profile(i+1,3) = mat(y,x) + br_height(5);  
        elseif ismember(mat2(y,x),[76])   %6th bridge
    profile(i+1,3) = mat(y,x) + br_height(6);  
        elseif ismember(mat2(y,x),[77])   %7th bridge
    profile(i+1,3) = mat(y,x) + br_height(7);  
        elseif ismember(mat2(y,x),[78])   %8th bridge
    profile(i+1,3) = mat(y,x) + br_height(8);  
        elseif ismember(mat2(y,x),[79])   %9th bridge
    profile(i+1,3) = mat(y,x) + br_height(9);  
        elseif ismember(mat2(y,x),[80])   %10th bridge
    profile(i+1,3) = mat(y,x) + br_height(10);  
        elseif ismember(mat2(y,x),[81])   %11th bridge
    profile(i+1,3) = mat(y,x) + br_height(11);  
        elseif ismember(mat2(y,x),[82])   %12th bridge
    profile(i+1,3) = mat(y,x) + br_height(12);  
        elseif ismember(mat2(y,x),[83])   %13th bridge
    profile(i+1,3) = mat(y,x) + br_height(13);  
        elseif ismember(mat2(y,x),[84])   %14th bridge
    profile(i+1,3) = mat(y,x) + br_height(14);  
        elseif ismember(mat2(y,x),[85])   %15th bridge
    profile(i+1,3) = mat(y,x) + br_height(15);   
    else
    profile(i+1,3) = mat(y,x);
    end
    profile(i + 1, 4) = dist([x y], [a(1) a(2)], res);
end

profile(end,3)=mat(y,x);               %To recover the true height of the receiving point

function d = dist(a, b, res)
% DIST determines the distance between two points
% Dist (a, b, res)
% INPUT:
% a - 1st point
% b - 2nd point
% res - resolution
%
% OUTPUT:
% d - distance between two points

d = sqrt((a - b) * (a - b)'*res(1)*res(2));
