# Line of sight between points in a logical matrix

10 views (last 30 days)

Show older comments

Hi all,

I am trying to work out the 'line of sight' between points in a logical matrix, i.e. if zeros represent floor space and ones represent walls, I would like to know all the points in the matrix than could be observed from a specific row,column location directly, without the line of sight crossing a wall.

I have tried using bwdist and bwdistgeodesic to solve this problem, my first thought was to take the difference between the maps produced using these which would give me a map of all pixels where the euclidean distance and distance traversing around walls are the same (directly visible). Alas this doesn't seem to be the case (see code below).

There must be an easy way to do this but I just can't seem to work it out...

Thanks for any help,

Rod.

emap = true(30,30);

emap = padarray(emap,[1 1],false,'both');

emap(1:15,16) = false;

spoint = [9,6];

figure

subplot(2,2,1)

imshow(emap)

daspect([1 1 1])

axis xy

title('Shape of interest')

subplot(2,2,2)

D1 = bwdistgeodesic(emap,spoint(2),spoint(1),'quasi-euclidean');

imagesc(D1)

daspect([1 1 1])

axis xy

c = caxis;

title('Geodesic distance')

subplot(2,2,3)

bw = false(size(emap));

bw(spoint(1),spoint(2)) = true;

D2 = bwdist(bw,'quasi-euclidean');

D2(~emap) = 0;

imagesc(D2)

daspect([1 1 1])

axis xy

caxis(c)

title('Euclidean distance')

subplot(2,2,4)

imagesc(D1-D2)

daspect([1 1 1])

axis xy

title('Difference')

##### 0 Comments

### Answers (2)

John BG
on 8 Apr 2018

Hi Roddy

1.- the map

clear all;close all;clc

S1=30;S2=30;

emap = true(S1,S2);

emap = padarray(emap,[1 1],false,'both');

emap(1:15,16) = false;

2.-

obstacle definition

K=[16*ones(15,1) [1:15]'];

3.-

the spotter

nx_spoint=9;

ny_spoint=6;

spoint=[nx_spoint,ny_spoint];

figure(1);imshow(emap);ax1=gca

4.- map points free of obstacles

[nx,ny]=find(emap==true);

5.- checking the start area free of obstacles has the correct points

A=zeros(size(emap)); %

hf2=figure(2);imshow(A);ax2=hf2.CurrentAxes

plot(ax2,ny,nx,'r*');

axis(ax2,'equal');

6.- spotter on the map

hold(ax1,'all');plot(ax1,spoint(1),spoint(2),'g.'); % 'LineSpec' Style-Marker-Color 'LineStyle'

hold(ax2,'all');plot(ax2,spoint(1),spoint(2),'g*'); % 'LineSpec' Style-Marker-Color 'LineStyle'

P=[nx ny]; % free space points, excluding obstacles

7.-

remove spotter point from area to sweep

L1=[];

for k=1:1:size(P,1)

if P(k,1)==nx_spoint && P(k,2)==ny_spoint

L1=[L1 k];

end

end

P(L1,:)=[];

8.-

defining outer perimeter: the fence along which the spotter is going to sweep along

left_side=[ones(1,S1)' [1:1:S1]']

top_side=[[2:1:S2-1]' S1*ones(1,S2-2)']

right_side=[S2*ones(1,S1)' [S1:-1:1]']

bottom_side=[[S2-1:-1:2]' ones(1,S2-2)']

Fence=[top_side;right_side;bottom_side;left_side];

9.-

check obstacle points correctly on map

plot(ax2,K(:,1),K(:,2),'c*')

for k2=1:1:size(Fence,1)

P0=Fence(k2,:);

nxP0=P0(1);nyP0=P0(2);

plot(ax2,nxP0,nyP0,'b*');

10.-

define one ray points, spotter - fence

Lx=floor(linspace(nx_spoint,nxP0,max(abs(nxP0-nx_spoint),abs(nyP0-ny_spoint))));

Ly=floor(linspace(ny_spoint,nyP0,max(abs(nxP0-nx_spoint),abs(nyP0-ny_spoint))));

L1=[Lx' Ly'];

11.-

intersect finds, if any, intersect points for 2D vectors, intersect.m attached along with this script.

[U,n1,n2]=intersectN(L1,K)

12.-

have to decide whether truncate ray or not

if isempty(U)

plot(ax2,Lx,Ly,'y.');drawnow; % ray without obstacle

% pause(.25)

else

if size(U,2)>1 % solve intersectN may return more than 1 point

U=U(1,:);

end

13.-

Following bank of ifs to make sure rays cover from all 4 quadrants

if spoint(1)<U(1)

Lx2obst_=[spoint(1):1:U(1)]';

end

if spoint(1)>U(1)

Lx2obst_=[spoint(1):-1:U(1)]';

end

if spoint(1)==U(1)

Lx2obst=spoint(1);

end

if spoint(2)<U(2)

Ly2obst_=[spoint(2):1:U(2)]';

end

if spoint(2)>U(2)

Ly2obst_=[spoint(2):-1:U(2)]';

end

if spoint(2)==U(2)

Ly2obst=spoint(2);

end

Lx2obst=floor(linspace(spoint(1),U(1),max(numel(Lx2obst_),numel(Ly2obst_))))

Ly2obst=floor(linspace(spoint(2),U(2),max(numel(Lx2obst_),numel(Ly2obst_))))

plot(ax2,Lx2obst',Ly2obst','y.') % ray up to obstacle

% pause(.25)

end

end

.

.

.

Note that the plot function shows the image reversed compared to the imshow of logical 2D maps that you start with, but the coordinates are all the same.

Roddy

f you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance for time and attention

John BG

##### 4 Comments

Walter Roberson
on 8 Apr 2018

Neuropragmatist
on 8 Apr 2018

##### 2 Comments

John BG
on 11 Apr 2018

Edited: John BG
on 11 Apr 2018

Hi Roddy

this is John BG jgb2012@sky.com

You have fixed the quadrants that looked as if spotter out of square, following is with spotter [52 25] which is ok because the observer is not supposed to be out of square

Yet,

.

.

the few pixels on the left of the ray for spotter on [25 25] are pixels that shouldn't see but apparently your answer claims to see beyond the cut ray. Is this ok to you?

.

.

With these simple answers we are supply it's understandable a +-1 pixels error away from the correct pixels, but +-2 pixels away, why is it?

2 pixels too many may be the difference between missing or hitting a pedestrian, don't you agree?

regards

John BG

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!