soccer jersey forums.silverfrost.com :: View topic - Intersection of two curves
forums.silverfrost.com Forum Index forums.silverfrost.com
Welcome to the Silverfrost forums
 
 FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

Intersection of two curves

 
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> General
View previous topic :: View next topic  
Author Message
arctica



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

PostPosted: Thu Dec 28, 2023 9:35 am    Post subject: Intersection of two curves Reply with quote

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
View user's profile Send private message
JohnCampbell



Joined: 16 Feb 2006
Posts: 2584
Location: Sydney

PostPosted: Thu Dec 28, 2023 12:30 pm    Post subject: Reply with quote

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
Back to top
View user's profile Send private message
arctica



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

PostPosted: Thu Dec 28, 2023 1:13 pm    Post subject: Reply with quote

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
Back to top
View user's profile Send private message
JohnCampbell



Joined: 16 Feb 2006
Posts: 2584
Location: Sydney

PostPosted: Sun Dec 31, 2023 3:41 am    Post subject: Reply with quote

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
View user's profile Send private message
arctica



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

PostPosted: Mon Jan 08, 2024 2:48 pm    Post subject: Reply with quote

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(Smile, intent(in) :: x1, y1, x2, y2
real, dimension(Smile, 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
Back to top
View user's profile Send private message
Kenneth_Smith



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

PostPosted: Tue Jan 09, 2024 10:47 am    Post subject: Reply with quote

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 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
Back to top
View user's profile Send private message Visit poster's website
mecej4



Joined: 31 Oct 2006
Posts: 1896

PostPosted: Tue Jan 09, 2024 3:35 pm    Post subject: Reply with quote

Praiseworthy, excellent graphical illustration!
Back to top
View user's profile Send private message
arctica



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

PostPosted: Wed Jan 10, 2024 8:45 am    Post subject: Reply with quote

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, 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.

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
View user's profile Send private message
Kenneth_Smith



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

PostPosted: Wed Jan 10, 2024 11:50 am    Post subject: Reply with quote

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
View user's profile Send private message Visit poster's website
arctica



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

PostPosted: Thu Jan 11, 2024 1:34 pm    Post subject: Reply with quote

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
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> General All times are GMT + 1 Hour
Page 1 of 1

 
Jump to:  
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