Find the intersection of a square pixel with a closed 2D curve

Dear MatLab Experts,

I would appreciate your help with the following problem.

I have a set of closed concave 2D curves of arbitrary shape. The coordinates of the polygonal curve are known.

I cover the areas inclosed by the curves with square pixels whose side is a parameter. Therefore the pixels may be smaller or larger.

The rim of the area (curve) must be covered by the square pixels as well. Please, see the attached pictures.

Question:

How can I find the pixels that intesect the curve as distinguished from the pixels totally contained inside the curve?

Thank you in advance for any suggestion and help.

Kind regards,

Maura E. M.

[darova]

DId you try just round coordinates of a curve?

a = 1.5;            % side of a pixel
x = round(x/a)*a;
y = round(y/a)*a;

[Fabio Freschi]

Have you got a particular data structure for the pixel?

My approache would be:

• check for each point of the grid if it is inside the curve (see function inpolygon)
• if a pixel has all four corners inside the curve, the pixel is inside (the implementation of this check depends on the data structure)

If you can share the data, I can try to be more specific

[Maura E. Monville]

Thank you.

I have attached the polygons corresponding to t he pictures I have attached earlier.

The pixels are squares whose side is an input parameter. So they can be bigger or smaller...

The challenge is to find the pixels (squares) that intersect the curve. Alternatively, all the pixels (squares) completely inside the closed curve.

The reason is that the pixels crossing the rim of the closed area must have twice the resolution of the pixelse inside the closed curve.

Regards and many thanks,

maura

Answer 1

[Maura E. Monville]

Dear All,

I have attached

• the main script "Pixelize_Collimator_Aperture_v4.m"
• the funcions called by the script "CollimatorGrid_v3.m", "curveintersect.m"
• the coordinates of the closed curve that causes the problem "26-Jul-2018_A170258" in text format
• the picture showing the problem on the attached curve "A170258.jpg"
• the function I downloaded from FEX that find curves intersections "curveintersect.m"

To run the code, please do the following:

1. run "Pixelize_Collimator_Aperture_v4.m"
2. upon promp select the folder containing the curve "26-Jul-2018_A170258.txt"
3. upon prompt click on the curve file name "26-Jul-2018_A170257.txt"
4. upon prompt enter 1 for the number of nodes (of the Grid) per unit square area. This data can be any positive number (not necessaily integer). If there are too many points per unit square area then it becomes hard to see what is going on.
5. the program will stop on line 155 of script "Pixelize_Collimator_Aprture_v4.m"

Given a closed 2D arbitrary curve, the program builds a grid over the curve bounding box. The goal is to approximate the area inclosed by the curve with square pixels whose side is sqrt(#nodes per square area). The goal is to keep the pixels that

1. are totally inside the closed area
2. intersect the curve on more than one point
3. double the resolution only of the pixels that intersect the curve in more than one
4. repeat steps 1. and 2. on the new smaller pixels in order to better approximate the area

The logic of my script is wrong as it would keep a pixel (yellow in the attached picture) that is totally outside the closed curve but only has the upper left corner on the curve. This corner coordinates are printed on the picture and coincide with one of the 4 points making up the pixel as printed by the script.

I wrongly toughts that MatLab function "inpolygon" would separate the points inside a curve from the points on the curve. Instead this is returns as "inside" a point which is on the curve.

I use the contributed function "curveintersect" as a double-check. Luckily, as it helped me detect my error.

********* Question *********

In order to find the pixels that intersect the curve on more than one pont may I totally rely on function "curveintersect"?

Actually, I still have to use fynction "inpolygon" to keep the pixels totally inside the closed area. Whereas I will generate 4 smaller pixels from each original pixel that intesect the curve.

Thank you in advance for any help and/or suggestion.

Kind regards,

maura

[darova]

My english is not good enough to do you wrote

The question is still not clear for me. You trying to find intersection of a grid and a curve?

Why don't just loop over all lines and use intersect?

[Matt J]

Maura's response moved here:

Hi,

Don't worry. Your English is good enough.

My goal is to cover a closed area, delimited by an arbitrary closed curve, with square pixels making sure also the points on area rim (curve) lie inside a pixel.

My logic is:

1. cover with square pixel the close curve bounding box. This is the grid
2. Get rid of all pixels which lie outside the area
3. Get rid of all pixels that only have one corner on the curve but lie outside the area
4. Keep the pixel that are completely inside the area
5. Double the resolution of pixels that have more than one intersection with the curve. That is pixels that lie partially inside and partislly outside the area. Each of these pixels will give rise to 4 pixels whose side is half as much the side of the orifinal pixel. For each of the new smalller pixels repeat steps 2, 3,4

The goal of doubling the resolution of borderline pixels is to refine the approximation of the closed area shape.

I use function inpolygon to check if the pixel is inside or outside the area

I use function curveintersect to check if the pixel intersects the curve or not.

Thank you

Regards,

Maura

[darova]

I have an idea based on triangulation:

Create grid - number each pixel (pixels have number of nodes)

Loop over all pixels. If inpolygon() returns 0 or 1 nodes inside - delete pixel (pixel is outside), 4 nodes - pixel is inside, 2-3 - create new 4 pixels and delete an old one (double resolution).

Repeat procedure if needed

[darova]

I wanted you to accept my answer =(

Answer 2

[darova]

Doesn't consider below case (use intersect() if it matters to you):

Only when nodes inside our outside (base on inpolygon())

I succeded

[Matt J]

Agreed. The general approach should be

1. Use inpolygon to find rectangles whose vertices are entirely inside the polygon. These can be removed from further investigation
2. Loop over the remaining rectangles (a greatly reduced number, presumably) and use polyshape.intersect() to determine whether there is an intersection or not.

Answer 3

[KSSV]

Use this file exchange:

https://in.mathworks.com/matlabcentral/fileexchange/22441-curve-intersections?focused=5165138&tab=function

[Maura E. Monville]

I tried this function. It crashes immediately.

Basically I have a closed polygonal made up of roughly 150 points and a closed square made up of 5 points. Obviously for either curve the star and the last point coincide as either curve is closed.

It seems that this function takes difference between pairs of points from each curve. It crashes as the number of points making up each curve is not the same.

[darova]

Can you please attach the script?

[KSSV]

This function is apt....it should give you what you want. Attach your code...let us see how you have used the function.

[Maura E. Monville]

Again, i tried InterX. It crashes apparently because I am trying to find the intersections between two polygonals made up of a different number of points. Please, take a look at the attached pictures.

Does InterX expect the two input curve to vave to SAME number of points?

In my case I have a suare pixel (4 points, 5 if I close it) and apolygonal with an arbitrary number of points >> 5. This seems to be the problem.

Any suggestion?

Thank you.

Best regards,

Maura E.

Answer 4

[Maura E. Monville]

There are all the triky cases when a pixel has only one corner on the curve. The other 3 nodes may be inside (--> keep the pixel) or outside (--> discard the pixel).

Furthermore, function "inpolygon does not discriminates between in and on nodes.

If a pixel has only 1 node on the curve then it can have 3 other nodes inside or 2 other nodes outside.

My code i along the line you suggest. I am trying to develop a function that test that checks and returns the relationship between the pixel and the curve so taht it can be called for the pixel generated on the boundingnox and also for the pixels created by doubling the resolution of the borderline ones.

I will post my code once I am done with debugging.

Someone might suggest some improvements.

Thank you

Answer 5

[DAdler]

Hope this will get you started in ~20 lines. It is along the lines of what you described. It can be slow if you use too many pixels, but it's a start. The code uses InterX.

Have you looked at the command polyarea for the area inside the curve?

Remarks:

• I tried this code with your data and it looks fine (to me).
• No pixel refinement is done here. Change N instead!
• You require more than 1 intersection for the boundary curve. This may be problematic if the curve is tangential to the pixel boundaries, thus leaving the boundary open. If you still want this, instead of ~isempty in definining the boundary, use size("same commands InterX etc",2)>1.

clc,clf, clear all, close all
% CURVE, CIRCLE OF RADIUS 3
t = linspace(0,2*pi);
Curve = [3*cos(t);3*sin(t)];

N = 21; % (N-1)X(N-1) PIXELS TO COVER THE DOMAIN

% DEFINE THE BOUNDING BOX AND PIXEL GRID COORDINATES 
Xmax = 4; Xmin = -4; 
[X,Y] = meshgrid(linspace(Xmin,Xmax,N));
VV = [X(:),Y(:)];

% PIXEL CONNECTIVITY LIST
[x,y] = meshgrid(1:N);
V = [x(:),y(:)];
j = setdiff(1:N^2-N,N:N:N^2-N)';
C = mat2cell([0,1,N+1,N,0]+j,ones(1,(N-1)^2),5);

% BOUNDARY PIXELS
boundary = find(cell2mat(cellfun(@(c)~isempty(InterX([VV(c,1),VV(c,2)]',Curve)),C,'UniformOutput',0)));

% INTERIOR PIXELS
interior = find(cell2mat(cellfun(@(c) all(inpolygon(VV(c,1),VV(c,2),Curve(1,:),Curve(2,:))),C,'UniformOutput',0)));
interior = setdiff(interior,boundary);

figure, hold on

% DRAW PIXELS
for i=boundary'
    fill(VV(C{i},1),VV(C{i},2),'k','edgecolor','c')  
end

for i=interior'
    fill(VV(C{i},1),VV(C{i},2),'m','edgecolor','c')  
end

% DRAW CURVE
plot(Curve(1,:),Curve(2,:),'r')
axis equal

[Maura E. Monville]

It is an excellent piece of code.

I am having a hard time at understanding how you isolate the borderline pixels. I am not that good at manipulating cells.

A problem may be the finer resolution taht is only required for the borderline pixels in order to beter approximate the curve.

There is no need to refine the resolution of pixels inside the curve. It would just slow down the performence without adding any benefit.

Thank you

[DAdler]

InterX is used to check whether each of the pixels (essentially a square) intersects (or is tangential) to the curve. Similarly inpolygon is used to determine whether all vertices of the pixel lie within the curve.

If you don't like cellfun (which maybe slower than the for loop) this

C = mat2cell([0,1,N+1,N,0]+j,ones(1,(N-1)^2),5);
find(cell2mat(cellfun(@(c)~isempty(InterX([VV(c,1),VV(c,2)]',Curve)),C,'UniformOutput',0)));
interior = find(cell2mat(cellfun(@(c) all(inpolygon(VV(c,1),VV(c,2),Curve(1,:),Curve(2,:))),C,'UniformOutput',0)));
interior = setdiff(interior,boundary);

should be the same as

C = [0,1,N+1,N,0]+j;
boundary = zeros((N-1)^2,1);
interior = zeros((N-1)^2,1);
for i=1:(N-1)^2
    if ~isempty(InterX([VV(C(i,:),1),VV(C(i,:),2)]',Curve))
       boundary(i) = i; 
    elseif all(inpolygon(VV(C(i,:),1),VV(C(i,1:4),2),Curve(1,1:4),Curve(2,:)))
       interior(i) = i;
    end
end
boundary(boundary==0)=[];
interior(interior==0)=[];

[Maura E. Monville]

I tried again to use InterX. Again it crashed:

I managed to attach a picture of the printed error and a picture of how I call InterX.

I simply call InterX passing vector sqx containing the x_coordinates of the square pixel, vector sqy containing the y_coordinates of the square pixel,

vector XC containing the x_coordinates of the closed polygonal, vector YC containing the y_coordinates of the closed polygonal.

I ask again the question I asked a week ago:

The square pixel has 4 point in all (actually 5 as it is a closed polygonal). The curve has an arbitrary number of points, for sure >> 5 and is closed.

• Do the surves input to InterX have to have the same number of points?
• Both the square pixel and the curve are closed. Which implies that the first and lthe ast point coincide. Is this s problem for InterX?

Please, help me get over this impass.

Thank you.

Best regards,

maura

[DAdler]

Most likely you are calling InterX incorrectly, i.e you calling it with Xc and Yc being column vectors. They should be row vectors to form a 2xN array,

No, the curves do not need to be of the same length!

It does not matter if the 2 points coincide to close the curve, I guess. The code tests if segments intersect.

Answer 6

[Maura E. Monville]

I agree that is very likely to be the way I called InterX. On the other hand, the providd documentation does not explicitly explains how InterX expects the input data. Nor it was clea to me from the only provided example.

At the moment I found function "curveintersect" to work sufficiently well for my need. I am saying "sufficiently" as I found out the intersection may be between the interpolater polygonal and the square pixel. I expected the intersection point to belomng to the closed curve. However, "curveintersect" returns intersection coordinates that do not belomg to the curve but are very closed to it.

Later on I will give InterX a try again. Not immediately though as this is only one of my many tasks.

I will keep you posted about the outcome.

Thank you.

Regards,

maura

Answer 7

[DAdler]

Maybe this would have resolved your issue if you provide XC and YC as column vectors?

P = InterX([sqx;sqy],[XC,YC]')

Answer 8

[Maura E. Monville]

I tried the following simple example using functions "curveintersect" and "IntersX".

If I plot the square and the rectangle I can easily tell that they share a segment.

I get two different answers. Please, see attached screenshots as I cannot paste anything inside this window.

Regards,

Maura