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

10 views (last 30 days)
Maura Monville on 17 Jun 2019
Answered: Maura Monville on 4 Jul 2019
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 on 18 Jun 2019
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 on 18 Jun 2019
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 Monville on 18 Jun 2019
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

Maura Monville on 22 Jun 2019
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

Show 1 older comment
Matt J on 23 Jun 2019
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 on 23 Jun 2019
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 on 1 Jul 2019
I wanted you to accept my answer =(

darova on 24 Jun 2019
Doesn't consider below case (use intersect() if it matters to you): Only when nodes inside our outside (base on inpolygon())
I succeded #### 1 Comment

Matt J on 24 Jun 2019
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.

KSSV on 18 Jun 2019

Show 1 older comment
darova on 21 Jun 2019
Can you please attach the script?
KSSV on 21 Jun 2019

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

Maura Monville on 26 Jun 2019
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.

Maura Monville on 23 Jun 2019
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

Edited: DAdler on 23 Jun 2019
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

Show 1 older comment
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 Monville on 26 Jun 2019
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.
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
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.

Maura Monville on 30 Jun 2019
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

Maybe this would have resolved your issue if you provide XC and YC as column vectors?
P = InterX([sqx;sqy],[XC,YC]')

Maura Monville on 4 Jul 2019
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