Line of sight between points in a logical matrix

Question

Neuropragmatist 2018-4-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/393352-line-of-sight-between-points-in-a-logical-matrix

编辑： John BG 2018-4-11

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John BG 2018-4-8

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/393352-line-of-sight-between-points-in-a-logical-matrix#answer_313984

在 MATLAB Online 中打开

Hi Roddy

this is John BG <mailto:jgb2012@sky.com jgb2012@sky.com>

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

<mailto:jgb2012@sky.com jgb2012@sky.com>

4 个评论
显示 2更早的评论隐藏 2更早的评论

Neuropragmatist 2018-4-8

Hi John,

I have added some new code as a comment to my own answer, mainly I changed the interpolation to nearest and added a tolerance limit for detecting wall crossings (although I don't think the tolerance is necessary if the interpolation is 'nearest'). On testing this seems to work for any point I try now.

I didn't know there was a rule about accepting your own answer - very often my questions go unanswered, but if I find a solution I like to come back and add the info for further people... I have unaccepted my answer if you would like to accept it instead.

Rod.

Walter Roberson 2018-4-8

You cannot accept your own answer to someone else's question. You can accept your own answer to your own question.

请先登录，再进行评论。

Answer 2

Neuropragmatist 2018-4-8

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/393352-line-of-sight-between-points-in-a-logical-matrix#answer_314024

在 MATLAB Online 中打开

I think I have found a relatively simple solution that should also be quite fast. I use the inbuilt function improfile which was mentioned by image analyst somewhere (I have lost which post it was!) to interpolate out the values between my starting position and any pixel in the map.

Next I say that if this interpolated vector contains any zeros (the value assigned to walls) the pixel is given a value representing 'not visible'.

This approach results in some missing pixels, I assume this is because the interpolation passes between wall pixels perfectly and the zeros are missed out. To rectify this I use bwareaopen to fill in these 'holes'.

Here is the code and a plot of the result:

% create environment
S1=30;
S2=30;
emap = true(S1,S2);
emap = padarray(emap,[1 1],false,'both');
emap(1:15,16) = false;
% define current position
nx_spoint=9;
ny_spoint=6;
spoint=[nx_spoint,ny_spoint];
% show matrix
figure(1);
subplot(1,3,1)
imshow(emap);
daspect([1 1 1])
axis xy
title('environment')
% create 'viewmap'
vmap = zeros(size(emap));
indx = find(emap);
for ii = 1:length(indx)
    [y,x] = ind2sub(size(emap),indx(ii));
      xi = [spoint(2),x];
      yi = [spoint(1),y];
      C = improfile(emap,xi,yi,'bilinear');
      if sum(~C)>0
          vmap(indx(ii)) = 5;
      else
          vmap(indx(ii)) = 10;        
      end
  end
% show resulting viewmap
subplot(1,3,2)
imagesc(vmap);
daspect([1 1 1])
axis xy
title('viewmap')
% close single pixel 'holes' in viewmap caused by lines passing between pixels without hitting them
vmap2 = vmap;
vmap2(vmap<10) = 0;
vmap2 = logical(vmap2);
vmap2 = bwareaopen(vmap2,2);
% show resulting viewmap
subplot(1,3,3)
imshow(vmap2);
daspect([1 1 1])
axis xy
title('final viewmap')

2 个评论
显示无隐藏无

Neuropragmatist 2018-4-8

在 MATLAB Online 中打开

As per John's testing, the correct code should actually be:

% create environment
S1=30;
S2=30;
emap = true(S1,S2);
emap = padarray(emap,[1 1],false,'both');
emap(1:15,16) = false;
% define current position
nx_spoint=25;
ny_spoint=5;
spoint=[nx_spoint,ny_spoint];
% show matrix
figure(1);
subplot(1,3,1)
imshow(emap);
hold on
plot(ny_spoint,nx_spoint,'kx')
daspect([1 1 1])
axis xy
title('environment')
% create 'viewmap'
vmap = zeros(size(emap));
indx = find(emap);
tolerance = 0.3;
for ii = 1:length(indx)
    [y,x] = ind2sub(size(emap),indx(ii));
      xi = [spoint(2),x];
      yi = [spoint(1),y];
      C = improfile(emap,xi,yi,'nearest');
      if sum(C<tolerance)>0
        vmap(indx(ii)) = 5;
      else
        vmap(indx(ii)) = 10;        
      end
  end
% show resulting viewmap
subplot(1,3,2)
imagesc(vmap);
hold on
plot(ny_spoint,nx_spoint,'kx')
daspect([1 1 1])
axis xy
title('viewmap')
% close single pixel 'holes' in viewmap caused by lines passing between pixels without hitting them
vmap2 = vmap;
vmap2(vmap<10) = 0;
vmap2 = logical(vmap2);
vmap2 = bwareaopen(vmap2,2);
% show resulting viewmap
subplot(1,3,3)
imshow(vmap2);
hold on
plot(ny_spoint,nx_spoint,'kx')
daspect([1 1 1])
axis xy
title('final viewmap')

John BG 2018-4-11

编辑：John BG 2018-4-11

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

jgb2012@sky.com

请先登录，再进行评论。

Line of sight between points in a logical matrix

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

2 个评论
显示无隐藏无

另请参阅

类别

标签

Community Treasure Hunt

Line of sight between points in a logical matrix

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

2 个评论 显示 无隐藏 无

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

2 个评论
显示无隐藏无