Announcement

**Mark Hunter** · 19 Sep 2009, 08:38 PM

Unfortunately I can’t edit the code I posted originally because the time limit has past.

The function “Astride” above will fail when one of the rays is vertical. This doesn’t seem to make any difference unless it occurs at the beginning, rather than at or near the end, of the closing-in process. It seems that the only case where it matters is when the Eye lies on the vertical line through the center of the circle, that is, when a = 0.

The algorithm is robust in that if a is extremely close to zero (but not zero) it will find the solution. For example I tried a = 10^-20 and it worked perfectly. Still, the algorithm fails if a = 0.

Here’s the program again, with that problem fixed.

Another point. As noted in the original, there’s no solution if Eye and Object are on opposite sides of the circle. They’re on opposite sides when the line segment between them intersects the circle (given they are both outside the circle). (Note that this is different from the line through them intersecting the circle.) This is not too hard to determine analytically, so one could avoid even trying to find a solution in that case. However, it’s kind of fun watching the algorithm fail when the Eye can’t see the Object, so I didn’t add that refinement.

Code:

' Compiler: PBCC
' Could use PBWIN by nixing the Console statement
' and using MsgBox for Waitkey$

' You have a circular mirror, the silver on its convex, bulging, side.
' Given an Object and an Eye, where on the mirror does the
' Eye see the Object?
'
' This point on the circle can be determined as follows.
'
' Start by finding the two points on the circle, call them L and R,
' where lines from the circle's center to the points E and O intersect
' the circle.  Then find the point between them, call it M, on the
' circle.
'
' Like the ancient Greeks did, think of an "eye-ray" going from
' the Eye to the mirror, reflecting so it hits the Object.  Real
' light rays simply backtrack eye-rays.

' Shoot an "eye-ray" from E to each of the three points L, M, and R,
' and consider the reflected rays.  A pair of them (from L&M or M&R)
' will straddle O. (Sometimes you'll have to consider fictional
' internal reflections.)
'
' Take the midpoint (on the circle) of this pair.  Again shoot the
' eye-rays and again choose the pair of points where the reflected
' rays straddle O.  Etc, dividing the circular arc in half each time.
' Thus you home in on the point on the circle which is the exact
' reflection point.  (Running the rays in reverse doesn't change
' the reflection point.)
'------------------------------------------------------------

#Compile Exe
#Dim All

Global a, b As Single       'Eye at (a,b)
Global c, d As Single       'Object at (c,d)

Global xxP As Single, yyP As Single  'head of reflected vector

%radius = 200               'radius of circle (cross-section of cylinder)
%radius2 = %radius*%radius  'radius squared
%dotradius = 3              'dot radius

%maxIterations = 20         'more than enough for tolerance (see Sub Found)
%delay0 = 250               'initial delay
%delay = 300                'delay between iterations
'------------------------------------------------------------

'length of vector
Function Norm(x As Single, y As Single) As Single
 Function = Sqr(x*x + y*y)
End Function
'------------------------------------------------------------

'normalize vector -- make it unit length
Sub Normalize(x As Single, y As Single)
 Local s As Single
 s = Sqr(x*x + y*y)
 x = x/s
 y = y/s
End Sub
'------------------------------------------------------------

'change length of vector to N
Sub Lengthen(x As Single, y As Single, N As Long)
 Local s As Single
 s = Sqr(x*x + y*y)
 x = x*N/s
 y = y*N/s
End Sub
'------------------------------------------------------------

' Given point (x0,y0) outside of circle
' Return projection onto circle (x,y)
' Note: x and y returned by reference
Sub Project(x0 As Single, y0 As Single, x As Single, y As Single)
 Local s As Single
 s = Norm(x0,y0)
 x = x0*%radius/s
 y = y0*%radius/s
End Sub
'------------------------------------------------------------

' Given two rays, one directed from (xs,ys) to (x1,y1) and beyond,
' the other from (xt,yt) to (x2,y2) and beyond,
' with the angle between them less than 180 degrees on the side facing
' the Object point (c,d).
' Return the result: Is the Object point (c,d) between these rays?
Function Astride(xs As Single, ys As Single, x1 As Single, y1 As Single, _
                 xt As Single, yt As Single, x2 As Single, y2 As Single) As Long
 Local ms, mt As Single                 'slopes
 Local xi, yi As Single                 'interesection of rays
 Local xx1, yy1, xx2, yy2 As Single     'translated to xi,yi
 Local cc, dd As Single                 'ditto
 Local xx, yy As Single                 'chord intersection
 Local d1, d2, d3 As Single             'distances
 'find back intersection of reflected eye rays
 'the slopes are ms and mt
 ' if ms = mt then no solution, but won't happen.
 'equation of first ray
 ' y = ms*(x - xs) + ys
 'second ray
 ' y = mt*(x - xt) + yt
 'Worry one of rays is vertical (can't both be)
 If Abs(x1 - xs) < 10^-20 Then          'x1 = xs
  xi = xs
  yi = ys
' elseif Abs(x2 - xt) < 10^-20 then      'x2 = xt
'  xi = xt
'  yi = (y1 - ys)*(xi - xs)/(x1 - xs) + ys
 Else
  ms = (y1 - ys)/(x1 - xs)
  mt = (y2 - yt)/(x2 - xt)
  xi = (ms*xs - mt*xt - ys + yt)/(ms - mt)
  yi = ms*(xi - xs) + ys
 End If
 'translate (xi,yi) to (0,0)
 xx1 = x1 - xi
 yy1 = y1 - yi
 xx2 = x2 - xi
 yy2 = y2 - yi
 cc = c - xi
 dd = d - yi
 'make first two the same length
 If Norm(xx1,yy1) > Norm(xx2,yy2) Then
  Lengthen xx1,yy1, Norm(xx2,yy2)
 Else
  Lengthen xx2,yy2, Norm(xx1,yy1)
 End If

'find intersection of (xi,yi)-(c,d) with (x1,y1)-(x2,y2)
'or rather
'(0,0)-(cc,dd) with (xx1,yy1)-(xx2,yy2)
' m = dd/cc
' mm = (yy2 - yy1)/(xx2 - xx1)
' xx = (m*cc - dd - mm*xx1 + yy1)/(m - mm)
' yy = m*(xx - cc) + dd
 If Abs(cc) > Abs(dd) Then          'worry about division by near 0
  xx = (xx2*yy1 - xx1*yy2)*cc/(dd*(xx2 - xx1) - cc*(yy2 - yy1))
  yy = dd*xx/cc
 Else
  yy = (xx2*yy1 - xx1*yy2)*dd/(dd*(xx2 - xx1) - cc*(yy2 - yy1))
  xx = cc*yy/dd
 End If
 'is xx,yy between (xx1,yy1) and (xx2,yy2) ?
 'same as: do distances add:  d1 + d2 = d3
 d1 = Norm(xx - xx1, yy - yy1)
 d2 = Norm(xx - xx2, yy - yy2)
 d3 = Norm(xx1 - xx2, yy1 - yy2)
 Function = (Abs(d1 + d2 - d3) < .1)
End Function
'------------------------------------------------------------

'This is called at the end of each iteration.
'No solution if Eye and Object are on different sides of circle.
'In that case reflection is internal instead of external
'Always give this function (x,y) = (xM,yM),
'on the circle near the point of reflection
Function NoSolution(x As Single, y As Single) As Long     'boolean
 Local xV,yV, xW,yW As Single
 Local xU,yU As Single
 'at this point
 'x,y is on circle, near the point of reflection
 'a,b is Eye
 'xxP,yyP  is NOT near Object (c,d)
 xV = a - x                    'turn into vectors
 yV = b - y
 xW = xxP - x
 yW = yyP - y
 Lengthen xV,yV, 1000          'make same length, say 1000
 Lengthen xW,yW, 1000
 xU = xV + xW                  'add them
 yU = yV + yW
 Lengthen xU,yU, %radius       'make length of circle radius
 xU = x + xU                   'move tail to (xM,yM)
 yU = y + yU
 'internal reflection if head towards center of circle (0,0)
 Function = (Abs(xU) < %radius*.5 And Abs(yU) < %radius*.5)
End Function
'------------------------------------------------------------

'Does the line from the reflection point (x,y) to (xxP,yyP)
'go through (within tolerance) the Object point (c,d) ?
Function Found(x As Single, y As Single) As Long
 Local AA, BB, CC As Single
 'equation of line in X,Y
 ' (yyP - y)*X - (xxP - x)*Y - (yyP - y)*x + (xxP - x)*y = 0
 AA = yyP - y
 BB = -xxP + x
 CC = -(yyP - y)*x + (xxP - x)*y
 Function = (Abs(AA*c + BB*d + CC)/Norm(AA,BB) <= .01)
End Function
'------------------------------------------------------------

' Given (x,y) on the circle,
' and hence the unit vector in the direction from the Eye to that point.
' Return the vector that is the reflection of this vector
' off the tangent at (x,y).
' Comments:
' Let V' = V reflected.
' Let u  = unit vector parallel to the tangent to circle at (x,y),
' either direction.
'  u = (y,-x) normalized = (y,-x)/%radius
' Let Vpara = V's component parallel to u.
' Let Vperp = V's component perpendicular to u.
' Then
'  V  = Vpara + Vperp
'  V' = Vpara - Vperp
' Thus
'  V' = 2Vpara - V
' Since Vpara = (V*u)u  where * is the dot product,
'  V' = 2(V*u)u - V
Sub Reflect(x As Single, y As Single, xR As Single, yR As Single)
 Local vx, vy As Single
 Local ux, uy As Single
 Local c As Single
 vx = x - a
 vy = y - b
 ux = y/%radius
 uy = -x/%radius
 c = 2*(vx*ux + vy*uy)
 xR = c*ux - vx
 yR = c*uy - vy
End Sub
'------------------------------------------------------------

' Given point on circle (x,y) and color,
' Shoot eye-ray at it, draw it and reflected ray.
Function Shoot(x As Single, y As Single, c As Long) As Long
 Local xP, yP As Single                 'reflection of an Eye-ray
 'The trouble with not using internal reflections is that
 'sometimes it would result in giving up too soon.
 Graphic Line(a,b)-(x,y), c
 Reflect x, y, xP, yP
 Normalize xP,yP
 xxP = x + xP*600                       'magnify so end beyond window
 yyP = y + yP*600
 Graphic Line(x,y)-(xxP,yyP), c
 Function = Found(x,y)
End Function
'------------------------------------------------------------

Sub DisplayFixedElements
 'Circular mirror
 Graphic Ellipse (-%radius, %radius) - (%radius, -%radius), %Red, %Red
 'Instructions
 Graphic Set Pos (-387,100-6):  Graphic Print "Press Esc to quit,"
 Graphic Set Pos (-387,100-23): Graphic Print "Enter or any other key"
 Graphic Set Pos (-387,100-40): Graphic Print "for new points."
 'Eye and Object
 Graphic Ellipse (a-%dotradius, b-%dotradius) - (a+%dotradius, b+%dotradius), %Blue, %Blue
 Graphic Ellipse (c-%dotradius, d-%dotradius) - (c+%dotradius, d+%dotradius), %Yellow, %Yellow
End Sub
'------------------------------------------------------------

Function PBMain
 Local winC As Dword                    'handle of window
 Local aa As Single, bb As Single       'projection of Eye
 Local cc As Single, dd As Single       'projection of Object
 Local x, y As Single                   'point on circle
 Local xL, yL As Single                 'point L at (xL,yL)
 Local xR, yR As Single                 'point R at (xR,yR)
 Local xM, yM As Single                 'midpoint M at (xM,yM)
 Local getkey As String                 'result of WaitKey$
 Local xxP_L, yyP_L As Single           'saved values of xxP, yyP
 Local xxP_R, yyP_R As Single
 Local xxP_M, yyP_M As Single
 Local count As Long                    'count iterations

 'set up 800x600 graphic window
 Graphic Window "", 0, 0, 800, 600 To winC
 Graphic Attach winC,0, ReDraw

 'origin (0,0) at middle, 100 down from top
 Graphic Scale(-%radius*2,100)-(%radius*2,-500)

 Console Set Loc 0,600                  'hide console
 Console Set Focus                      'but give it focus for keystrokes

 Randomize Timer
 Graphic Font "Times New Roman", 12, 0

Start:                                   'main loop

 Do
  a = Rnd(-400+10,400-10)                'random Eye point
  b = Rnd(-500+10,100-10)
 Loop Until a*a + b*b > %radius2 + 10    'outside of circle

 Do
  c = Rnd(-400+10,400-10)                'random Object point
  d = Rnd(-500+10,100-10)
 Loop Until c*c + d*d > %radius2 + 10    'outside of circle

 DisplayFixedElements
 Graphic ReDraw


 Sleep %delay0

 'determine starting points L and M by
 'projecting points Eye and Object onto circle
 Project a, b, aa, bb
 xL = aa
 yL = bb
 Project c, d, cc, dd
 xR = cc
 yR = dd

 For count = 1 To %maxIterations
  Graphic Clear
  DisplayFixedElements

  'Given the two points L = (xL,yL) and M = (xR,yR) on the circle
  'find the midpoint between them (x,y) on the circle
  Project xL + xR, yL + yR, xM, yM       'returns (xM,yM)

  'Now have three points on circle L, M, and R.

  If Shoot(xL,yL,%Green) GoTo LastRay    'shoot ray to L
  xxP_L = xxP                            'save (xxP,yyP)
  yyP_L = yyP

  If Shoot(xM,yM,%Magenta) GoTo LastRay  'shoot ray to M
  xxP_M = xxP                            'save (xxP,yyP)
  yyP_M = yyP

  If Astride(xL,yL, xxP_L,yyP_L, xM,yM, xxP_M,yyP_M) Then   'between L and M ?
   xR = xM
   yR = yM
  Else                                   'try R and M
   If Shoot(xR,yR,%Gray) Then GoTo Lastray  'shoot ray to R
   Graphic Set Pos (0,-140)
   If Astride(xR,yR, xxP_R,yyP_R, xM,yM, xxP_M,yyP_M) Then  'between R and M ?
    xL = xM
    yL = yM
   Else
    Graphic Set Pos (-20,-140)           'should never get here
    Graphic Print "Failed"
    Graphic ReDraw
    GoTo Done
   End If
  End If

  Graphic ReDraw
  If InStat Then GoTo Done
  Sleep %delay/count        'less delay between iterations as we progress
 Next

LastRay:
 If NoSolution(xM,yM) Then       'Eye and Object on different sides of circle?
  Graphic Clear
  DisplayFixedElements
  Graphic Set Pos (-35,-140)
  Graphic Print " No solution "
  Graphic ReDraw
  GoTo Done
 End If

 Graphic Clear
 DisplayFixedElements
 Shoot xM,yM, %White
 Graphic ReDraw

Done:
 getkey = WaitKey$
 Graphic Clear
 If getkey$ <> Chr$(27) Then GoTo Start
End Function

Announcement

Forum Guidelines

Ray Tracing for a Circular Mirror

PBCC Ray Tracing for a Circular Mirror

Comment