1 /* $Id: curves.c,v 1.4 2004-12-21 11:58:14 btb Exp $ */
3 THE COMPUTER CODE CONTAINED HEREIN IS THE SOLE PROPERTY OF PARALLAX
4 SOFTWARE CORPORATION ("PARALLAX"). PARALLAX, IN DISTRIBUTING THE CODE TO
5 END-USERS, AND SUBJECT TO ALL OF THE TERMS AND CONDITIONS HEREIN, GRANTS A
6 ROYALTY-FREE, PERPETUAL LICENSE TO SUCH END-USERS FOR USE BY SUCH END-USERS
7 IN USING, DISPLAYING, AND CREATING DERIVATIVE WORKS THEREOF, SO LONG AS
8 SUCH USE, DISPLAY OR CREATION IS FOR NON-COMMERCIAL, ROYALTY OR REVENUE
9 FREE PURPOSES. IN NO EVENT SHALL THE END-USER USE THE COMPUTER CODE
10 CONTAINED HEREIN FOR REVENUE-BEARING PURPOSES. THE END-USER UNDERSTANDS
11 AND AGREES TO THE TERMS HEREIN AND ACCEPTS THE SAME BY USE OF THIS FILE.
12 COPYRIGHT 1993-1998 PARALLAX SOFTWARE CORPORATION. ALL RIGHTS RESERVED.
17 * curve generation stuff
22 static char rcsid[] = "$Id: curves.c,v 1.4 2004-12-21 11:58:14 btb Exp $";
48 #define ONE_OVER_SQRT2 F1_0 * 0.707106781
53 segment *OriginalMarkedSeg;
55 int OriginalMarkedSide;
56 segment *CurveSegs[MAX_SEGMENTS];
58 const fix Mh[4][4] = { { 2*F1_0, -2*F1_0, 1*F1_0, 1*F1_0 },
59 {-3*F1_0, 3*F1_0, -2*F1_0, -1*F1_0 },
60 { 0*F1_0, 0*F1_0, 1*F1_0, 0*F1_0 },
61 { 1*F1_0, 0*F1_0, 0*F1_0, 0*F1_0 } };
63 void generate_banked_curve(fix maxscale, vms_equation coeffs);
65 void create_curve(vms_vector *p1, vms_vector *p4, vms_vector *r1, vms_vector *r4, vms_equation *coeffs) {
66 // Q(t) = (2t^3 - 3t^2 + 1) p1 + (-2t^3 + 3t^2) p4 + (t~3 - 2t^2 + t) r1 + (t^3 - t^2 ) r4
68 coeffs->n.x3 = fixmul(2*F1_0,p1->x) - fixmul(2*F1_0,p4->x) + r1->x + r4->x;
69 coeffs->n.x2 = fixmul(-3*F1_0,p1->x) + fixmul(3*F1_0,p4->x) - fixmul(2*F1_0,r1->x) - fixmul(1*F1_0,r4->x);
72 coeffs->n.y3 = fixmul(2*F1_0,p1->y) - fixmul(2*F1_0,p4->y) + r1->y + r4->y;
73 coeffs->n.y2 = fixmul(-3*F1_0,p1->y) + fixmul(3*F1_0,p4->y) - fixmul(2*F1_0,r1->y) - fixmul(1*F1_0,r4->y);
76 coeffs->n.z3 = fixmul(2*F1_0,p1->z) - fixmul(2*F1_0,p4->z) + r1->z + r4->z;
77 coeffs->n.z2 = fixmul(-3*F1_0,p1->z) + fixmul(3*F1_0,p4->z) - fixmul(2*F1_0,r1->z) - fixmul(1*F1_0,r4->z);
83 vms_vector evaluate_curve(vms_equation *coeffs, int degree, fix t) {
87 if (degree!=3) printf("ERROR: for Hermite Curves degree must be 3\n");
89 t2 = fixmul(t,t); t3 = fixmul(t2,t);
91 coord.x = fixmul(coeffs->n.x3,t3) + fixmul(coeffs->n.x2,t2) + fixmul(coeffs->n.x1,t) + coeffs->n.x0;
92 coord.y = fixmul(coeffs->n.y3,t3) + fixmul(coeffs->n.y2,t2) + fixmul(coeffs->n.y1,t) + coeffs->n.y0;
93 coord.z = fixmul(coeffs->n.z3,t3) + fixmul(coeffs->n.z2,t2) + fixmul(coeffs->n.z1,t) + coeffs->n.z0;
99 fix curve_dist(vms_equation *coeffs, int degree, fix t0, vms_vector *p0, fix dist) {
103 if (degree!=3) printf("ERROR: for Hermite Curves degree must be 3\n");
105 for (t=t0;t<1*F1_0;t+=0.001*F1_0) {
106 coord = evaluate_curve(coeffs, 3, t);
107 diff = dist - vm_vec_dist(&coord, p0);
108 if (diff<ACCURACY) //&&(diff>-ACCURACY))
115 void curve_dir(vms_equation *coeffs, int degree, fix t0, vms_vector *dir) {
118 if (degree!=3) printf("ERROR: for Hermite Curves degree must be 3\n");
122 dir->x = fixmul(3*F1_0,fixmul(coeffs->n.x3,t2)) + fixmul(2*F1_0,fixmul(coeffs->n.x2,t0)) + coeffs->n.x1;
123 dir->y = fixmul(3*F1_0,fixmul(coeffs->n.y3,t2)) + fixmul(2*F1_0,fixmul(coeffs->n.y2,t0)) + coeffs->n.y1;
124 dir->z = fixmul(3*F1_0,fixmul(coeffs->n.z3,t2)) + fixmul(2*F1_0,fixmul(coeffs->n.z2,t0)) + coeffs->n.z1;
125 vm_vec_normalize( dir );
129 void plot_parametric(vms_equation *coeffs, fix min_t, fix max_t, fix del_t) {
130 vms_vector coord, dcoord;
134 gr_box( 75, 40, 325, 290 );
135 gr_box( 75, 310, 325, 560 );
136 gr_box( 475, 310, 725, 560 );
137 //gr_pal_fade_in( grd_curscreen->pal );
139 for (t=min_t;t<max_t-del_t;t+=del_t) {
142 coord = evaluate_curve(coeffs, 3, t);
143 dcoord = evaluate_curve(coeffs, 3, dt);
146 gr_line ( 75*F1_0 + coord.x, 290*F1_0 - coord.z, 75*F1_0 + dcoord.x, 290*F1_0 - dcoord.z );
148 gr_line ( 75*F1_0 + coord.x, 560*F1_0 - coord.y, 75*F1_0 + dcoord.x, 560*F1_0 - dcoord.y );
150 gr_line ( 475*F1_0 + coord.z, 560*F1_0 - coord.y, 475*F1_0 + dcoord.z, 560*F1_0 - dcoord.y );
157 vms_vector *vm_vec_interp(vms_vector *result, vms_vector *v0, vms_vector *v1, fix scale) {
160 vm_vec_sub(&tvec, v1, v0);
161 vm_vec_scale_add(result, v0, &tvec, scale);
162 vm_vec_normalize(result);
166 vms_vector p1, p4, r1, r4;
167 vms_vector r4t, r1save;
169 int generate_curve( fix r1scale, fix r4scale ) {
170 vms_vector vec_dir, tvec;
171 vms_vector coord,prev_point;
173 fix enddist, nextdist;
176 fixang rangle, uangle;
178 compute_center_point_on_side( &p1, Cursegp, Curside );
182 extract_right_vector_from_segment(Cursegp, &r1);
183 vm_vec_scale( &r1, -F1_0 );
186 extract_up_vector_from_segment(Cursegp, &r1);
189 extract_right_vector_from_segment(Cursegp, &r1);
192 extract_up_vector_from_segment(Cursegp, &r1);
193 vm_vec_scale( &r1, -F1_0 );
196 extract_forward_vector_from_segment(Cursegp, &r1);
199 extract_forward_vector_from_segment(Cursegp, &r1);
200 vm_vec_scale( &r1, -F1_0 );
204 compute_center_point_on_side( &p4, Markedsegp, Markedside );
206 switch( Markedside ) {
208 extract_right_vector_from_segment(Markedsegp, &r4);
209 extract_up_vector_from_segment(Markedsegp, &r4t);
212 extract_up_vector_from_segment(Markedsegp, &r4);
213 vm_vec_scale( &r4, -F1_0 );
214 extract_forward_vector_from_segment(Markedsegp, &r4t);
215 vm_vec_scale( &r4t, -F1_0 );
218 extract_right_vector_from_segment(Markedsegp, &r4);
219 vm_vec_scale( &r4, -F1_0 );
220 extract_up_vector_from_segment(Markedsegp, &r4t);
223 extract_up_vector_from_segment(Markedsegp, &r4);
224 extract_forward_vector_from_segment(Markedsegp, &r4t);
227 extract_forward_vector_from_segment(Markedsegp, &r4);
228 vm_vec_scale( &r4, -F1_0 );
229 extract_up_vector_from_segment(Markedsegp, &r4t);
232 extract_forward_vector_from_segment(Markedsegp, &r4);
233 extract_up_vector_from_segment(Markedsegp, &r4t);
239 vm_vec_scale(&r1,r1scale);
240 vm_vec_scale(&r4,r4scale);
242 create_curve( &p1, &p4, &r1, &r4, &coeffs );
243 OriginalSeg = Cursegp;
244 OriginalMarkedSeg = Markedsegp;
245 OriginalSide = Curside;
246 OriginalMarkedSide = Markedside;
248 coord = prev_point = p1;
252 enddist = F1_0; nextdist = 0;
253 while ( enddist > fixmul( nextdist, 1.5*F1_0 )) {
254 vms_matrix rotmat,rotmat2;
260 extract_forward_vector_from_segment(Cursegp, &tvec);
261 nextdist = vm_vec_mag(&tvec); // nextdist := distance to next point
262 t = curve_dist(&coeffs, 3, t, &prev_point, nextdist); // t = argument at which function is forward vector magnitude units away from prev_point (in 3-space, not along curve)
263 coord = evaluate_curve(&coeffs, 3, t); // coord := point about forward vector magnitude units away from prev_point
264 enddist = vm_vec_dist(&coord, &p4); // enddist := distance from current to end point, vec_dir used as a temporary variable
265 //vm_vec_normalize(vm_vec_sub(&vec_dir, &coord, &prev_point));
266 vm_vec_normalized_dir(&vec_dir, &coord, &prev_point);
267 if (!med_attach_segment( Cursegp, &New_segment, Curside, AttachSide )) {
268 med_extract_matrix_from_segment( Cursegp,&rotmat ); // rotmat := matrix describing orientation of Cursegp
269 vm_vec_rotate(&tdest,&vec_dir,&rotmat); // tdest := vec_dir in reference frame of Cursegp
272 vm_vector_2_matrix(&rotmat2,&vec_dir,NULL,NULL);
273 // mprintf(0, "[ [%6.2f %6.2f %6.2f]", f2fl(rotmat2.m1), f2fl(rotmat2.m2), f2fl(rotmat2.m3));
274 // mprintf(0, " [%6.2f %6.2f %6.2f]", f2fl(rotmat2.m4), f2fl(rotmat2.m5), f2fl(rotmat2.m6));
275 // mprintf(0, " [%6.2f %6.2f %6.2f] ]\n", f2fl(rotmat2.m7), f2fl(rotmat2.m8), f2fl(rotmat2.m9));
277 med_rotate_segment( Cursegp, &rotmat2 );
279 Curside = Side_opposite[AttachSide];
281 CurveSegs[CurveNumSegs]=Cursegp;
286 extract_up_vector_from_segment( Cursegp,&tvec );
287 uangle = vm_vec_delta_ang( &tvec, &r4t, &r4 );
288 if (uangle >= F1_0 * 1/8) uangle -= F1_0 * 1/4;
289 if (uangle >= F1_0 * 1/8) uangle -= F1_0 * 1/4;
290 if (uangle <= -F1_0 * 1/8) uangle += F1_0 * 1/4;
291 if (uangle <= -F1_0 * 1/8) uangle += F1_0 * 1/4;
292 extract_right_vector_from_segment( Cursegp,&tvec );
293 rangle = vm_vec_delta_ang( &tvec, &r4t, &r4 );
294 if (rangle >= F1_0/8) rangle -= F1_0/4;
295 if (rangle >= F1_0/8) rangle -= F1_0/4;
296 if (rangle <= -F1_0/8) rangle += F1_0/4;
297 if (rangle <= -F1_0/8) rangle += F1_0/4;
299 if ((uangle != 0) && (rangle != 0)) {
300 maxscale = CurveNumSegs*F1_0;
301 // mprintf(0, "Banked Curve Generation.. %f.\n", f2fl(maxscale));
302 generate_banked_curve(maxscale, coeffs);
306 med_form_bridge_segment( Cursegp, Side_opposite[AttachSide], Markedsegp, Markedside );
307 CurveSegs[CurveNumSegs] = &Segments[ Markedsegp->children[Markedside] ];
311 Cursegp = OriginalSeg;
312 Curside = OriginalSide;
314 med_create_new_segment_from_cursegp();
316 //warn_if_concave_segments();
318 if (CurveNumSegs) return 1;
322 void generate_banked_curve(fix maxscale, vms_equation coeffs) {
323 vms_vector vec_dir, tvec, b4r4t;
324 vms_vector coord,prev_point;
325 fix enddist, nextdist;
327 fixang rangle, uangle, angle, scaled_ang=0;
332 extract_up_vector_from_segment( Cursegp,&b4r4t );
333 uangle = vm_vec_delta_ang( &b4r4t, &r4t, &r4 );
334 if (uangle >= F1_0 * 1/8) uangle -= F1_0 * 1/4;
335 if (uangle >= F1_0 * 1/8) uangle -= F1_0 * 1/4;
336 if (uangle <= -F1_0 * 1/8) uangle += F1_0 * 1/4;
337 if (uangle <= -F1_0 * 1/8) uangle += F1_0 * 1/4;
338 // mprintf(0, "up angle %f\n", f2fl(uangle)*360);
340 extract_right_vector_from_segment( Cursegp,&b4r4t );
341 rangle = vm_vec_delta_ang( &b4r4t, &r4t, &r4 );
342 if (rangle >= F1_0/8) rangle -= F1_0/4;
343 if (rangle >= F1_0/8) rangle -= F1_0/4;
344 if (rangle <= -F1_0/8) rangle += F1_0/4;
345 if (rangle <= -F1_0/8) rangle += F1_0/4;
346 // mprintf(0, "right angle %f\n", f2fl(rangle)*360);
349 if (abs(rangle) < abs(uangle)) angle = rangle;
353 coord = prev_point = p1;
355 #define MAGIC_NUM 0.707*F1_0
358 scaled_ang = fixdiv(angle,fixmul(maxscale,MAGIC_NUM));
359 mprintf((0, "scaled angle = %f\n", f2fl(scaled_ang)));
364 enddist = F1_0; nextdist = 0;
365 while ( enddist > fixmul( nextdist, 1.5*F1_0 )) {
366 vms_matrix rotmat,rotmat2;
372 extract_forward_vector_from_segment(Cursegp, &tvec);
373 nextdist = vm_vec_mag(&tvec); // nextdist := distance to next point
374 t = curve_dist(&coeffs, 3, t, &prev_point, nextdist); // t = argument at which function is forward vector magnitude units away from prev_point (in 3-space, not along curve)
375 coord = evaluate_curve(&coeffs, 3, t); // coord := point about forward vector magnitude units away from prev_point
376 enddist = vm_vec_dist(&coord, &p4); // enddist := distance from current to end point, vec_dir used as a temporary variable
377 //vm_vec_normalize(vm_vec_sub(&vec_dir, &coord, &prev_point));
378 vm_vec_normalized_dir(&vec_dir, &coord, &prev_point);
379 if (!med_attach_segment( Cursegp, &New_segment, Curside, AttachSide )) {
380 med_extract_matrix_from_segment( Cursegp,&rotmat ); // rotmat := matrix describing orientation of Cursegp
381 vm_vec_rotate(&tdest,&vec_dir,&rotmat); // tdest := vec_dir in reference frame of Cursegp
383 vm_vec_ang_2_matrix(&rotmat2,&vec_dir,scaled_ang);
384 // mprintf((0, "[ [%6.2f %6.2f %6.2f]", f2fl(rotmat2.m1), f2fl(rotmat2.m2), f2fl(rotmat2.m3)));
385 // mprintf((0, " [%6.2f %6.2f %6.2f]", f2fl(rotmat2.m4), f2fl(rotmat2.m5), f2fl(rotmat2.m6)));
386 // mprintf((0, " [%6.2f %6.2f %6.2f] ]\n", f2fl(rotmat2.m7), f2fl(rotmat2.m8), f2fl(rotmat2.m9)));
388 med_rotate_segment( Cursegp, &rotmat2 );
390 Curside = Side_opposite[AttachSide];
392 CurveSegs[CurveNumSegs]=Cursegp;
400 void delete_curve() {
403 for (i=0; i<CurveNumSegs; i++) {
404 // mprintf((0, "[%d] %d\n", i, CurveSegs[i]->segnum ));
405 if (CurveSegs[i]->segnum != -1)
406 med_delete_segment(CurveSegs[i]);
408 Markedsegp = OriginalMarkedSeg;
409 Markedside = OriginalMarkedSide;
410 Cursegp = OriginalSeg;
411 Curside = OriginalSide;
412 med_create_new_segment_from_cursegp();
414 // mprintf((0, "Num_segments %d\n", Num_segments));
417 //warn_if_concave_segments();
428 vms_vector test, test2, tvec;
435 printf("Enter p1 (x,y,z): ");
436 scanf("%f %f %f", &x, &y, &z);
437 p1.x = x*F1_0; p1.y = y*F1_0; p1.z = z*F1_0;
438 printf("Enter p4 (x,y,z): ");
439 scanf("%f %f %f", &x, &y, &z);
440 p4.x = x*F1_0; p4.y = y*F1_0; p4.z = z*F1_0;
441 printf("Enter r1 <x,y,z>: ");
442 scanf("%f %f %f", &x, &y, &z);
443 r1.x = x*F1_0; r1.y = y*F1_0; r1.z = z*F1_0;
444 printf("Enter r4 <x,y,z>: ");
445 scanf("%f %f %f", &x, &y, &z);
446 r4.x = x*F1_0; r4.y = y*F1_0; r4.z = z*F1_0;
448 create_curve( &p1, &p4, &r1, &r4, &coeffs );
451 printf("x [%6.3f %6.3f %6.3f %6.3f]\n", f2fl(coeffs.n.x3), f2fl(coeffs.n.x2), f2fl(coeffs.n.x1), f2fl(coeffs.n.x0));
452 printf(" y [%6.3f %6.3f %6.3f %6.3f]\n", f2fl(coeffs.n.y3), f2fl(coeffs.n.y2), f2fl(coeffs.n.y1), f2fl(coeffs.n.y0));
453 printf(" z [%6.3f %6.3f %6.3f %6.3f]\n", f2fl(coeffs.n.z3), f2fl(coeffs.n.z2), f2fl(coeffs.n.z1), f2fl(coeffs.n.z0));
455 printf("\nChecking direction vectors.\n");
457 for (t=0*F1_0;t<1*F1_0;t+=0.1*F1_0) {
458 curve_dir(&coeffs, 3, t, &test);
459 printf(" t = %.3f dir = <%6.3f, %6.3f, %6.3f >\n", f2fl(t), f2fl(test.x), f2fl(test.y), f2fl(test.z) );
462 printf("\nChecking distance function.\n");
463 printf("Enter a distance: ");
466 printf("Enter a (0<t<1) value: ");
470 gr_init(15); // 800x600 mode
471 plot_parametric(&coeffs, 0*F1_0, 1*F1_0, 0.05*F1_0);
473 test = evaluate_curve(&coeffs, 3, t0);
474 t = curve_dist(&coeffs, 3, t0, &test, distance);
475 test2 = evaluate_curve(&coeffs, 3, t);
477 dist = vm_vec_mag(vm_vec_sub(&tvec, &test, &test2));
481 gr_rect( 74+f2fl(test.x), 289-f2fl(test.z), 76+f2fl(test.x), 291-f2fl(test.z) );
482 gr_rect( 74+f2fl(test.x), 559-f2fl(test.y), 76+f2fl(test.x), 561-f2fl(test.y) );
483 gr_rect( 474+f2fl(test.z), 559-f2fl(test.y), 476+f2fl(test.z), 561-f2fl(test.y) );
485 gr_rect( 74+f2fl(test2.x), 289-f2fl(test2.z), 76+f2fl(test2.x), 291-f2fl(test2.z) );
486 gr_rect( 74+f2fl(test2.x), 559-f2fl(test2.y), 76+f2fl(test2.x), 561-f2fl(test2.y) );
487 gr_rect( 474+f2fl(test2.z), 559-f2fl(test2.y), 476+f2fl(test2.z), 561-f2fl(test2.y) );
492 if (key == KEY_ESC) break;
493 else key = key_getch();
499 printf("From t=%.3f to t=1.000, ", f2fl(t0));
500 printf("two points separated by the distance %.3f\n do not exist on this curve.\n", x);
503 printf("\nThe distance between points at:\n");
504 printf(" t0 = %.3f ( %6.3f,%6.3f,%6.3f ) and\n", f2fl(t0), f2fl(test.x), f2fl(test.y), f2fl(test.z));
505 printf(" t = %.3f ( %6.3f,%6.3f,%6.3f ) is:\n", f2fl(t), f2fl(test2.x), f2fl(test2.y), f2fl(test2.z));
506 printf(" expected: %.3f\n", x);
507 printf(" actual : %.3f\n", f2fl(dist) );