2D matrix viewshed

Hi everyone,

I have a binary matrix which represents the floorplan of a building, given a position in this floorplan I would like to know all of the positions visible by direct line of sight. At the moment I am implementing the viewshed function by passing by binary matrix as an elevation map, i.e.:

floorplan = padarray(false(150,150),[1 1],true,'both'); % 1s = walls, 0s = space
floorplan(75,1:75) = true;

Z = floorplan.*100;
R = georefcells([0 1],[0 1],size(floorplan));

ind = 3000;

[rx,cx] = ind2sub(size(floorplan),ind);        
xnow = (rx-1)./(size(floorplan,1)-1);
ynow = (cx-1)./(size(floorplan,2)-1);
los = viewshed(Z,R,xnow,ynow,50,0,"AGL","AGL",1e2);

figure
subplot(1,2,1)
imshow(floorplan); hold on;
plot(cx,rx,'r+')
axis on xy
title('floorplan')

subplot(1,2,2)
imshow(los); hold on;
plot(cx,rx,'r+')
axis on xy
title('line of sight')

However, the viewshed function is a bit overpowered for this purpose because it is designed to work with 3D elevation data. To increase speed I was wondering if there is a faster way to achieve what I want? The documentation for viewshed doesn't specify how it works, so I'm not sure if I can simplify it for a 2D situation?

Thanks for any help.

Answer 1

[Maneet Kaur Bagga]

Hi,

As per my understanding you want to compute the line of sight without using the "viewshed" function.

One of the possible ways to do this is to use "Bresenham's line algorithm". Using this you can determine which pixels should be drawn to form a straight line between two points on a grid. It can be adapted to check for visibility by stopping the line drawing process when a wall (a 1 in the matrix) is encountered.

Please refer to the code snippet below using "Bresenham's line algorithm" to check the visibility for a single ray. The function checks if there is a clear line of sight between the start and end positions within the floorplan.

function isVisible = checkVisibility(floorplan, startRow, startCol, endRow, endCol)
    % Initialize visibility as false
    isVisible = false;

    % Calculate deltas
    deltaRow = abs(endRow - startRow);
    deltaCol = abs(endCol - startCol);

    % Calculate error
    error = 0;
    
    % Determine the direction of the steps
    if startCol < endCol
        stepCol = 1;
    else
        stepCol = -1;
    end
    if startRow < endRow
        stepRow = 1;
    else
        stepRow = -1;
    end

    % Initialize current position
    row = startRow;
    col = startCol;

    % Choose the dominant direction
    if deltaCol > deltaRow
        % Horizontal dominant
        for c = startCol:stepCol:endCol
            % Check for a wall
            if floorplan(row, col)
                return;
            end
            error = error + deltaRow;
            if 2 * error >= deltaCol
                row = row + stepRow;
                error = error - deltaCol;
            end
            col = col + stepCol;
        end
    else
        % Vertical dominant
        for r = startRow:stepRow:endRow
            % Check for a wall
            if floorplan(row, col)
                return;
            end
            error = error + deltaCol;
            if 2 * error >= deltaRow
                col = col + stepCol;
                error = error - deltaRow;
            end
            row = row + stepRow;
        end
    end

    % If we reach here, the path is clear
    isVisible = true;
end

Please refer to the File Exchange link below for better understanding of the implementation: https://in.mathworks.com/matlabcentral/fileexchange/25544-line-drawing-by-bresenham-algorithm

Hope this helps!

[Neurolab]

Hi Maneet,

This is a really nice piece of code, thank you. It is much faster than viewshed, however, there are some errors when a piece of wall is 'hidden' behind a neighbor, presumably because the algorithm tries to take the same path to both and encounters a wall. For example, see the black pixels along the top of the right graph here:

The wall positions at the top are marked as not visible (black) from the position of the red cross, because they are obscured by the closer wall positions, they need to be marked as visible for my purposes. I'm not sure the best way to correct for this?

[Maneet Kaur Bagga]

Hi

To solve this, one approach is to slightly extend the line beyond the point where it hits the wall, just enough to include the edge of the wall in the visibility check. This can be done by adjusting the stopping condition of the algorithm to continue the line to the next point after hitting a wall, ensuring that the edge of the wall is marked as visible.

Here's an updated version of the "checkVisibility" function that includes this adjustment:

function isVisible = checkVisibility(floorplan, startRow, startCol, endRow, endCol)
    % Initialize visibility as false
    isVisible = false;

    % Calculate deltas
    deltaRow = abs(endRow - startRow);
    deltaCol = abs(endCol - startCol);

    % Calculate error
    error = 0;
    
    % Determine the direction of the steps
    stepCol = sign(endCol - startCol);
    stepRow = sign(endRow - startRow);

    % Initialize current position
    row = startRow;
    col = startCol;

    % Choose the dominant direction
    if deltaCol > deltaRow
        % Horizontal dominant
        while col ~= endCol
            if floorplan(row, col)
                % If we hit a wall, check the next position before stopping
                if row ~= endRow
                    isVisible = ~floorplan(row + stepRow, col);
                else
                    isVisible = true;
                end
                return;
            end
            error = error + deltaRow;
            if 2 * error >= deltaCol
                row = row + stepRow;
                error = error - deltaCol;
            end
            col = col + stepCol;
        end
    else
        % Vertical dominant
        while row ~= endRow
            if floorplan(row, col)
                % If we hit a wall, check the next position before stopping
                if col ~= endCol
                    isVisible = ~floorplan(row, col + stepCol);
                else
                    isVisible = true;
                end
                return;
            end
            error = error + deltaCol;
            if 2 * error >= deltaRow
                col = col + stepCol;
                error = error - deltaRow;
            end
            row = row + stepRow;
        end
    end

    % If we reach here, the path is clear
    isVisible = true;
end

Hope this helps!