Skip to content
Browse files

Split eval.js into several multiple files

  • Loading branch information...
1 parent 88bba48 commit 1b82cc99b15d40c4f6bf22537fef54f34bd8c4cc @pboyer committed Jan 21, 2014
Showing with 4,235 additions and 4,199 deletions.
  1. +8 −1 Gruntfile.js
  2. +1,207 −1,197 build/verb.js
  3. +1,208 −1,198 build/verbEval.js
  4. +0 −1,803 src/eval/eval.js
  5. +519 −0 src/eval/geom.js
  6. +1,023 −0 src/eval/intersect.js
  7. +270 −0 src/eval/tesselate.js
View
9 Gruntfile.js
@@ -17,8 +17,15 @@ module.exports = function(grunt) {
'src/geom/base/NurbsSurface.js',
'src/geom/*.js',
'src/intersect/*.js',
+ 'src/eval/intersect.js',
+ 'src/eval/tesselate.js',
+ 'src/eval/geom.js',
'src/eval/eval.js' ],
- 'build/<%= pkg.name %>Eval.js': ['src/eval/header.js', 'src/eval/eval.js']
+ 'build/<%= pkg.name %>Eval.js': [ 'src/eval/header.js',
+ 'src/eval/intersect.js',
+ 'src/eval/tesselate.js',
+ 'src/eval/geom.js',
+ 'src/eval/eval.js']
}
}
},
View
2,404 build/verb.js
1,207 additions, 1,197 deletions not shown because the diff is too large. Please use a local Git client to view these changes.
View
2,406 build/verbEval.js
1,208 additions, 1,198 deletions not shown because the diff is too large. Please use a local Git client to view these changes.
View
1,803 src/eval/eval.js
@@ -1,1808 +1,5 @@
-// ###verb.eval
-// This defines verb's core geometry library which is called by the current Engine.
-//
-// ####get_sweep1_surface( profile_knots, profile_degree, profile_control_points, profile_weights, rail_knots, rail_degree, rail_control_points, rail_weights )
-//
-// Generate the control points, weights, and knots of an elliptical arc
-//
-// **params**
-// + *Array*, the center
-// + *Array*, the xaxis
-// + *Array*, orthogonal yaxis
-// + *Number*, xradius of the ellipse arc
-// + *Number*, yradius of the ellipse arc
-// + *Number*, start angle of the ellipse arc, between 0 and 2pi, where 0 points at the xaxis
-// + *Number*, end angle of the arc, between 0 and 2pi, greater than the start angle
-//
-// **returns**
-// + *Object*, an object with the following properties: control_points, weights, knots, degree
-//
-verb.eval.nurbs.get_sweep1_surface = function( profile_knots, profile_degree, profile_control_points, profile_weights, rail_knots, rail_degree, rail_control_points, rail_weights ) {
-
- // for each point on rail, move all of the points
- var homo_rail = verb.eval.nurbs.homogenize_1d( rail_control_points, rail_weights )
- , rail_start = verb.eval.nurbs.rational_curve_point( rail_degree, rail_knots, homo_rail, 0 )
- , span = 1.0 / rail_control_points.length
- , control_points = []
- , weights = [];
-
- for (var i = 0; i < rail_control_points.length; i++ ){
-
- // evaluate the point on the curve, subtracting it from the first point
- var rail_point = verb.eval.nurbs.rational_curve_point( rail_degree, rail_knots, homo_rail, i * span )
- , rail_offset = numeric.sub( rail_point, rail_start )
- , row_control_points = []
- , row_weights = [];
-
- for (var j = 0; j < profile_control_points.length; j++ ){
-
- row_control_points.push( numeric.add(rail_offset, profile_control_points[j] ) );
- row_weights.push( profile_weights[j] * rail_weights[i] );
-
- }
-
- control_points.push( row_control_points);
- weights.push( row_weights );
- }
-
- return {"knots_u": rail_knots,
- "knots_v": profile_knots,
- "control_points": control_points,
- "degree_u": rail_degree,
- "degree_v": profile_degree,
- "weights": weights };
-
-}
-
-//
-// ####get_ellipse_arc( center, xaxis, yaxis, xradius, yradius, start_angle, end_angle )
-//
-// Generate the control points, weights, and knots of an elliptical arc
-//
-// **params**
-// + *Array*, the center
-// + *Array*, the xaxis
-// + *Array*, orthogonal yaxis
-// + *Number*, xradius of the ellipse arc
-// + *Number*, yradius of the ellipse arc
-// + *Number*, start angle of the ellipse arc, between 0 and 2pi, where 0 points at the xaxis
-// + *Number*, end angle of the arc, between 0 and 2pi, greater than the start angle
-//
-// **returns**
-// + *Object*, an object with the following properties: control_points, weights, knots, degree
-//
-
-verb.eval.nurbs.get_ellipse_arc = function( center, xaxis, yaxis, xradius, yradius, start_angle, end_angle ) {
-
- // if the end angle is less than the start angle, do a circle
- if (end_angle < start_angle) end_angle = 2 * Math.PI + start_angle;
-
- var theta = end_angle - start_angle
- , narcs = 0;
-
- // how many arcs?
- if (theta <= Math.PI / 2) {
- narcs = 1;
- } else {
- if (theta <= Math.PI){
- narcs = 2;
- } else if (theta <= 3 * Math.PI / 2){
- narcs = 3;
- } else {
- narcs = 4;
- }
- }
-
- var dtheta = theta / narcs
- , n = 2 * narcs
- , w1 = Math.cos( dtheta / 2)
- , P0 = numeric.add( center, numeric.mul( xradius, Math.cos(start_angle), xaxis), numeric.mul( yradius, Math.sin(start_angle), yaxis ) )
- , T0 = numeric.sub( numeric.mul( Math.cos(start_angle), yaxis ), numeric.mul( Math.sin(start_angle), xaxis) )
- , Pw = verb.eval.nurbs.zeros_1d( narcs * 2 )
- , U = verb.eval.nurbs.zeros_1d( 2 *narcs + 3 )
- , index = 0
- , angle = start_angle
- , W = verb.eval.nurbs.zeros_1d( narcs * 2 );
-
- Pw[0] = P0;
- W[0] = 1;
-
- for (var i = 1; i <= narcs; i++){
-
- angle += dtheta;
- var P2 = numeric.add( center, numeric.mul( xradius, Math.cos(angle), xaxis), numeric.mul( yradius, Math.sin(angle), yaxis ) )
-
- W[index+2] = 1;
- Pw[index+2] = P2;
-
- var T2 = numeric.sub( numeric.mul( Math.cos(angle), yaxis ), numeric.mul( Math.sin(angle), xaxis) )
-
- var params = verb.eval.geom.intersect_rays(P0, numeric.mul( 1 / numeric.norm2(T0), T0), P2, numeric.mul( 1 / numeric.norm2(T2), T2));
- var P1 = numeric.add( P0, numeric.mul(T0, params[0]));
-
- W[index+1] = w1;
- Pw[index+1] = P1;
-
- index += 2;
-
- if (i < narcs){
- P0 = P2;
- T0 = T2;
- }
- }
-
- var j = 2 * narcs + 1;
-
- for (var i = 0; i < 3; i++){
- U[i] = 0.0;
- U[i+j] = 1.0;
- }
-
- switch (narcs){
- case 1: break;
- case 2: U[3] = U[4] = 0.5; break;
- case 3: U[3] = U[4] = 1/3;
- U[5] = U[6] = 2/3; break;
- case 4: U[3] = U[4] = 0.25;
- U[5] = U[6] = 0.5;
- U[7] = U[8] = 0.75; break;
- }
-
- return {knots: U, control_points: Pw, degree: 2, weights: W };
-
-}
-
-//
-// ####get_sphere_surface( center, axis, xaxis, radius )
-//
-// Generate the control points, weights, and knots of a sphere
-//
-// **params**
-// + *Array*, the center of the sphere
-// + *Array*, normalized axis of sphere
-// + *Array*, vector perpendicular to axis of sphere, starting the rotation of the sphere
-// + *Number*, radius of the sphere
-//
-// **returns**
-// + *Object*, an object with the following properties: control_points, weights, knots_u, knots_v, degree_u, degree_v
-//
-
-verb.eval.nurbs.get_sphere_surface = function( center, axis, xaxis, radius ){
-
- var arc = verb.eval.nurbs.get_arc(center, numeric.mul(axis, -1), xaxis, radius, 0, Math.PI );
-
- return verb.eval.nurbs.get_revolved_surface( center, axis, 2 * Math.PI, arc.knots, arc.degree, arc.control_points, arc.weights );
-
-}
-
-
-//
-// ####get_polyline_curve( pts )
-//
-// Generate the control points, weights, and knots of a polyline curve
-//
-// **params**
-// + *Array*, array of points in curve
-//
-// **returns**
-// + *Object*, an object with the following properties: control_points, weights, knots, degree
-//
-
-verb.eval.nurbs.get_polyline_curve = function( pts ){
-
- var num_spans = pts.length - 1
- , span = 1.0 / num_spans
- , knots = [0,0];
-
- for (var i = 1; i < num_spans; i++){
- knots.push(i * span);
- }
-
- knots.push(1);
- knots.push(1);
-
- var weights = [];
-
- for (var i = 0; i < pts.length; i++){
- weights.push(1);
- }
-
- return {
- "knots": knots,
- "control_points": pts.slice(0),
- "degree": 1,
- "weights": weights
- };
-
-}
-
-//
-// ####get_4pt_surface( p1, p2, p3, p4 )
-//
-// Generate the control points, weights, and knots of a surface define by 3 points
-//
-// **params**
-// + *Array*, first point in counter-clockwise form
-// + *Array*, second point in counter-clockwise form
-// + *Array*, third point in counter-clockwise form
-// + *Array*, forth point in counter-clockwise form
-//
-// **returns**
-// + *Object*, an object with the following properties: control_points, weights, knots_u, knots_v, degree_u, degree_v
-//
-
-verb.eval.nurbs.get_4pt_surface = function( p1, p2, p3, p4 ){
-
- return {"knots_u": [0,0,1,1],
- "knots_v": [0,0,1,1],
- "control_points": [ [p1, p4], [p2, p3] ],
- "degree_u": 1,
- "degree_v": 1,
- "weights": [ [1, 1], [1, 1] ] };
-
-}
-
-//
-// ####get_cylinder_surface( axis, xaxis, base, height, radius )
-//
-// Generate the control points, weights, and knots of a cylinder
-//
-// **params**
-// + *Array*, normalized axis of cylinder
-// + *Array*, xaxis in plane of cylinder
-// + *Array*, position of base of cylinder
-// + *Number*, height from base to top
-// + *Number*, radius of the cylinder
-//
-// **returns**
-// + *Object*, an object with the following properties: control_points, weights, knots_u, knots_v, degree_u, degree_v
-//
-
-verb.eval.nurbs.get_cylinder_surface = function( axis, xaxis, base, height, radius ){
-
- var yaxis = crossprod( axis, xaxis )
- , angle = 2 * Math.PI
- , circ = verb.eval.nurbs.get_arc( base, xaxis, yaxis, radius, 0, 2 * Math.PI );
-
- return verb.eval.nurbs.get_extruded_surface( axis, height, circ.knots, circ.degree, circ.control_points, circ.weights );
-
-}
-
-//
-// ####get_cone_surface( axis, xaxis, base, height, radius )
-//
-// Generate the control points, weights, and knots of a cone
-//
-// **params**
-// + *Array*, normalized axis of cone
-// + *Array*, position of base of cone
-// + *Number*, height from base to tip
-// + *Number*, radius at the base of the cone
-//
-// **returns**
-// + *Object*, an object with the following properties: control_points, weights, knots, degree
-//
-
-verb.eval.nurbs.get_cone_surface = function( axis, xaxis, base, height, radius ){
-
- var angle = 2 * Math.PI
- , prof_degree = 1
- , prof_ctrl_pts = [ numeric.add( base, numeric.mul( height, axis ) ), numeric.add( base, numeric.mul( radius, xaxis ) )]
- , prof_knots = [0,0,1,1]
- , prof_weights = [1,1];
-
- return verb.eval.nurbs.get_revolved_surface(base, axis, angle, prof_knots, prof_degree, prof_ctrl_pts, prof_weights);
-
-}
-
-//
-// ####get_extruded_surface( axis, length, prof_knots, prof_degree, prof_control_points, prof_weights)
-//
-// Generate the control points, weights, and knots of an extruded surface
-//
-// **params**
-// + *Array*, axis of the extrusion
-// + *Array*, length of the extrusion
-// + *Number*, degree of the profile
-// + *Number*, control points of the profile
-// + *Number*, weights of the profile
-//
-// **returns**
-// + *Object*, an object with the following properties: control_points, weights, knots, degree
-//
-
-verb.eval.nurbs.get_extruded_surface = function( axis, length, prof_knots, prof_degree, prof_control_points, prof_weights){
-
- var control_points = verb.eval.nurbs.zeros_2d( 2, prof_control_points.length )
- , weights = verb.eval.nurbs.zeros_2d( 2, prof_control_points.length );
-
- // original control points
- for (var j = 0; j < prof_control_points.length; j++){
- control_points[0][j] = prof_control_points[j];
- weights[0][j] = prof_weights[j];
- }
-
- // build translated control points
- var translation = numeric.mul(axis, length);
-
- for (var j = 0; j < prof_control_points.length; j++){
- control_points[1][j] = numeric.add( translation, prof_control_points[j] );
- weights[1][j] = prof_weights[j];
- }
-
- // return all parameters
- return {"knots_u": [0,0,1,1],
- "knots_v": prof_knots,
- "control_points": control_points,
- "degree_u": 1,
- "degree_v": prof_degree,
- "weights": weights };
-}
-
-//
-// ####get_revolved_surface( center, axis, theta, prof_knots, prof_degree, prof_control_points, prof_weights)
-//
-// Generate the control points, weights, and knots of a revolved surface
-// (Corresponds to Algorithm A7.1 from Piegl & Tiller)
-//
-// **params**
-// + *Array*, center of the rotation axis
-// + *Array*, axis of the rotation axis
-// + *Number*, angle to revolve around axis
-// + *Number*, degree of the generatrix
-// + *Number*, control points of the generatrix
-// + *Number*, weights of the generatrix
-//
-// **returns**
-// + *Object*, an object with the following properties: control_points, weights, knots, degree
-//
-
-// helper method
-
-function crossprod(u,v) {
- return [u[1]*v[2]-u[2]*v[1],u[2]*v[0]-u[0]*v[2],u[0]*v[1]-u[1]*v[0]];
-}
-
-verb.eval.nurbs.get_revolved_surface = function( center, axis, theta, prof_knots, prof_degree, prof_control_points, prof_weights){
-
- var narcs, knots_u, control_points, weights;
-
- if (theta <= Math.PI / 2) { // less than 90
- narcs = 1;
- knots_u = verb.eval.nurbs.zeros_1d( 6 + 2 * (narcs-1) );
- } else {
- if (theta <= Math.PI){ // between 90 and 180
- narcs = 2;
- knots_u = verb.eval.nurbs.zeros_1d( 6 + 2 * (narcs-1) );
- knots_u[3]= knots_u[4] = 0.5;
- } else if (theta <= 3 * Math.PI / 2){ // between 180 and 270
- narcs = 3;
- knots_u = verb.eval.nurbs.zeros_1d( 6 + 2 * (narcs-1) );
- knots_u[3]= knots_u[4] = 1/3;
- knots_u[5]= knots_u[6] = 2/3;
- } else { // between 270 and 360
- narcs = 4;
- knots_u = verb.eval.nurbs.zeros_1d( 6 + 2 * (narcs-1) );
- knots_u[3]= knots_u[4] = 1/4;
- knots_u[5]= knots_u[6] = 1/2;
- knots_u[7]= knots_u[8] = 3/4;
- }
- }
-
- var dtheta = theta / narcs // divide the interval into several points
- , j = 3 + 2 * (narcs-1);
-
- // initialize the start and end knots
- // keep in mind that we only return the knot vector for the
-
- for (var i = 0; i < 3; j++, i++){
- knots_u[i] = 0.0;
- knots_u[j] = 1.0;
- }
-
- // do some initialization
- var n = 2 * narcs
- , wm = Math.cos( dtheta/2.0 )
- , angle = 0.0
- , sines = verb.eval.nurbs.zeros_1d( narcs + 1)
- , cosines = verb.eval.nurbs.zeros_1d( narcs + 1)
- , control_points = verb.eval.nurbs.zeros_2d( 2*narcs + 1, prof_control_points.length )
- , weights = verb.eval.nurbs.zeros_2d( 2*narcs + 1, prof_control_points.length );
-
- // initialize the sines and cosines
- for (var i = 1; i <= narcs; i++){
- angle += dtheta;
- cosines[i] = Math.cos(angle);
- sines[i] = Math.sin(angle);
- }
-
- // for each pt in the generatrix
- // i.e. for each row of the 2d knot vectors
- for (j = 0; j < prof_control_points.length; j++){
-
- // get the closest point of the generatrix point on the axis
- var O = verb.eval.geom.closest_point_on_ray(prof_control_points[j], center, axis)
- // X is the vector from the axis to generatrix control pt
- , X = numeric.sub( prof_control_points[j], O )
- // radius at that height
- , r = numeric.norm2(X)
- // Y is perpendicular to X and axis, and complete the coordinate system
- , Y = crossprod(axis,X);
-
- if ( r > verb.EPSILON ){
- X = numeric.mul( 1 / r, X);
- Y = numeric.mul( 1 / r, Y);
- }
-
- // the first row of control_points and weights is just the generatrix
- control_points[0][j] = prof_control_points[j];
- var P0 = prof_control_points[j];
- weights[0][j] = prof_weights[j];
-
- // store T0 as the Y vector
- var T0 = Y
- , index = 0
- , angle = 0.0;
-
- // proceed around the circle
- for (var i = 1; i <= narcs; i++){
-
- // O + r * cos(theta) * X + r * sin(theta) * Y
- // rotated generatrix pt
- var P2 = r == 0 ? O : numeric.add( O, numeric.mul( r, cosines[i], X), numeric.mul( r, sines[i], Y) );
-
- control_points[index+2][j] = P2;
- weights[index+2][j] = prof_weights[j];
-
- // construct the vector tangent to the rotation
- var T2 = numeric.sub( numeric.mul( cosines[i], Y), numeric.mul(sines[i], X));
-
- // construct the next control pt
- if (r == 0){
- control_points[index+1][j] = O;
- } else {
- var params = verb.eval.geom.intersect_rays(P0, numeric.mul( 1 / numeric.norm2(T0), T0), P2, numeric.mul( 1 / numeric.norm2(T2), T2));
- var P1 = numeric.add( P0, numeric.mul(T0, params[0]));
-
- control_points[index+1][j] = P1;
- }
-
- weights[index+1][j] = wm * prof_weights[j];
-
- index += 2;
-
- if (i < narcs)
- {
- P0 = P2;
- T0 = T2;
- }
-
- }
- }
-
- // store all of the parameters
- return {"knots_u": knots_u,
- "knots_v": prof_knots,
- "control_points": control_points,
- "degree_u": 2,
- "degree_v": prof_degree,
- "weights": weights };
-
-}
-
-//
-// ####get_arc( center, xaxis, yaxis, radius, start_angle, end_angle )
-//
-// Generate the control points, weights, and knots of an arbitrary arc
-// (Corresponds to Algorithm A7.1 from Piegl & Tiller)
-//
-// **params**
-// + *Array*, the center of the arc
-// + *Array*, the xaxis of the arc
-// + *Array*, orthogonal yaxis of the arc
-// + *Number*, radius of the arc
-// + *Number*, start angle of the arc, between 0 and 2pi
-// + *Number*, end angle of the arc, between 0 and 2pi, greater than the start angle
-//
-// **returns**
-// + *Object*, an object with the following properties: control_points, weights, knots, degree
-//
-
-verb.eval.nurbs.get_arc = function( center, xaxis, yaxis, radius, start_angle, end_angle ) {
-
- return verb.eval.nurbs.get_ellipse_arc( center, xaxis, yaxis, radius, radius, start_angle, end_angle );
-
-}
-
-//
-// ####intersect_meshes_by_aabb( points1, tris1, uvs1, points2, tris2, uvs2 )
-//
-// Intersect two meshes via aabb intersection
-//
-// **params**
-// + *Number*, integer degree of surface in u direction
-// + *Array*, array of nondecreasing knot values in u direction
-// + *Number*, integer degree of surface in v direction
-// + *Array*, array of nondecreasing knot values in v direction
-// + *Array*, 3d array of control points, top to bottom is increasing u direction, left to right is increasing v direction,
-// and where each control point is an array of length (dim+1)
-// + *Number*, u parameter at which to evaluate the surface point
-// + *Number*, v parameter at which to evaluate the surface point
-//
-// **returns**
-// + *Array*, a point represented by an array of length (dim)
-//
-
-verb.eval.mesh.intersect_meshes_by_aabb = function( points1, tris1, uvs1, points2, tris2, uvs2 ) {
-
- // build aabb for each mesh
- var tri_indices1 = _.range(tris1.length)
- , tri_indices2 = _.range(tris2.length)
- , aabb1 = verb.eval.mesh.make_mesh_aabb_tree( points1, tris1, tri_indices1 )
- , aabb2 = verb.eval.mesh.make_mesh_aabb_tree( points2, tris2, tri_indices2 )
-
- // intersect and get the pairs of triangle intersctions
- , intersection_pairs = verb.eval.mesh.intersect_aabb_tree( points1, tris1, points2, tris2, aabb1, aabb2 );
-
- // get the segments of the intersection crv with uvs
-
- // sort the intersection segments
-
- // TODO
-
-}
-
-//
-// ####intersect_tris( points1, tri1, uvs1, points2, tri2, uvs2 )
-//
-// Intersect two triangles
-//
-// **params**
-// + *Array*, array of length 3 arrays of numbers representing the points of mesh1
-// + *Array*, array of length 3 arrays of number representing the triangles of mesh1
-// + *Array*, array of length 3 arrays of numbers representing the points of mesh2
-// + *Array*, array of length 3 arrays of number representing the triangles of mesh2
-//
-// **returns**
-// + *Array*, a point represented by an array of length (dim)
-//
-
-verb.eval.geom.intersect_tris = function( points1, tri1, uvs1, points2, tri2, uvs2 ) {
-
- var seg1a = [ points1[ tr1[0] ], points1[ tr1[1] ] ]
- , seg1b = [ points1[ tr1[1] ], points1[ tr1[2] ] ]
- , seg1c = [ points1[ tr1[2] ], points1[ tr1[0] ] ]
- , seg2a = [ points2[ tr2[0] ], points2[ tr2[1] ] ]
- , seg2b = [ points2[ tr2[1] ], points2[ tr2[2] ] ]
- , seg2c = [ points2[ tr2[2] ], points2[ tr2[0] ] ]
- , segs1 = [ seg1a, seg1b, seg1c ]
- , segs2 = [ seg2a, seg2b, seg2c ]
- , int_results = []
- , tri2norm = verb.eval.geom.get_tri_norm(points2, tri2)
- , pt2 = points2[ tr2[0] ];
-
- for (var i = 0; i < 3; i++){
-
- var result = verb.eval.geom.intersect_segment_with_plane( segs1[i][0], segs2[i][1], pt2, tri2norm );
-
- if ( result.intersects ){
- int_results.push( result );
- }
-
- }
-
- // if you don't have 2 intersections you do not have an intersection,
- // 0 would mean a glancing intersection
- // 3 means we don't have a triangle
- if ( int_results.length !== 2 ){
- return null;
- }
-
- // what portions of the segment lie within triangle 2
-
- // intersect edges of triangle 2 with the segment, obtaining the "inner" triangle
- var seg = [int_results[0].point, int_results[1].point ]
- , seg_int_results = [];
-
- for (var i = 0; i < 3; i++){
- var seg_int_result = verb.eval.geom.intersect_segments( seg[0], seg[1], seg, b1, tol );
- if ( seg_int_result ){
- seg_int_results.push( seg_int_result );
- }
- }
-
- // all intersections should be with uv's
-
- if ( seg_int_results.length === 0 ) {
-
- // tri1 is intersecting and the intersection segment
- // is inside tri2
-
- // TODO
-
- // return the two outer points
-
- } else if ( seg_int_results.length === 1 ) {
-
- // tri1 is intersecting and the intersection segment
- // is partially inside tri2
-
- // TODO
-
- // return the end point of seg that is INSIDE tri2 and the intersection
-
- } else if ( seg_int_results.length === 2 ) {
-
- // tri1 is intersecting and the intersection segment's
- // end points are outside of tri2
-
- // TODO
-
- // return the two seg_int_results
-
- }
-
-}
-
-//
-// ####intersect_rational_curve_surface_by_aabb( degree_u, knots_u, degree_v, knots_v, homo_control_points, degree_crv, knots_crv, homo_control_points_crv, sample_tol, tol )
-//
-// Get the intersection of a NURBS curve and a NURBS surface by axis-aligned bounding box intersection and refinement
-//
-// **params**
-// + *Number*, integer degree of surface in u direction
-// + *Array*, array of nondecreasing knot values in u direction
-// + *Number*, integer degree of surface in v direction
-// + *Array*, array of nondecreasing knot values in v direction
-// + *Array*, 3d array of homogeneous control points, top to bottom is increasing u direction, left to right is increasing v direction,
-// and where each control point is an array of length (dim+1)
-// + *Number*, integer degree of curve
-// + *Array*, array of nondecreasing knot values
-// + *Array*, 2d array of homogeneous control points, where each control point is an array of length (dim+1)
-// and form (wi*pi, wi)
-// + *Number*, the sample tolerance of the curve
-// + *Number*, tolerance for the curve intersection
-// + *Number*, integer number of divisions of the surface in the U direction for initial approximation (placeholder until adaptive tess of surfaces)
-// + *Number*, integer number of divisions of the surface in the V direction for initial approximation (placeholder until adaptive tess of surfaces)
-//
-// **returns**
-// + *Array*, array of intersection objects, each holding:
-// - a "point" property where intersections took place
-// - a "p" the parameter on the curve
-// - a "uv" the parameter on the surface
-// - a "face" the index of the face where the intersection took place
-//
-
-verb.eval.nurbs.intersect_rational_curve_surface_by_aabb_refine = function( degree_u, knots_u, degree_v, knots_v, homo_control_points_srf, degree_crv, knots_crv, homo_control_points_crv, sample_tol, tol, divs_u, divs_v ) {
-
- // get the approximate intersections
- var ints = verb.eval.nurbs.intersect_rational_curve_surface_by_aabb( degree_u, knots_u, degree_v, knots_v, homo_control_points_srf, degree_crv, knots_crv, homo_control_points_crv, sample_tol, tol, divs_u, divs_v );
-
- // refine them
- return ints.map(function( inter ){
-
- // get intersection params
- var start_params = [inter.p, inter.uv[0], inter.uv[1] ]
-
- // refine the parameters
- , refined_params = verb.eval.nurbs.refine_rational_curve_surface_intersection( degree_u, knots_u, degree_v, knots_v, homo_control_points_srf, degree_crv, knots_crv, homo_control_points_crv, start_params );
-
- // update the inter object
- inter.p = refined_params[0];
- inter.uv[0] = refined_params[1];
- inter.uv[1] = refined_params[2];
- inter.distance = refined_params[3];
- delete inter.face;
-
- return inter;
-
- });
-
-}
-
-//
-// ####refine_rational_curve_surface_intersection( degree_u, knots_u, degree_v, knots_v, homo_control_points_srf, degree_crv, knots_crv, homo_control_points_crv, start_params )
-//
-// Refine an intersection pair for a surface and curve given an initial guess. This is an unconstrained minimization,
-// so the caller is responsible for providing a very good initial guess.
-//
-// **params**
-// + *Number*, integer degree of surface in u direction
-// + *Array*, array of nondecreasing knot values in u direction
-// + *Number*, integer degree of surface in v direction
-// + *Array*, array of nondecreasing knot values in v direction
-// + *Array*, 3d array of homogeneous control points, top to bottom is increasing u direction, left to right is increasing v direction,
-// and where each control point is an array of length (dim+1)
-// + *Number*, integer degree of curve
-// + *Array*, array of nondecreasing knot values
-// + *Array*, 2d array of homogeneous control points, where each control point is an array of length (dim+1)
-// and form (wi*pi, wi)
-// + *Array*, array of initial parameter values [ u_crv, u_srf, v_srf ]
-//
-// **returns**
-// + *Array*, a length 3 array containing the [ u_crv, u_srf, v_srf, final_distance ]
-//
-
-verb.eval.nurbs.refine_rational_curve_surface_intersection = function( degree_u, knots_u, degree_v, knots_v, homo_control_points_srf, degree_crv, knots_crv, homo_control_points_crv, start_params ) {
-
- var objective = function(x) {
-
- var p1 = verb.eval.nurbs.rational_curve_point(degree_crv, knots_crv, homo_control_points_crv, x[0])
- , p2 = verb.eval.nurbs.rational_surface_point( degree_u, knots_u, degree_v, knots_v, homo_control_points_srf, x[1], x[2] )
- , p1_p2 = numeric.sub(p1, p2);
-
- return numeric.dot(p1_p2, p1_p2);
- }
-
- var sol_obj = numeric.uncmin( objective, start_params);
- return sol_obj.solution.concat( sol_obj.f );
-
-}
-
-
-
-//
-// ####intersect_rational_curve_surface_by_aabb( degree_u, knots_u, degree_v, knots_v, homo_control_points, degree_crv, knots_crv, homo_control_points_crv, sample_tol, tol )
-//
-// Approximate the intersection of two nurbs surface by axis-aligned bounding box intersection.
-//
-// **params**
-// + *Number*, integer degree of surface in u direction
-// + *Array*, array of nondecreasing knot values in u direction
-// + *Number*, integer degree of surface in v direction
-// + *Array*, array of nondecreasing knot values in v direction
-// + *Array*, 3d array of homogeneous control points, top to bottom is increasing u direction, left to right is increasing v direction,
-// and where each control point is an array of length (dim+1)
-// + *Number*, integer degree of curve
-// + *Array*, array of nondecreasing knot values
-// + *Array*, 2d array of homogeneous control points, where each control point is an array of length (dim+1)
-// and form (wi*pi, wi)
-// + *Array*, array of initial parameter values [ u_crv, u_srf, v_srf ]
-// + *Number*, the sample tolerance of the curve
-// + *Number*, tolerance for the curve intersection
-// + *Number*, integer number of divisions of the surface in the U direction
-// + *Number*, integer number of divisions of the surface in the V direction
-//
-// **returns**
-// + *Array*, array of intersection objects, each holding:
-// - a "point" property where intersections took place
-// - a "p" the parameter on the polyline
-// - a "uv" the parameter on the mesh
-// - a "face" the index of the face where the intersection took place
-//
-
-verb.eval.nurbs.intersect_rational_curve_surface_by_aabb = function( degree_u, knots_u, degree_v, knots_v, homo_control_points_srf, degree_crv, knots_crv, homo_control_points_crv, sample_tol, tol, divs_u, divs_v ) {
-
- // tesselate the curve
- var crv = verb.eval.nurbs.rational_curve_adaptive_sample( degree_crv, knots_crv, homo_control_points_crv, sample_tol, true)
-
- // tesselate the surface
- , mesh = verb.eval.nurbs.tesselate_rational_surface_naive( degree_u, knots_u, degree_v, knots_v, homo_control_points_srf, divs_u, divs_v )
-
- // separate parameters from points in the polyline (params are the first index in the array)
- , u1 = crv.map( function(el) { return el[0]; })
- , p1 = crv.map( function(el) { return el.slice(1) })
-
- // perform intersection
- , res = verb.eval.nurbs.intersect_parametric_polyline_mesh_by_aabb( p1, u1, mesh, _.range(mesh.faces.length), tol );
-
- // eliminate duplicate intersections
- return verb.unique( res, function(a, b){
- return numeric.norm2( numeric.sub( a.point, b.point ) ) < tol && Math.abs( a.p - b.p ) < tol && numeric.norm2( numeric.sub( a.uv, b.uv ) ) < tol
- });
-
-}
-
-//
-// ####intersect_parametric_polyline_mesh_by_aabb( crv_points, crv_param_points, mesh, included_faces, tol )
-//
-// Approximate the intersection of a polyline and mesh while maintaining parameter information
-//
-// **params**
-// + *Array*, array of 3d points on the curve
-// + *Array*, array of parameters corresponding to the parameters on the curve
-// + *Object*, a triangular mesh with a "faces" attribute and "points" attribute
-// + *Array*, an array of indices, representing the faces to include in the intersection operation
-// + *Number*, tolerance for the intersection
-//
-// **returns**
-// + *Array*, array of intersection objects (with potential duplicates ) each holding:
-// - a "point" property where intersections took place
-// - a "p" the parameter on the polyline
-// - a "uv" the parameter on the mesh
-// - a "face" the index of the face where the intersection took place
-//
-
-verb.eval.nurbs.intersect_parametric_polyline_mesh_by_aabb = function( crv_points, crv_param_points, mesh, included_faces, tol ) {
-
- // check if two bounding boxes intersect
- var pl_bb = new verb.geom.BoundingBox( crv_points )
- , mesh_bb = verb.eval.mesh.make_mesh_aabb( mesh.points, mesh.faces, included_faces );
-
- // if bounding boxes do not intersect, return empty array
- if ( !pl_bb.intersects( mesh_bb, tol ) ) {
- return [];
- }
-
- if ( crv_points.length === 2 && included_faces.length === 1 ){
-
- // intersect segment and triangle
-
- var inter = verb.eval.geom.intersect_segment_with_tri( crv_points[0], crv_points[1], mesh.points, mesh.faces[ included_faces[0] ] );
-
- if ( inter != null ){
-
- // map the parameters of the segment to the parametric space of the entire polyline
- var p = inter.p * ( crv_param_points[1]-crv_param_points[0] ) + crv_param_points[0];
-
- // map the parameters of the single triangle to the entire parametric space of the triangle
- var index_v0 = mesh.faces[ included_faces ][0]
- , index_v1 = mesh.faces[ included_faces ][1]
- , index_v2 = mesh.faces[ included_faces ][2]
- , uv_v0 = mesh.uvs[ index_v0 ]
- , uv_v1 = mesh.uvs[ index_v1 ]
- , uv_v2 = mesh.uvs[ index_v2 ]
- , uv_s_diff = numeric.sub( uv_v1, uv_v0 )
- , uv_t_diff = numeric.sub( uv_v2, uv_v0 )
- , uv = numeric.add( uv_v0, numeric.mul( inter.s, uv_s_diff ), numeric.mul( inter.t, uv_t_diff ) );
-
- // a pair representing the param on the polyline and the param on the mesh
- return [ { point: inter.point, p: p, uv: uv, face: included_faces[0] } ];
-
- }
-
- } else if ( included_faces.length === 1 ) {
-
- // intersect triangle and polyline
-
- // divide polyline in half, rightside includes the pivot
- var crv_points_a = verb.left( crv_points )
- , crv_points_b = verb.rightWithPivot( crv_points )
- , crv_param_points_a = verb.left( crv_param_points )
- , crv_param_points_b = verb.rightWithPivot( crv_param_points );
-
- return verb.eval.nurbs.intersect_parametric_polyline_mesh_by_aabb( crv_points_a, crv_param_points_a, mesh, included_faces, tol )
- .concat( verb.eval.nurbs.intersect_parametric_polyline_mesh_by_aabb( crv_points_b, crv_param_points_b, mesh, included_faces, tol ) );
-
-
- } else if ( crv_points.length === 2 ) {
-
- // intersect mesh >2 faces and line
-
- // divide mesh in "half" by first sorting then dividing array in half
- var sorted_included_faces = verb.eval.mesh.sort_tris_on_longest_axis( mesh_bb, mesh.points, mesh.faces, included_faces )
- , included_faces_a = verb.left( sorted_included_faces )
- , included_faces_b = verb.right( sorted_included_faces );
-
- return verb.eval.nurbs.intersect_parametric_polyline_mesh_by_aabb( crv_points, crv_param_points, mesh, included_faces_a, tol )
- .concat( verb.eval.nurbs.intersect_parametric_polyline_mesh_by_aabb( crv_points, crv_param_points, mesh, included_faces_b, tol ));
-
-
- } else {
-
- // intersect mesh with >2 faces and polyline
-
- // divide mesh in "half"
- var sorted_included_faces = verb.eval.mesh.sort_tris_on_longest_axis( mesh_bb, mesh.points, mesh.faces, included_faces )
- , included_faces_a = verb.left( sorted_included_faces )
- , included_faces_b = verb.right( sorted_included_faces );
-
- // divide polyline in half, rightside includes the pivot
- var crv_points_a = verb.left( crv_points )
- , crv_points_b = verb.rightWithPivot( crv_points )
- , crv_param_points_a = verb.left( crv_param_points )
- , crv_param_points_b = verb.rightWithPivot( crv_param_points );
-
- return verb.eval.nurbs.intersect_parametric_polyline_mesh_by_aabb( crv_points_a, crv_param_points_a, mesh, included_faces_a, tol )
- .concat( verb.eval.nurbs.intersect_parametric_polyline_mesh_by_aabb( crv_points_a, crv_param_points_a, mesh, included_faces_b, tol ) )
- .concat( verb.eval.nurbs.intersect_parametric_polyline_mesh_by_aabb( crv_points_b, crv_param_points_b, mesh, included_faces_a, tol ) )
- .concat( verb.eval.nurbs.intersect_parametric_polyline_mesh_by_aabb( crv_points_b, crv_param_points_b, mesh, included_faces_b, tol ) );
-
- }
-
- return [];
-
-}
-
-//
-// ####intersect_segment_with_tri( p1, p0, points, tri )
-//
-// Intersect segment with triangle (from http://geomalgorithms.com/a06-_intersect-2.html)
-//
-// **params**
-// + *Array*, array of length 3 representing first point of the segment
-// + *Array*, array of length 3 representing second point of the segment
-// + *Array*, array of length 3 arrays representing the points of the triangle
-// + *Array*, array of length 3 containing int indices in the array of points, this allows passing a full mesh
-//
-// **returns**
-// + *Object*, an object with an "intersects" property that is true or false and if true, a
-// "s" property giving the param on u, and "t" is the property on v, where u is the
-// axis from v0 to v1, and v is v0 to v1, a "point" property
-// where the intersection took place, and "p" property representing the parameter along the segment
-//
-
-verb.eval.geom.intersect_segment_with_tri = function( p0, p1, points, tri ) {
-
- var v0 = points[ tri[0] ]
- , v1 = points[ tri[1] ]
- , v2 = points[ tri[2] ]
- , u = numeric.sub( v1, v0 )
- , v = numeric.sub( v2, v0 )
- , n = numeric.cross( u, v );
-
- var dir = numeric.sub( p1, p0 )
- , w0 = numeric.sub( p0, v0 )
- , a = -numeric.dot( n, w0 )
- , b = numeric.dot( n, dir )
-
- // is ray is parallel to triangle plane?
- if ( Math.abs( b ) < verb.EPSILON ){
- return null;
- }
-
- var r = a / b;
-
- // segment goes away from triangle or is beyond segment
- if ( r < 0 || r > 1 ){
- return null;
- }
-
- // get proposed intersection
- var pt = numeric.add( p0, numeric.mul( r, dir ) );
-
- // is I inside T?
- var uv = numeric.dot(u,v)
- , uu = numeric.dot(u,u)
- , vv = numeric.dot(v,v)
- , w = numeric.sub( pt, v0 )
- , wu = numeric.dot( w, u )
- , wv = numeric.dot( w, v )
- , denom = uv * uv - uu * vv
- , s = ( uv * wv - vv * wu ) / denom
- , t = ( uv * wu - uu * wv ) / denom;
-
- if (s > 1.0 + verb.EPSILON || t > 1.0 + verb.EPSILON || t < -verb.EPSILON || s < -verb.EPSILON || s + t > 1.0 + verb.EPSILON){
- return null;
- }
-
- return { point: pt, s: s, t: t, p: r };
-
-}
-
-//
-// ####intersect_segment_with_plane( p0, p1, v0, n )
-//
-// Intersect ray/segment with plane (from http://geomalgorithms.com/a06-_intersect-2.html)
-//
-// If intersecting a ray, the param needs to be between 0 and 1 and the caller is responsible
-// for making that check
-//
-// **params**
-// + *Array*, array of length 3 representing first point of the segment
-// + *Array*, array of length 3 representing second point of the segment
-// + *Array*, array of length 3 representing an origin point on the plane
-// + *Array*, array of length 3 representing the normal of the plane
-//
-// **returns**
-// null or an object with a p property representing the param on the segment
-//
-
-verb.eval.geom.intersect_segment_with_plane = function( p0, p1, v0, n ) {
-
- var denom = numeric.dot( n, numeric.sub(p0,p1) );
-
- // parallel case
- if ( abs( denom ) < EPSILON ) {
- return null;
- }
-
- var numer = numeric.dot( n, numeric.sub(v0,p0) );
-
- return { p: numer / denom };
-
-}
-
-//
-// ####intersect_aabb_trees( points1, tris1, points2, tris2, aabb_tree1, aabb_tree2 )
-//
-// Intersect two aabb trees - a recursive function
-//
-// **params**
-// + *Array*, array of length 3 arrays of numbers representing the points of mesh1
-// + *Array*, array of length 3 arrays of number representing the triangles of mesh1
-// + *Array*, array of length 3 arrays of numbers representing the points of mesh2
-// + *Array*, array of length 3 arrays of number representing the triangles of mesh2
-// + *Object*, nested object representing the aabb tree of the first mesh
-// + *Object*, nested object representing the aabb tree of the second mesh
-//
-// **returns**
-// + *Array*, a list of pairs of triangle indices for mesh1 and mesh2 that are intersecting
-//
-
-verb.eval.geom.intersect_aabb_trees = function( points1, tris1, points2, tris2, aabb_tree1, aabb_tree2 ) {
-
- var intersects = aabb_tree1.bounding_box.intersects( aabb_tree2.bounding_box );
-
- if (!intersects){
- return [];
- }
-
- if (aabb_tree1.children.length === 0 && aabb_tree2.children.length === 0){
-
- return [ [aabb_tree1.triangle, aabb_tree2.triangle ] ];
-
- } else if (aabb_tree1.children.length === 0 && aabb_tree2.children.length != 0){
-
- return verb.eval.geom.intersect_aabb_trees( points1, tris1, points2, tris2, aabb_tree1, aabb_tree2.children[0] )
- .concat( verb.eval.geom.intersect_aabb_trees( points1, tris1, points2, tris2, aabb_tree1, aabb_tree2.children[1] ) );
-
- } else if (aabb_tree1.children.length != 0 && aabb_tree2.children.length === 0){
-
- return verb.eval.geom.intersect_aabb_trees( points1, tris1, points2, tris2, aabb_tree1.children[0], aabb_tree2 )
- .concat( verb.eval.geom.intersect_aabb_trees( points1, tris1, points2, tris2, aabb_tree1.children[1], aabb_tree2 ) );
-
- } else if (aabb_tree1.children.length != 0 && aabb_tree2.children.length != 0){
-
- return verb.eval.geom.intersect_aabb_trees( points1, tris1, points2, tris2, aabb_tree1.children[0], aabb_tree2.children[0] )
- .concat( verb.eval.geom.intersect_aabb_trees( points1, tris1, points2, tris2, aabb_tree1.children[0], aabb_tree2.children[1] ) )
- .concat( verb.eval.geom.intersect_aabb_trees( points1, tris1, points2, tris2, aabb_tree1.children[1], aabb_tree2.children[0] ) )
- .concat( verb.eval.geom.intersect_aabb_trees( points1, tris1, points2, tris2, aabb_tree1.children[1], aabb_tree2.children[1] ) );
-
- }
-
-}
-
-//
-// ####make_mesh_aabb_tree( points, tris, tri_indices )
-//
-// Make tree of axis aligned bounding boxes
-//
-// **params**
-// + *Array*, array of length 3 arrays of numbers representing the points
-// + *Array*, array of length 3 arrays of number representing the triangles
-// + *Array*, array of numbers representing the relevant triangles to use to form aabb
-//
-// **returns**
-// + *Array*, a point represented by an array of length (dim)
-//
-
-verb.eval.mesh.make_mesh_aabb_tree = function( points, tris, tri_indices ) {
-
- // build bb
- var aabb = { bounding_box: verb.eval.mesh.make_mesh_aabb( points, tris, tri_indices ),
- children: [] };
-
- // if only one ele, terminate recursion and store the triangles
- if (tri_indices.length === 1){
- aabb.triangle = tri_indices[0];
- return aabb;
- }
-
- // sort triangles in sub mesh
- var sorted_tri_indices = verb.eval.mesh.sort_tris_on_longest_axis( aabb.bounding_box, points, tris, tri_indices )
- , tri_indices_a = sorted_tri_indices.slice( 0, Math.floor( sorted_tri_indices.length / 2 ) )
- , tri_indices_b = sorted_tri_indices.slice( Math.floor( sorted_tri_indices.length / 2 ), sorted_tri_indices.length );
-
- // recurse
- aabb.children = [ verb.eval.mesh.make_mesh_aabb_tree(points, tris, tri_indices_a),
- verb.eval.mesh.make_mesh_aabb_tree(points, tris, tri_indices_b) ];
-
- // return result
- return aabb;
-
-}
-
-//
-// ####make_mesh_aabb( points, tris, tri_indices )
-//
-// Form axis-aligned bounding box from triangles of mesh
-//
-// **params**
-// + *Array*, array of length 3 arrays of numbers representing the points
-// + *Array*, array of length 3 arrays of number representing the triangles
-// + *Array*, array of numbers representing the relevant triangles
-//
-// **returns**
-// + *Array*, a point represented by an array of length (dim)
-//
-
-verb.eval.mesh.make_mesh_aabb = function( points, tris, tri_indices ) {
-
- var bb = new verb.geom.BoundingBox();
-
- tri_indices.forEach(function(x){
- bb.add( points[ tris[ x ][0] ] );
- bb.add( points[ tris[ x ][1] ] );
- bb.add( points[ tris[ x ][2] ] );
- });
-
- return bb;
-
-}
-
-//
-// ####sort_tris_on_longest_axis( container_bb, points, tris, tri_indices )
-//
-// Sort triangles on longest axis
-//
-// **params**
-// + *Number*, integer degree of surface in u direction
-// + *Array*, array of nondecreasing knot values in u direction
-// + *Number*, integer degree of surface in v direction
-// + *Array*, array of nondecreasing knot values in v direction
-//
-// **returns**
-// + *Array*, a point represented by an array of length (dim)
-//
-
-verb.eval.mesh.sort_tris_on_longest_axis = function( container_bb, points, tris, tri_indices ) {
-
- var long_axis = container_bb.get_longest_axis();
-
- var axis_position_map = [];
- for (var i = tri_indices.length - 1; i >= 0; i--) {
-
- var tri_i = tri_indices[i],
- tri_min = verb.eval.mesh.get_min_coordinate_on_axis( points, tris[ tri_i ], long_axis );
-
- axis_position_map.push( [ tri_min, tri_i ] );
-
- }
-
- axis_position_map.sort(function(a,b) { return a[0] > b[0] } );
-
- var sorted_tri_indices = [];
- for (var i = 0, l = axis_position_map.length; i < l; i++) {
- sorted_tri_indices.push( axis_position_map[i][1] );
- }
-
- return sorted_tri_indices;
-
-}
-
-//
-// ####get_min_coordinate_on_axis( points, tri, axis )
-//
-// Get min coordinate on axis
-//
-// **params**
-// + *Array*, array of length 3 arrays of numbers representing the points
-// + *Array*, length 3 array of point indices for the triangle
-//
-// **returns**
-// + *Number*, a point represented by an array of length 3
-//
-
-verb.eval.mesh.get_min_coordinate_on_axis = function( points, tri, axis ) {
-
- var axis_coords = [];
-
- for (var i = 0; i < 3; i++){
- axis_coords.push( points[ tri[i] ][ axis ] );
- }
-
- return Math.min.apply(Math, axis_coords);
-};
-
-//
-// ####get_tri_centroid( points, tri )
-//
-// Get triangle centroid
-//
-// **params**
-// + *Array*, array of length 3 arrays of numbers representing the points
-// + *Array*, length 3 array of point indices for the triangle
-//
-// **returns**
-// + *Array*, a point represented by an array of length 3
-//
-
-verb.eval.geom.get_tri_centroid = function( points, tri ) {
-
- var centroid = [0,0,0];
-
- for (var i = 0; i < 3; i++){
- for (var j = 0; j < 3; j++){
- centroid[j] += points[ tri[i] ][j];
- }
- }
-
- for (var i = 0; i < 3; i++){
- centroid[i] /= 3;
- }
-
- return centroid;
-
-};
-
-//
-// ####get_tri_norm( points, tri )
-//
-// Get triangle normal
-//
-// **params**
-// + *Array*, array of length 3 arrays of numbers representing the points
-// + *Array*, length 3 array of point indices for the triangle
-//
-// **returns**
-// + *Array*, a normal vector represented by an array of length 3
-//
-
-verb.eval.geom.get_tri_norm = function( points, tri ) {
-
- var v0 = points[ tri[0] ]
- , v1 = points[ tri[1] ]
- , v2 = points[ tri[2] ]
- , u = numeric.sub( v1, v0 )
- , v = numeric.sub( v2, v0 )
- , n = numeric.cross( u, v );
-
- return numeric.mul( 1 / numeric.norm2( n ), n );
-
-};
-
-//
-// ####tesselate_rational_surface_naive( degree_u, knots_u, degree_v, knots_v, homo_control_points, divs_u, divs_v )
-//
-// Tesselate a nurbs surface
-//
-// **params**
-// + *Number*, integer degree of surface in u direction
-// + *Array*, array of nondecreasing knot values in u direction
-// + *Number*, integer degree of surface in v direction
-// + *Array*, array of nondecreasing knot values in v direction
-// + *Array*, 3d array of control points, top to bottom is increasing u direction, left to right is increasing v direction,
-// and where each control point is an array of length (dim+1)
-//
-// **returns**
-// + *Array*, first element of array is an array of positions, second element are 3-tuple of triangle windings, third element is the
- // uvs
-
-verb.eval.nurbs.tesselate_rational_surface_naive = function( degree_u, knots_u, degree_v, knots_v, homo_control_points, divs_u, divs_v ) {
-
- if ( divs_u < 1 ) {
- divs_u = 1;
- }
-
- if ( divs_v < 1 ) {
- divs_v = 1;
- }
-
- var span_u = 1 / divs_u,
- span_v = 1 / divs_v;
-
- var points = [];
- var uvs = [];
- var normals = [];
-
- for (var i = 0; i < divs_u + 1; i++) {
- for (var j = 0; j < divs_v + 1; j++) {
-
- var pt_u = i * span_u,
- pt_v = j * span_v;
-
- uvs.push( [pt_u, pt_v] );
-
- var derivs = verb.eval.nurbs.rational_surface_derivs( degree_u, knots_u, degree_v, knots_v, homo_control_points, 1, pt_u, pt_v );
- var pt = derivs[0][0];
-
- points.push( pt );
-
- var normal = numeric.cross( derivs[0][1], derivs[1][0] );
- normals.push( normal );
-
- }
- }
-
- var faces = [];
-
- for (var i = 0; i < divs_u ; i++) {
- for (var j = 0; j < divs_v ; j++) {
-
- var a_i = i * (divs_v + 1) + j,
- b_i = (i + 1) * (divs_v + 1) + j,
- c_i = b_i + 1,
- d_i = a_i + 1,
- abc = [a_i, b_i, c_i],
- acd = [a_i, c_i, d_i];
-
- faces.push(abc);
- faces.push(acd);
-
- }
- }
-
- return { points: points, faces : faces, uvs: uvs, normals: normals };
-
-}
-
-//
-// ####intersect_rational_curves_by_aabb_refine( degree1, knots1, homo_control_points1, degree2, knots2, homo_control_points2, sample_tol, tol )
-//
-// Approximate the intersection of two nurbs surface by axis-aligned bounding box intersection and then refine all solutions.
-//
-// **params**
-// + *Number*, integer degree of curve1
-// + *Array*, array of nondecreasing knot values for curve 1
-// + *Array*, 2d array of homogeneous control points, where each control point is an array of length (dim+1) and form (wi*pi, wi) for curve 1
-// + *Number*, integer degree of curve2
-// + *Array*, array of nondecreasing knot values for curve 2
-// + *Array*, 2d array of homogeneous control points, where each control point is an array of length (dim+1) and form (wi*pi, wi) for curve 2
-// + *Number*, tolerance for the intersection
-//
-// **returns**
-// + *Array*, a 2d array specifying the intersections on u params of intersections on curve 1 and cruve 2
-//
-
-verb.eval.nurbs.intersect_rational_curves_by_aabb_refine = function( degree1, knots1, homo_control_points1, degree2, knots2, homo_control_points2, sample_tol, tol ) {
-
- var ints = verb.eval.nurbs.intersect_rational_curves_by_aabb( degree1, knots1, homo_control_points1, degree2, knots2, homo_control_points2, sample_tol, tol );
-
- return ints.map(function(start_params){
- return verb.eval.nurbs.refine_rational_curve_intersection( degree1, knots1, homo_control_points1, degree2, knots2, homo_control_points2, start_params )
- });
-
-}
-
-
-//
-// ####rational_curve_curve_bb_intersect_refine( degree1, knots1, control_points1, degree2, knots2, control_points2, start_params )
-//
-// Refine an intersection pair for two curves given an initial guess. This is an unconstrained minimization,
-// so the caller is responsible for providing a very good initial guess.
-//
-// **params**
-// + *Number*, integer degree of curve1
-// + *Array*, array of nondecreasing knot values for curve 1
-// + *Array*, 2d array of homogeneous control points, where each control point is an array of length (dim+1)
- // and form (wi*pi, wi) for curve 1
-// + *Number*, integer degree of curve2
-// + *Array*, array of nondecreasing knot values for curve 2
-// + *Array*, 2d array of homogeneous control points, where each control point is an array of length (dim+1)
- // and form (wi*pi, wi) for curve 2
-// + *Array*, length 2 array with first param guess in first position and second param guess in second position
-//
-// **returns**
-// + *Array*, a length 3 array containing the [ distance// distance, u1, u2 ]
-//
-
-verb.eval.nurbs.refine_rational_curve_intersection = function( degree1, knots1, control_points1, degree2, knots2, control_points2, start_params ) {
-
- var objective = function(x) {
-
- var p1 = verb.eval.nurbs.rational_curve_point(degree1, knots1, control_points1, x[0])
- , p2 = verb.eval.nurbs.rational_curve_point(degree2, knots2, control_points2, x[1])
- , p1_p2 = numeric.sub(p1, p2);
-
- return numeric.dot(p1_p2, p1_p2);
- }
-
- var sol_obj = numeric.uncmin( objective, start_params);
- return sol_obj.solution.concat( sol_obj.f );
-
-}
-
-//
-// ####intersect_rational_curves_by_aabb( degree1, knots1, homo_control_points1, degree2, knots2, homo_control_points2, sample_tol, tol )
-//
-// Approximate the intersection of two nurbs surface by axis-aligned bounding box intersection.
-//
-// **params**
-// + *Number*, integer degree of curve1
-// + *Array*, array of nondecreasing knot values for curve 1
-// + *Array*, 2d array of homogeneous control points, where each control point is an array of length (dim+1) and form (wi*pi, wi) for curve 1
-// + *Number*, integer degree of curve2
-// + *Array*, array of nondecreasing knot values for curve 2
-// + *Array*, 2d array of homogeneous control points, where each control point is an array of length (dim+1) and form (wi*pi, wi) for curve 2
-// + *Number*, tolerance for the intersection
-//
-// **returns**
-// + *Array*, array of parameter pairs representing the intersection of the two parameteric polylines
-//
-
-verb.eval.nurbs.intersect_rational_curves_by_aabb = function( degree1, knots1, homo_control_points1, degree2, knots2, homo_control_points2, sample_tol, tol ) {
-
- var up1 = verb.eval.nurbs.rational_curve_adaptive_sample( degree1, knots1, homo_control_points1, sample_tol, true)
- , up2 = verb.eval.nurbs.rational_curve_adaptive_sample( degree2, knots2, homo_control_points2, sample_tol, true)
- , u1 = up1.map( function(el) { return el[0]; })
- , u2 = up2.map( function(el) { return el[0]; })
- , p1 = up1.map( function(el) { return el.slice(1) })
- , p2 = up2.map( function(el) { return el.slice(1) });
-
- return verb.eval.nurbs.intersect_parametric_polylines_by_aabb( p1, p2, u1, u2, tol );
-
-}
-
-//
-// ####intersect_parametric_polylines_by_aabb( p1, p2, u1, u2, tol )
-//
-// Intersect two polyline curves, keeping track of parameterization on each
-//
-// **params**
-// + *Array*, array of point values for curve 1
-// + *Array*, array of parameter values for curve 1, same length as first arg
-// + *Array*, array of point values for curve 2
-// + *Array*, array of parameter values for curve 2, same length as third arg
-// + *Number*, tolerance for the intersection
-//
-// **returns**
-// + *Array*, array of parameter pairs representing the intersection of the two parameteric polylines
-//
-
-verb.eval.nurbs.intersect_parametric_polylines_by_aabb = function( p1, p2, u1, u2, tol ) {
-
- var bb1 = new verb.geom.BoundingBox(p1)
- , bb2 = new verb.geom.BoundingBox(p2);
-
- if ( !bb1.intersects(bb2, tol) ) {
- return [];
- }
-
- if (p1.length === 2 && p2.length === 2 ){
-
- var inter = verb.eval.geom.intersect_segments(p1[0],p1[1], p2[0], p2[1], tol);
-
- if ( inter != null ){
-
- // map the parameters of the segment to the parametric space of the entire polyline
- inter[0][0] = inter[0][0] * ( u1[1]-u1[0] ) + u1[0];
- inter[1][0] = inter[1][0] * ( u2[1]-u2[0] ) + u2[0];
-
- return [ [ inter[0][0], inter[1][0] ] ];
-
- }
-
- } else if (p1.length === 2) {
-
- var p2_mid = Math.ceil( p2.length / 2 ),
- p2_a = p2.slice( 0, p2_mid ),
- p2_b = p2.slice( p2_mid-1 ),
- u2_a = u2.slice( 0, p2_mid ),
- u2_b = u2.slice( p2_mid-1 );
-
- return verb.eval.nurbs.intersect_parametric_polylines_by_aabb(p1, p2_a, u1, u2_a, tol)
- .concat( verb.eval.nurbs.intersect_parametric_polylines_by_aabb(p1, p2_b, u1, u2_b, tol) );
-
- } else if (p2.length === 2) {
-
- var p1_mid = Math.ceil( p1.length / 2 ),
- p1_a = p1.slice( 0, p1_mid ),
- p1_b = p1.slice( p1_mid-1 ),
- u1_a = u1.slice( 0, p1_mid ),
- u1_b = u1.slice( p1_mid-1 );
-
- return verb.eval.nurbs.intersect_parametric_polylines_by_aabb(p1_a, p2, u1_a, u2, tol)
- .concat( verb.eval.nurbs.intersect_parametric_polylines_by_aabb(p1_b, p2, u1_b, u2, tol) );
-
- } else {
-
- var p1_mid = Math.ceil( p1.length / 2 ),
- p1_a = p1.slice( 0, p1_mid ),
- p1_b = p1.slice( p1_mid-1 ),
- u1_a = u1.slice( 0, p1_mid ),
- u1_b = u1.slice( p1_mid-1 ),
-
- p2_mid = Math.ceil( p2.length / 2 ),
- p2_a = p2.slice( 0, p2_mid ),
- p2_b = p2.slice( p2_mid-1 ),
- u2_a = u2.slice( 0, p2_mid ),
- u2_b = u2.slice( p2_mid-1 );
-
- return verb.eval.nurbs.intersect_parametric_polylines_by_aabb(p1_a, p2_a, u1_a, u2_a, tol)
- .concat( verb.eval.nurbs.intersect_parametric_polylines_by_aabb(p1_a, p2_b, u1_a, u2_b, tol) )
- .concat( verb.eval.nurbs.intersect_parametric_polylines_by_aabb(p1_b, p2_a, u1_b, u2_a, tol) )
- .concat( verb.eval.nurbs.intersect_parametric_polylines_by_aabb(p1_b, p2_b, u1_b, u2_b, tol) );
-
- }
-
- return [];
-
-}
-
-//
-// ####intersect_segments( a0, a1, b0, b1, tol )
-//
-// Find the closest parameter on two rays, see http://geomalgorithms.com/a07-_distance.html
-//
-// **params**
-// + *Array*, first point on a
-// + *Array*, second point on a
-// + *Array*, first point on b
-// + *Array*, second point on b
-// + *Number*, tolerance for the intersection
-//
-// **returns**
-// + *Array*, a 2d array specifying the intersections on u params of intersections on curve 1 and cruve 2
-//
-
-verb.eval.geom.intersect_segments = function( a0, a1, b0, b1, tol ) {
-
- // get axis and length of segments
- var a1ma0 = numeric.sub(a1, a0),
- aN = Math.sqrt( numeric.dot(a1ma0, a1ma0) ),
- a = numeric.mul( 1/ aN, a1ma0 ),
- b1mb0 = numeric.sub(b1, b0),
- bN = Math.sqrt( numeric.dot(b1mb0, b1mb0) ),
- b = numeric.mul( 1 / bN, b1mb0 ),
- int_params = verb.eval.geom.intersect_rays(a0, a, b0, b);
-
- if ( int_params != null ) {
-
- var u1 = Math.min( Math.max( 0, int_params[0] / aN ), 1.0),
- u2 = Math.min( Math.max( 0, int_params[1] / bN ), 1.0),
- int_a = numeric.add( numeric.mul( u1, a1ma0 ), a0 ),
- int_b = numeric.add( numeric.mul( u2, b1mb0 ), b0 ),
- dist = numeric.norm2Squared( numeric.sub(int_a, int_b) );
-
- if ( dist < tol*tol ) {
- return [ [u1].concat(int_a), [u2].concat(int_b) ] ;
- }
-
- }
-
- return null;
-
- }
-
-//
-// ####closest_point_on_ray( pt, o, r )
-//
-// Find the closest point on a ray
-//
-// **params**
-// + *Array*, point to project
-// + *Array*, origin for ray
-// + *Array*, direction of ray 1, assumed normalized
-//
-// **returns**
-// + *Array*, [param, pt]
-//
-
-verb.eval.geom.closest_point_on_ray = function( pt, o, r ) {
-
- var o2pt = numeric.sub(pt,o)
- , do2ptr = numeric.dot(o2pt, r)
- , proj = numeric.add(o, numeric.mul(do2ptr, r));
-
- return proj;
-
- }
-
-//
-// ####intersect_rays( a0, a, b0, b )
-//
-// Find the closest parameter on two rays, see http://geomalgorithms.com/a07-_distance.html
-//
-// **params**
-// + *Array*, origin for ray 1
-// + *Array*, direction of ray 1, assumed normalized
-// + *Array*, origin for ray 1
-// + *Array*, direction of ray 1, assumed normalized
-//
-// **returns**
-// + *Array*, a 2d array specifying the intersections on u params of intersections on curve 1 and curve 2
-//
-
-verb.eval.geom.intersect_rays = function( a0, a, b0, b ) {
-
- var dab = numeric.dot( a, b ),
- dab0 = numeric.dot( a, b0 ),
- daa0 = numeric.dot( a, a0 ),
- dbb0 = numeric.dot( b, b0 ),
- dba0 = numeric.dot( b, a0 ),
- daa = numeric.dot( a, a ),
- dbb = numeric.dot( b, b ),
- div = daa*dbb - dab*dab;
-
- // parallel case
- if ( Math.abs( div ) < verb.EPSILON ) {
- return null;
- }
-
- var num = dab * (dab0-daa0) - daa * (dbb0-dba0),
- w = num / div,
- t = (dab0 - daa0 + w * dab)/daa;
-
- return [t, w];
-
- }
-
-
-//
-// ####rational_curve_regular_sample( degree, knots, control_points, num_samples [, include_u] )
-//
-// Sample a NURBS curve assuming parameterization 0 to 1, corresponds to http://ariel.chronotext.org/dd/defigueiredo93adaptive.pdf
-//
-// **params**
-// + *Number*, integer degree
-// + *Array*, array of nondecreasing knot values
-// + *Array*, 1d array of homogeneous control points, where each control point is an array of length (dim+1) and form (wi*pi, wi)
-// + *Number*, integer number of samples
-//
-// **returns**
-// + *Array*, an array of points, prepended by the point param
-//
-
-verb.eval.nurbs.rational_curve_regular_sample = function( degree, knots, control_points, num_samples, include_u ) {
-
- return verb.eval.nurbs.rational_curve_regular_sample_range( degree, knots, control_points, 0, 1.0, num_samples, include_u);
-
-}
-
-//
-// ####rational_curve_regular_sample_range( degree, knots, control_points, start_u, end_u, num_samples, include_u )
-//
-// Sample a NURBS curve assuming parameterization 0 to 1, corresponds to http://ariel.chronotext.org/dd/defigueiredo93adaptive.pdf
-//
-// **params**
-// + *Number*, integer degree
-// + *Array*, array of nondecreasing knot values
-// + *Array*, 1d array of homogeneous control points, where each control point is an array of length (dim+1) and form (wi*pi, wi)
-// + *Number*, start parameter for sampling
-// + *Number*, end parameter for sampling
-// + *Number*, integer number of samples
-// + *Boolean*, whether to prefix the point with the parameter
-//
-// **returns**
-// + *Array*, an dictionary of parameter - point pairs
-//
-
-verb.eval.nurbs.rational_curve_regular_sample_range = function( degree, knots, control_points, start_u, end_u, num_samples, include_u ) {
-
- if (num_samples < 1){
- num_samples = 2;
- }
-
- var p = [],
- span = (end_u - start_u) / (num_samples - 1),
- u = 0;
-
- for (var i = 0; i < num_samples; i++){
-
- u = start_u + span * i;
- if ( include_u ){
- p.push( [u].concat( verb.eval.nurbs.rational_curve_point(degree, knots, control_points, u) ) );
- } else {
- p.push( verb.eval.nurbs.rational_curve_point(degree, knots, control_points, u) );
- }
-
- }
-
- return p;
-
-}
-
-//
-// ####rational_curve_adaptive_sample( degree, knots, control_points, tol, include_u )
-//
-// Sample a NURBS curve assuming parameterization 0 to 1, corresponds to http://ariel.chronotext.org/dd/defigueiredo93adaptive.pdf
-//
-// **params**
-// + *Number*, integer degree
-// + *Array*, array of nondecreasing knot values
-// + *Array*, 2d array of homogeneous control points, where each control point is an array of length (dim+1)
-// and form (wi*pi, wi)
-// + *Number*, tolerance for the adaptive scheme
-// + *Boolean*, whether to prefix the point with the parameter
-//
-// **returns**
-// + *Array*, an array of dim + 1 length where the first element is the param where it was sampled and the remaining the pt
-//
-
-verb.eval.nurbs.rational_curve_adaptive_sample = function( degree, knots, control_points, tol, include_u ) {
-
- // if degree is 1, just return the dehomogenized control points
- if (degree === 1){
- if ( !include_u ) {
- return control_points.map( verb.eval.nurbs.dehomogenize );
- } else {
- // the first element of each array is the parameter
- return control_points.map(function(x, i){
- return [ knots[i+1] ].concat( verb.eval.nurbs.dehomogenize( x ) );
- });
- }
- }
-
- return verb.eval.nurbs.rational_curve_adaptive_sample_range( degree, knots, control_points, 0, 1.0, tol, include_u );
-
-}
-
-//
-// ####rational_curve_adaptive_sample_range( degree, knots, control_points, start_u, end_u, tol, include_u )
-//
-// Sample a NURBS curve at 3 points, facilitating adaptive sampling
-//
-// **params**
-// + *Number*, integer degree
-// + *Array*, array of nondecreasing knot values
-// + *Array*, 2d array of homogeneous control points, where each control point is an array of length (dim+1)
-// and form (wi*pi, wi)
-// + *Number*, start parameter for sampling
-// + *Number*, end parameter for sampling
-// + *Boolean*, whether to prefix the point with the parameter
-//
-// **returns**
-// + *Array*, an array of dim + 1 length where the first element is the param where it was sampled and the remaining the pt
-//
-
-verb.eval.nurbs.rational_curve_adaptive_sample_range = function( degree, knots, control_points, start_u, end_u, tol, include_u ) {
-
- // sample curve at three pts
- var p1 = verb.eval.nurbs.rational_curve_point(degree, knots, control_points, start_u),
- p3 = verb.eval.nurbs.rational_curve_point(degree, knots, control_points, end_u),
- t = 0.5 + 0.2 * Math.random(),
- mid_u = start_u + (end_u - start_u) * t,
- p2 = verb.eval.nurbs.rational_curve_point(degree, knots, control_points, mid_u);
-
- // if the two end control points are coincident, the three point test will always return 0, let's split the curve
- var diff = numeric.sub( p1, p3);
- var diff2 = numeric.sub( p1, p2);
-
- // the first condition checks if the curve makes up a loop, if so, we will need to continue evaluation
- if ( ( numeric.dot( diff, diff ) < tol && numeric.dot( diff2, diff2 ) > tol ) || !verb.eval.nurbs.three_points_are_flat( p1, p2, p3, tol ) ) {
-
- // get the exact middle
- var exact_mid_u = start_u + (end_u - start_u) * 0.5;
-
- // recurse on the two halves
- var left_pts = verb.eval.nurbs.rational_curve_adaptive_sample_range( degree, knots, control_points, start_u, exact_mid_u, tol, include_u )
- , right_pts = verb.eval.nurbs.rational_curve_adaptive_sample_range( degree, knots, control_points, exact_mid_u, end_u, tol, include_u );
-
- // concatenate the two
- return left_pts.slice(0, -1).concat(right_pts);
-
- } else {
-
- if (include_u){
- return [ [ start_u ].concat(p1) , [end_u].concat(p3) ];
- } else {
- return [ p1, p3 ];
- }
-
- }
-}
-
-
-
-//
-// ####three_points_are_flat( p1, p2, p3, tol )
-//
-// Determine if three points form a straight line within a given tolerance for their 2 * squared area
-//
-// * p2
-// / \
-// / \
-// / \
-// / \
-// * p1 ---- * p3
-//
-// The area metric is 2 * the squared norm of the cross product of two edges, requiring no square roots and no divisions
-//
-// **params**
-// + *Array*, p1
-// + *Array*, p2
-// + *Array*, p3
-// + *Number*, The tolerance for whether the three points form a line
-//
-// **returns**
-// + *Number*, Whether the triangle passes the test
-//
-verb.eval.nurbs.three_points_are_flat = function( p1, p2, p3, tol ) {
-
-
- // find the area of the triangle without using a square root
- var p2mp1 = numeric.sub( p2, p1 )
- , p3mp1 = numeric.sub( p3, p1 )
- , norm = crossprod( p2mp1, p3mp1 )
- , area = numeric.dot( norm, norm );
-
- return area < tol;
-
-}
-
//
// ####curve_knot_insert( degree, knots, control_points, u, s, r )
//
View
519 src/eval/geom.js
@@ -0,0 +1,519 @@
+
+// ###verb.eval
+// This defines verb's core geometry library which is called by the current Engine.
+//
+// ####get_sweep1_surface( profile_knots, profile_degree, profile_control_points, profile_weights, rail_knots, rail_degree, rail_control_points, rail_weights )
+//
+// Generate the control points, weights, and knots of an elliptical arc
+//
+// **params**
+// + *Array*, the center
+// + *Array*, the xaxis
+// + *Array*, orthogonal yaxis
+// + *Number*, xradius of the ellipse arc
+// + *Number*, yradius of the ellipse arc
+// + *Number*, start angle of the ellipse arc, between 0 and 2pi, where 0 points at the xaxis
+// + *Number*, end angle of the arc, between 0 and 2pi, greater than the start angle
+//
+// **returns**
+// + *Object*, an object with the following properties: control_points, weights, knots, degree
+//
+verb.eval.nurbs.get_sweep1_surface = function( profile_knots, profile_degree, profile_control_points, profile_weights, rail_knots, rail_degree, rail_control_points, rail_weights ) {
+
+ // for each point on rail, move all of the points
+ var homo_rail = verb.eval.nurbs.homogenize_1d( rail_control_points, rail_weights )
+ , rail_start = verb.eval.nurbs.rational_curve_point( rail_degree, rail_knots, homo_rail, 0 )
+ , span = 1.0 / rail_control_points.length
+ , control_points = []
+ , weights = [];
+
+ for (var i = 0; i < rail_control_points.length; i++ ){
+
+ // evaluate the point on the curve, subtracting it from the first point
+ var rail_point = verb.eval.nurbs.rational_curve_point( rail_degree, rail_knots, homo_rail, i * span )
+ , rail_offset = numeric.sub( rail_point, rail_start )
+ , row_control_points = []
+ , row_weights = [];
+
+ for (var j = 0; j < profile_control_points.length; j++ ){
+
+ row_control_points.push( numeric.add(rail_offset, profile_control_points[j] ) );
+ row_weights.push( profile_weights[j] * rail_weights[i] );
+
+ }
+
+ control_points.push( row_control_points);
+ weights.push( row_weights );
+ }
+
+ return {"knots_u": rail_knots,
+ "knots_v": profile_knots,
+ "control_points": control_points,
+ "degree_u": rail_degree,
+ "degree_v": profile_degree,
+ "weights": weights };
+
+}
+
+//
+// ####get_ellipse_arc( center, xaxis, yaxis, xradius, yradius, start_angle, end_angle )
+//
+// Generate the control points, weights, and knots of an elliptical arc
+//
+// **params**
+// + *Array*, the center
+// + *Array*, the xaxis
+// + *Array*, orthogonal yaxis
+// + *Number*, xradius of the ellipse arc
+// + *Number*, yradius of the ellipse arc
+// + *Number*, start angle of the ellipse arc, between 0 and 2pi, where 0 points at the xaxis
+// + *Number*, end angle of the arc, between 0 and 2pi, greater than the start angle
+//
+// **returns**
+// + *Object*, an object with the following properties: control_points, weights, knots, degree
+//
+
+verb.eval.nurbs.get_ellipse_arc = function( center, xaxis, yaxis, xradius, yradius, start_angle, end_angle ) {
+
+ // if the end angle is less than the start angle, do a circle
+ if (end_angle < start_angle) end_angle = 2 * Math.PI + start_angle;
+
+ var theta = end_angle - start_angle
+ , narcs = 0;
+
+ // how many arcs?
+ if (theta <= Math.PI / 2) {
+ narcs = 1;
+ } else {
+ if (theta <= Math.PI){
+ narcs = 2;
+ } else if (theta <= 3 * Math.PI / 2){
+ narcs = 3;
+ } else {
+ narcs = 4;
+ }
+ }
+
+ var dtheta = theta / narcs
+ , n = 2 * narcs
+ , w1 = Math.cos( dtheta / 2)
+ , P0 = numeric.add( center, numeric.mul( xradius, Math.cos(start_angle), xaxis), numeric.mul( yradius, Math.sin(start_angle), yaxis ) )
+ , T0 = numeric.sub( numeric.mul( Math.cos(start_angle), yaxis ), numeric.mul( Math.sin(start_angle), xaxis) )
+ , Pw = verb.eval.nurbs.zeros_1d( narcs * 2 )
+ , U = verb.eval.nurbs.zeros_1d( 2 *narcs + 3 )
+ , index = 0
+ , angle = start_angle
+ , W = verb.eval.nurbs.zeros_1d( narcs * 2 );
+
+ Pw[0] = P0;
+ W[0] = 1;
+
+ for (var i = 1; i <= narcs; i++){
+
+ angle += dtheta;
+ var P2 = numeric.add( center, numeric.mul( xradius, Math.cos(angle), xaxis), numeric.mul( yradius, Math.sin(angle), yaxis ) )
+
+ W[index+2] = 1;
+ Pw[index+2] = P2;
+
+ var T2 = numeric.sub( numeric.mul( Math.cos(angle), yaxis ), numeric.mul( Math.sin(angle), xaxis) )
+
+ var params = verb.eval.geom.intersect_rays(P0, numeric.mul( 1 / numeric.norm2(T0), T0), P2, numeric.mul( 1 / numeric.norm2(T2), T2));
+ var P1 = numeric.add( P0, numeric.mul(T0, params[0]));
+
+ W[index+1] = w1;
+ Pw[index+1] = P1;
+
+ index += 2;
+
+ if (i < narcs){
+ P0 = P2;
+ T0 = T2;
+ }
+ }
+
+ var j = 2 * narcs + 1;
+
+ for (var i = 0; i < 3; i++){
+ U[i] = 0.0;
+ U[i+j] = 1.0;
+ }
+
+ switch (narcs){
+ case 1: break;
+ case 2: U[3] = U[4] = 0.5; break;
+ case 3: U[3] = U[4] = 1/3;
+ U[5] = U[6] = 2/3; break;
+ case 4: U[3] = U[4] = 0.25;
+ U[5] = U[6] = 0.5;
+ U[7] = U[8] = 0.75; break;
+ }
+
+ return {knots: U, control_points: Pw, degree: 2, weights: W };
+
+}
+
+//
+// ####get_sphere_surface( center, axis, xaxis, radius )
+//
+// Generate the control points, weights, and knots of a sphere
+//
+// **params**
+// + *Array*, the center of the sphere
+// + *Array*, normalized axis of sphere
+// + *Array*, vector perpendicular to axis of sphere, starting the rotation of the sphere
+// + *Number*, radius of the sphere
+//
+// **returns**
+// + *Object*, an object with the following properties: control_points, weights, knots_u, knots_v, degree_u, degree_v
+//
+
+verb.eval.nurbs.get_sphere_surface = function( center, axis, xaxis, radius ){
+
+ var arc = verb.eval.nurbs.get_arc(center, numeric.mul(axis, -1), xaxis, radius, 0, Math.PI );
+
+ return verb.eval.nurbs.get_revolved_surface( center, axis, 2 * Math.PI, arc.knots, arc.degree, arc.control_points, arc.weights );
+
+}
+
+
+//
+// ####get_polyline_curve( pts )
+//
+// Generate the control points, weights, and knots of a polyline curve
+//
+// **params**
+// + *Array*, array of points in curve
+//
+// **returns**
+// + *Object*, an object with the following properties: control_points, weights, knots, degree
+//
+
+verb.eval.nurbs.get_polyline_curve = function( pts ){
+
+ var num_spans = pts.length - 1
+ , span = 1.0 / num_spans
+ , knots = [0,0];
+
+ for (var i = 1; i < num_spans; i++){
+ knots.push(i * span);
+ }
+
+ knots.push(1);
+ knots.push(1);
+
+ var weights = [];
+
+ for (var i = 0; i < pts.length; i++){
+ weights.push(1);
+ }
+
+ return {
+ "knots": knots,
+ "control_points": pts.slice(0),
+ "degree": 1,
+ "weights": weights
+ };
+
+}
+
+//
+// ####get_4pt_surface( p1, p2, p3, p4 )
+//
+// Generate the control points, weights, and knots of a surface define by 3 points
+//
+// **params**
+// + *Array*, first point in counter-clockwise form
+// + *Array*, second point in counter-clockwise form
+// + *Array*, third point in counter-clockwise form
+// + *Array*, forth point in counter-clockwise form
+//
+// **returns**
+// + *Object*, an object with the following properties: control_points, weights, knots_u, knots_v, degree_u, degree_v
+//
+
+verb.eval.nurbs.get_4pt_surface = function( p1, p2, p3, p4 ){
+
+ return {"knots_u": [0,0,1,1],
+ "knots_v": [0,0,1,1],
+ "control_points": [ [p1, p4], [p2, p3] ],
+ "degree_u": 1,
+ "degree_v": 1,
+ "weights": [ [1, 1], [1, 1] ] };
+
+}
+
+//
+// ####get_cylinder_surface( axis, xaxis, base, height, radius )
+//
+// Generate the control points, weights, and knots of a cylinder
+//
+// **params**
+// + *Array*, normalized axis of cylinder
+// + *Array*, xaxis in plane of cylinder
+// + *Array*, position of base of cylinder
+// + *Number*, height from base to top
+// + *Number*, radius of the cylinder
+//
+// **returns**
+// + *Object*, an object with the following properties: control_points, weights, knots_u, knots_v, degree_u, degree_v
+//
+
+verb.eval.nurbs.get_cylinder_surface = function( axis, xaxis, base, height, radius ){
+
+ var yaxis = crossprod( axis, xaxis )
+ , angle = 2 * Math.PI
+ , circ = verb.eval.nurbs.get_arc( base, xaxis, yaxis, radius, 0, 2 * Math.PI );
+
+ return verb.eval.nurbs.get_extruded_surface( axis, height, circ.knots, circ.degree, circ.control_points, circ.weights );
+
+}
+
+//
+// ####get_cone_surface( axis, xaxis, base, height, radius )
+//
+// Generate the control points, weights, and knots of a cone
+//
+// **params**
+// + *Array*, normalized axis of cone
+// + *Array*, position of base of cone
+// + *Number*, height from base to tip
+// + *Number*, radius at the base of the cone
+//
+// **returns**
+// + *Object*, an object with the following properties: control_points, weights, knots, degree
+//
+
+verb.eval.nurbs.get_cone_surface = function( axis, xaxis, base, height, radius ){
+
+ var angle = 2 * Math.PI
+ , prof_degree = 1
+ , prof_ctrl_pts = [ numeric.add( base, numeric.mul( height, axis ) ), numeric.add( base, numeric.mul( radius, xaxis ) )]
+ , prof_knots = [0,0,1,1]
+ , prof_weights = [1,1];
+
+ return verb.eval.nurbs.get_revolved_surface(base, axis, angle, prof_knots, prof_degree, prof_ctrl_pts, prof_weights);
+
+}
+
+//
+// ####get_extruded_surface( axis, length, prof_knots, prof_degree, prof_control_points, prof_weights)
+//
+// Generate the control points, weights, and knots of an extruded surface
+//
+// **params**
+// + *Array*, axis of the extrusion
+// + *Array*, length of the extrusion
+// + *Number*, degree of the profile
+// + *Number*, control points of the profile
+// + *Number*, weights of the profile
+//
+// **returns**
+// + *Object*, an object with the following properties: control_points, weights, knots, degree
+//
+
+verb.eval.nurbs.get_extruded_surface = function( axis, length, prof_knots, prof_degree, prof_control_points, prof_weights){
+
+ var control_points = verb.eval.nurbs.zeros_2d( 2, prof_control_points.length )
+ , weights = verb.eval.nurbs.zeros_2d( 2, prof_control_points.length );
+
+ // original control points
+ for (var j = 0; j < prof_control_points.length; j++){
+ control_points[0][j] = prof_control_points[j];
+ weights[0][j] = prof_weights[j];
+ }
+
+ // build translated control points
+ var translation = numeric.mul(axis, length);
+
+ for (var j = 0; j < prof_control_points.length; j++){
+ control_points[1][j] = numeric.add( translation, prof_control_points[j] );
+ weights[1][j] = prof_weights[j];
+ }
+
+ // return all parameters
+ return {"knots_u": [0,0,1,1],
+ "knots_v": prof_knots,
+ "control_points": control_points,
+ "degree_u": 1,
+ "degree_v": prof_degree,
+ "weights": weights };
+}
+
+//
+// ####get_revolved_surface( center, axis, theta, prof_knots, prof_degree, prof_control_points, prof_weights)
+//
+// Generate the control points, weights, and knots of a revolved surface
+// (Corresponds to Algorithm A7.1 from Piegl & Tiller)
+//
+// **params**
+// + *Array*, center of the rotation axis
+// + *Array*, axis of the rotation axis
+// + *Number*, angle to revolve around axis
+// + *Number*, degree of the generatrix
+// + *Number*, control points of the generatrix
+// + *Number*, weights of the generatrix
+//
+// **returns**
+// + *Object*, an object with the following properties: control_points, weights, knots, degree
+//
+
+// helper method
+
+function crossprod(u,v) {
+ return [u[1]*v[2]-u[2]*v[1],u[2]*v[0]-u[0]*v[2],u[0]*v[1]-u[1]*v[0]];
+}
+
+verb.eval.nurbs.get_revolved_surface = function( center, axis, theta, prof_knots, prof_degree, prof_control_points, prof_weights){
+
+ var narcs, knots_u, control_points, weights;
+
+ if (theta <= Math.PI / 2) { // less than 90
+ narcs = 1;
+ knots_u = verb.eval.nurbs.zeros_1d( 6 + 2 * (narcs-1) );
+ } else {
+ if (theta <= Math.PI){ // between 90 and 180
+ narcs = 2;
+ knots_u = verb.eval.nurbs.zeros_1d( 6 + 2 * (narcs-1) );
+ knots_u[3]= knots_u[4] = 0.5;
+ } else if (theta <= 3 * Math.PI / 2){ // between 180 and 270
+ narcs = 3;
+ knots_u = verb.eval.nurbs.zeros_1d( 6 + 2 * (narcs-1) );
+ knots_u[3]= knots_u[4] = 1/3;
+ knots_u[5]= knots_u[6] = 2/3;
+ } else { // between 270 and 360
+ narcs = 4;
+ knots_u = verb.eval.nurbs.zeros_1d( 6 + 2 * (narcs-1) );
+ knots_u[3]= knots_u[4] = 1/4;
+ knots_u[5]= knots_u[6] = 1/2;
+ knots_u[7]= knots_u[8] = 3/4;
+ }
+ }
+
+ var dtheta = theta / narcs // divide the interval into several points
+ , j = 3 + 2 * (narcs-1);
+
+ // initialize the start and end knots
+ // keep in mind that we only return the knot vector for the
+
+ for (var i = 0; i < 3; j++, i++){
+ knots_u[i] = 0.0;
+ knots_u[j] = 1.0;
+ }
+
+ // do some initialization
+ var n = 2 * narcs
+ , wm = Math.cos( dtheta/2.0 )
+ , angle = 0.0
+ , sines = verb.eval.nurbs.zeros_1d( narcs + 1)
+ , cosines = verb.eval.nurbs.zeros_1d( narcs + 1)
+ , control_points = verb.eval.nurbs.zeros_2d( 2*narcs + 1, prof_control_points.length )
+ , weights = verb.eval.nurbs.zeros_2d( 2*narcs + 1, prof_control_points.length );
+
+ // initialize the sines and cosines
+ for (var i = 1; i <= narcs; i++){
+ angle += dtheta;
+ cosines[i] = Math.cos(angle);
+ sines[i] = Math.sin(angle);
+ }
+
+ // for each pt in the generatrix
+ // i.e. for each row of the 2d knot vectors
+ for (j = 0; j < prof_control_points.length; j++){
+
+ // get the closest point of the generatrix point on the axis
+ var O = verb.eval.geom.closest_point_on_ray(prof_control_points[j], center, axis)
+ // X is the vector from the axis to generatrix control pt
+ , X = numeric.sub( prof_control_points[j], O )
+ // radius at that height
+ , r = numeric.norm2(X)
+ // Y is perpendicular to X and axis, and complete the coordinate system
+ , Y = crossprod(axis,X);
+
+ if ( r > verb.EPSILON ){
+ X = numeric.mul( 1 / r, X);
+ Y = numeric.mul( 1 / r, Y);
+ }
+
+ // the first row of control_points and weights is just the generatrix
+ control_points[0][j] = prof_control_points[j];
+ var P0 = prof_control_points[j];
+ weights[0][j] = prof_weights[j];
+
+ // store T0 as the Y vector
+ var T0 = Y
+ , index = 0
+ , angle = 0.0;
+
+ // proceed around the circle
+ for (var i = 1; i <= narcs; i++){
+
+ // O + r * cos(theta) * X + r * sin(theta) * Y
+ // rotated generatrix pt
+ var P2 = r == 0 ? O : numeric.add( O, numeric.mul( r, cosines[i], X), numeric.mul( r, sines[i], Y) );
+
+ control_points[index+2][j] = P2;
+ weights[index+2][j] = prof_weights[j];
+
+ // construct the vector tangent to the rotation
+ var T2 = numeric.sub( numeric.mul( cosines[i], Y), numeric.mul(sines[i], X));
+
+ // construct the next control pt
+ if (r == 0){
+ control_points[index+1][j] = O;
+ } else {
+ var params = verb.eval.geom.intersect_rays(P0, numeric.mul( 1 / numeric.norm2(T0), T0), P2, numeric.mul( 1 / numeric.norm2(T2), T2));
+ var P1 = numeric.add( P0, numeric.mul(T0, params[0]));
+
+ control_points[index+1][j] = P1;
+ }
+
+ weights[index+1][j] = wm * prof_weights[j];
+
+ index += 2;
+
+ if (i < narcs)
+ {
+ P0 = P2;
+ T0 = T2;
+ }
+
+ }
+ }
+
+ // store all of the parameters
+ return {"knots_u": knots_u,
+ "knots_v": prof_knots,
+ "control_points": control_points,
+ "degree_u": 2,
+ "degree_v": prof_degree,
+ "weights": weights };
+
+}
+
+
+
+//
+// ####get_arc( center, xaxis, yaxis, radius, start_angle, end_angle )
+//
+// Generate the control points, weights, and knots of an arbitrary arc
+// (Corresponds to Algorithm A7.1 from Piegl & Tiller)
+//
+// **params**
+// + *Array*, the center of the arc
+// + *Array*, the xaxis of the arc
+// + *Array*, orthogonal yaxis of the arc
+// + *Number*, radius of the arc
+// + *Number*, start angle of the arc, between 0 and 2pi
+// + *Number*, end angle of the arc, between 0 and 2pi, greater than the start angle
+//
+// **returns**
+// + *Object*, an object with the following properties: control_points, weights, knots, degree
+//
+
+verb.eval.nurbs.get_arc = function( center, xaxis, yaxis, radius, start_angle, end_angle ) {
+
+ return verb.eval.nurbs.get_ellipse_arc( center, xaxis, yaxis, radius, radius, start_angle, end_angle );
+
+}
+
View
1,023 src/eval/intersect.js
@@ -0,0 +1,1023 @@
+
+
+//
+// ####intersect_rational_curve_surface_by_aabb( degree_u, knots_u, degree_v, knots_v, homo_control_points, degree_crv, knots_crv, homo_control_points_crv, sample_tol, tol )
+//
+// Get the intersection of a NURBS curve and a NURBS surface by axis-aligned bounding box intersection and refinement
+//
+// **params**
+// + *Number*, integer degree of surface in u direction
+// + *Array*, array of nondecreasing knot values in u direction
+// + *Number*, integer degree of surface in v direction
+// + *Array*, array of nondecreasing knot values in v direction
+// + *Array*, 3d array of homogeneous control points, top to bottom is increasing u direction, left to right is increasing v direction,
+// and where each control point is an array of length (dim+1)
+// + *Number*, integer degree of curve
+// + *Array*, array of nondecreasing knot values
+// + *Array*, 2d array of homogeneous control points, where each control point is an array of length (dim+1)
+// and form (wi*pi, wi)
+// + *Number*, the sample tolerance of the curve
+// + *Number*, tolerance for the curve intersection
+// + *Number*, integer number of divisions of the surface in the U direction for initial approximation (placeholder until adaptive tess of surfaces)
+// + *Number*, integer number of divisions of the surface in the V direction for initial approximation (placeholder until adaptive tess of surfaces)
+//
+// **returns**
+// + *Array*, array of intersection objects, each holding:
+// - a "point" property where intersections took place
+// - a "p" the parameter on the curve
+// - a "uv" the parameter on the surface
+// - a "face" the index of the face where the intersection took place
+//
+
+verb.eval.nurbs.intersect_rational_curve_surface_by_aabb_refine = function( degree_u, knots_u, degree_v, knots_v, homo_control_points_srf, degree_crv, knots_crv, homo_control_points_crv, sample_tol, tol, divs_u, divs_v ) {
+
+ // get the approximate intersections
+ var ints = verb.eval.nurbs.intersect_rational_curve_surface_by_aabb( degree_u, knots_u, degree_v, knots_v, homo_control_points_srf, degree_crv, knots_crv, homo_control_points_crv, sample_tol, tol, divs_u, divs_v );
+
+ // refine them
+ return ints.map(function( inter ){
+
+ // get intersection params
+ var start_params = [inter.p, inter.uv[0], inter.uv[1] ]
+
+ // refine the parameters
+ , refined_params = verb.eval.nurbs.refine_rational_curve_surface_intersection( degree_u, knots_u, degree_v, knots_v, homo_control_points_srf, degree_crv, knots_crv, homo_control_points_crv, start_params );
+
+ // update the inter object
+ inter.p = refined_params[0]; </