QB64.org Forum

Active Forums => Programs => Topic started by: STxAxTIC on January 31, 2020, 12:38:04 pm

Title: A Two-Roads Problem
Post by: STxAxTIC on January 31, 2020, 12:38:04 pm
Hi folks. So I was pondering a statement made by OldMoses about bplus's work on tangents https://www.qb64.org/forum/index.php?topic=2132.0 (https://www.qb64.org/forum/index.php?topic=2132.0):

Quote
Nice. I think I'll file this one away as a possible basis for an orbital insertion routine.

So I decided to look at this. Along the way, I solved a toy problem that I think some of you would enjoy... as a challenge first if you aren't sick of them. I can imagine bplus and Steve looking at this a few different ways, I do wonder would OldMoses might come up with.

So the challenge is like this: Consider two curved roads or train tracks, modeled by parabolas to make it easy:

Code: QB64: [Select]
  1. FUNCTION ya (x)
  2.     ya = -x ^ 2 / 4 - 1
  3.  
  4. FUNCTION yb (x)
  5.     yb = x ^ 2 / 4 + 1

All that does is define two arc-shaped roads that don't intersect. So far so good?

Now consider the two points (-2,-2) and (1,1.25). Each of these points is on one of the roads. The job is to connect the two roads at the two points, but smoothly.

The code below produces the screenshot attached. I made a crude guess that seems to work for the first curve, but misses the target on the second curve. How do we connect the curves with a curvy line that joins each road smoothly?

Code: QB64: [Select]
  1.  
  2. CALL cline(0, 0, _WIDTH / 2, 0, 7)
  3. CALL cline(0, 0, -_WIDTH / 2, 0, 7)
  4. CALL cline(0, 0, 0, _HEIGHT / 2, 7)
  5. CALL cline(0, 0, 0, -_HEIGHT / 2, 7)
  6.  
  7. zoom = 75
  8. CALL ccircle(-2 * zoom, ya(-2) * zoom, 10, 14)
  9. CALL ccircle(1 * zoom, yb(1) * zoom, 10, 15)
  10. FOR x = -2.5 TO 2 STEP .01
  11.     CALL ccircle(x * zoom, ya(x) * zoom, 1, 14)
  12.     CALL ccircle(x * zoom, yb(x) * zoom, 1, 15)
  13.     CALL ccircle(x * zoom, yc(x) * zoom, 1, 5)
  14.  
  15.  
  16. FUNCTION ya (x)
  17.     ya = -x ^ 2 / 4 - 1
  18.  
  19. FUNCTION yb (x)
  20.     yb = x ^ 2 / 4 + 1
  21.  
  22. FUNCTION yc (x)
  23.     yc = x ' What should this be in order to smoothly join each curve?
  24.  
  25. SUB cline (x1, y1, x2, y2, col)
  26.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
  27.  
  28. SUB ccircle (x1, y1, r, col)
  29.     CIRCLE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), r, col
  30.  

Title: Re: A Two-Roads Problem
Post by: bplus on January 31, 2020, 12:52:56 pm
Well I think I've solved the problem of what a "cured" road is, it is one without a v in it. ;-))
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on January 31, 2020, 12:55:08 pm
(Damn you shoulda quoted cause I fixed that right away.)
Title: Re: A Two-Roads Problem
Post by: SMcNeill on January 31, 2020, 01:20:35 pm
Here’s an idea:   Draw a straight line from the top circle, and a straight line from the bottom circle.  Then choose a midpoint between those lines and create a circle. 

To illustrate: In your picture above, you’d have what looks like the end of a paperclip slid between the two points, with the “U” being off to the left side of the two arcs.

Would such a solution as the above be considered valid for smoothly connecting the two points?  It seems like a very simple solution to the problem at hand.
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on January 31, 2020, 03:28:03 pm
Nice thinking Steve, but theres a problem with that proposed curve in that it isnt smooth. In a car driving analogy, note that the steering wheel would have to jerk *instantly* when you transition from the line to the circle. A real road/car would forbid an instant jump like that, so you need to enter/exit the circle more gradually than a brutal tangent line. This is especially evident for smaller circles but true for all circles. That doesnt mean a part of the solution cant be line-like or circle-like though.

Said another way, you want a curve to join at the points, yes. You also want the slope of the curve to match the slope of the road where it joins, check. Here's the subtle part: the slope of the slope of the curve must also match the slope of the slope of the road at the intersection points. This guarantees a jerk-free transition.
Title: Re: A Two-Roads Problem
Post by: OldMoses on January 31, 2020, 04:43:15 pm
OMG, did I start something I can't finish?

I'm thinking something along the lines of a cube function, but getting it to fit is a bit beyond my pay grade.

Code: QB64: [Select]
  1. A& = _NEWIMAGE(650, 650, 32)
  2. _TITLE "Graph"
  3.  
  4. WINDOW (-15, 15)-(15, -15)
  5. g = -16
  6.     g = g + 1
  7.     LINE (-15, g)-(15, g), &H7F7F7F7F '                         horizontal grid lines
  8.     LINE (g, -15)-(g, 15), &H7F7F7F7F '
  9. LOOP UNTIL g = 15
  10. LINE (0, -15)-(0, 15), &HFF0000FF, BF
  11. LINE (-15, 0)-(15, 0), &HFF0000FF, BF
  12. FOR x = -15 TO 15 STEP .01
  13.     y = (.5 * x) ^ 3
  14.     PSET (x, y)
  15.  
  16.  
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on February 01, 2020, 12:01:01 am
I'll post the answer as a hint, but this is by no means a solution.

Code: QB64: [Select]
  1.  
  2. CALL cline(0, 0, _WIDTH / 2, 0, 7)
  3. CALL cline(0, 0, -_WIDTH / 2, 0, 7)
  4. CALL cline(0, 0, 0, _HEIGHT / 2, 7)
  5. CALL cline(0, 0, 0, -_HEIGHT / 2, 7)
  6.  
  7. zoom = 75
  8. CALL ccircle(-2 * zoom, ya(-2) * zoom, 10, 14)
  9. CALL ccircle(1 * zoom, yb(1) * zoom, 10, 15)
  10. FOR x = -2.5 TO 2 STEP .01
  11.     CALL ccircle(x * zoom, ya(x) * zoom, 1, 14)
  12.     CALL ccircle(x * zoom, yb(x) * zoom, 1, 15)
  13.     CALL ccircle(x * zoom, yc(x) * zoom, 1, 5)
  14.  
  15.  
  16. FUNCTION ya (x)
  17.     ya = -x ^ 2 / 4 - 1
  18.  
  19. FUNCTION yb (x)
  20.     yb = x ^ 2 / 4 + 1
  21.  
  22. FUNCTION yc (x)
  23.     yc = 0.5123456790123457 + 1.1419753086419753 * x - 0.40895061728395077 * x ^ 2 - 0.12345679012345677 * x ^ 3 + 0.09413580246913578 * x ^ 4 + 0.03395061728395061 * x ^ 5
  24.  
  25. SUB cline (x1, y1, x2, y2, col)
  26.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
  27.  
  28. SUB ccircle (x1, y1, r, col)
  29.     CIRCLE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), r, col
Title: Re: A Two-Roads Problem
Post by: SMcNeill on February 02, 2020, 01:05:35 pm
Code: [Select]
SCREEN _NEWIMAGE(640, 480, 32)
WINDOW (-3.2, -2.4)-(3.2, 2.4)
FOR x = -3.2 TO 3.2 STEP .01 'The two arcs
    PSET (x, ya(x)), -1
    PSET (x, yb(x)), -1
NEXT

LowPoint = 2.4: HighPoint = -2.4
FOR y = 2.4 TO 0 STEP -0.01
    FOR x = -3.2 TO 3.2 STEP 0.01
        IF POINT(x, y) <> 0 AND y < LowPoint THEN LowPoint = y
    NEXT
NEXT
FOR y = -2.4 TO 0 STEP 0.01
    FOR x = -3.2 TO 3.2 STEP 0.01
        IF POINT(x, y) <> 0 AND y > HighPoint THEN HighPoint = y
    NEXT
NEXT

Midpoint = (HighPoint - LowPoint) / 2

CIRCLE (0, LowPoint + Midpoint / 2), Midpoint / 2
CIRCLE (0, HighPoint - Midpoint / 2), Midpoint / 2






CIRCLE (-2, -2), .1, &HFFFF0000 'one point
CIRCLE (1, 1.25), .1, &HFFFF0000 'second point

SLEEP

FUNCTION ya (x)
    ya = -x ^ 2 / 4 - 1
END FUNCTION

FUNCTION yb (x)
    yb = x ^ 2 / 4 + 1
END FUNCTION

Start at the bottom left, follow the arc to your first point and board the train to the apex and get on the circle.  Continue to the next circle and get on it.  Exit the circle at the apex of the top arc and continue on to your destination at the second point.
Title: Re: A Two-Roads Problem
Post by: OldMoses on February 02, 2020, 01:12:32 pm
Bravo! That's what I call thinking out of the box.
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on February 02, 2020, 01:13:32 pm
*facepalm*

Not because it's wrong (you are supposed to intersect at the points), but I explained this more than once (not shown here). You could have described this without the circles ya know... Isn't that almost what you did above?

Make sure your method can start anywhere on curve one, and end anywhere on curve 2.

OldMoses: You really want your orbital simulation doing stuff like that? Lol
Title: Re: A Two-Roads Problem
Post by: SMcNeill on February 02, 2020, 01:26:06 pm
*facepalm*

Not because it's wrong (you are supposed to intersect at the points), but I explained this more than once (not shown here). You could have described this without the circles ya know... Isn't that almost what you did above?

Make sure your method can start anywhere on curve one, and end anywhere on curve 2.

Why wouldn't it work with any spot on curve one and end anywhere on curve 2?  Unless you want to go from the same side to the same side?  In that case, you just build a third circle to connect the other two and you use it as a direction changer.
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on February 02, 2020, 01:34:56 pm
I see the construction you're going for - not to shit on it. But I kinda got rid of circles already if the point wasn't too subtle. I will try to reiterate why your path still isn't smooth.

The instantaneous curvature of the curve your're on is equal to one divided the radius of an imaginary circle (somewhere in space) whose edge instantaneously coincides with your curve. Straight lines have zero curvature, as a circle of infinite radius parked very far away would be needed to match its edge to a straight line.

Okay. So when you jump from a straight line to a curve, you *instantly* change curvature from 0 to 1/r. Can't have instant jumps. It gets even worse when you go from one circle to another circle that bends the other way, because your curvature jump is 1/r1 + 1/r2. We need a path that is truly smooth.

Hint: The actual answer is above. See if you can get close.
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on February 02, 2020, 01:58:57 pm
Here's the same problem with less unusable road:
Title: Re: A Two-Roads Problem
Post by: SMcNeill on February 02, 2020, 02:02:29 pm
Here's the same problem with less unusable road:

And here, we just use our warp drive, fold dimensional space in two and stick a pencil in it, and go instantly from point A to point B at whatever heading, direction, and speed which we want to arrive at...  :P
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on February 02, 2020, 02:21:39 pm
Okay DarkStar :-)
Title: Re: A Two-Roads Problem
Post by: bplus on February 02, 2020, 02:37:09 pm
Why wouldn't it work with any spot on curve one and end anywhere on curve 2?  Unless you want to go from the same side to the same side?  In that case, you just build a third circle to connect the other two and you use it as a direction changer.

RE: Steve's reply #10 with 3 circles
Hey it works! There is no jump in slope, they are always matched at tangent points before changing curves PLUS you get from any point on the curve to any other without building a special road for any 2 points. Did anyone ask for the shortest curve?

I tried to fit a chunk of a parabola from one point to the other but couldn't get the numbers to fit, maybe a chunk of ellipse? Wasted 2 days ;(
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on February 02, 2020, 02:44:03 pm
The shortest curve would be mostly straight with just a little curvature at each intersection point.

bplus: Matching just the slope is insufficient. You need to match the slope of the slope to experience no unjustified accelerations.
Title: Re: A Two-Roads Problem
Post by: bplus on February 02, 2020, 04:16:34 pm
Well how smooth is this?
Code: QB64: [Select]
  1. CONST xmax = 600, ymax = 600
  2. SCREEN _NEWIMAGE(xmax, ymax, 32)
  3. WINDOW (-10, 10)-(10, -10)
  4. TYPE XY
  5.     x AS SINGLE
  6.     y AS SINGLE
  7.  
  8. REDIM arr(19) AS XY
  9.  
  10. FOR x = -10 TO 10 STEP .1
  11.     IF x > -2.5 AND x <= -1.5 THEN 'collect points around jump from red
  12.         arr(i).x = x
  13.         arr(i).y = ya(x)
  14.         i = i + 1
  15.     END IF
  16.     IF x >= .5 AND x < 1.5 THEN 'collect points around jump to blue
  17.         arr(i).x = x
  18.         arr(i).y = yb(x)
  19.         i = i + 1
  20.     END IF
  21.     IF x > -10 THEN
  22.         LINE (x - 1, ya(x - 1))-STEP(1, ya(x) - ya(x - 1)), &HFFFF0000
  23.         LINE (x - 1, yb(x - 1))-STEP(1, yb(x) - yb(x - 1)), &HFF0000FF
  24.         'LINE (x - 1, soln(x - 1))-STEP(1, soln(x) - soln(x - 1)), &HFF00FF00 'test soln
  25.     END IF
  26. 'connect (-2,ya(-2)) with (1, yb(1))
  27. CIRCLE (-2, ya(-2)), .1: CIRCLE (1, yb(1)), .1
  28. PRINT "yadx(-2)  ="; yadx(-2); "  ybdx(1) ="; ybdx(1) '1 , .5
  29. LINE (-1, -1)-(-3, -3) 'slope over  -2, -2
  30. LINE (0, .75)-(2, 1.75) 'slope over 1, 1.25
  31. SMOOTH arr(), 200, 100
  32. FOR i = 0 TO 200
  33.     CIRCLE (arr(i).x, arr(i).y), .05, &HFF008800
  34.  
  35. FUNCTION ya (x)
  36.     ya = -(x ^ 2) / 4 - 1
  37.  
  38. FUNCTION yadx (x)
  39.     yadx = -x / 2
  40.  
  41. FUNCTION yb (x)
  42.     yb = x ^ 2 / 4 + 1
  43.  
  44. FUNCTION ybdx (x)
  45.     ybdx = x / 2
  46.  
  47. FUNCTION soln (x) 'failed to eyeball in the right curve fit
  48.     soln = -((x - 4) ^ 2) / 12 + 2
  49.  
  50. '======================= Feature SUB =======================================================================
  51. ' This code takes a dynamic points array and adds and modifies points to smooth out the data,
  52. ' to be used as Toolbox SUB. b+ 2020-01-24 adapted and modified from:
  53. ' Curve smoother by STxAxTIC https://www.qb64.org/forum/index.php?topic=184.msg963#msg963
  54. SUB SMOOTH (arr() AS XY, targetPoints AS INTEGER, smoothIterations AS INTEGER)
  55.     'TYPE XY
  56.     '    x AS SINGLE
  57.     '    y AS SINGLE
  58.     'END TYPE
  59.     ' targetPoints is the number of points to be in finished smoothed out array
  60.     ' smoothIterations is number of times to try and round out corners
  61.  
  62.     DIM rad2Max, kmax, k, numPoints, xfac, yfac, rad2, j
  63.     numPoints = UBOUND(arr)
  64.     REDIM _PRESERVE arr(0 TO targetPoints) AS XY
  65.     REDIM temp(0 TO targetPoints) AS XY
  66.     DO
  67.         '
  68.         ' Determine the pair of neighboring points that have the greatest separation of all pairs.
  69.         '
  70.         rad2Max = -1
  71.         kmax = -1
  72.         FOR k = 1 TO numPoints - 1
  73.             xfac = arr(k).x - arr(k + 1).x
  74.             yfac = arr(k).y - arr(k + 1).y
  75.             rad2 = xfac ^ 2 + yfac ^ 2
  76.             IF rad2 > rad2Max THEN
  77.                 kmax = k
  78.                 rad2Max = rad2
  79.             END IF
  80.         NEXT
  81.         '
  82.         ' Starting next to kmax, create a `gap' by shifting all other points by one index.
  83.         '
  84.         FOR j = numPoints TO kmax + 1 STEP -1
  85.             arr(j + 1).x = arr(j).x
  86.             arr(j + 1).y = arr(j).y
  87.         NEXT
  88.  
  89.         '
  90.         ' Fill the gap with a new point whose position is determined by the average of its neighbors.
  91.         '
  92.         arr(kmax + 1).x = .5 * (arr(kmax).x + arr(kmax + 2).x)
  93.         arr(kmax + 1).y = .5 * (arr(kmax).y + arr(kmax + 2).y)
  94.  
  95.         numPoints = numPoints + 1
  96.     LOOP UNTIL (numPoints = targetPoints)
  97.     '
  98.     ' At this stage, the curve still has all of its sharp edges. Use a `relaxation method' to smooth.
  99.     ' The new position of a point is equal to the average position of its neighboring points.
  100.     '
  101.     FOR j = 1 TO smoothIterations
  102.         FOR k = 2 TO numPoints - 1
  103.             temp(k).x = .5 * (arr(k - 1).x + arr(k + 1).x)
  104.             temp(k).y = .5 * (arr(k - 1).y + arr(k + 1).y)
  105.         NEXT
  106.         FOR k = 2 TO numPoints - 1
  107.             arr(k).x = temp(k).x
  108.             arr(k).y = temp(k).y
  109.         NEXT
  110.     NEXT
  111.  
  112.  
  113.  

  [ This attachment cannot be displayed inline in 'Print Page' view ]  
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on February 02, 2020, 04:21:29 pm
This is why I love you bplus.

That's as close to a perfectly acceptable solution as you can get. I bet it overlays the one I provided perfectly. Bravo, sir. You always find a way!

And hey, where'd you get that smoothing code? I literally never saw that coming.

XD

EDIT:

Here's your code with my solution plotted as well:

Code: QB64: [Select]
  1. CONST xmax = 600, ymax = 600
  2. SCREEN _NEWIMAGE(xmax, ymax, 32)
  3. WINDOW (-10, 10)-(10, -10)
  4. TYPE XY
  5.     x AS SINGLE
  6.     y AS SINGLE
  7.  
  8. REDIM arr(19) AS XY
  9.  
  10. FOR x = -10 TO 10 STEP .1
  11.     IF x > -2.5 AND x <= -1.5 THEN 'collect points around jump from red
  12.         arr(i).x = x
  13.         arr(i).y = ya(x)
  14.         i = i + 1
  15.     END IF
  16.     IF x >= .5 AND x < 1.5 THEN 'collect points around jump to blue
  17.         arr(i).x = x
  18.         arr(i).y = yb(x)
  19.         i = i + 1
  20.     END IF
  21.     IF x > -10 THEN
  22.         LINE (x - 1, ya(x - 1))-STEP(1, ya(x) - ya(x - 1)), &HFFFF0000
  23.         LINE (x - 1, yb(x - 1))-STEP(1, yb(x) - yb(x - 1)), &HFF0000FF
  24.         LINE (x - 1, soln(x - 1))-STEP(1, soln(x) - soln(x - 1)), &HFF00FF00 'test soln
  25.     END IF
  26. 'connect (-2,ya(-2)) with (1, yb(1))
  27. CIRCLE (-2, ya(-2)), .1: CIRCLE (1, yb(1)), .1
  28. PRINT "yadx(-2)  ="; yadx(-2); "  ybdx(1) ="; ybdx(1) '1 , .5
  29. LINE (-1, -1)-(-3, -3) 'slope over  -2, -2
  30. LINE (0, .75)-(2, 1.75) 'slope over 1, 1.25
  31. SMOOTH arr(), 200, 100
  32. FOR i = 0 TO 200
  33.     CIRCLE (arr(i).x, arr(i).y), .05, &HFF008800
  34.  
  35. FUNCTION ya (x)
  36.     ya = -(x ^ 2) / 4 - 1
  37.  
  38. FUNCTION yadx (x)
  39.     yadx = -x / 2
  40.  
  41. FUNCTION yb (x)
  42.     yb = x ^ 2 / 4 + 1
  43.  
  44. FUNCTION ybdx (x)
  45.     ybdx = x / 2
  46.  
  47. FUNCTION soln (x) 'failed to eyeball in the right curve fit
  48.     yc = 0.5123456790123457 + 1.1419753086419753 * x - 0.40895061728395077 * x ^ 2 - 0.12345679012345677 * x ^ 3 + 0.09413580246913578 * x ^ 4 + 0.03395061728395061 * x ^ 5
  49.     soln = yc
  50.  
  51. '======================= Feature SUB =======================================================================
  52. ' This code takes a dynamic points array and adds and modifies points to smooth out the data,
  53. ' to be used as Toolbox SUB. b+ 2020-01-24 adapted and modified from:
  54. ' Curve smoother by STxAxTIC https://www.qb64.org/forum/index.php?topic=184.msg963#msg963
  55. SUB SMOOTH (arr() AS XY, targetPoints AS INTEGER, smoothIterations AS INTEGER)
  56.     'TYPE XY
  57.     '    x AS SINGLE
  58.     '    y AS SINGLE
  59.     'END TYPE
  60.     ' targetPoints is the number of points to be in finished smoothed out array
  61.     ' smoothIterations is number of times to try and round out corners
  62.  
  63.     DIM rad2Max, kmax, k, numPoints, xfac, yfac, rad2, j
  64.     numPoints = UBOUND(arr)
  65.     REDIM _PRESERVE arr(0 TO targetPoints) AS XY
  66.     REDIM temp(0 TO targetPoints) AS XY
  67.     DO
  68.         '
  69.         ' Determine the pair of neighboring points that have the greatest separation of all pairs.
  70.         '
  71.         rad2Max = -1
  72.         kmax = -1
  73.         FOR k = 1 TO numPoints - 1
  74.             xfac = arr(k).x - arr(k + 1).x
  75.             yfac = arr(k).y - arr(k + 1).y
  76.             rad2 = xfac ^ 2 + yfac ^ 2
  77.             IF rad2 > rad2Max THEN
  78.                 kmax = k
  79.                 rad2Max = rad2
  80.             END IF
  81.         NEXT
  82.         '
  83.         ' Starting next to kmax, create a `gap' by shifting all other points by one index.
  84.         '
  85.         FOR j = numPoints TO kmax + 1 STEP -1
  86.             arr(j + 1).x = arr(j).x
  87.             arr(j + 1).y = arr(j).y
  88.         NEXT
  89.  
  90.         '
  91.         ' Fill the gap with a new point whose position is determined by the average of its neighbors.
  92.         '
  93.         arr(kmax + 1).x = .5 * (arr(kmax).x + arr(kmax + 2).x)
  94.         arr(kmax + 1).y = .5 * (arr(kmax).y + arr(kmax + 2).y)
  95.  
  96.         numPoints = numPoints + 1
  97.     LOOP UNTIL (numPoints = targetPoints)
  98.     '
  99.     ' At this stage, the curve still has all of its sharp edges. Use a `relaxation method' to smooth.
  100.     ' The new position of a point is equal to the average position of its neighboring points.
  101.     '
  102.     FOR j = 1 TO smoothIterations
  103.         FOR k = 2 TO numPoints - 1
  104.             temp(k).x = .5 * (arr(k - 1).x + arr(k + 1).x)
  105.             temp(k).y = .5 * (arr(k - 1).y + arr(k + 1).y)
  106.         NEXT
  107.         FOR k = 2 TO numPoints - 1
  108.             arr(k).x = temp(k).x
  109.             arr(k).y = temp(k).y
  110.         NEXT
  111.     NEXT
  112.  
Title: Re: A Two-Roads Problem
Post by: bplus on February 02, 2020, 04:36:38 pm
LOL STxAxTIC, how could you not like your own work?

I made a close up or zoom in:
Code: QB64: [Select]
  1. CONST xmax = 600, ymax = 600
  2. SCREEN _NEWIMAGE(xmax, ymax, 32)
  3. WINDOW (-3, 3)-(3, -3)
  4. TYPE XY
  5.     x AS SINGLE
  6.     y AS SINGLE
  7. REDIM arr(19) AS XY
  8. FOR x = -3 TO 3 STEP .1
  9.     IF x > -2.5 AND x <= -1.5 THEN 'collect points around jump from red
  10.         arr(i).x = x
  11.         arr(i).y = ya(x)
  12.         i = i + 1
  13.     END IF
  14.     IF x >= .5 AND x < 1.5 THEN 'collect points around jump to blue
  15.         arr(i).x = x
  16.         arr(i).y = yb(x)
  17.         i = i + 1
  18.     END IF
  19.     IF x > -3 THEN
  20.         LINE (x - .1, ya(x - .1))-(x, ya(x)), &HFFFF0000
  21.         LINE (x - .1, yb(x - .1))-(x, yb(x)), &HFF0000FF
  22.     END IF
  23. 'connect (-2,ya(-2)) with (1, yb(1))
  24. CIRCLE (-2, ya(-2)), .1: CIRCLE (1, yb(1)), .1
  25. PRINT "yadx(-2)  ="; yadx(-2); "  ybdx(1) ="; ybdx(1) '1 , .5
  26. LINE (-1, -1)-(-3, -3) 'slope over  -2, -2
  27. LINE (0, .75)-(2, 1.75) 'slope over 1, 1.25
  28. SMOOTH arr(), 200, 100
  29. FOR i = 1 TO 200
  30.     LINE (arr(i - 1).x, arr(i - 1).y)-(arr(i).x, arr(i).y), &HFFFFFF00
  31.  
  32. FUNCTION ya (x)
  33.     ya = -(x ^ 2) / 4 - 1
  34.  
  35. FUNCTION yadx (x)
  36.     yadx = -x / 2
  37.  
  38. FUNCTION yb (x)
  39.     yb = x ^ 2 / 4 + 1
  40.  
  41. FUNCTION ybdx (x)
  42.     ybdx = x / 2
  43.  
  44. FUNCTION soln (x) 'failed to eyeball in the right curve fit
  45.     soln = -((x - 4) ^ 2) / 12 + 2
  46.  
  47. '======================= Feature SUB =======================================================================
  48. ' This code takes a dynamic points array and adds and modifies points to smooth out the data,
  49. ' to be used as Toolbox SUB. b+ 2020-01-24 adapted and modified from:
  50. ' Curve smoother by STxAxTIC https://www.qb64.org/forum/index.php?topic=184.msg963#msg963
  51. SUB SMOOTH (arr() AS XY, targetPoints AS INTEGER, smoothIterations AS INTEGER)
  52.     'TYPE XY
  53.     '    x AS SINGLE
  54.     '    y AS SINGLE
  55.     'END TYPE
  56.     ' targetPoints is the number of points to be in finished smoothed out array
  57.     ' smoothIterations is number of times to try and round out corners
  58.  
  59.     DIM rad2Max, kmax, k, numPoints, xfac, yfac, rad2, j
  60.     numPoints = UBOUND(arr)
  61.     REDIM _PRESERVE arr(0 TO targetPoints) AS XY
  62.     REDIM temp(0 TO targetPoints) AS XY
  63.     DO
  64.         '
  65.         ' Determine the pair of neighboring points that have the greatest separation of all pairs.
  66.         '
  67.         rad2Max = -1
  68.         kmax = -1
  69.         FOR k = 1 TO numPoints - 1
  70.             xfac = arr(k).x - arr(k + 1).x
  71.             yfac = arr(k).y - arr(k + 1).y
  72.             rad2 = xfac ^ 2 + yfac ^ 2
  73.             IF rad2 > rad2Max THEN
  74.                 kmax = k
  75.                 rad2Max = rad2
  76.             END IF
  77.         NEXT
  78.         '
  79.         ' Starting next to kmax, create a `gap' by shifting all other points by one index.
  80.         '
  81.         FOR j = numPoints TO kmax + 1 STEP -1
  82.             arr(j + 1).x = arr(j).x
  83.             arr(j + 1).y = arr(j).y
  84.         NEXT
  85.  
  86.         '
  87.         ' Fill the gap with a new point whose position is determined by the average of its neighbors.
  88.         '
  89.         arr(kmax + 1).x = .5 * (arr(kmax).x + arr(kmax + 2).x)
  90.         arr(kmax + 1).y = .5 * (arr(kmax).y + arr(kmax + 2).y)
  91.  
  92.         numPoints = numPoints + 1
  93.     LOOP UNTIL (numPoints = targetPoints)
  94.     '
  95.     ' At this stage, the curve still has all of its sharp edges. Use a `relaxation method' to smooth.
  96.     ' The new position of a point is equal to the average position of its neighboring points.
  97.     '
  98.     FOR j = 1 TO smoothIterations
  99.         FOR k = 2 TO numPoints - 1
  100.             temp(k).x = .5 * (arr(k - 1).x + arr(k + 1).x)
  101.             temp(k).y = .5 * (arr(k - 1).y + arr(k + 1).y)
  102.         NEXT
  103.         FOR k = 2 TO numPoints - 1
  104.             arr(k).x = temp(k).x
  105.             arr(k).y = temp(k).y
  106.         NEXT
  107.     NEXT
  108.  
  109.  
  110.  

Title: Re: A Two-Roads Problem
Post by: STxAxTIC on February 07, 2020, 12:52:14 am
Alright, sorry this took so long - here is the long awaited (by nobody by now I bet) solution to the two-roads problem. It brushes over calculus and linear algebra without apology. Scroll to the end for a picture.

Ah, and my numerical solution came from sxript code:

Code: QB64: [Select]
  1. include(`../test/linalg.txt'):
  2.  
  3. let(q1,-2):
  4. let(f0,-2):
  5. let(f1,1):
  6. let(f2,-1/2):
  7.  
  8. let(q2,1):
  9. let(g0,5/4):
  10. let(g1,1/2):
  11. let(g2,1/2):
  12.  
  13. let(d,
  14.   <
  15.   <1,[q1],[q1]^2,[q1]^3,[q1]^4,[q1]^5,    [f0]>,
  16.   <1,[q2],[q2]^2,[q2]^3,[q2]^4,[q2]^5,    [g0]>,
  17.   <0,1,2*[q1],3*[q1]^2,4*[q1]^3,5*[q1]^4, [f1]>,
  18.   <0,1,2*[q2],3*[q2]^2,4*[q2]^3,5*[q2]^4, [g1]>,
  19.   <0,0,2,3*2*[q1],4*3*[q1]^2,5*4*[q1]^3,  [f2]/2>,
  20.   <0,0,2,3*2*[q2],4*3*[q2]^2,5*4*[q2]^3,  [g2]/2>
  21.   >
  22. ):
  23.  
  24. let(r,syslinsolve([d])):
  25.  
  26. print_let(r,smooth(<for(<k,1,len([r]),1>,{elem(elem([r],[k]),len([r])+1),})>)) ,\n:
  27.  
  28. func(f6,{for(<k,1,6,1>,{elem([r],[k])*[x]^([k]-1)})}):
  29.  
  30. print_plot(f6,[q1]-1,[q2]+1,.1,60,40,1)
  31.  
Title: Re: A Two-Roads Problem
Post by: SMcNeill on February 07, 2020, 01:03:15 am
Syntax error on Line 1.  Your QB64 code is broken.
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on February 07, 2020, 01:23:36 am
Alright, let me make it fair for the QB64-minded.

I did this in QB64 by the following steps.

1) Create a sub-language strong enough to do math. That project is here:

http://barnes.x10host.com/sxript/index.html (http://barnes.x10host.com/sxript/index.html)

2) Compile Sxript.exe from Sxript.bas as I have done, and then drag+drop the NON QB64 SOURCE CODE, BUT CODE FOR A QB64 PROGRAM NONETHELESS onto Sxript.exe:

Quote
include(`../test/linalg.txt'):

let(q1,-2):
let(f0,-2):
let(f1,1):
let(f2,-1/2):

let(q2,1):
let(g0,5/4):
let(g1,1/2):
let(g2,1/2):

let(d,
  <
  <1,[q1],[q1]^2,[q1]^3,[q1]^4,[q1]^5,    [f0]>,
  <1,[q2],[q2]^2,[q2]^3,[q2]^4,[q2]^5,    [g0]>,
  <0,1,2*[q1],3*[q1]^2,4*[q1]^3,5*[q1]^4, [f1]>,
  <0,1,2*[q2],3*[q2]^2,4*[q2]^3,5*[q2]^4, [g1]>,
  <0,0,2,3*2*[q1],4*3*[q1]^2,5*4*[q1]^3,  [f2]/2>,
  <0,0,2,3*2*[q2],4*3*[q2]^2,5*4*[q2]^3,  [g2]/2>
  >
):

let(r,syslinsolve([d])):

print_let(r,smooth(<for(<k,1,len([r]),1>,{elem(elem([r],[k]),len([r])+1),})>)) ,\n:

func(f6,{for(<k,1,6,1>,{elem([r],[k])*
  • ^([k]-1)})}):


print_plot(f6,[q1]-1,[q2]+1,.1,60,40,1)


...to generate the screenshot below. Contains the answer and an ASCII plot. Have a nice morning.

Edit: I can actually see the quote box did some autoformatting that garbled my function, but for the sake of not triggering anyone again, I'll leave it as a quote and you can scroll back up to the function.

Title: Re: A Two-Roads Problem
Post by: _vince on February 07, 2020, 04:06:27 pm
Code: [Select]
defsng a-z

const sw = 800
const sh = 600
declare SUB ccircle (x1, y1, r, col)
declare SUB cline (x1, y1, x2, y2, col)
declare FUNCTION yc (x)
declare FUNCTION yb (x)
declare FUNCTION ya (x)
   
SCREEN 12
 
CALL cline(0, 0, sw / 2, 0, 7)
CALL cline(0, 0, -sw / 2, 0, 7)
CALL cline(0, 0, 0, sh / 2, 7)
CALL cline(0, 0, 0, -sh / 2, 7)
 
zoom = 75



CALL ccircle(-2 * zoom, ya(-2) * zoom, 10, 14)
CALL ccircle(1 * zoom, yb(1) * zoom, 10, 15)
FOR x = -2.5 TO 2 STEP .01
    CALL ccircle(x * zoom, ya(x) * zoom, 1, 14)
    CALL ccircle(x * zoom, yb(x) * zoom, 1, 15)
    'CALL ccircle(x * zoom, yc(x) * zoom, 1, 5)
   
    a = -2/27
    b = -7/36
    c = 10/9
    d = 11/27
    y = a*x*x*x + b*x*x + c*x + d
    pset (x*zoom + sw/2, sh/2 - y*zoom)
NEXT
 
sleep
system
END
 
FUNCTION ya (x)
    ya = -x ^ 2 / 4 - 1
END FUNCTION
 
FUNCTION yb (x)
    yb = x ^ 2 / 4 + 1
END FUNCTION
 
FUNCTION yc (x)
    yc = x ' What should this be in order to smoothly join each curve?
END FUNCTION
 
SUB cline (x1, y1, x2, y2, col)
    LINE (sw / 2 + x1, -y1 + sh / 2)-(sw / 2 + x2, -y2 + sh / 2), col
END SUB
 
SUB ccircle (x1, y1, r, col)
    CIRCLE (sw / 2 + x1, -y1 + sh / 2), r, col
END SUB
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on February 07, 2020, 08:30:47 pm
Yes! This ^ is excellent. Best answer goes to _vince. If you aren't satisfied by a lack of explanation of why this answer is so good, see Two Roads.pdf above and imagine chopping the series at the O(3) term instead of O(5).
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on February 26, 2020, 09:28:17 pm
So the boys were talking about MS paint in Discord and I got thinking about the one feature in PAINT that is a little mathy - and that's the auto-curve-drawing button. Thinking about how it acted, I had kindof a Eureka moment and realize it solves a spline calculation, which is exactly the solution to this problem. Check out the screenshot. I made the bold line using 3 clicks in paint, which solves the problem graphically. Try it yourself and you can kindof see what the math is doing.
Title: Re: A Two-Roads Problem
Post by: _vince on February 27, 2020, 07:59:19 pm
I think the MS paint is a simple cubic bezier curve.

Here is an interactive n-order bezier curve based on the iterative definition:  https://en.wikipedia.org/wiki/B%C3%A9zier_curve#Explicit_definition (https://en.wikipedia.org/wiki/B%C3%A9zier_curve#Explicit_definition). You can change "STEP 0.01" on line 70ish to something smaller for a finer curve

Code: QB64: [Select]
  1. deflng a-z
  2. const sw = 800
  3. const sh = 600
  4. pi = 4*atn(1)
  5. declare sub circlef(x, y, r, c)
  6.  
  7. n = 8
  8. dim x(n), y(n)
  9.  
  10. x(0) = 0: y(0) = -200
  11. x(1) = -130: y(1) = -130
  12. x(2) = -200: y(2) = 0
  13. x(3) = -30: y(3) = 30
  14. x(4) = 0: y(4) = 200
  15. x(5) = 230: y(5) = 230
  16. x(6) = 200: y(6) = 0
  17. x(7) = 130: y(7) = -30
  18. 'x(n) = x(0)
  19. 'y(n) = y(0)
  20.  
  21. 'for i=0 to n-1
  22. '       x(i) = 200*cos(2*pi*i/n)
  23. '       y(i) = 200*sin(2*pi*i/n)
  24. 'next
  25.  
  26. 'screenres sw, sh, 32
  27. screen _newimage(sw, sh, 32)
  28. redraw = -1
  29. drag = 0
  30.         'm = getmouse(mx, my, mw, mb)
  31.         do
  32.                 mx = _mousex
  33.                 my = _mousey
  34.                 mb = -_mousebutton(1)
  35.         loop while _mouseinput
  36.        
  37.         if drag then
  38.                 if mb = 1 then
  39.                         x(ii) = mx - sw/2
  40.                         y(ii) = sh/2 - my
  41.                 else
  42.                         drag = 0
  43.                 end if
  44.         else
  45.                 dmin = 1e10
  46.                 for i=0 to n - 1
  47.                         dx = (mx - sw/2) - x(i)
  48.                         dy = (sh/2 - my) - y(i)
  49.                         d = (dx*dx + dy*dy)
  50.                         if d < dmin then
  51.                                 dmin = d
  52.                                 ii = i
  53.                         end if
  54.                 next
  55.  
  56.                 if mb = 1 then
  57.                         omx = mx
  58.                         omy = my
  59.                         drag = -1
  60.                 end if
  61.         end if
  62.  
  63.  
  64.         redraw = -1
  65.  
  66.         if redraw then
  67.                 'screenlock
  68.                 line (0,0)-(sw,sh),_rgb(0,0,0),bf
  69.                 locate 1,1: ? m, mx, my, mw, mb
  70.  
  71.  
  72.  
  73.                 dim t as double, bx as double, by as double
  74.                 pset (sw/2 + x(0), sh/2 - y(0)), _rgb(0,255,0)
  75.                 for t=0 to 1 step 0.01
  76.                         bx = 0
  77.                         by = 0
  78.  
  79.                         for i=0 to n-1
  80.                                 bin = 1
  81.                                 for j=1 to i
  82.                                         bin = bin*(n - j)/j
  83.                                 next
  84.  
  85.                                 bx = bx + bin*((1 - t)^(n - 1 - i))*(t^i)*x(i)
  86.                                 by = by + bin*((1 - t)^(n - 1 - i))*(t^i)*y(i)
  87.                         next
  88.  
  89.                         line -(sw/2 + bx, sh/2 - by), _rgb(0,255,0)
  90.                 next
  91.  
  92.  
  93.  
  94.  
  95.                 'x(n) = x(0)
  96.                 'y(n) = y(0)
  97.                 x0 = sw/2 + x(0)
  98.                 y0 = sh/2 - y(0)
  99.                 circlef x0, y0, 8, _rgb(55,55,55)
  100.                 for i=1 to n-1
  101.                         x1 = sw/2 + x(i)
  102.                         y1 = sh/2 - y(i)
  103.                         line (x0, y0)-(x1, y1), _rgb(55,55,55)
  104.                         'line -(sw/2 + cx, sh/2 - cy)
  105.                         circlef x1, y1, 8, _rgb(55,55,55)
  106.                         x0 = x1
  107.                         y0 = y1
  108.                 next
  109.                 circlef sw/2 + x(ii), sh/2 - y(ii), 8, _rgb(255,0,0)
  110.  
  111.                 'screenunlock
  112.                 'screensync
  113.  
  114.                 redraw = 0
  115.  
  116.                 _display
  117.         end if
  118.  
  119.  
  120. sub circlef(x, y, r, c)
  121.         x0 = r
  122.         y0 = 0
  123.         e = -r
  124.  
  125.         do while y0 < x0
  126.                 if e <=0 then
  127.                         y0 = y0 + 1
  128.                         line (x - x0, y + y0)-(x + x0, y + y0), c, bf
  129.                         line (x - x0, y - y0)-(x + x0, y - y0), c, bf
  130.                         e = e + 2*y0
  131.                 else
  132.                         line (x - y0, y - x0)-(x + y0, y - x0), c, bf
  133.                         line (x - y0, y + x0)-(x + y0, y + x0), c, bf
  134.                         x0 = x0 - 1
  135.                         e = e - 2*x0
  136.                 end if
  137.         loop
  138.         line (x - r, y)-(x + r, y), c, bf
  139.  

I might program in the ability to interactively add more points, I'll edit this post if so
Title: Re: A Two-Roads Problem
Post by: FellippeHeitor on February 27, 2020, 08:07:18 pm
Here is an interactive n-order bezier curve based on the iterative definition:

Wow, _vince.
Title: Re: A Two-Roads Problem
Post by: _vince on February 27, 2020, 11:02:44 pm
Thanks, Fil

Nice thinking Steve, but theres a problem with that proposed curve in that it isnt smooth. In a car driving analogy, note that the steering wheel would have to jerk *instantly* when you transition from the line to the circle. A real road/car would forbid an instant jump like that, so you need to enter/exit the circle more gradually than a brutal tangent line. This is especially evident for smaller circles but true for all circles. That doesnt mean a part of the solution cant be line-like or circle-like though.

Said another way, you want a curve to join at the points, yes. You also want the slope of the curve to match the slope of the road where it joins, check. Here's the subtle part: the slope of the slope of the curve must also match the slope of the slope of the road at the intersection points. This guarantees a jerk-free transition.
Funny you trying to say derivative without saying derivative and giving away the solution (to who?). The fact that there isn't a closed expression for G means physics math will remain junk math.
Title: Re: A Two-Roads Problem
Post by: STxAxTIC on February 28, 2020, 05:38:02 am
Quote
The fact that there isn't a closed expression for G means physics math will remain junk math.

Wait what?