Welcome to our forums
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
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
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
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 )
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
My first attempt:
winapp
program intersect_between_curves
USE clrwin
integer, parameter :: n=301
integer :: i, j, n_intersects
real :: x1(n), y1(n), x2(n), y2(n), th1, a, dth1, xi, yi
real, allocatable :: xip(:), yip(:)
logical :: intersects
! generate x-y data arrays, x1,y1 defines a Fermant spiral
a = 1.0 ; th1 = 0.0 ; dth1 = 4.0*atan(1.0)/25.0
do i = 1, n, 1
x1(i) = a * sqrt(th1) * cos(th1) ; y1(i) = a * sqrt(th1) * sin(th1) ; th1 = th1 + dth1
end do
! x2,y2 is the same spiral displaced by +1 in the x and y directions
x2 = x1 + 1.0 ; y2 = y1 + 1.0
! test each possible pair of line segments from the two curves for intersection
open(unit=10,status='scratch',form='unformatted')
n_intersects = 0
do i = 1, n-1, 1
do j = 1, n-1, 1
call LIP(x1(i), y1(i), x1(i+1), y1(i+1), x2(j), y2(j), x2(j+1), y2(j+1), intersects, xi, yi)
if (intersects) then
write(10) xi, yi
n_intersects = n_intersects + 1
end if
end do
end do
! copy intersection points from scratch file to arrays xip and yip
if (allocated(xip)) deallocate(xip)
if (allocated(yip)) deallocate(yip)
allocate(xip(n_intersects),yip(n_intersects))
rewind 10
do i = 1, n_intersects, 1
read(10) xip(i), yip(i)
write(6,'('(',F6.3,',',F6.3,')')') xip(i), yip(i)
end do
close(unit=10)
! configure plot
call winop@('%pl[native,n_graphs=3,independent,frame,gridlines,smoothing=4,x_array,x_min=-6]')
call winop@('%pl[link=lines,colour=red,symbol=0,link=lines,colour=blue,symbol=0,link=none,colour=black,symbol=6]')
i = winio@('%fn[consolas]%ts%bf&',1.5d0)
i = winio@('%pl',900,900,[n,n,n_intersects],dble(x1),dble(y1),dble(x2),dble(y2),dble(xip),dble(yip))
end program intersect_between_curves
subroutine LIP(x1, y1, x2, y2, x3, y3, x4, y4, intersects, xi, yi)
! LIP 'Line Intersection Point'
! Line 1 is from point (x1,y1) to point (x2,y2). Line 2 is from point (x3,y3) to point (x4,y4)
! If lines 1 and 2 intersect, intersects is set to .TRUE. and the coordinates of the intersection
! are returned as (xi,yi), otherwise intersects is set to .FALSE.
real, intent(in) :: x1, y1, x2, y2, x3, y3, x4, y4
logical, intent(out) :: intersects
real, intent(out) :: xi, yi
real :: tol = epsilon(1.0), det, a, b
intersects = .false.
det = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4)
! If det is zero lines are parallel
if (abs(det) .lt. tol) return
! calculate intersection point
a = (x1 * y2 - y1 * x2)
b = (x3 * y4 - y3 * x4)
xi = (a * (x3 - x4) - (x1 - x2) * b) / det
yi = (a * (y3 - y4) - (y1 - y2) * b) / det
! check if intersection point lies on both line segments
intersects = ((xi .ge. min(x1, x2) .and. xi .le. max(x1, x2)) .and. (xi .ge. min(x3, x4) .and. xi .le. max(x3, x4)) .and. &
(yi .ge. min(y1, y2) .and. yi .le. max(y1, y2)) .and. (yi .ge. min(y3, y4) .and. yi .le. max(y3, y4)))
end subroutine LIP
Praiseworthy, excellent graphical illustration!
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:
! test each possible pair of line segments from the two curves for intersection
open(unit=10,status='scratch',form='unformatted')
n_intersects = 0
do i = 1, k-1, 1
do j = 1, k-1, 1
call LIP (vec(i), real_part(i), vec(i+1), real_part(i+1), vec(j),&
imaginary_part(j), vec(j+1), imaginary_part(j+1), intersects, xi, yi)
if (intersects) then
write(10) xi, yi
n_intersects = n_intersects + 1
end if
end do
end do
! copy intersection points from scratch file to arrays xip and yip
if (allocated(xip)) deallocate(xip)
if (allocated(yip)) deallocate(yip)
allocate(xip(n_intersects),yip(n_intersects))
rewind 10
do i = 1, n_intersects, 1
read(10) xip(i), yip(i)
where (yip .lt. 0.01) yip=0.0
write(6,'('(',F6.3,',',F6.3,')')') xip(i), yip(i)
end do
close(unit=10)
Assed in a tolerance aspect to narrow things down to get the non-trivial zero locations.
Tried to add the xip, yip and symbol syntax as you showed in the following and it would compile but either make a blank plot or give a run-time error.
iw = winio@('%mn[Exit]&','Exit')
iw = winio@('%mn[Export]&',export_cb)
iw = winio@('%fn[Times New Roman]%ts%bf&',1.5d0)
call winop@('%pl[native,x_array,frame,gridlines,width=1.5,link=curves,colour=red,colour=blue]')
call winop@('%pl[xaxis='Real(t)',yaxis='Zeta(s)']')
call winop@('%pl[Title='Riemann Zeta function (s = 0.5 + i*t)']')
call winop@('%pl[dx=10, dy=1]')
call winop@('%pl[x_min=0,x_max=100]')
call winop@('%pl[y_min=-3,y_max=5]')
iw = winio@('%pl&',1400,800,k,dble(vec),dble(real_part),dble(imaginary_part))
iw = winio@('')
Adding the reference to the symbol in winop and xip and yip in winio as you demonstrated did not seem to work for some reason?
Leaving the code to just show the real and imaginary plot is fine and I get the computed intersections.
Lester
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].
i = winio@('%pl',900,900,[n,n,n_intersects],dble(x1),dble(y1),dble(x2),dble(y2),dble(xip),dble(yip))
With the %pl option [x_array] and WITHOUT [independent]. This means that each of the three curves share a common x coordinate array. Thus there are four coordinate arrays passed as arguments to the winio call with %pl (i.e. x_common, y1, y2, y3). In this case the number of points in each curve is the same, and this is passed as a scalar.
You should also include the option [n_graphs=3].
Hi Kenneth
Managed to get the plot working, just followed you original syntax:
iw = winio@('%mn[Exit]&','Exit')
iw = winio@('%mn[Export]&',export_cb)
iw = winio@('%fn[Times New Roman]%ts%bf&',1.5d0)
!call winop@('%pl[native,x_array,frame,gridlines,width=1.5,link=curves,colour=red,colour=blue]')
call winop@('%pl[native,n_graphs=3,independent,frame,gridlines,smoothing=4,x_array]')
call winop@('%pl[link=curves,colour=red,symbol=0,link=curves,colour=blue,symbol=0,link=none,colour=black,symbol=6]')
call winop@('%pl[xaxis='Real(t)',yaxis='Zeta(s)']')
call winop@('%pl[Title='Riemann Zeta function (s = 0.5 + i*t)']')
call winop@('%pl[dx=10, dy=1]')
call winop@('%pl[x_min=0,x_max=100]')
call winop@('%pl[y_min=-3,y_max=5]')
!iw = winio@('%pl&',1400,800,k,dble(vec),dble(real_part),dble(imaginary_part))
iw = winio@('%pl',1400,800,[k,k,n_intersects],dble(vec),dble(real_part),dble(vec),dble(imaginary_part),dble(xip),dble(yip))
iw = winio@('')
real_part and imaginary_part have the same number of point, but n_intersects has fewer, so using independent as an option worked.
Thanks for the help
Lester