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