WebSVN - shark - Rev 55 - /shark/branches/xen/ports/mesa/src-glu/nurbsutl.c

/* $Id: nurbsutl.c,v 1.1 2003-02-28 11:42:07 pj Exp $ */

/*
* Mesa 3-D graphics library
* Version: 3.3
* Copyright (C) 1995-2000 Brian Paul
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Library General Public
* License as published by the Free Software Foundation; either
* version 2 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Library General Public License for more details.
*
* You should have received a copy of the GNU Library General Public
* License along with this library; if not, write to the Free
* Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/

/*
* NURBS implementation written by Bogdan Sikorski (bogdan@cira.it)
* See README2 for more info.
*/

#ifdef PC_HEADER
#include "all.h"
#else
#include <math.h>
#include <stdlib.h>
#include "gluP.h"
#include "nurbs.h"
#endif

GLenum test_knot(GLint nknots, GLfloat * knot, GLint order)
{
GLsizei i;
GLint knot_mult;
GLfloat tmp_knot;

tmp_knot = knot[0];
knot_mult = 1;
for (i = 1; i < nknots; i++) {
if (knot[i] < tmp_knot)
return GLU_NURBS_ERROR4;
if (fabs(tmp_knot - knot[i]) > EPSILON) {
if (knot_mult > order)
return GLU_NURBS_ERROR5;
knot_mult = 1;
tmp_knot = knot[i];
}
else
++knot_mult;
}
return GLU_NO_ERROR;
}

static int
/* qsort function */
#if defined(WIN32) && !defined(OPENSTEP)
__cdecl
#endif
knot_sort(const void *a, const void *b)
{
GLfloat x, y;

x = *((GLfloat *) a);
y = *((GLfloat *) b);
if (fabs(x - y) < EPSILON)
return 0;
if (x > y)
return 1;
return -1;
}

/* insert into dest knot all values within the valid range from src knot */
/* that do not appear in dest */
void
collect_unified_knot(knot_str_type * dest, knot_str_type * src,
GLfloat maximal_min_knot, GLfloat minimal_max_knot)
{
GLfloat *src_knot, *dest_knot;
GLint src_t_min, src_t_max, dest_t_min, dest_t_max;
GLint src_nknots, dest_nknots;
GLint i, j, k, new_cnt;
GLboolean not_found_flag;

src_knot = src->unified_knot;
dest_knot = dest->unified_knot;
src_t_min = src->t_min;
src_t_max = src->t_max;
dest_t_min = dest->t_min;
dest_t_max = dest->t_max;
src_nknots = src->unified_nknots;
dest_nknots = dest->unified_nknots;

k = new_cnt = dest_nknots;
for (i = src_t_min; i <= src_t_max; i++)
if (src_knot[i] - maximal_min_knot > -EPSILON &&
src_knot[i] - minimal_max_knot < EPSILON) {
not_found_flag = GL_TRUE;
for (j = dest_t_min; j <= dest_t_max; j++)
if (fabs(dest_knot[j] - src_knot[i]) < EPSILON) {
not_found_flag = GL_FALSE;
break;
}
if (not_found_flag) {
/* knot from src is not in dest - add this knot to dest */
dest_knot[k++] = src_knot[i];
++new_cnt;
++(dest->t_max); /* the valid range widens */
++(dest->delta_nknots); /* increment the extra knot value counter */
}
}
dest->unified_nknots = new_cnt;
qsort((void *) dest_knot, (size_t) new_cnt, (size_t) sizeof(GLfloat),
&knot_sort);
}

/* basing on the new common knot range for all attributes set */
/* t_min and t_max values for each knot - they will be used later on */
/* by explode_knot() and calc_new_ctrl_pts */
static void
set_new_t_min_t_max(knot_str_type * geom_knot, knot_str_type * color_knot,
knot_str_type * normal_knot, knot_str_type * texture_knot,
GLfloat maximal_min_knot, GLfloat minimal_max_knot)
{
GLuint t_min = 0, t_max = 0, cnt = 0;

if (minimal_max_knot - maximal_min_knot < EPSILON) {
/* knot common range empty */
geom_knot->t_min = geom_knot->t_max = 0;
color_knot->t_min = color_knot->t_max = 0;
normal_knot->t_min = normal_knot->t_max = 0;
texture_knot->t_min = texture_knot->t_max = 0;
}
else {
if (geom_knot->unified_knot != NULL) {
cnt = geom_knot->unified_nknots;
for (t_min = 0; t_min < cnt; t_min++)
if (fabs((geom_knot->unified_knot)[t_min] - maximal_min_knot) <
EPSILON) break;
for (t_max = cnt - 1; t_max; t_max--)
if (fabs((geom_knot->unified_knot)[t_max] - minimal_max_knot) <
EPSILON) break;
}
else if (geom_knot->nknots) {
cnt = geom_knot->nknots;
for (t_min = 0; t_min < cnt; t_min++)
if (fabs((geom_knot->knot)[t_min] - maximal_min_knot) < EPSILON)
break;
for (t_max = cnt - 1; t_max; t_max--)
if (fabs((geom_knot->knot)[t_max] - minimal_max_knot) < EPSILON)
break;
}
geom_knot->t_min = t_min;
geom_knot->t_max = t_max;
if (color_knot->unified_knot != NULL) {
cnt = color_knot->unified_nknots;
for (t_min = 0; t_min < cnt; t_min++)
if (fabs((color_knot->unified_knot)[t_min] - maximal_min_knot) <
EPSILON) break;
for (t_max = cnt - 1; t_max; t_max--)
if (fabs((color_knot->unified_knot)[t_max] - minimal_max_knot) <
EPSILON) break;
color_knot->t_min = t_min;
color_knot->t_max = t_max;
}
if (normal_knot->unified_knot != NULL) {
cnt = normal_knot->unified_nknots;
for (t_min = 0; t_min < cnt; t_min++)
if (fabs((normal_knot->unified_knot)[t_min] - maximal_min_knot) <
EPSILON) break;
for (t_max = cnt - 1; t_max; t_max--)
if (fabs((normal_knot->unified_knot)[t_max] - minimal_max_knot) <
EPSILON) break;
normal_knot->t_min = t_min;
normal_knot->t_max = t_max;
}
if (texture_knot->unified_knot != NULL) {
cnt = texture_knot->unified_nknots;
for (t_min = 0; t_min < cnt; t_min++)
if (fabs((texture_knot->unified_knot)[t_min] - maximal_min_knot)
< EPSILON)
break;
for (t_max = cnt - 1; t_max; t_max--)
if (fabs((texture_knot->unified_knot)[t_max] - minimal_max_knot)
< EPSILON)
break;
texture_knot->t_min = t_min;
texture_knot->t_max = t_max;
}
}
}

/* modify all knot valid ranges in such a way that all have the same */
/* range, common to all knots */
/* do this by knot insertion */
GLenum
select_knot_working_range(GLUnurbsObj * nobj, knot_str_type * geom_knot,
knot_str_type * color_knot,
knot_str_type * normal_knot,
knot_str_type * texture_knot)
{
GLint max_nknots;
GLfloat maximal_min_knot, minimal_max_knot;
GLint i;

/* find the maximum modified knot length */
max_nknots = geom_knot->nknots;
if (color_knot->unified_knot)
max_nknots += color_knot->nknots;
if (normal_knot->unified_knot)
max_nknots += normal_knot->nknots;
if (texture_knot->unified_knot)
max_nknots += texture_knot->nknots;
maximal_min_knot = (geom_knot->knot)[geom_knot->t_min];
minimal_max_knot = (geom_knot->knot)[geom_knot->t_max];
/* any attirb data ? */
if (max_nknots != geom_knot->nknots) {
/* allocate space for the unified knots */
if ((geom_knot->unified_knot =
(GLfloat *) malloc(sizeof(GLfloat) * max_nknots)) == NULL) {
call_user_error(nobj, GLU_OUT_OF_MEMORY);
return GLU_ERROR;
}
/* copy the original knot to the unified one */
geom_knot->unified_nknots = geom_knot->nknots;
for (i = 0; i < geom_knot->nknots; i++)
(geom_knot->unified_knot)[i] = (geom_knot->knot)[i];
if (color_knot->unified_knot) {
if ((color_knot->knot)[color_knot->t_min] - maximal_min_knot >
EPSILON)
maximal_min_knot = (color_knot->knot)[color_knot->t_min];
if (minimal_max_knot - (color_knot->knot)[color_knot->t_max] >
EPSILON)
minimal_max_knot = (color_knot->knot)[color_knot->t_max];
if ((color_knot->unified_knot =
(GLfloat *) malloc(sizeof(GLfloat) * max_nknots)) == NULL) {
free(geom_knot->unified_knot);
call_user_error(nobj, GLU_OUT_OF_MEMORY);
return GLU_ERROR;
}
/* copy the original knot to the unified one */
color_knot->unified_nknots = color_knot->nknots;
for (i = 0; i < color_knot->nknots; i++)
(color_knot->unified_knot)[i] = (color_knot->knot)[i];
}
if (normal_knot->unified_knot) {
if ((normal_knot->knot)[normal_knot->t_min] - maximal_min_knot >
EPSILON)
maximal_min_knot = (normal_knot->knot)[normal_knot->t_min];
if (minimal_max_knot - (normal_knot->knot)[normal_knot->t_max] >
EPSILON)
minimal_max_knot = (normal_knot->knot)[normal_knot->t_max];
if ((normal_knot->unified_knot =
(GLfloat *) malloc(sizeof(GLfloat) * max_nknots)) == NULL) {
free(geom_knot->unified_knot);
free(color_knot->unified_knot);
call_user_error(nobj, GLU_OUT_OF_MEMORY);
return GLU_ERROR;
}
/* copy the original knot to the unified one */
normal_knot->unified_nknots = normal_knot->nknots;
for (i = 0; i < normal_knot->nknots; i++)
(normal_knot->unified_knot)[i] = (normal_knot->knot)[i];
}
if (texture_knot->unified_knot) {
if ((texture_knot->knot)[texture_knot->t_min] - maximal_min_knot >
EPSILON)
maximal_min_knot = (texture_knot->knot)[texture_knot->t_min];
if (minimal_max_knot - (texture_knot->knot)[texture_knot->t_max] >
EPSILON)
minimal_max_knot = (texture_knot->knot)[texture_knot->t_max];
if ((texture_knot->unified_knot =
(GLfloat *) malloc(sizeof(GLfloat) * max_nknots)) == NULL) {
free(geom_knot->unified_knot);
free(color_knot->unified_knot);
free(normal_knot->unified_knot);
call_user_error(nobj, GLU_OUT_OF_MEMORY);
return GLU_ERROR;
}
/* copy the original knot to the unified one */
texture_knot->unified_nknots = texture_knot->nknots;
for (i = 0; i < texture_knot->nknots; i++)
(texture_knot->unified_knot)[i] = (texture_knot->knot)[i];
}
/* work on the geometry knot with all additional knot values */
/* appearing in attirbutive knots */
if (minimal_max_knot - maximal_min_knot < EPSILON) {
/* empty working range */
geom_knot->unified_nknots = 0;
color_knot->unified_nknots = 0;
normal_knot->unified_nknots = 0;
texture_knot->unified_nknots = 0;
}
else {
if (color_knot->unified_knot)
collect_unified_knot(geom_knot, color_knot, maximal_min_knot,
minimal_max_knot);
if (normal_knot->unified_knot)
collect_unified_knot(geom_knot, normal_knot, maximal_min_knot,
minimal_max_knot);
if (texture_knot->unified_knot)
collect_unified_knot(geom_knot, texture_knot, maximal_min_knot,
minimal_max_knot);
/* since we have now built the "unified" geometry knot */
/* add same knot values to all attributive knots */
if (color_knot->unified_knot)
collect_unified_knot(color_knot, geom_knot, maximal_min_knot,
minimal_max_knot);
if (normal_knot->unified_knot)
collect_unified_knot(normal_knot, geom_knot, maximal_min_knot,
minimal_max_knot);
if (texture_knot->unified_knot)
collect_unified_knot(texture_knot, geom_knot, maximal_min_knot,
minimal_max_knot);
}
}
set_new_t_min_t_max(geom_knot, color_knot, normal_knot, texture_knot,
maximal_min_knot, minimal_max_knot);
return GLU_NO_ERROR;
}

void
free_unified_knots(knot_str_type * geom_knot, knot_str_type * color_knot,
knot_str_type * normal_knot, knot_str_type * texture_knot)
{
if (geom_knot->unified_knot)
free(geom_knot->unified_knot);
if (color_knot->unified_knot)
free(color_knot->unified_knot);
if (normal_knot->unified_knot)
free(normal_knot->unified_knot);
if (texture_knot->unified_knot)
free(texture_knot->unified_knot);
}

GLenum explode_knot(knot_str_type * the_knot)
{
GLfloat *knot, *new_knot;
GLint nknots, n_new_knots = 0;
GLint t_min, t_max;
GLint ord;
GLsizei i, j, k;
GLfloat tmp_float;

if (the_knot->unified_knot) {
knot = the_knot->unified_knot;
nknots = the_knot->unified_nknots;
}
else {
knot = the_knot->knot;
nknots = the_knot->nknots;
}
ord = the_knot->order;
t_min = the_knot->t_min;
t_max = the_knot->t_max;

for (i = t_min; i <= t_max;) {
tmp_float = knot[i];
for (j = 0; j < ord && (i + j) <= t_max; j++)
if (fabs(tmp_float - knot[i + j]) > EPSILON)
break;
n_new_knots += ord - j;
i += j;
}
/* alloc space for new_knot */
if (
(new_knot =
(GLfloat *) malloc(sizeof(GLfloat) * (nknots + n_new_knots + 1))) == NULL) {
return GLU_OUT_OF_MEMORY;
}
/* fill in new knot */
for (j = 0; j < t_min; j++)
new_knot[j] = knot[j];
for (i = j; i <= t_max; i++) {
tmp_float = knot[i];
for (k = 0; k < ord; k++) {
new_knot[j++] = knot[i];
if (tmp_float == knot[i + 1])
i++;
}
}
for (i = t_max + 1; i < (int) nknots; i++)
new_knot[j++] = knot[i];
/* fill in the knot structure */
the_knot->new_knot = new_knot;
the_knot->delta_nknots += n_new_knots;
the_knot->t_max += n_new_knots;
return GLU_NO_ERROR;
}

GLenum calc_alphas(knot_str_type * the_knot)
{
GLfloat tmp_float;
int i, j, k, m, n;
int order;
GLfloat *alpha, *alpha_new, *tmp_alpha;
GLfloat denom;
GLfloat *knot, *new_knot;

knot = the_knot->knot;
order = the_knot->order;
new_knot = the_knot->new_knot;
n = the_knot->nknots - the_knot->order;
m = n + the_knot->delta_nknots;
if ((alpha = (GLfloat *) malloc(sizeof(GLfloat) * n * m)) == NULL) {
return GLU_OUT_OF_MEMORY;
}
if ((alpha_new = (GLfloat *) malloc(sizeof(GLfloat) * n * m)) == NULL) {
free(alpha);
return GLU_OUT_OF_MEMORY;
}
for (j = 0; j < m; j++) {
for (i = 0; i < n; i++) {
if ((knot[i] <= new_knot[j]) && (new_knot[j] < knot[i + 1]))
tmp_float = 1.0;
else
tmp_float = 0.0;
alpha[i + j * n] = tmp_float;
}
}
for (k = 1; k < order; k++) {
for (j = 0; j < m; j++)
for (i = 0; i < n; i++) {
denom = knot[i + k] - knot[i];
if (fabs(denom) < EPSILON)
tmp_float = 0.0;
else
tmp_float = (new_knot[j + k] - knot[i]) / denom *
alpha[i + j * n];
denom = knot[i + k + 1] - knot[i + 1];
if (fabs(denom) > EPSILON)
tmp_float += (knot[i + k + 1] - new_knot[j + k]) / denom *
alpha[(i + 1) + j * n];
alpha_new[i + j * n] = tmp_float;
}
tmp_alpha = alpha_new;
alpha_new = alpha;
alpha = tmp_alpha;
}
the_knot->alpha = alpha;
free(alpha_new);
return GLU_NO_ERROR;
}

GLenum
calc_new_ctrl_pts(GLfloat * ctrl, GLint stride, knot_str_type * the_knot,
GLint dim, GLfloat ** new_ctrl, GLint * ncontrol)
{
GLsizei i, j, k, l, m, n;
GLsizei index1, index2;
GLfloat *alpha;
GLfloat *new_knot;

new_knot = the_knot->new_knot;
n = the_knot->nknots - the_knot->order;
alpha = the_knot->alpha;

m = the_knot->t_max + 1 - the_knot->t_min - the_knot->order;
k = the_knot->t_min;
/* allocate space for new control points */
if ((*new_ctrl = (GLfloat *) malloc(sizeof(GLfloat) * dim * m)) == NULL) {
return GLU_OUT_OF_MEMORY;
}
for (j = 0; j < m; j++) {
for (l = 0; l < dim; l++)
(*new_ctrl)[j * dim + l] = 0.0;
for (i = 0; i < n; i++) {
index1 = i + (j + k) * n;
index2 = i * stride;
for (l = 0; l < dim; l++)
(*new_ctrl)[j * dim + l] += alpha[index1] * ctrl[index2 + l];
}
}
*ncontrol = (GLint) m;
return GLU_NO_ERROR;
}

static GLint
calc_factor(GLfloat * pts, GLint order, GLint indx, GLint stride,
GLfloat tolerance, GLint dim)
{
GLdouble model[16], proj[16];
GLint viewport[4];
GLdouble x, y, z, w, winx1, winy1, winz, winx2, winy2;
GLint i;
GLdouble len, dx, dy;

glGetDoublev(GL_MODELVIEW_MATRIX, model);
glGetDoublev(GL_PROJECTION_MATRIX, proj);
glGetIntegerv(GL_VIEWPORT, viewport);
if (dim == 4) {
w = (GLdouble) pts[indx + 3];
x = (GLdouble) pts[indx] / w;
y = (GLdouble) pts[indx + 1] / w;
z = (GLdouble) pts[indx + 2] / w;
gluProject(x, y, z, model, proj, viewport, &winx1, &winy1, &winz);
len = 0.0;
for (i = 1; i < order; i++) {
w = (GLdouble) pts[indx + i * stride + 3];
x = (GLdouble) pts[indx + i * stride] / w;
y = (GLdouble) pts[indx + i * stride + 1] / w;
z = (GLdouble) pts[indx + i * stride + 2] / w;
if (gluProject
(x, y, z, model, proj, viewport, &winx2, &winy2, &winz)) {
dx = winx2 - winx1;
dy = winy2 - winy1;
len += sqrt(dx * dx + dy * dy);
}
winx1 = winx2;
winy1 = winy2;
}
}
else {
x = (GLdouble) pts[indx];
y = (GLdouble) pts[indx + 1];
if (dim == 2)
z = 0.0;
else
z = (GLdouble) pts[indx + 2];
gluProject(x, y, z, model, proj, viewport, &winx1, &winy1, &winz);
len = 0.0;
for (i = 1; i < order; i++) {
x = (GLdouble) pts[indx + i * stride];
y = (GLdouble) pts[indx + i * stride + 1];
if (dim == 2)
z = 0.0;
else
z = (GLdouble) pts[indx + i * stride + 2];
if (gluProject
(x, y, z, model, proj, viewport, &winx2, &winy2, &winz)) {
dx = winx2 - winx1;
dy = winy2 - winy1;
len += sqrt(dx * dx + dy * dy);
}
winx1 = winx2;
winy1 = winy2;
}
}
len /= tolerance;
return ((GLint) len + 1);
}

/* we can't use the Mesa evaluators - no way to get the point coords */
/* so we use our own Bezier point calculus routines */
/* because I'm lazy, I reuse the ones from eval.c */

static void
bezier_curve(GLfloat * cp, GLfloat * out, GLfloat t,
GLuint dim, GLuint order, GLint offset)
{
GLfloat s, powert;
GLuint i, k, bincoeff;

if (order >= 2) {
bincoeff = order - 1;
s = 1.0 - t;

for (k = 0; k < dim; k++)
out[k] = s * cp[k] + bincoeff * t * cp[offset + k];

for (i = 2, cp += 2 * offset, powert = t * t; i < order;
i++, powert *= t, cp += offset) {
bincoeff *= order - i;
bincoeff /= i;

for (k = 0; k < dim; k++)
out[k] = s * out[k] + bincoeff * powert * cp[k];
}
}
else { /* order=1 -> constant curve */

for (k = 0; k < dim; k++)
out[k] = cp[k];
}
}

static GLint
calc_parametric_factor(GLfloat * pts, GLint order, GLint indx, GLint stride,
GLfloat tolerance, GLint dim)
{
GLdouble model[16], proj[16];
GLint viewport[4];
GLdouble x, y, z, w, x1, y1, z1, x2, y2, z2, x3, y3, z3;
GLint i;
GLint P;
GLfloat bez_pt[4];
GLdouble len = 0.0, tmp, z_med;

P = 2 * (order + 2);
glGetDoublev(GL_MODELVIEW_MATRIX, model);
glGetDoublev(GL_PROJECTION_MATRIX, proj);
glGetIntegerv(GL_VIEWPORT, viewport);
z_med = (viewport[2] + viewport[3]) * 0.5;
switch (dim) {
case 4:
for (i = 1; i < P; i++) {
bezier_curve(pts + indx, bez_pt, (GLfloat) i / (GLfloat) P, 4,
order, stride);
w = (GLdouble) bez_pt[3];
x = (GLdouble) bez_pt[0] / w;
y = (GLdouble) bez_pt[1] / w;
z = (GLdouble) bez_pt[2] / w;
gluProject(x, y, z, model, proj, viewport, &x3, &y3, &z3);
z3 *= z_med;
bezier_curve(pts + indx, bez_pt, (GLfloat) (i - 1) / (GLfloat) P, 4,
order, stride);
w = (GLdouble) bez_pt[3];
x = (GLdouble) bez_pt[0] / w;
y = (GLdouble) bez_pt[1] / w;
z = (GLdouble) bez_pt[2] / w;
gluProject(x, y, z, model, proj, viewport, &x1, &y1, &z1);
z1 *= z_med;
bezier_curve(pts + indx, bez_pt, (GLfloat) (i + 1) / (GLfloat) P, 4,
order, stride);
w = (GLdouble) bez_pt[3];
x = (GLdouble) bez_pt[0] / w;
y = (GLdouble) bez_pt[1] / w;
z = (GLdouble) bez_pt[2] / w;
gluProject(x, y, z, model, proj, viewport, &x2, &y2, &z2);
z2 *= z_med;
/* calc distance between point (x3,y3,z3) and line segment */
/* <x1,y1,z1><x2,y2,z2> */
x = x2 - x1;
y = y2 - y1;
z = z2 - z1;
tmp = sqrt(x * x + y * y + z * z);
x /= tmp;
y /= tmp;
z /= tmp;
tmp = x3 * x + y3 * y + z3 * z - x1 * x - y1 * y - z1 * z;
x = x1 + x * tmp - x3;
y = y1 + y * tmp - y3;
z = z1 + z * tmp - z3;
tmp = sqrt(x * x + y * y + z * z);
if (tmp > len)
len = tmp;
}
break;
case 3:
for (i = 1; i < P; i++) {
bezier_curve(pts + indx, bez_pt, (GLfloat) i / (GLfloat) P, 3,
order, stride);
x = (GLdouble) bez_pt[0];
y = (GLdouble) bez_pt[1];
z = (GLdouble) bez_pt[2];
gluProject(x, y, z, model, proj, viewport, &x3, &y3, &z3);
z3 *= z_med;
bezier_curve(pts + indx, bez_pt, (GLfloat) (i - 1) / (GLfloat) P, 3,
order, stride);
x = (GLdouble) bez_pt[0];
y = (GLdouble) bez_pt[1];
z = (GLdouble) bez_pt[2];
gluProject(x, y, z, model, proj, viewport, &x1, &y1, &z1);
z1 *= z_med;
bezier_curve(pts + indx, bez_pt, (GLfloat) (i + 1) / (GLfloat) P, 3,
order, stride);
x = (GLdouble) bez_pt[0];
y = (GLdouble) bez_pt[1];
z = (GLdouble) bez_pt[2];
gluProject(x, y, z, model, proj, viewport, &x2, &y2, &z2);
z2 *= z_med;
/* calc distance between point (x3,y3,z3) and line segment */
/* <x1,y1,z1><x2,y2,z2> */
x = x2 - x1;
y = y2 - y1;
z = z2 - z1;
tmp = sqrt(x * x + y * y + z * z);
x /= tmp;
y /= tmp;
z /= tmp;
tmp = x3 * x + y3 * y + z3 * z - x1 * x - y1 * y - z1 * z;
x = x1 + x * tmp - x3;
y = y1 + y * tmp - y3;
z = z1 + z * tmp - z3;
tmp = sqrt(x * x + y * y + z * z);
if (tmp > len)
len = tmp;
}
break;
case 2:
for (i = 1; i < P; i++) {
bezier_curve(pts + indx, bez_pt, (GLfloat) i / (GLfloat) P, 2,
order, stride);
x = (GLdouble) bez_pt[0];
y = (GLdouble) bez_pt[1];
z = 0.0;
gluProject(x, y, z, model, proj, viewport, &x3, &y3, &z3);
z3 *= z_med;
bezier_curve(pts + indx, bez_pt, (GLfloat) (i - 1) / (GLfloat) P, 2,
order, stride);
x = (GLdouble) bez_pt[0];
y = (GLdouble) bez_pt[1];
z = 0.0;
gluProject(x, y, z, model, proj, viewport, &x1, &y1, &z1);
z1 *= z_med;
bezier_curve(pts + indx, bez_pt, (GLfloat) (i + 1) / (GLfloat) P, 2,
order, stride);
x = (GLdouble) bez_pt[0];
y = (GLdouble) bez_pt[1];
z = 0.0;
gluProject(x, y, z, model, proj, viewport, &x2, &y2, &z2);
z2 *= z_med;
/* calc distance between point (x3,y3,z3) and line segment */
/* <x1,y1,z1><x2,y2,z2> */
x = x2 - x1;
y = y2 - y1;
z = z2 - z1;
tmp = sqrt(x * x + y * y + z * z);
x /= tmp;
y /= tmp;
z /= tmp;
tmp = x3 * x + y3 * y + z3 * z - x1 * x - y1 * y - z1 * z;
x = x1 + x * tmp - x3;
y = y1 + y * tmp - y3;
z = z1 + z * tmp - z3;
tmp = sqrt(x * x + y * y + z * z);
if (tmp > len)
len = tmp;
}
break;

}
if (len < tolerance)
return (order);
else
return (GLint) (sqrt(len / tolerance) * (order + 2) + 1);
}

static GLenum
calc_sampling_3D(new_ctrl_type * new_ctrl, GLfloat tolerance, GLint dim,
GLint uorder, GLint vorder, GLint ** ufactors,
GLint ** vfactors)
{
GLfloat *ctrl;
GLint tmp_factor1, tmp_factor2;
GLint ufactor_cnt, vfactor_cnt;
GLint offset1, offset2, offset3;
GLint i, j;

ufactor_cnt = new_ctrl->s_bezier_cnt;
vfactor_cnt = new_ctrl->t_bezier_cnt;
if ((*ufactors = (GLint *) malloc(sizeof(GLint) * ufactor_cnt * 3))
== NULL) {
return GLU_OUT_OF_MEMORY;
}
if ((*vfactors = (GLint *) malloc(sizeof(GLint) * vfactor_cnt * 3))
== NULL) {
free(*ufactors);
return GLU_OUT_OF_MEMORY;
}
ctrl = new_ctrl->geom_ctrl;
offset1 = new_ctrl->geom_t_stride * vorder;
offset2 = new_ctrl->geom_s_stride * uorder;
for (j = 0; j < vfactor_cnt; j++) {
*(*vfactors + j * 3 + 1) = tmp_factor1 = calc_factor(ctrl, vorder,
j * offset1, dim,
tolerance, dim);
/* loop ufactor_cnt-1 times */
for (i = 1; i < ufactor_cnt; i++) {
tmp_factor2 = calc_factor(ctrl, vorder,
j * offset1 + i * offset2, dim, tolerance,
dim);
if (tmp_factor2 > tmp_factor1)
tmp_factor1 = tmp_factor2;
}
/* last time for the opposite edge */
*(*vfactors + j * 3 + 2) = tmp_factor2 = calc_factor(ctrl, vorder,
j * offset1 +
i * offset2 -
new_ctrl->
geom_s_stride, dim,
tolerance, dim);
if (tmp_factor2 > tmp_factor1)
*(*vfactors + j * 3) = tmp_factor2;
else
*(*vfactors + j * 3) = tmp_factor1;
}
offset3 = new_ctrl->geom_s_stride;
offset2 = new_ctrl->geom_s_stride * uorder;
for (j = 0; j < ufactor_cnt; j++) {
*(*ufactors + j * 3 + 1) = tmp_factor1 = calc_factor(ctrl, uorder,
j * offset2,
offset3, tolerance,
dim);
/* loop vfactor_cnt-1 times */
for (i = 1; i < vfactor_cnt; i++) {
tmp_factor2 = calc_factor(ctrl, uorder,
j * offset2 + i * offset1, offset3,
tolerance, dim);
if (tmp_factor2 > tmp_factor1)
tmp_factor1 = tmp_factor2;
}
/* last time for the opposite edge */
*(*ufactors + j * 3 + 2) = tmp_factor2 = calc_factor(ctrl, uorder,
j * offset2 +
i * offset1 -
new_ctrl->
geom_t_stride,
offset3, tolerance,
dim);
if (tmp_factor2 > tmp_factor1)
*(*ufactors + j * 3) = tmp_factor2;
else
*(*ufactors + j * 3) = tmp_factor1;
}
return GL_NO_ERROR;
}

static GLenum
calc_sampling_param_3D(new_ctrl_type * new_ctrl, GLfloat tolerance, GLint dim,
GLint uorder, GLint vorder, GLint ** ufactors,
GLint ** vfactors)
{
GLfloat *ctrl;
GLint tmp_factor1, tmp_factor2;
GLint ufactor_cnt, vfactor_cnt;
GLint offset1, offset2, offset3;
GLint i, j;

ufactor_cnt = new_ctrl->s_bezier_cnt;
vfactor_cnt = new_ctrl->t_bezier_cnt;
if ((*ufactors = (GLint *) malloc(sizeof(GLint) * ufactor_cnt * 3))
== NULL) {
return GLU_OUT_OF_MEMORY;
}
if ((*vfactors = (GLint *) malloc(sizeof(GLint) * vfactor_cnt * 3))
== NULL) {
free(*ufactors);
return GLU_OUT_OF_MEMORY;
}
ctrl = new_ctrl->geom_ctrl;
offset1 = new_ctrl->geom_t_stride * vorder;
offset2 = new_ctrl->geom_s_stride * uorder;
for (j = 0; j < vfactor_cnt; j++) {
*(*vfactors + j * 3 + 1) = tmp_factor1 =
calc_parametric_factor(ctrl, vorder, j * offset1, dim, tolerance,
dim);
/* loop ufactor_cnt-1 times */
for (i = 1; i < ufactor_cnt; i++) {
tmp_factor2 = calc_parametric_factor(ctrl, vorder,
j * offset1 + i * offset2, dim,
tolerance, dim);
if (tmp_factor2 > tmp_factor1)
tmp_factor1 = tmp_factor2;
}
/* last time for the opposite edge */
*(*vfactors + j * 3 + 2) = tmp_factor2 =
calc_parametric_factor(ctrl, vorder,
j * offset1 + i * offset2 -
new_ctrl->geom_s_stride, dim, tolerance, dim);
if (tmp_factor2 > tmp_factor1)
*(*vfactors + j * 3) = tmp_factor2;
else
*(*vfactors + j * 3) = tmp_factor1;
}
offset3 = new_ctrl->geom_s_stride;
offset2 = new_ctrl->geom_s_stride * uorder;
for (j = 0; j < ufactor_cnt; j++) {
*(*ufactors + j * 3 + 1) = tmp_factor1 =
calc_parametric_factor(ctrl, uorder, j * offset2, offset3, tolerance,
dim);
/* loop vfactor_cnt-1 times */
for (i = 1; i < vfactor_cnt; i++) {
tmp_factor2 = calc_parametric_factor(ctrl, uorder,
j * offset2 + i * offset1,
offset3, tolerance, dim);
if (tmp_factor2 > tmp_factor1)
tmp_factor1 = tmp_factor2;
}
/* last time for the opposite edge */
*(*ufactors + j * 3 + 2) = tmp_factor2 =
calc_parametric_factor(ctrl, uorder,
j * offset2 + i * offset1 -
new_ctrl->geom_t_stride, offset3, tolerance,
dim);
if (tmp_factor2 > tmp_factor1)
*(*ufactors + j * 3) = tmp_factor2;
else
*(*ufactors + j * 3) = tmp_factor1;
}
return GL_NO_ERROR;
}

static GLenum
calc_sampling_2D(GLfloat * ctrl, GLint cnt, GLint order,
GLfloat tolerance, GLint dim, GLint ** factors)
{
GLint factor_cnt;
GLint tmp_factor;
GLint offset;
GLint i;

factor_cnt = cnt / order;
if ((*factors = (GLint *) malloc(sizeof(GLint) * factor_cnt)) == NULL) {
return GLU_OUT_OF_MEMORY;
}
offset = order * dim;
for (i = 0; i < factor_cnt; i++) {
tmp_factor = calc_factor(ctrl, order, i * offset, dim, tolerance, dim);
if (tmp_factor == 0)
(*factors)[i] = 1;
else
(*factors)[i] = tmp_factor;
}
return GL_NO_ERROR;
}

static void
set_sampling_and_culling(GLUnurbsObj * nobj)
{
if (nobj->auto_load_matrix == GL_FALSE) {
GLint i;
GLfloat m[4];

glPushAttrib((GLbitfield) (GL_VIEWPORT_BIT | GL_TRANSFORM_BIT));
for (i = 0; i < 4; i++)
m[i] = nobj->sampling_matrices.viewport[i];
glViewport(m[0], m[1], m[2], m[3]);
glMatrixMode(GL_PROJECTION);
glPushMatrix();
glLoadMatrixf(nobj->sampling_matrices.proj);
glMatrixMode(GL_MODELVIEW);
glPushMatrix();
glLoadMatrixf(nobj->sampling_matrices.model);
}
}

static void
revert_sampling_and_culling(GLUnurbsObj * nobj)
{
if (nobj->auto_load_matrix == GL_FALSE) {
glMatrixMode(GL_MODELVIEW);
glPopMatrix();
glMatrixMode(GL_PROJECTION);
glPopMatrix();
glPopAttrib();
}
}

GLenum
glu_do_sampling_3D(GLUnurbsObj * nobj, new_ctrl_type * new_ctrl,
GLint ** sfactors, GLint ** tfactors)
{
GLint dim;
GLenum err;

*sfactors = NULL;
*tfactors = NULL;
dim = nobj->surface.geom.dim;
set_sampling_and_culling(nobj);
if ((err = calc_sampling_3D(new_ctrl, nobj->sampling_tolerance, dim,
nobj->surface.geom.sorder,
nobj->surface.geom.torder, sfactors,
tfactors)) == GLU_ERROR) {
revert_sampling_and_culling(nobj);
call_user_error(nobj, err);
return GLU_ERROR;
}
revert_sampling_and_culling(nobj);
return GLU_NO_ERROR;
}

GLenum
glu_do_sampling_uv(GLUnurbsObj * nobj, new_ctrl_type * new_ctrl,
GLint ** sfactors, GLint ** tfactors)
{
GLint s_cnt, t_cnt, i;
GLint u_steps, v_steps;

s_cnt = new_ctrl->s_bezier_cnt;
t_cnt = new_ctrl->t_bezier_cnt;
*sfactors = NULL;
*tfactors = NULL;
if ((*sfactors = (GLint *) malloc(sizeof(GLint) * s_cnt * 3))
== NULL) {
return GLU_OUT_OF_MEMORY;
}
if ((*tfactors = (GLint *) malloc(sizeof(GLint) * t_cnt * 3))
== NULL) {
free(*sfactors);
return GLU_OUT_OF_MEMORY;
}
u_steps = nobj->u_step;
v_steps = nobj->v_step;
for (i = 0; i < s_cnt; i++) {
*(*sfactors + i * 3) = u_steps;
*(*sfactors + i * 3 + 1) = u_steps;
*(*sfactors + i * 3 + 2) = u_steps;
}
for (i = 0; i < t_cnt; i++) {
*(*tfactors + i * 3) = v_steps;
*(*tfactors + i * 3 + 1) = v_steps;
*(*tfactors + i * 3 + 2) = v_steps;
}
return GLU_NO_ERROR;
}

GLenum
glu_do_sampling_param_3D(GLUnurbsObj * nobj, new_ctrl_type * new_ctrl,
GLint ** sfactors, GLint ** tfactors)
{
GLint dim;
GLenum err;

*sfactors = NULL;
*tfactors = NULL;
dim = nobj->surface.geom.dim;
set_sampling_and_culling(nobj);
if (
(err =
calc_sampling_param_3D(new_ctrl, nobj->parametric_tolerance, dim,
nobj->surface.geom.sorder,
nobj->surface.geom.torder, sfactors,
tfactors)) == GLU_ERROR) {
revert_sampling_and_culling(nobj);
call_user_error(nobj, err);
return GLU_ERROR;
}
revert_sampling_and_culling(nobj);
return GLU_NO_ERROR;
}

static GLenum
glu_do_sampling_2D(GLUnurbsObj * nobj, GLfloat * ctrl, GLint cnt, GLint order,
GLint dim, GLint ** factors)
{
GLenum err;

set_sampling_and_culling(nobj);
err = calc_sampling_2D(ctrl, cnt, order, nobj->sampling_tolerance, dim,
factors);
revert_sampling_and_culling(nobj);
return err;
}

static GLenum
glu_do_sampling_u(GLUnurbsObj * nobj, GLfloat * ctrl, GLint cnt, GLint order,
GLint dim, GLint ** factors)
{
GLint i;
GLint u_steps;

cnt /= order;
if ((*factors = (GLint *) malloc(sizeof(GLint) * cnt))
== NULL) {
return GLU_OUT_OF_MEMORY;
}
u_steps = nobj->u_step;
for (i = 0; i < cnt; i++)
(*factors)[i] = u_steps;
return GLU_NO_ERROR;
}

static GLenum
glu_do_sampling_param_2D(GLUnurbsObj * nobj, GLfloat * ctrl, GLint cnt,
GLint order, GLint dim, GLint ** factors)
{
GLint i;
GLint u_steps;
GLfloat tolerance;

set_sampling_and_culling(nobj);
tolerance = nobj->parametric_tolerance;
cnt /= order;
if ((*factors = (GLint *) malloc(sizeof(GLint) * cnt))
== NULL) {
revert_sampling_and_culling(nobj);
return GLU_OUT_OF_MEMORY;
}
u_steps = nobj->u_step;
for (i = 0; i < cnt; i++) {
(*factors)[i] = calc_parametric_factor(ctrl, order, 0,
dim, tolerance, dim);

}
revert_sampling_and_culling(nobj);
return GLU_NO_ERROR;
}

GLenum
glu_do_sampling_crv(GLUnurbsObj * nobj, GLfloat * ctrl, GLint cnt,
GLint order, GLint dim, GLint ** factors)
{
GLenum err;

*factors = NULL;
switch (nobj->sampling_method) {
case GLU_PATH_LENGTH:
if ((err = glu_do_sampling_2D(nobj, ctrl, cnt, order, dim, factors)) !=
GLU_NO_ERROR) {
call_user_error(nobj, err);
return GLU_ERROR;
}
break;
case GLU_DOMAIN_DISTANCE:
if ((err = glu_do_sampling_u(nobj, ctrl, cnt, order, dim, factors)) !=
GLU_NO_ERROR) {
call_user_error(nobj, err);
return GLU_ERROR;
}
break;
case GLU_PARAMETRIC_ERROR:
if (
(err =
glu_do_sampling_param_2D(nobj, ctrl, cnt, order, dim,
factors)) != GLU_NO_ERROR) {
call_user_error(nobj, err);
return GLU_ERROR;
}
break;
default:
abort();
}

return GLU_NO_ERROR;
}

/* TODO - i don't like this culling - this one just tests if at least one */
/* ctrl point lies within the viewport . Also the point_in_viewport() */
/* should be included in the fnctions for efficiency reasons */

static GLboolean
point_in_viewport(GLfloat * pt, GLint dim)
{
GLdouble model[16], proj[16];
GLint viewport[4];
GLdouble x, y, z, w, winx, winy, winz;

glGetDoublev(GL_MODELVIEW_MATRIX, model);
glGetDoublev(GL_PROJECTION_MATRIX, proj);
glGetIntegerv(GL_VIEWPORT, viewport);
if (dim == 3) {
x = (GLdouble) pt[0];
y = (GLdouble) pt[1];
z = (GLdouble) pt[2];
gluProject(x, y, z, model, proj, viewport, &winx, &winy, &winz);
}
else {
w = (GLdouble) pt[3];
x = (GLdouble) pt[0] / w;
y = (GLdouble) pt[1] / w;
z = (GLdouble) pt[2] / w;
gluProject(x, y, z, model, proj, viewport, &winx, &winy, &winz);
}
if ((GLint) winx >= viewport[0] && (GLint) winx < viewport[2] &&
(GLint) winy >= viewport[1] && (GLint) winy < viewport[3])
return GL_TRUE;
return GL_FALSE;
}

GLboolean
fine_culling_test_3D(GLUnurbsObj * nobj, GLfloat * pts, GLint s_cnt,
GLint t_cnt, GLint s_stride, GLint t_stride, GLint dim)
{
GLint i, j;

if (nobj->culling == GL_FALSE)
return GL_FALSE;
set_sampling_and_culling(nobj);

if (dim == 3) {
for (i = 0; i < s_cnt; i++)
for (j = 0; j < t_cnt; j++)
if (point_in_viewport(pts + i * s_stride + j * t_stride, dim)) {
revert_sampling_and_culling(nobj);
return GL_FALSE;
}
}
else {
for (i = 0; i < s_cnt; i++)
for (j = 0; j < t_cnt; j++)
if (point_in_viewport(pts + i * s_stride + j * t_stride, dim)) {
revert_sampling_and_culling(nobj);
return GL_FALSE;
}
}
revert_sampling_and_culling(nobj);
return GL_TRUE;
}

/*GLboolean
fine_culling_test_3D(GLUnurbsObj *nobj,GLfloat *pts,GLint s_cnt,GLint t_cnt,
GLint s_stride,GLint t_stride, GLint dim)
{
GLint visible_cnt;
GLfloat feedback_buffer[5];
GLsizei buffer_size;
GLint i,j;

if(nobj->culling==GL_FALSE)
return GL_FALSE;
buffer_size=5;
set_sampling_and_culling(nobj);

glFeedbackBuffer(buffer_size,GL_2D,feedback_buffer);
glRenderMode(GL_FEEDBACK);
if(dim==3)
{
for(i=0;i<s_cnt;i++)
{
glBegin(GL_LINE_LOOP);
for(j=0;j<t_cnt;j++)
glVertex3fv(pts+i*s_stride+j*t_stride);
glEnd();
}
for(j=0;j<t_cnt;j++)
{
glBegin(GL_LINE_LOOP);
for(i=0;i<s_cnt;i++)
glVertex3fv(pts+i*s_stride+j*t_stride);
glEnd();
}
}
else
{
for(i=0;i<s_cnt;i++)
{
glBegin(GL_LINE_LOOP);
for(j=0;j<t_cnt;j++)
glVertex4fv(pts+i*s_stride+j*t_stride);
glEnd();
}
for(j=0;j<t_cnt;j++)
{
glBegin(GL_LINE_LOOP);
for(i=0;i<s_cnt;i++)
glVertex4fv(pts+i*s_stride+j*t_stride);
glEnd();
}
}
visible_cnt=glRenderMode(GL_RENDER);

revert_sampling_and_culling(nobj);
return (GLboolean)(visible_cnt==0);
}*/

GLboolean
fine_culling_test_2D(GLUnurbsObj * nobj, GLfloat * pts, GLint cnt,
GLint stride, GLint dim)
{
GLint i;

if (nobj->culling == GL_FALSE)
return GL_FALSE;
set_sampling_and_culling(nobj);

if (dim == 3) {
for (i = 0; i < cnt; i++)
if (point_in_viewport(pts + i * stride, dim)) {
revert_sampling_and_culling(nobj);
return GL_FALSE;
}
}
else {
for (i = 0; i < cnt; i++)
if (point_in_viewport(pts + i * stride, dim)) {
revert_sampling_and_culling(nobj);
return GL_FALSE;
}
}
revert_sampling_and_culling(nobj);
return GL_TRUE;
}

/*GLboolean
fine_culling_test_2D(GLUnurbsObj *nobj,GLfloat *pts,GLint cnt,
GLint stride, GLint dim)
{
GLint visible_cnt;
GLfloat feedback_buffer[5];
GLsizei buffer_size;
GLint i;

if(nobj->culling==GL_FALSE)
return GL_FALSE;
buffer_size=5;
set_sampling_and_culling(nobj);

glFeedbackBuffer(buffer_size,GL_2D,feedback_buffer);
glRenderMode(GL_FEEDBACK);
glBegin(GL_LINE_LOOP);
if(dim==3)
{
for(i=0;i<cnt;i++)
glVertex3fv(pts+i*stride);
}
else
{
for(i=0;i<cnt;i++)
glVertex4fv(pts+i*stride);
}
glEnd();
visible_cnt=glRenderMode(GL_RENDER);

revert_sampling_and_culling(nobj);
return (GLboolean)(visible_cnt==0);
}*/

Subversion Repositories shark

(root)/shark/branches/xen/ports/mesa/src-glu/nurbsutl.c - Rev 55