function img = drawLine(img, p1, p2, line_color, line_width)
% The points belonging to the line are interpolated. Their intensity
% is inverse proportional to the distance to middle of the line.
% Example: drawLine(img, [10 20], [19 90], 1, 3);

[image_height image_width colorPlanes]=size(img);
if(nargin < 3)
    disp('Insufficient number of input params');
else
    if(nargin == 3)
        if(colorPlanes > 1)
            line_color = [1 1 1];
        else
            line_color = 1;
        end
        line_width = 1;
    end
end

startpoint = p1;
endpoint   = p2;
if(p1(2) > p2(2))
    startpoint = p2;
    endpoint   = p1;
end

if(p1(2) == p2(2))
    for x = max(1, floor(p1(2))):min(image_width, floor(p1(2))+line_width-1)
        if(startpoint(1) < endpoint(1))
            for y = startpoint(1):min(endpoint(1), image_height)
                if(colorPlanes > 1)
                    img(y,x,:) = line_color;
                else
                   img(y,x) = line_color;
                end
            end
        else
            for y = endpoint(1):min(startpoint(1), image_height)
                if(colorPlanes > 1)
                    img(y,x,:) = line_color;
                else
                   img(y,x) = line_color;
                end
            end        
        end
    end
else
    m = (startpoint(1)-endpoint(1))/(startpoint(2)-endpoint(2));
    n = (endpoint(1)*startpoint(2)-endpoint(2)*startpoint(1))/(startpoint(2)-endpoint(2));
    normalizingFactor = sqrt(1+m*m);
    for x = max(1, startpoint(2)):min(endpoint(2), image_width)
        y = m*x+n;
        if((y > 0) && (y < image_height))
            y1 = floor(y-line_width);
            if(y1 < 1)
                y1 = 1;
            end
            y2 = floor(y)+line_width+1;
            if(y2 > image_height)
                y2 = image_height;
            end
            constantSum = m*x+n;
            for ylocal=y1:y2             
                dist = abs((constantSum-ylocal)/normalizingFactor);
                if(dist < line_width/2.0)
                    if(colorPlanes > 1)
                        img(ylocal,x, :) = line_color*(line_width-2*dist)/line_width;
                    else
                        img(ylocal,x) = line_color*(line_width-2*dist)/line_width;
                    end
                end
            end
        end
    end

end