00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "stdafx.h"
00027 #include "shapes/nurbs.h"
00028 #include "shapes/trianglemesh.h"
00029 #include "paramset.h"
00030 #include "texture.h"
00031
00032
00033 static int KnotOffset(const float *knot, int order, int np, float t) {
00034 int firstKnot = order - 1;
00035
00036 int knotOffset = firstKnot;
00037 while (t > knot[knotOffset+1])
00038 ++knotOffset;
00039 Assert(knotOffset < np);
00040 Assert(t >= knot[knotOffset] && t <= knot[knotOffset + 1]);
00041 return knotOffset;
00042 }
00043
00044
00045
00046
00047
00048 struct Homogeneous3 {
00049 Homogeneous3() { x = y = z = w = 0.; }
00050 Homogeneous3(float xx, float yy, float zz, float ww) {
00051 x = xx; y = yy; z = zz; w = ww;
00052 }
00053
00054
00055 float x, y, z, w;
00056 };
00057
00058
00059 static Homogeneous3
00060 NURBSEvaluate(int order, const float *knot, const Homogeneous3 *cp, int np,
00061 int cpStride, float t, Vector *deriv = NULL) {
00062
00063 float alpha;
00064
00065 int knotOffset = KnotOffset(knot, order, np, t);
00066 knot += knotOffset;
00067
00068 int cpOffset = knotOffset - order + 1;
00069 Assert(cpOffset >= 0 && cpOffset < np);
00070
00071 Homogeneous3 *cpWork = ALLOCA(Homogeneous3, order);
00072 for (int i = 0; i < order; ++i)
00073 cpWork[i] = cp[(cpOffset+i) * cpStride];
00074
00075 for (int i = 0; i < order - 2; ++i)
00076 for (int j = 0; j < order - 1 - i; ++j) {
00077 alpha = (knot[1 + j] - t) /
00078 (knot[1 + j] - knot[j + 2 - order + i]);
00079 Assert(alpha >= 0. && alpha <= 1.);
00080
00081 cpWork[j].x = cpWork[j].x * alpha + cpWork[j+1].x * (1 - alpha);
00082 cpWork[j].y = cpWork[j].y * alpha + cpWork[j+1].y * (1 - alpha);
00083 cpWork[j].z = cpWork[j].z * alpha + cpWork[j+1].z * (1 - alpha);
00084 cpWork[j].w = cpWork[j].w * alpha + cpWork[j+1].w * (1 - alpha);
00085 }
00086
00087 alpha = (knot[1] - t) / (knot[1] - knot[0]);
00088 Assert(alpha >= 0. && alpha <= 1.);
00089
00090 Homogeneous3 val(cpWork[0].x * alpha + cpWork[1].x * (1 - alpha),
00091 cpWork[0].y * alpha + cpWork[1].y * (1 - alpha),
00092 cpWork[0].z * alpha + cpWork[1].z * (1 - alpha),
00093 cpWork[0].w * alpha + cpWork[1].w * (1 - alpha));
00094
00095 if (deriv) {
00096 float factor = (order - 1) / (knot[1] - knot[0]);
00097 Homogeneous3 delta((cpWork[1].x - cpWork[0].x) * factor,
00098 (cpWork[1].y - cpWork[0].y) * factor,
00099 (cpWork[1].z - cpWork[0].z) * factor,
00100 (cpWork[1].w - cpWork[0].w) * factor);
00101
00102 deriv->x = delta.x / val.w - (val.x * delta.w / (val.w * val.w));
00103 deriv->y = delta.y / val.w - (val.y * delta.w / (val.w * val.w));
00104 deriv->z = delta.z / val.w - (val.z * delta.w / (val.w * val.w));
00105 }
00106
00107 return val;
00108 }
00109
00110
00111
00112 static Point
00113 NURBSEvaluateSurface(int uOrder, const float *uKnot, int ucp, float u,
00114 int vOrder, const float *vKnot, int vcp, float v,
00115 const Homogeneous3 *cp, Vector *dPdu, Vector *dPdv) {
00116 Homogeneous3 *iso = ALLOCA(Homogeneous3, max(uOrder, vOrder));
00117
00118 int uOffset = KnotOffset(uKnot, uOrder, ucp, u);
00119 int uFirstCp = uOffset - uOrder + 1;
00120 Assert(uFirstCp >= 0 && uFirstCp + uOrder - 1 < ucp);
00121
00122 for (int i = 0; i < uOrder; ++i)
00123 iso[i] = NURBSEvaluate(vOrder, vKnot, &cp[uFirstCp + i], vcp,
00124 ucp, v);
00125
00126 int vOffset = KnotOffset(vKnot, vOrder, vcp, v);
00127 int vFirstCp = vOffset - vOrder + 1;
00128 Assert(vFirstCp >= 0 && vFirstCp + vOrder - 1 < vcp);
00129
00130 Homogeneous3 P = NURBSEvaluate(uOrder, uKnot, iso - uFirstCp, ucp,
00131 1, u, dPdu);
00132
00133 if (dPdv) {
00134 for (int i = 0; i < vOrder; ++i)
00135 iso[i] = NURBSEvaluate(uOrder, uKnot, &cp[(vFirstCp+i)*ucp], ucp,
00136 1, u);
00137 (void)NURBSEvaluate(vOrder, vKnot, iso - vFirstCp, vcp, 1, v, dPdv);
00138 }
00139 return Point(P.x/P.w, P.y/P.w, P.z/P.w);;
00140 }
00141
00142
00143
00144
00145
00146 NURBS::NURBS(const Transform *o2w, const Transform *w2o,
00147 bool ro, int numu, int uo, const float *uk,
00148 float u0, float u1, int numv, int vo, const float *vk,
00149 float v0, float v1, const float *p, bool homogeneous)
00150 : Shape(o2w, w2o, ro) {
00151 nu = numu; uorder = uo;
00152 umin = u0; umax = u1;
00153 nv = numv; vorder = vo;
00154 vmin = v0; vmax = v1;
00155 isHomogeneous = homogeneous;
00156 if (isHomogeneous) {
00157 P = new float[4*nu*nv];
00158 memcpy(P, p, 4*nu*nv*sizeof(float));
00159 } else {
00160 P = new float[3*nu*nv];
00161 memcpy(P, p, 3*nu*nv*sizeof(float));
00162 }
00163 uknot = new float[nu + uorder];
00164 memcpy(uknot, uk, (nu + uorder) * sizeof(float));
00165 vknot = new float[nv + vorder];
00166 memcpy(vknot, vk, (nv + vorder) * sizeof(float));
00167 }
00168
00169
00170 NURBS::~NURBS() {
00171 delete[] P;
00172 delete[] uknot;
00173 delete[] vknot;
00174 }
00175
00176
00177 BBox NURBS::ObjectBound() const {
00178 if (!isHomogeneous) {
00179
00180 float *pp = P;
00181 BBox bound;
00182 for (int i = 0; i < nu*nv; ++i, pp += 3)
00183 bound = Union(bound, Point(pp[0], pp[1], pp[2]));
00184 return bound;
00185 } else {
00186
00187 float *pp = P;
00188 BBox bound;
00189 for (int i = 0; i < nu*nv; ++i, pp += 4)
00190 bound = Union(bound, Point(pp[0] / pp[3], pp[1] / pp[3], pp[2] / pp[3]));
00191 return bound;
00192 }
00193 }
00194
00195
00196 BBox NURBS::WorldBound() const {
00197 if (!isHomogeneous) {
00198
00199 float *pp = P;
00200 BBox bound;
00201 for (int i = 0; i < nu*nv; ++i, pp += 3) {
00202 Point pt = (*ObjectToWorld)(Point(pp[0], pp[1], pp[2]));
00203 bound = Union(bound, pt);
00204 }
00205 return bound;
00206 } else {
00207
00208 float *pp = P;
00209 BBox bound;
00210 for (int i = 0; i < nu*nv; ++i, pp += 4) {
00211 Point pt = (*ObjectToWorld)(Point(pp[0]/pp[3],
00212 pp[1]/pp[3], pp[2]/pp[3]));
00213 bound = Union(bound, pt);
00214 }
00215 return bound;
00216 }
00217 }
00218
00219
00220
00221 void NURBS::Refine(vector<Reference<Shape> > &refined) const {
00222
00223 int diceu = 30, dicev = 30;
00224 float *ueval = new float[diceu];
00225 float *veval = new float[dicev];
00226 Point *evalPs = new Point[diceu*dicev];
00227 Normal *evalNs = new Normal[diceu*dicev];
00228 int i;
00229 for (i = 0; i < diceu; ++i)
00230 ueval[i] = Lerp((float)i / (float)(diceu-1), umin, umax);
00231 for (i = 0; i < dicev; ++i)
00232 veval[i] = Lerp((float)i / (float)(dicev-1), vmin, vmax);
00233
00234 memset(evalPs, 0, diceu*dicev*sizeof(Point));
00235 memset(evalNs, 0, diceu*dicev*sizeof(Point));
00236 float *uvs = new float[2*diceu*dicev];
00237
00238 Homogeneous3 *Pw = (Homogeneous3 *)P;
00239 if (!isHomogeneous) {
00240 Pw = new Homogeneous3[nu * nv];
00241 for (int i = 0; i < nu*nv; ++i) {
00242 Pw[i].x = P[3*i];
00243 Pw[i].y = P[3*i+1];
00244 Pw[i].z = P[3*i+2];
00245 Pw[i].w = 1.;
00246 }
00247 }
00248 for (int v = 0; v < dicev; ++v) {
00249 for (int u = 0; u < diceu; ++u) {
00250 uvs[2*(v*diceu+u)] = ueval[u];
00251 uvs[2*(v*diceu+u)+1] = veval[v];
00252
00253 Vector dPdu, dPdv;
00254 Point pt = NURBSEvaluateSurface(uorder, uknot, nu, ueval[u],
00255 vorder, vknot, nv, veval[v], Pw, &dPdu, &dPdv);
00256 evalPs[v*diceu + u].x = pt.x;
00257 evalPs[v*diceu + u].y = pt.y;
00258 evalPs[v*diceu + u].z = pt.z;
00259 evalNs[v*diceu + u] = Normal(Normalize(Cross(dPdu, dPdv)));
00260 }
00261 }
00262
00263 int nTris = 2*(diceu-1)*(dicev-1);
00264 int *vertices = new int[3 * nTris];
00265 int *vertp = vertices;
00266
00267 for (int v = 0; v < dicev-1; ++v) {
00268 for (int u = 0; u < diceu-1; ++u) {
00269 #define VN(u,v) ((v)*diceu+(u))
00270 *vertp++ = VN(u, v);
00271 *vertp++ = VN(u+1, v);
00272 *vertp++ = VN(u+1, v+1);
00273
00274 *vertp++ = VN(u, v);
00275 *vertp++ = VN(u+1, v+1);
00276 *vertp++ = VN(u, v+1);
00277 #undef VN
00278 }
00279 }
00280 int nVerts = diceu*dicev;
00281 ParamSet paramSet;
00282 paramSet.AddInt("indices", vertices, 3*nTris);
00283 paramSet.AddPoint("P", evalPs, nVerts);
00284 paramSet.AddFloat("uv", uvs, 2 * nVerts);
00285 paramSet.AddNormal("N", evalNs, nVerts);
00286 refined.push_back(CreateTriangleMeshShape(ObjectToWorld, WorldToObject,
00287 ReverseOrientation, paramSet));
00288
00289 if (Pw != (Homogeneous3 *)P) delete[] Pw;
00290 delete[] uvs;
00291 delete[] ueval;
00292 delete[] veval;
00293 delete[] evalPs;
00294 delete[] evalNs;
00295 delete[] vertices;
00296 }
00297
00298
00299
00300 NURBS *CreateNURBSShape(const Transform *o2w, const Transform *w2o,
00301 bool ReverseOrientation, const ParamSet ¶ms) {
00302 int nu = params.FindOneInt("nu", -1);
00303 int uorder = params.FindOneInt("uorder", -1);
00304 int nuknots, nvknots;
00305 const float *uknots = params.FindFloat("uknots", &nuknots);
00306 Assert(nu != -1 && uorder != -1 && uknots != NULL);
00307 Assert(nuknots == nu + uorder);
00308 float u0 = params.FindOneFloat("u0", uknots[uorder-1]);
00309 float u1 = params.FindOneFloat("u1", uknots[nu]);
00310
00311 int nv = params.FindOneInt("nv", -1);
00312 int vorder = params.FindOneInt("vorder", -1);
00313 const float *vknots = params.FindFloat("vknots", &nvknots);
00314 Assert(nv != -1 && vorder != -1 && vknots != NULL);
00315 Assert(nvknots == nv + vorder);
00316 float v0 = params.FindOneFloat("v0", vknots[vorder-1]);
00317 float v1 = params.FindOneFloat("v1", vknots[nv]);
00318
00319 bool isHomogeneous = false;
00320 int npts;
00321 const float *P = (const float *)params.FindPoint("P", &npts);
00322 if (!P) {
00323 P = params.FindFloat("Pw", &npts);
00324 npts /= 4;
00325 if (!P) return NULL;
00326 isHomogeneous = true;
00327 }
00328 Assert(P);
00329 Assert(npts == nu*nv);
00330 return new NURBS(o2w, w2o, ReverseOrientation, nu, uorder, uknots, u0, u1,
00331 nv, vorder, vknots, v0, v1, (float *)P,
00332 isHomogeneous);
00333 }
00334
00335