reliable method for getting all intersection points

[trevoravant]

Hello,

First of all, thank you so much for this document. This is the best Metapost reference I've found, and I've been using it regularly!

I'm not sure if this is the right place to ask this question, but I have been having some trouble determining a reliable method to find all of the intersection points between two paths. I have implemented both methods in Section 9.4.2, but they seem to fail sometimes. Below I have code with a simple example demonstrating what I'm talking about:

beginfig(1);

% first method in Section 9.4.2 of Drawing with Metapost
vardef get_all_intersections_one(suffix A,B) =
	save P,n,times;
	pair P[], times;
	numeric n;
	n=0;
	for i = 1 upto length(A):
		for j = 1 upto length(B):
			times := subpath (i-1,i) of A intersectiontimes subpath (j-1,j) of B;
			if xpart times > -1:
				P[incr n] = 1/2[point xpart times of subpath (i-1,i) of A,
				point ypart times of subpath (j-1,j) of B];
				draw P[n] withcolor red;
			fi
		endfor
	endfor
enddef;

% second method in Section 9.4.2 of Drawing with Metapost
vardef get_all_intersections_two(suffix A,B) =
	save P,n;
	pair P[];
	numeric n;
	n=0;
	path R;
	R := A; % take a copy of A
	forever:
		R := R cutbefore B; % snip where we cross B
		exitif length cuttings = 0; % stop if nothing was cut
		P[incr n] = point 0 of R; % capture the point
		R := subpath (epsilon, infinity) of R; % nudge along
		draw P[n] withcolor red;
	endfor
enddef;

% define a circular path "A"
pair d[];
d0 = (0,0);
d1 = (8,6);
path A;
A := d0..d1..d0;
draw A;

% create a bunch of horizontal paths "B", and find their intersections with "A"
path B;
for y=-3 upto 9:
	B := (-2,y)--(10,y);
	draw B withcolor blue;

	% to use method 1 or 2, uncomment one of the lines below
	get_all_intersections_one(A,B);
	%get_all_intersections_two(A,B);
endfor

endfig;

[thruston]

Thanks for opening this issue -- it is hard to define a completely general approach. It is always possible to find edge cases that fail in one way or another -- part of my reason for giving two different approaches was that one might work better in some situations than others, and you might like to work out how to adapt them to what you are actually trying to achieve. Let me look at your examples and see if I can find something more robust.

Do you have a particular use case for finding all intersections?

[thruston]

As it says in the discussion in §9.4.2, the first technique only works well if the "unit subpaths" are short -- it is only ever going to find one intersection for each unit subpath. Since you have rather long subpaths in your path A the first technique is not appropriate.

So you should try the second approach -- however there is a small bug in it. Try changing epsilon to eps and it should work a bit better.

[trevoravant]

Thanks for the quick response!

My use case is actually trying to find the centroid of blob-like objects. For example:

z0 = (0,0);
z1 = (60,40);
z2 = (40,90);
z3 = (10,70);
z4 = (30,50);
path blob;
blob = z0..z1..z2..z3..z4..cycle;

To find the centroid, I am approximating each blob as a series of horizontal rectangles. It is the same idea as using a rectangular approximation in a Riemann Sum. So my goal is to use the intersection points to define the ends of rectangles.

And thanks for the suggestions. I tried using both approaches, but unfortunately they both seem to miss some intersection points. I also just tried changing epsilon to eps but unfortunately it didn't seem to make a difference in that example

[thruston]

I've moved this to a discussion...
There is an easier way to find the centroid -- just take the average of your points..

pair c; c = 1/5 (z0 + z1 + z2 + z3 + z4);  

Admittedly this ignores any Bezier curves between the points, but it might be a good enough approximation.

[trevoravant]

I recently realized that, for my needs, the best way to determine the centroid of a closed path is to approximate the path as a polygon with many vertices. Then I can use the polygon centroid formula.

I have also written my own function to get all intersection points, which I will post here soon.

[trevoravant]

Ok, here is my own function to find all intersection points:

% get the length of an array
% from https://github.com/thruston/metapost-sort-routines/blob/master/sorts.mp
vardef array_length(suffix A) = 
	save i;
	numeric i;
	i := -1;
	forever:
		i := i+1;
		exitif not known A[i];
	endfor
	i
enddef;

% get all intersections of paths "A" and "B", and save the intersection points
% in an array of pairs "P"
vardef get_all_intersections(expr A,B)(suffix P) =
	save L,R,t,t_start,t_end,a,a_delta,n;
	path L,R;
	a_delta = .01;

	(t,whatever) = A intersectiontimes B;
	if not (t=-1):
		% save intersection point
		n = array_length(P);
		P[n] = point t of A;

		% determine endpoints of left and right segments
		a = arclength subpath (0,t) of A;
		t_end = arctime a-a_delta of A;
		t_start = arctime a+a_delta of A;
			
		% apply function to left segment
		if t_end>0:
			L = subpath (0, t_end) of A;
			get_all_intersections(L,B)(P);
		fi

		% apply function to right segment
		if t_start<length A:
			R = subpath (t_start, length A) of A;
			get_all_intersections(R,B)(P);
		fi
	fi
enddef;

% simple example showing usage
path A,B;
pair P[];
A = fullcircle scaled 20;
B = unitsquare scaled 20;
get_all_intersections(A,B)(P);
n = array_length(P);
draw A;
draw B;
for i=0 upto n-1: draw P[i] withcolor red; endfor

The function operates recursively. It finds an intersection point, splits the first path into two subpaths around the intersection, and then calls the function on each of the subpaths.

To ensure that the subpaths do not include the intersection at the beginning/end of the path, I have defined the variable t_delta. This value is defined as the t parameter corresponding to an arc length of .001 on path A. I had more success defining t_delta this way rather than as a constant value such as epsilon. But someone may want to modify it if isn't working for their particular case.

example from my original post:

more complex example:

[thruston]

Thanks for this -- I will review it and try to work your ideas into the document.

[trevoravant]

Ok great, thank you. I also just edited the function I posted to improve the bisection procedure.