Project 2: 12/10

Sushovan De

Intersection of a line and a bezier curve

Please evaluate this entire notebook by running the menu command Evaluation > Evaluate Notebook so that all the figures display properly.

due: 9/28

Use subdivision to find the intersection(s) between a Bezier curve and a straight line. As an example, take the curve with control points

b={{0,0,0},{1,0,0},{1,0,1},{1,1,1},{1,1,0},{0,1,0}};

Find 3D straight lines intersecting this curve 0, 1, 2 times. 

Use tolerances of a) 10E -1 and b) 10E -3

Your output: six plots of successive minmax boxes, with the result(s) highlighted. Also plot the straight lines.

Solution

To start off, we define the De Casteljau' s algorithm, and plot the Bezier curve to make sure everything is in order

pts = {{0,0,0},{1,0,0},{1,0,1},{1,1,1},{1,1,0},{0,1,0}};dca[points_,r_,i_,t_]:=If [r==0,points[[i+1]],(1-t)*dca[points,r-1,i,t]+t*dca[points,r-1,i+1,t]];
n = Length[pts]-1;

curve=ParametricPlot3D[dca[pts,n,0,t],{t,0,1},PlotStyle->{Thick, Green}];
poly=Graphics3D[{Thick,Blue,Line[pts]}];
points=Graphics3D[{PointSize[Large],Red,Point[pts]}];
Show[{poly,curve, points}]

The following function, FindMinMaxBox, takes a control polygon and returns the end-points of the diagonal of the the min-max box of the polygon

FindMinMaxBox[ctrlpoly_] := 
Module[{temp}, temp = Transpose[ctrlpoly]; 
{{Min[temp[[1]]], Min[temp[[2]]], Min[temp[[3]]]}, {Max[temp[[1]]], Max[temp[[2]]], Max[temp[[3]]]}}]

The following function, Interpolate, divides the values y1 and y2 in the same ratio that xp divides x1 and x2

Interpolate[y1_, y2_, x1_, x2_, xp_]:=
y1+ (((xp-x1))/((x2-x1)))(y2-y1)

The following function, IsIntersectingX, takes the following arguments
plane: the endpoints of the diagonal of a rectangle that represents the plane. The direction of the plane is specified by the x parameter
line: the endpoints of the line segment to test intersection with.
x: this is the direction in which the normal to the plane is pointing. 1 specifies x-direction, 2 specifies y-direction, 3 specifies z-direction
y, z: specifies the other directions
and returns True if the line intersects the plane, and False otherwise.

IsIntersectingX[plane_, line_, x_, y_, z_]:=
Module[{intpty, intptz},If[ (line[[1]][[x]] - plane[[1]][[x]])* (line[[2]][[x]] - plane[[1]][[x]])<0,
(intpty=Interpolate[line[[1]][[y]], line[[2]][[y]], line[[1]][[x]], line[[2]][[x]], plane[[1]][[x]]];
intptz=Interpolate[line[[1]][[z]], line[[2]][[z]], line[[1]][[x]], line[[2]][[x]], plane[[1]][[x]]];
If[((intpty - plane[[1]][[y]])*(intpty - plane[[2]][[y]])<0) && 
((intptz - plane[[1]][[z]]) * (intptz - plane[[2]][[z]])<0), True, False]),
False]
];

MyPlane1= {{0, 1, -1}, {0, -1, 1}};
MyLine1 = {{-1, 1, -1}, {1, -1, 1000}};
IsIntersectingX[MyPlane1, MyLine1, 1, 2, 3]

The following function, IsIntersectingBox, takes the following arguments
box: the endpoints of the diagonal of the box with which to test the line segment
line: the endpoints of the line segment to test intersection with.

It calls IsIntersectingX function for each of the six sides of the box, and returns immediately if it finds a hit.

IsIntersectingBox[box_, line_]:=
Module[{reply, planept1, planept2},

(* check the back plane *)
planept1=box[[1]];
planept2=box[[2]]; planept2[[3]]=box[[1]][[3]];
reply = IsIntersectingX[{planept1, planept2}, line, 3, 1, 2];
If[reply, Return[True]];

(*check the left plane*)
planept1=box[[1]];
planept2 = box[[2]]; planept2[[1]]=box[[1]][[1]];
reply = IsIntersectingX[{planept1, planept2}, line, 1, 2, 3];
If[reply, Return[True]];

(*check the bottom plane*)
planept1=box[[1]];
planept2 = box[[2]]; planept2[[2]]=box[[1]][[2]];
reply = IsIntersectingX[{planept1, planept2}, line, 2, 3, 1];
If[reply, Return[True]];

(*check the top plane*)
planept2 = box[[2]]; 
planept1=box[[1]]; planept1[[2]] = box[[2]][[2]];
reply = IsIntersectingX[{planept1, planept2}, line, 2, 3, 1];
If[reply, Return[True]];

(*check the right plane*)
planept2 = box[[2]]; 
planept1=box[[1]]; planept1[[1]] = box[[2]][[1]];
reply = IsIntersectingX[{planept1, planept2}, line, 1, 2, 3];
If[reply, Return[True]];

(*check the front plane*)
planept2 = box[[2]]; 
planept1=box[[1]]; planept1[[3]] = box[[2]][[3]];
reply = IsIntersectingX[{planept1, planept2}, line, 3, 1, 2];
If[reply, Return[True]];

Return[False];
];

The following block contains test cases to check if the IsIntersectingBox function is working correctly.

MyBox1 = {{-1,-1,-1}, {1, 1, 1}};
MyLine1= {{-10.5, 0, 0}, {-9.5, 0, 0}};
IsIntersectingBox[MyBox1, MyLine1] (* False *)
MyLine1= {{-1.5, 0, 0}, {-0.5, 0, 0}};
IsIntersectingBox[MyBox1, MyLine1] (* True *)
MyLine1= {{0, 0.5, 0}, {0, 1.5, 0}};
IsIntersectingBox[MyBox1, MyLine1] (* True *)
MyLine1= {{-1.5, -1.5, 0}, {-0.5, 0.5, 0}};
IsIntersectingBox[MyBox1, MyLine1] (* True *)
MyLine1= {{-1.5, -1.5, 0}, {-2.5, -3.5, 0}};
IsIntersectingBox[MyBox1, MyLine1] (* False *)
MyLine1={{0.8, -0.2, 1.2}, {-0.5, 0.5, -0.5}};
IsIntersectingBox[MyBox1, MyLine1] (* True *)
grline1 = Graphics3D[Line[MyLine1]];
grcube = Graphics3D[{EdgeForm[{Thick,Blue}],FaceForm[{Pink,Opacity[0.2]}],Cuboid[MyBox1]}];
Show[{grline1, grcube}, Axes-> True, AxesLabel -> {x,y, z}]

Subdivide takes a control polygon and returns two control polygons by subdividing them at parameter value t.

Subdivide[ctrlpoly_, t_] :=
{Table[dca[ctrlpoly, r, 0, t], {r, 0, n}],
Table[dca[ctrlpoly, r, n-r, t], {r, n, 0, -1}]};

Let's check if Subdivide is working properly. The previous curve has been divided into a red curve and a green curve by subdividing them at t = 0.4.

tt=0.4;curve=ParametricPlot3D[dca[(Subdivide[pts, tt])[[1]],n,0,t],{t,0,1},PlotStyle->{Thick, Green}];
curve2=ParametricPlot3D[dca[(Subdivide[pts, tt])[[2]],n,0,t],{t,0,1},PlotStyle->{Thick, Red}];
poly=Graphics3D[{Thick,Blue,Line[pts]}];
points=Graphics3D[{PointSize[Large],Red,Point[pts]}];
Show[{poly,curve, curve2, points}]

The FindIntersectionPointsRec function is the recursive function that repeatedly refines the min-max box of the bezier till it is smaller than the threshold.

Arguments:
ctrlpoly: the current control polygon
line: the line to test intersection with
tolerance: the maximum length of the diagonal of the minmax box at which to stop recursion
depth: the current depth of recursion

Returns:
{True, Point of intersection} if a point of intersection was found
{False, {}} if a point was not found

Sows:
In addition to returning the point of intersection, if a point of intersection was found, this function also Sows the current bounding box with the current recursion depth, which can be later read back using the Reap function.

Set the variable debug to True to print messages about current box size, function depth etc.

debug=False;

FindIntersectionPointsRec[ctrlpoly_, line_, tolerance_, depth_]:=
Module[{boundbox, leftpoly, rightpoly, leftsuccess, leftintpt, rightsuccess, rightintpt, reply},
(*Print["FindIntersectionPointsRec ", ctrlpoly, line, tolerance];*)
boundbox = FindMinMaxBox[ctrlpoly];
If[debug, Print["[", depth, "] SizeOfBox -> ", Norm[boundbox[[1]]-boundbox[[2]]]]];
If[debug, Print["[",depth, "] Is Intersecting Box -> ", IsIntersectingBox[boundbox, line]]];
If[!IsIntersectingBox[boundbox, line], Return[{False, {}}]];
If[Norm[boundbox[[1]]-boundbox[[2]]]<tolerance, 
Sow[{depth, boundbox}]; 
Return[{True, {(((boundbox[[1]]+boundbox[[2]]))/2)}}]];
{leftpoly, rightpoly} = Subdivide[ctrlpoly, 0.5];
If[debug, Print["[", depth,"] Going Left"]];
{leftsuccess, leftintpt}=FindIntersectionPointsRec[leftpoly, line, tolerance, depth+1];
If[debug,Print["[", depth, "] Going Right"]];
{rightsuccess, rightintpt}=FindIntersectionPointsRec[rightpoly, line, tolerance, depth+1];
reply = List[];
If[leftsuccess,If[debug, Print["[", depth, "] LeftIntPt: ",leftintpt]];Sow[{depth, boundbox}];reply=Join[reply, leftintpt]]; 
If[rightsuccess,If[debug,Print["[", depth, "] RighttIntPt: ",rightintpt]];Sow[{depth, boundbox}];reply=Join[reply, rightintpt]];
Return[{leftsuccess||rightsuccess, reply}];
];

The FindIntersectionPoints function is master function that calls the recursive worker function with the correct depth.

Arguments:
ctrlpoly: the control polygon
line: the line to test intersection with
tolerance: the maximum length of the diagonal of the minmax box at which to stop recursion

FindIntersectionPoints[ctrlpoly_, line_, tolerance_]:=
Module[{a, b},
{a, b}=FindIntersectionPointsRec[ctrlpoly, line, tolerance, 1];
If[a, Return[b], Return[List[]]];
];

The para of code below tests the case where there are two intersections points (one at t=0.2, the other at t=0.7) with tolerance of 0.1

ptoncurve1=dca[pts, n, 0, 0.2];
ptoncurve2=dca[pts, n, 0, 0.7];
GetPointOnLine[pt1_, pt2_, t_]:=(1-t)pt1 + t pt2;
myline2 = {GetPointOnLine[ptoncurve1, ptoncurve2, -0.8], GetPointOnLine[ptoncurve1, ptoncurve2, 1.8]};
temp=Reap[FindIntersectionPoints[pts,myline2, .1]];
ptOfIntersection = temp[[1]];
intermediateboxes = If[Length[temp[[2]]]>0, temp[[2]][[1]], {}];curve=ParametricPlot3D[dca[pts,n,0,t],{t,0,1},PlotStyle->{Thick, Green}];
poly=Graphics3D[{Thin,Blue,Line[pts]}];
grline1=Graphics3D[{Thick, Blue, Line[myline2]}];
grcube = Graphics3D[{EdgeForm[{Thin,Blue}],FaceForm[{Pink,Opacity[0.2]}],Cuboid[FindMinMaxBox[pts]]}];
grsmallcube=Graphics3D[Table[Cuboid[{i-0.02, i+0.02}],{i, ptOfIntersection}]];
Manipulate[Show[{poly,curve, grline1,Graphics3D[{EdgeForm[{Thin,Red}],FaceForm[{Opacity[0.1]}],Table[Cuboid[i[[2]]],{i,Select[Sort[intermediateboxes], #[[1]]==depth&]}]}],   grsmallcube}, Axes->True, AxesLabel->{x, y, z}], {depth, 1, 6, 1}]

The Figure above shows the intersection point of a line which intersects the bezier curve at two points. Move the slider to view the bounding boxes at various depths from 1 to 6. Since a dot was nearly invisible, I represented the point of intersection with tiny cuboids, which is more visible.

ptoncurve3=dca[pts, n, 0, 0.2];
ptoncurve4={0.8, 1, 0.8};
GetPointOnLine[pt1_, pt2_, t_]:=(1-t)pt1 + t pt2;
myline3 = {GetPointOnLine[ptoncurve3, ptoncurve4, -0.8], GetPointOnLine[ptoncurve3, ptoncurve4, 1.4]};
temp2=Reap[FindIntersectionPoints[pts,myline3, .1]];
ptOfIntersection2 = temp2[[1]];
intermediateboxes2 = If[Length[temp2[[2]]]>0, temp2[[2]][[1]], {}];
grline2=Graphics3D[{Thick, Blue, Line[myline3]}];
grsmallcube2=Graphics3D[Table[Cuboid[{i-0.02, i+0.02}],{i, ptOfIntersection2}]];
Manipulate[Show[{poly,curve, grline2,Graphics3D[{EdgeForm[{Thin,Red}],FaceForm[{Opacity[0.1]}],Table[Cuboid[i[[2]]],{i,Select[Sort[intermediateboxes2], #[[1]]==depth&]}]}],   grsmallcube2}, Axes->True, AxesLabel->{x, y, z}], {depth, 1, 6, 1}]

The Figure above shows the intersection point of a line which intersects the bezier curve at one point (t=0.2). Move the slider to view the bounding boxes at various depths from 1 to 6. Since a dot was nearly invisible, I represented the point of intersection with a tiny cuboid, which is more visible. The threshold is still 0.1

ptoncurve5=dca[pts, n, 0, 0.2];
ptoncurve6=dca[pts, n, 0, 0.7];
GetPointOnLine[pt1_, pt2_, t_]:=(1-t)pt1 + t pt2;
myline4 = {GetPointOnLine[ptoncurve5, ptoncurve6, -0.8], GetPointOnLine[ptoncurve5, ptoncurve6, 1.4]};
temp3=Reap[FindIntersectionPoints[pts,myline4, .001]];
ptOfIntersection3 = temp3[[1]];
intermediateboxes3 = If[Length[temp3[[2]]]>0, temp3[[2]][[1]], {}];
grline3=Graphics3D[{Thick, Blue, Line[myline4]}];
grsmallcube3=Graphics3D[Table[Cuboid[{i-0.01, i+0.01}],{i, ptOfIntersection3}]];
Manipulate[Show[{poly,curve, grline3,Graphics3D[{EdgeForm[{Thin,Red}],FaceForm[{Opacity[0.1]}],Table[Cuboid[i[[2]]],{i,Select[Sort[intermediateboxes3], #[[1]]==depth&]}]}],   grsmallcube3}, Axes->True, AxesLabel->{x, y, z}], {depth, 1, 13, 1}]

Two point intersection with threshold set to 0.001. The bounding boxes go to level 13 in this case.

ptoncurve7=dca[pts, n, 0, 0.2];
ptoncurve8={0.8, 1.0, 0.8};
GetPointOnLine[pt1_, pt2_, t_]:=(1-t)pt1 + t pt2;
myline5 = {GetPointOnLine[ptoncurve7, ptoncurve8, -0.8], GetPointOnLine[ptoncurve7, ptoncurve8, 1.4]};
temp4=Reap[FindIntersectionPoints[pts,myline5, .001]];
ptOfIntersection4 = temp4[[1]];
intermediateboxes4 = If[Length[temp4[[2]]]>0, temp4[[2]][[1]], {}];
grline4=Graphics3D[{Thick, Blue, Line[myline5]}];
grsmallcube4=Graphics3D[Table[Cuboid[{i-0.01, i+0.01}],{i, ptOfIntersection4}]];
Manipulate[Show[{poly,curve, grline4,Graphics3D[{EdgeForm[{Thin,Red}],FaceForm[{Opacity[0.1]}],Table[Cuboid[i[[2]]],{i,Select[Sort[intermediateboxes4], #[[1]]==depth&]}]}],   grsmallcube4}, Axes->True, AxesLabel->{x, y, z}], {depth, 1, 13, 1}]

One point intersection with threshold set to 0.001. The bounding boxes go to level 13 in this case.

ptoncurve9={1.2, 1.2, 1.2};
ptoncurve10={0.8, 1.0, 0.8};
GetPointOnLine[pt1_, pt2_, t_]:=(1-t)pt1 + t pt2;
myline6 = {GetPointOnLine[ptoncurve9, ptoncurve10, -0.8], GetPointOnLine[ptoncurve9, ptoncurve10, 1.4]};
temp5=Reap[FindIntersectionPoints[pts,myline6, .001]];
ptOfIntersection5 = temp5[[1]];
intermediateboxes5 = If[Length[temp5[[2]]]>0, temp5[[2]][[1]], {}];
grline5=Graphics3D[{Thick, Blue, Line[myline6]}];
grsmallcube5=Graphics3D[Table[Cuboid[{i-0.01, i+0.01}],{i, ptOfIntersection5}]];
Manipulate[Show[{poly,curve, grline5,Graphics3D[{EdgeForm[{Thin,Red}],FaceForm[{Opacity[0.1]}],Table[Cuboid[i[[2]]],{i,Select[Sort[intermediateboxes5], #[[1]]==depth&]}]}],   grsmallcube5}, Axes->True, AxesLabel->{x, y, z}], {depth, 1, 13, 1}]

Zero point intersection with threshold set to 0.001. No bounding boxes or intersection points appear.

temp5