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