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 "shape.h"
00026 #include "paramset.h"
00027 #include "dynload.h"
00028 #include "texture.h"
00029 #include <set>
00030 #include <map>
00031 using std::set;
00032 using std::map;
00033
00034 #define NEXT(i) (((i)+1)%3)
00035 #define PREV(i) (((i)+2)%3)
00036
00037 struct SDFace;
00038 struct SDFace;
00039 struct SDVertex {
00040
00041 SDVertex(Point pt = Point(0,0,0))
00042 : P(pt), startFace(NULL), child(NULL),
00043 regular(false), boundary(false) {
00044 }
00045
00046 int valence();
00047 void oneRing(Point *P);
00048 Point P;
00049 SDFace *startFace;
00050 SDVertex *child;
00051 bool regular, boundary;
00052 };
00053 struct SDFace {
00054
00055 SDFace() {
00056 int i;
00057 for (i = 0; i < 3; ++i) {
00058 v[i] = NULL;
00059 f[i] = NULL;
00060 }
00061 for (i = 0; i < 4; ++i)
00062 children[i] = NULL;
00063 }
00064
00065 int vnum(SDVertex *vert) const {
00066 for (int i = 0; i < 3; ++i)
00067 if (v[i] == vert) return i;
00068 Severe("Basic logic error in SDFace::vnum()");
00069 return -1;
00070 }
00071 SDFace *nextFace(SDVertex *vert) {
00072 return f[vnum(vert)];
00073 }
00074 SDFace *prevFace(SDVertex *vert) {
00075 return f[PREV(vnum(vert))];
00076 }
00077 SDVertex *nextVert(SDVertex *vert) {
00078 return v[NEXT(vnum(vert))];
00079 }
00080 SDVertex *prevVert(SDVertex *vert) {
00081 return v[PREV(vnum(vert))];
00082 }
00083 SDVertex *otherVert(SDVertex *v0, SDVertex *v1) {
00084 for (int i = 0; i < 3; ++i)
00085 if (v[i] != v0 && v[i] != v1)
00086 return v[i];
00087 Severe("Basic logic error in SDVertex::otherVert()");
00088 return NULL;
00089 }
00090 SDVertex *v[3];
00091 SDFace *f[3];
00092 SDFace *children[4];
00093 };
00094 struct SDEdge {
00095
00096 SDEdge(SDVertex *v0 = NULL, SDVertex *v1 = NULL) {
00097 v[0] = min(v0, v1);
00098 v[1] = max(v0, v1);
00099 f[0] = f[1] = NULL;
00100 f0edgeNum = -1;
00101 }
00102
00103 bool operator<(const SDEdge &e2) const {
00104 if (v[0] == e2.v[0]) return v[1] < e2.v[1];
00105 return v[0] < e2.v[0];
00106 }
00107 SDVertex *v[2];
00108 SDFace *f[2];
00109 int f0edgeNum;
00110 };
00111
00112 class LoopSubdiv : public Shape {
00113 public:
00114
00115 LoopSubdiv(const Transform &o2w, bool ro,
00116 int nt, int nv, const int *vi,
00117 const Point *P, int nlevels);
00118 ~LoopSubdiv();
00119 bool CanIntersect() const;
00120 void Refine(vector<Reference<Shape> > &refined) const;
00121 BBox ObjectBound() const;
00122 BBox WorldBound() const;
00123 private:
00124
00125 static float beta(int valence) {
00126 if (valence == 3) return 3.f/16.f;
00127 else return 3.f / (8.f * valence);
00128 }
00129 static Point weightOneRing(SDVertex *vert, float beta);
00130 static Point weightBoundary(SDVertex *vert, float beta);
00131 static float gamma(int valence) {
00132 return 1.f / (valence + 3.f / (8.f * beta(valence)));
00133 }
00134
00135 int nLevels;
00136 vector<SDVertex *> vertices;
00137 vector<SDFace *> faces;
00138 };
00139
00140 inline int SDVertex::valence() {
00141 SDFace *f = startFace;
00142 if (!boundary) {
00143
00144 int nf = 1;
00145 while ((f = f->nextFace(this)) != startFace)
00146 ++nf;
00147 return nf;
00148 }
00149 else {
00150
00151 int nf = 1;
00152 while ((f = f->nextFace(this)) != NULL)
00153 ++nf;
00154 f = startFace;
00155 while ((f = f->prevFace(this)) != NULL)
00156 ++nf;
00157 return nf+1;
00158 }
00159 }
00160
00161 LoopSubdiv::LoopSubdiv(const Transform &o2w, bool ro,
00162 int nfaces, int nvertices,
00163 const int *vertexIndices,
00164 const Point *P, int nl)
00165 : Shape(o2w, ro) {
00166 nLevels = nl;
00167
00168 int i;
00169 SDVertex *verts = new SDVertex[nvertices];
00170 for (i = 0; i < nvertices; ++i) {
00171 verts[i] = SDVertex(P[i]);
00172 vertices.push_back(&verts[i]);
00173 }
00174 SDFace *fs = new SDFace[nfaces];
00175 for (i = 0; i < nfaces; ++i)
00176 faces.push_back(&fs[i]);
00177
00178 const int *vp = vertexIndices;
00179 for (i = 0; i < nfaces; ++i) {
00180 SDFace *f = faces[i];
00181 for (int j = 0; j < 3; ++j) {
00182 SDVertex *v = vertices[vp[j]];
00183 f->v[j] = v;
00184 v->startFace = f;
00185 }
00186 vp += 3;
00187 }
00188
00189 set<SDEdge> edges;
00190 for (i = 0; i < nfaces; ++i) {
00191 SDFace *f = faces[i];
00192 for (int edgeNum = 0; edgeNum < 3; ++edgeNum) {
00193
00194 int v0 = edgeNum, v1 = NEXT(edgeNum);
00195 SDEdge e(f->v[v0], f->v[v1]);
00196 if (edges.find(e) == edges.end()) {
00197
00198 e.f[0] = f;
00199 e.f0edgeNum = edgeNum;
00200 edges.insert(e);
00201 }
00202 else {
00203
00204 e = *edges.find(e);
00205 e.f[0]->f[e.f0edgeNum] = f;
00206 f->f[edgeNum] = e.f[0];
00207 edges.erase(e);
00208 }
00209 }
00210 }
00211
00212 for (i = 0; i < nvertices; ++i) {
00213 SDVertex *v = vertices[i];
00214 SDFace *f = v->startFace;
00215 do {
00216 f = f->nextFace(v);
00217 } while (f && f != v->startFace);
00218 v->boundary = (f == NULL);
00219 if (!v->boundary && v->valence() == 6)
00220 v->regular = true;
00221 else if (v->boundary && v->valence() == 4)
00222 v->regular = true;
00223 else
00224 v->regular = false;
00225 }
00226 }
00227 LoopSubdiv::~LoopSubdiv() {
00228 delete[] vertices[0];
00229 delete[] faces[0];
00230 }
00231 BBox LoopSubdiv::ObjectBound() const {
00232 BBox b;
00233 for (u_int i = 0; i < vertices.size(); i++)
00234 b = Union(b, vertices[i]->P);
00235 return b;
00236 }
00237 BBox LoopSubdiv::WorldBound() const {
00238 BBox b;
00239 for (u_int i = 0; i < vertices.size(); i++)
00240 b = Union(b, ObjectToWorld(vertices[i]->P));
00241 return b;
00242 }
00243 bool LoopSubdiv::CanIntersect() const {
00244 return false;
00245 }
00246 void
00247 LoopSubdiv::Refine(vector<Reference<Shape> > &refined)
00248 const {
00249 vector<SDFace *> f = faces;
00250 vector<SDVertex *> v = vertices;
00251 ObjectArena<SDVertex> vertexArena;
00252 ObjectArena<SDFace> faceArena;
00253 for (int i = 0; i < nLevels; ++i) {
00254
00255 vector<SDFace *> newFaces;
00256 vector<SDVertex *> newVertices;
00257
00258 for (u_int j = 0; j < v.size(); ++j) {
00259 v[j]->child = new (vertexArena) SDVertex;
00260 v[j]->child->regular = v[j]->regular;
00261 v[j]->child->boundary = v[j]->boundary;
00262 newVertices.push_back(v[j]->child);
00263 }
00264 for (u_int j = 0; j < f.size(); ++j)
00265 for (int k = 0; k < 4; ++k) {
00266 f[j]->children[k] = new (faceArena) SDFace;
00267 newFaces.push_back(f[j]->children[k]);
00268 }
00269
00270
00271 for (u_int j = 0; j < v.size(); ++j) {
00272 if (!v[j]->boundary) {
00273
00274 if (v[j]->regular)
00275 v[j]->child->P = weightOneRing(v[j], 1.f/16.f);
00276 else
00277 v[j]->child->P =
00278 weightOneRing(v[j], beta(v[j]->valence()));
00279 }
00280 else {
00281
00282 v[j]->child->P = weightBoundary(v[j], 1.f/8.f);
00283 }
00284 }
00285
00286 map<SDEdge, SDVertex *> edgeVerts;
00287 for (u_int j = 0; j < f.size(); ++j) {
00288 SDFace *face = f[j];
00289 for (int k = 0; k < 3; ++k) {
00290
00291 SDEdge edge(face->v[k], face->v[NEXT(k)]);
00292 SDVertex *vert = edgeVerts[edge];
00293 if (!vert) {
00294
00295 vert = new (vertexArena) SDVertex;
00296 newVertices.push_back(vert);
00297 vert->regular = true;
00298 vert->boundary = (face->f[k] == NULL);
00299 vert->startFace = face->children[3];
00300
00301 if (vert->boundary) {
00302 vert->P = 0.5f * edge.v[0]->P;
00303 vert->P += 0.5f * edge.v[1]->P;
00304 }
00305 else {
00306 vert->P = 3.f/8.f * edge.v[0]->P;
00307 vert->P += 3.f/8.f * edge.v[1]->P;
00308 vert->P += 1.f/8.f *
00309 face->otherVert(edge.v[0], edge.v[1])->P;
00310 vert->P += 1.f/8.f *
00311 face->f[k]->otherVert(edge.v[0], edge.v[1])->P;
00312 }
00313 edgeVerts[edge] = vert;
00314 }
00315 }
00316 }
00317
00318
00319 for (u_int j = 0; j < v.size(); ++j) {
00320 SDVertex *vert = v[j];
00321 int vertNum = vert->startFace->vnum(vert);
00322 vert->child->startFace =
00323 vert->startFace->children[vertNum];
00324 }
00325
00326 for (u_int j = 0; j < f.size(); ++j) {
00327 SDFace *face = f[j];
00328 for (int k = 0; k < 3; ++k) {
00329
00330 face->children[3]->f[k] = face->children[NEXT(k)];
00331 face->children[k]->f[NEXT(k)] = face->children[3];
00332
00333 SDFace *f2 = face->f[k];
00334 face->children[k]->f[k] =
00335 f2 ? f2->children[f2->vnum(face->v[k])] : NULL;
00336 f2 = face->f[PREV(k)];
00337 face->children[k]->f[PREV(k)] =
00338 f2 ? f2->children[f2->vnum(face->v[k])] : NULL;
00339 }
00340 }
00341
00342 for (u_int j = 0; j < f.size(); ++j) {
00343 SDFace *face = f[j];
00344 for (int k = 0; k < 3; ++k) {
00345
00346 face->children[k]->v[k] = face->v[k]->child;
00347
00348 SDVertex *vert =
00349 edgeVerts[SDEdge(face->v[k], face->v[NEXT(k)])];
00350 face->children[k]->v[NEXT(k)] = vert;
00351 face->children[NEXT(k)]->v[k] = vert;
00352 face->children[3]->v[k] = vert;
00353 }
00354 }
00355
00356 f = newFaces;
00357 v = newVertices;
00358 }
00359
00360 Point *Plimit = new Point[v.size()];
00361 for (u_int i = 0; i < v.size(); ++i) {
00362 if (v[i]->boundary)
00363 Plimit[i] =
00364 weightBoundary(v[i], 1.f/5.f);
00365 else
00366 Plimit[i] =
00367 weightOneRing(v[i], gamma(v[i]->valence()));
00368 }
00369 for (u_int i = 0; i < v.size(); ++i)
00370 v[i]->P = Plimit[i];
00371
00372 vector<Normal> Ns;
00373 Ns.reserve(v.size());
00374 int ringSize = 16;
00375 Point *Pring = new Point[ringSize];
00376 for (u_int i = 0; i < v.size(); ++i) {
00377 SDVertex *vert = v[i];
00378 Vector S(0,0,0), T(0,0,0);
00379 int valence = vert->valence();
00380 if (valence > ringSize) {
00381 ringSize = valence;
00382 delete[] Pring;
00383 Pring = new Point[ringSize];
00384 }
00385 vert->oneRing(Pring);
00386
00387 if (!vert->boundary) {
00388
00389 for (int k = 0; k < valence; ++k) {
00390 S += cosf(2.f*M_PI*k/valence) * Vector(Pring[k]);
00391 T += sinf(2.f*M_PI*k/valence) * Vector(Pring[k]);
00392 }
00393 }
00394 else {
00395
00396 S = Pring[valence-1] - Pring[0];
00397 if (valence == 2)
00398 T = Vector(Pring[0] + Pring[1] - 2 * vert->P);
00399 else if (valence == 3)
00400 T = Pring[1] - vert->P;
00401 else if (valence == 4)
00402 T = Vector(-1*Pring[0] + 2*Pring[1] + 2*Pring[2] +
00403 -1*Pring[3] + -2*vert->P);
00404 else {
00405 float theta = M_PI / float(valence-1);
00406 T = Vector(sinf(theta) * (Pring[0] + Pring[valence-1]));
00407 for (int k = 1; k < valence-1; ++k) {
00408 float wt = (2*cosf(theta) - 2) * sinf((k) * theta);
00409 T += Vector(wt * Pring[k]);
00410 }
00411 T = -T;
00412 }
00413 }
00414 Ns.push_back(Normal(Cross(S, T)));
00415 }
00416
00417 u_int ntris = u_int(f.size());
00418 int *verts = new int[3*ntris];
00419 int *vp = verts;
00420 u_int totVerts = u_int(v.size());
00421 map<SDVertex *, int> usedVerts;
00422 for (u_int i = 0; i < totVerts; ++i)
00423 usedVerts[v[i]] = i;
00424 for (u_int i = 0; i < ntris; ++i) {
00425 for (int j = 0; j < 3; ++j) {
00426 *vp = usedVerts[f[i]->v[j]];
00427 ++vp;
00428 }
00429 }
00430 ParamSet paramSet;
00431 paramSet.AddInt("indices", verts, 3*ntris);
00432 paramSet.AddPoint("P", Plimit, totVerts);
00433 paramSet.AddNormal("N", &Ns[0], int(Ns.size()));
00434 refined.push_back(MakeShape("trianglemesh", ObjectToWorld,
00435 reverseOrientation, paramSet));
00436 delete[] verts;
00437 delete[] Plimit;
00438 }
00439 Point LoopSubdiv::weightOneRing(SDVertex *vert, float beta) {
00440
00441 int valence = vert->valence();
00442 Point *Pring = (Point *)alloca(valence * sizeof(Point));
00443 vert->oneRing(Pring);
00444 Point P = (1 - valence * beta) * vert->P;
00445 for (int i = 0; i < valence; ++i)
00446 P += beta * Pring[i];
00447 return P;
00448 }
00449 void SDVertex::oneRing(Point *P) {
00450 if (!boundary) {
00451
00452 SDFace *face = startFace;
00453 do {
00454 *P++ = face->nextVert(this)->P;
00455 face = face->nextFace(this);
00456 } while (face != startFace);
00457 }
00458 else {
00459
00460 SDFace *face = startFace, *f2;
00461 while ((f2 = face->nextFace(this)) != NULL)
00462 face = f2;
00463 *P++ = face->nextVert(this)->P;
00464 do {
00465 *P++ = face->prevVert(this)->P;
00466 face = face->prevFace(this);
00467 } while (face != NULL);
00468 }
00469 }
00470 Point LoopSubdiv::weightBoundary(SDVertex *vert,
00471 float beta) {
00472
00473 int valence = vert->valence();
00474 Point *Pring = (Point *)alloca(valence * sizeof(Point));
00475 vert->oneRing(Pring);
00476 Point P = (1-2*beta) * vert->P;
00477 P += beta * Pring[0];
00478 P += beta * Pring[valence-1];
00479 return P;
00480 }
00481 extern "C" DLLEXPORT Shape *CreateShape(const Transform &o2w,
00482 bool reverseOrientation, const ParamSet ¶ms) {
00483 int nlevels = params.FindOneInt("nlevels", 3);
00484 int nps, nIndices;
00485 const int *vi = params.FindInt("indices", &nIndices);
00486 const Point *P = params.FindPoint("P", &nps);
00487 if (!vi || !P) return NULL;
00488
00489
00490 string scheme = params.FindOneString("scheme", "loop");
00491
00492 return new LoopSubdiv(o2w, reverseOrientation, nIndices/3, nps,
00493 vi, P, nlevels);
00494 }