forums.silverfrost.com

[arctica]

Hello,

Is there a way to calculate the intersections of two different curves given the data of two arrays, (x1, y1) and (x2, y2)?

Thanks
Lester

[JohnCampbell]

Hi Lester,

The problem does not look too difficult.
A practical solution would depend on how the curves are defined.

I am understanding you have two arrays > ie tables of points,
like:
real :: array1(2,n1) and array2(2,n2).

This could imply the curves are approximated by straight lines between the list of points, ie a straight between array1(:,i) and array1(:.i+1).

A solution could be to find the two straight lines:-
between array1(:,i) and array1(:.i+1) and
between array2(:,j) and array2(:.j+1),
that do actually cross within their space.

This is equivalent to finding i,j so that the straight lines do cross.
Basically,
do i = 1,n1-1
do j = 1,n2-1 ;
is this i,j a solution to the intersection of 2 straight lines (high school math) exit
end do
end do

Alternatively, should the lines be defined by functions f1(x) and f2(x),
then the solution could simply be to find the zero of function f1(x) - f2(x)

You could start there and get an approach that generates a solution, then look for a better way ?

My apologies if I have misunderstood your question.

John

[arctica]

Hi John

Thanks for the ideas to test. Yes the curves are defined as arrays of x,y values and not functions.

I would think it logical to look for a small region to test an approximation, where a tolerance on (x1-x2, y1-y2) is small enough. I do have a Scilab code that does this, so may try to convert that as a start.

Will keep you posted if I crack this one!

Lester

[JohnCampbell]

I thought of a different approach to this problem.

My clearwin+ graphics is based on developing a "hidden line" removal approach when drawing 3D oblects made up of lines and flat surfaces. This is monitored with two 2D arrays of pixel colours and pixel depths.

A different approach to the curve intersection problem could be to draw the two curves on the screen pixel arrays and note which pixels are allocated to the first curve line elements and then for the second curve.
For the second curve, it can be noted when the 2 curves use the same pixel, but if not then scan the pixel array for the closest pixels of the two curves.

A different approach for a more general solution.

( I have always been amazed how fast the pixel based hidden line removal approach is. I developed this in about 1995 and increasing processor and memory speeds have just improved the performance ! Even after configuring for 4K screens, although I previously used larger virtual screens when the physical screen resolution was too small for high resolution plots )

[arctica]

Just for fun this is what chatGPT came up with, posted as normal text as it kept getting truncated via [code]!

program FindCurveIntersections
implicit none
integer, parameter :: num_points = 1000 ! Number of points defining each curve
real, dimension(num_points) :: x1, y1, x2, y2, intersections_x, intersections_y
integer :: i, num_intersections

! Generate random data for curve 1
call generate_random_data(num_points, x1, y1)

! Generate random data for curve 2
call generate_random_data(num_points, x2, y2)

! Find intersections between the curves
num_intersections = find_intersection(x1, y1, x2, y2, num_points, intersections_x, intersections_y)

! Output the intersections found
print *, "Number of intersections:", num_intersections
if (num_intersections > 0) then
print *, "Intersections (x, y):"
do i = 1, num_intersections
print *, intersections_x(i), intersections_y(i)
end do
else
print *, "No intersections found."
end if
end program FindCurveIntersections

subroutine generate_random_data(num_points, x, y)
implicit none
integer, intent(in) :: num_points
real, dimension(num_points), intent(out) :: x, y
integer :: i

do i = 1, num_points
call random_number(x(i))
call random_number(y(i))
end do
end subroutine

real function find_intersection(x1, y1, x2, y2, num_points, intersections_x, intersections_y)
implicit none
real, dimension(, intent(in) :: x1, y1, x2, y2
real, dimension(, intent(out) :: intersections_x, intersections_y
integer, intent(in) :: num_points
integer :: i, j, count
real :: tolerance
tolerance = 1.0e-6 ! Define a tolerance for comparing points

count = 0

do i = 1, num_points - 1
do j = 1, num_points - 1
! Check if line segments (i, i+1) from curve 1 and (j, j+1) from curve 2 intersect
if (do_line_segments_intersect(x1(i), y1(i), x1(i+1), y1(i+1), x2(j), y2(j), x2(j+1), y2(j+1), tolerance)) then
count = count + 1
call interpolate_intersection(x1(i), y1(i), x1(i+1), y1(i+1), x2(j), y2(j), x2(j+1), y2(j+1), &
intersections_x(count), intersections_y(count))
end if
end do
end do

find_intersection = count ! Return the number of intersections found
end function

logical function do_line_segments_intersect(x1, y1, x2, y2, x3, y3, x4, y4, tolerance)
implicit none
real, intent(in) :: x1, y1, x2, y2, x3, y3, x4, y4, tolerance
real :: dx1, dy1, dx2, dy2, delta
real :: s, t

dx1 = x2 - x1
dy1 = y2 - y1
dx2 = x4 - x3
dy2 = y4 - y3

delta = dx1 * dy2 - dy1 * dx2

if (abs(delta) < tolerance) then
do_line_segments_intersect = .false. ! Lines are parallel or coincident
else
s = (dx1 * (y1 - y3) - dy1 * (x1 - x3)) / delta
t = (dx2 * (y1 - y3) - dy2 * (x1 - x3)) / delta

if (s >= 0.0 .and. s <= 1.0 .and. t >= 0.0 .and. t <= 1.0) then
do_line_segments_intersect = .true. ! Segments intersect
else
do_line_segments_intersect = .false. ! Segments do not intersect
end if
end if
end function

subroutine interpolate_intersection(x1, y1, x2, y2, x3, y3, x4, y4, intersection_x, intersection_y)
implicit none
real, intent(in) :: x1, y1, x2, y2, x3, y3, x4, y4
real, intent(out) :: intersection_x, intersection_y
real :: dx1, dy1, dx2, dy2, s

dx1 = x2 - x1
dy1 = y2 - y1
dx2 = x4 - x3
dy2 = y4 - y3

s = ((x3 - x1) * dy1 * dx2 + y1 * dx1 * dx2 - y3 * dx2 * dx1) / (dy2 * dx1 - dy1 * dx2)

intersection_x = x3 + s * dx2
intersection_y = y3 + s * dy2
end subroutine

Lester

[Kenneth_Smith]

My first attempt:

[mecej4]

Praiseworthy, excellent graphical illustration!

[arctica]

Thanks for the suggestions Kenneth, the test code worked great.

Added the routines to my existing code and it works fine to compute the intersections, but having issues getting it plotting symbols at the intersections.

This all works fine:

[Kenneth_Smith]

Lester,

I suspect your problem is related to the %pl options [n_graphs=*], [x_array] and [independent].

In my example I used all of these. This means that for each curve in the plot an array of values corresponding to the x coordinates and the y coordinates of each curve must be provided in the call to %pl. Thus there are six coordinate arrays passed as arguments to the winio call with %pl (i.e. x1, y1, x2, y2, x3, y3). In this case the number of points in each curve is passed as an array. These can be seen below for [n_graphs=3,x_array,independent].

[arctica]

Hi Kenneth

Managed to get the plot working, just followed you original syntax: