
forums.silverfrost.com Welcome to the Silverfrost forums

View previous topic :: View next topic 
Author 
Message 
arctica
Joined: 10 Sep 2006 Posts: 105 Location: United Kingdom

Posted: Thu Dec 28, 2023 9:35 am Post subject: Intersection of two curves 


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 

Back to top 


JohnCampbell
Joined: 16 Feb 2006 Posts: 2580 Location: Sydney

Posted: Thu Dec 28, 2023 12:30 pm Post subject: 


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,n11
do j = 1,n21 ;
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 

Back to top 


arctica
Joined: 10 Sep 2006 Posts: 105 Location: United Kingdom

Posted: Thu Dec 28, 2023 1:13 pm Post subject: 


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 (x1x2, y1y2) 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 

Back to top 


JohnCampbell
Joined: 16 Feb 2006 Posts: 2580 Location: Sydney

Posted: Sun Dec 31, 2023 3:41 am Post subject: 


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 ) 

Back to top 


arctica
Joined: 10 Sep 2006 Posts: 105 Location: United Kingdom

Posted: Mon Jan 08, 2024 2:48 pm Post subject: 


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.0e6 ! 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 

Back to top 


Kenneth_Smith
Joined: 18 May 2012 Posts: 709 Location: Hamilton, Lanarkshire, Scotland.

Posted: Tue Jan 09, 2024 10:47 am Post subject: 


My first attempt:
Code:  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 xy 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, n1, 1
do j = 1, n1, 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 


Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 1896

Posted: Tue Jan 09, 2024 3:35 pm Post subject: 


Praiseworthy, excellent graphical illustration! 

Back to top 


arctica
Joined: 10 Sep 2006 Posts: 105 Location: United Kingdom

Posted: Wed Jan 10, 2024 8:45 am Post subject: 


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:
Code: 
! 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, k1, 1
do j = 1, k1, 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 nontrivial 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 runtime error.
Code: 
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 

Back to top 


Kenneth_Smith
Joined: 18 May 2012 Posts: 709 Location: Hamilton, Lanarkshire, Scotland.

Posted: Wed Jan 10, 2024 11:50 am Post subject: 


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].
Code:  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]. 

Back to top 


arctica
Joined: 10 Sep 2006 Posts: 105 Location: United Kingdom

Posted: Thu Jan 11, 2024 1:34 pm Post subject: 


Hi Kenneth
Managed to get the plot working, just followed you original syntax:
Code: 
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 

Back to top 




You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum

Powered by phpBB © 2001, 2005 phpBB Group
