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