source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/vrippack-0.31/src/march/slices.cc @ 37

Last change on this file since 37 was 37, checked in by (none), 14 years ago

Added original make3d
File size: 16.5 KB
Line 
1/*
2
3Name:         slices.cc
4
5Coded:        Paul Ning
6
7Modified by:  Brian Curless
8              Computer Graphics Laboratory
9              Stanford University
10
11Comment:      Processes all selected slices.
12
13Copyright (1997) The Board of Trustees of the Leland Stanford Junior
14University. Except for commercial resale, lease, license or other
15commercial transactions, permission is hereby given to use, copy,
16modify this software for academic purposes only.  No part of this
17software or any derivatives thereof may be used in the production of
18computer models for resale or for use in a commercial
19product. STANFORD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND
20CONCERNING THIS SOFTWARE.  No support is implied or provided.
21
22*/
23
24
25#include "Linear.h"
26#include <signal.h>
27#include <stdio.h>
28#include <string.h>
29#include "mc.h"
30#include <unistd.h>
31#include "ply.h"
32#include "OccGridRLE.h"
33#include "SectionRLE.h"
34#include "mc_more.h"
35#include <assert.h>
36
37//#define USE_VALUE_WEIGHT_PRODUCT
38
39#define GET_MORE_CUBE_TRIS
40
41struct PlyVertex {
42    float x, y, z;
43    float nx, ny, nz;
44    float confidence;
45    uchar realData;
46};
47
48struct PlyTri {
49    unsigned char nverts;
50    int *verts;
51};
52
53static char *elem_names[] = { 
54  "vertex", "face",
55};
56
57static PlyProperty vert_props[] =  {
58    {"x", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0},
59    {"y", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0},
60    {"z", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0},
61    {"confidence", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0},
62    {"real_data", PLY_UCHAR, PLY_UCHAR, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0},
63    {"nx", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0},
64    {"ny", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0},
65    {"nz", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0},
66};
67
68static PlyProperty face_props[] = { 
69  {"vertex_indices", PLY_INT, PLY_INT, 0, 1, PLY_UCHAR, PLY_UCHAR, 0},
70};
71
72
73int loadCube(Cube cube, SectionRLE **section, int index, int xx, int yy);
74void computeGradients(SectionRLE **section, int level);
75static void UpdateGeometry(int j, int k, SectionRLE **rleSection, 
76                           TriangleVertex *vertexList, uchar EdgeTableIndex);
77static void writePlyFile(char *filename, float tx, float ty, float tz);
78static Vertex *findVertex(TriangleVertex *vert, TriSet *set);
79static void loadSection(SectionRLE *rleSection, 
80                        OccGridRLE *occGrid, int sliceNum);
81static void copyElement(Point *sect, SectionElement *elem);
82static Vertex *addVertex(TriangleVertex *inVert, int j, int k, 
83                         SectionRLE **rleSection);
84int loadCubeAllowPartiallyValid(Cube cube, SectionRLE **section, 
85                                int index, int xx, int yy);
86
87
88static int numVerts;
89static int numTris;
90static ChunkAllocator *triChunker;
91static ChunkAllocator *vertChunker;
92
93
94void 
95DoSlicesOccRLE()
96{
97  SectionRLE *rleSection[5], *temp;
98  int i,j,k; /* indices within sections - relative to i0,j_0,k0 */
99  int slices,rows,cols; /* dimensions of stack of sections */
100  int iter;
101  Cube cube;
102  unsigned char EdgeTableIndex = 0x00;
103  int xx, yy, zz, ii;
104  int runType;
105  RunLength length;
106  SectionElement *element;   
107
108  numVerts = 0;
109  numTris = 0;
110
111  triChunker = new ChunkAllocator(CHUNK_SIZE);
112  vertChunker = new ChunkAllocator(CHUNK_SIZE);
113
114 
115
116  /*
117   * allocate memory for four sections
118   */
119  slices = i1 - i0 + 3;
120  rows = j_1 - j_0 + 3;
121  cols = k1 - k0 + 3;
122
123  if (UseValueWeightProduct) {
124     signal(SIGFPE, SIG_IGN);
125  }
126
127  for (i=0;i<5;i++)
128    rleSection[i] = new SectionRLE(rows,cols,CHUNK_SIZE);
129
130  /*
131   * read initial three slices (into section[1-3])
132   */
133  for (iter=0;iter<3;iter++) {
134    /*
135     * read appropriate regions of slice file
136     */
137      loadSection(rleSection[iter+1], occGrid, FirstSliceFileNumber+i0-1+iter);
138  }
139
140
141  /*
142   * compute normals in section[2]
143   */
144
145  computeGradients(rleSection, 2);
146
147  /*
148   * Main Loop - Read New Slice
149   */
150
151  fflush(stdout);
152
153  for (i=1;i<slices-2;i++) {
154
155    /* rotate section memories to make room for new slice */
156      temp = rleSection[0];
157      rleSection[0] = rleSection[1];
158      rleSection[1] = rleSection[2];
159      rleSection[2] = rleSection[3];
160      rleSection[3] = temp;
161
162      /*
163       * read appropriate regions of slice file
164       */
165     
166      loadSection(rleSection[3], occGrid, FirstSliceFileNumber+i0+i+1);     
167     
168      computeGradients(rleSection, 2);
169     
170      /* Process the cubes in raster order */
171      for (yy = 1; yy < rleSection[1]->ydim-2; yy++) {
172          rleSection[1]->setScanline(yy);
173          xx = 0;
174          while (TRUE) {
175
176              length = rleSection[1]->getNextLength();
177              runType = rleSection[1]->getRunType(&length);
178
179              if (runType == SectionRLE::END_OF_RUN)
180                  break;
181         
182              if (runType == SectionRLE::CONSTANT_DATA) {
183                  xx += length;
184              }
185              else {
186                  for (ii=0; ii<length; ii++, xx++) {
187                      if (xx == 0 || xx >= rleSection[2]->xdim-2)
188                          continue;
189
190#ifdef GET_MORE_CUBE_TRIS
191                      if (!loadCubeAllowPartiallyValid(cube, rleSection, 
192                                                       1, xx, yy))
193                          continue;
194#else
195                      if (!loadCube(cube, rleSection, 1, xx, yy))
196                          continue;
197#endif
198
199                      /* dispatch to DoCube with absolute i,j,k position */
200                      TriangleVertex *vertexList = 
201                          NewDoCube(cube,i0+i,j_0+xx,k0+yy,&EdgeTableIndex);
202                 
203                      if (TheEdgeTable[EdgeTableIndex].Ntriangles > 0)
204                          UpdateGeometry(xx, yy, rleSection, vertexList, 
205                                     EdgeTableIndex);
206                  }
207              }
208          }
209      }
210      printf("\rProcess slice %d of %d.  Triangle count = %d", 
211             i+3, slices+1, TotalTriangles);
212      fflush(stdout);
213     
214  } /* for (i) */
215
216  printf("\n\n");
217  fflush(stdout);
218
219
220  printf("Writing out triangles...\n");
221  writePlyFile(outfile, 
222               occGrid->origin[0] - occGrid->resolution, 
223               occGrid->origin[1] - occGrid->resolution, 
224               occGrid->origin[2] - occGrid->resolution);
225
226  printf("Done.\n\n");
227
228} /* DoSlices */
229
230
231static void
232UpdateGeometry(int j, int k, SectionRLE **rleSection, 
233               TriangleVertex *vertexList, uchar EdgeTableIndex)
234{
235    int n;
236    Triple triangle;
237    Tri *tri;
238    TriSet *set;
239
240    set = &rleSection[2]->getElement(j,k)->set;
241    set->ntris = 0;
242    TriangleVertex *v0, *v1, *v2;
243    for (n=0;n<TheEdgeTable[EdgeTableIndex].Ntriangles;n++) {
244        triangle = (TheEdgeTable[EdgeTableIndex].TriangleList)[n];
245
246        v0 = &vertexList[triangle.A];
247        v1 = &vertexList[triangle.B];
248        v2 = &vertexList[triangle.C];
249
250#ifdef GET_MORE_CUBE_TRIS
251        if (!v0->valid || !v1->valid || !v2->valid) {
252           continue;
253        }
254#endif
255
256        // Don't add degenerate triangles
257        if (((v0->x == v1->x) && (v0->y == v1->y) && (v0->z == v1->z)) ||
258            ((v0->x == v2->x) && (v0->y == v2->y) && (v0->z == v2->z)) ||
259            ((v2->x == v1->x) && (v2->y == v1->y) && (v2->z == v1->z))) {
260
261            continue;
262        }
263       
264        tri = (Tri *)triChunker->alloc(sizeof(Tri));
265        tri->verts[0] = addVertex(v0,j,k,rleSection);
266        tri->verts[1] = addVertex(v1,j,k,rleSection);
267        tri->verts[2] = addVertex(v2,j,k,rleSection);
268        numTris++;
269
270        set->tris[set->ntris] = tri;
271        set->ntris++;
272    }
273}
274
275
276static Vertex *
277addVertex(TriangleVertex *inVert, int j, int k, SectionRLE **rleSection)
278{
279    Vertex *outVert, *newVert, *oldVert, temp;
280    TriSet *set;
281    SectionElement *elem;
282
283    // We should only need to go to (j,k), not (j+1,k+1), but
284    // since we're dropping triangles with low confidence, we
285    // need to account for the possiblity that vertices below and
286    // behind may be dropped - really should just need to consider
287    // the (j+1,k+1) in the slice below...
288
289    for (int zz = 1; zz <= 2; zz++) {
290        for (int yy = k-1; yy <= k+1; yy++) {
291            for (int xx = j-1; xx <= j+1; xx++) {
292                elem = rleSection[zz]->getElement(xx,yy);
293                set = &elem->set;
294                oldVert = findVertex(inVert, set);
295                if (oldVert != NULL)
296                    return oldVert;
297
298            }
299        }
300    }
301
302    outVert = (Vertex *)vertChunker->alloc(sizeof(Vertex));
303    outVert->x = inVert->x;
304    outVert->y = inVert->y;
305    outVert->z = inVert->z;
306    outVert->nx = inVert->nx;
307    outVert->ny = inVert->ny;
308    outVert->nz = inVert->nz;
309    outVert->confidence = inVert->confidence;
310    outVert->realData = inVert->realData;
311    outVert->index = numVerts;
312
313    numVerts++;
314
315    return outVert;
316}
317
318
319static Vertex*
320findVertex(TriangleVertex *vert, TriSet *set)
321{
322    int i, j;
323    Vertex *oldVert;
324
325    for (i = 0; i < set->ntris; i++) {
326        for (j = 0; j < 3; j++) {
327            oldVert = set->tris[i]->verts[j];
328#if 0
329            if (fabs(oldVert->x - vert->x) < CULL_TRIANGLE_FACTOR*dx &&
330                fabs(oldVert->y - vert->y) < CULL_TRIANGLE_FACTOR*dy &&
331                fabs(oldVert->z - vert->z) < CULL_TRIANGLE_FACTOR*dz)
332                return oldVert;
333#else
334            if (oldVert->x == vert->x &&
335                oldVert->y == vert->y &&
336                oldVert->z == vert->z)
337                return oldVert;
338#endif
339        }
340    }
341
342    return NULL;
343}
344
345
346
347
348static void
349loadSection(SectionRLE *rleSection, OccGridRLE *occGrid, int sliceNum)
350{
351    RunLength length;
352    int runType, xx;
353    SectionElement sectElement;
354    OccElement *occElement;
355    ushort totalWeight;
356
357    rleSection->reset();
358
359    for (int yy = 0; yy < occGrid->ydim; yy++) {
360        occGrid->setScanline(yy, sliceNum);
361        rleSection->allocNewRun(yy);
362        xx = 0;
363        while (TRUE) {
364            length = occGrid->getNextLength();
365            rleSection->putNextLength(length);
366            runType = rleSection->getRunType(&length);
367
368            if (runType == SectionRLE::END_OF_RUN) {
369                break;
370            }
371
372            if (runType == SectionRLE::CONSTANT_DATA) {
373                occElement = occGrid->getNextElement();
374
375                // This line was not present, causing compiler
376                // warnings up until 6/5/06.  (Brian Curless)
377                totalWeight = occElement->totalWeight & 
378                   ~OccGridRLE::FALSE_DATA_BIT;
379
380                if (UseValueWeightProduct) {
381                   if (totalWeight > 0) {
382                      sectElement.density = (occElement->value/256.0 - 128)*
383                         occElement->totalWeight + 128;
384                   } else {
385                      sectElement.density = occElement->value/256.0;
386                   }
387                }
388                else {
389                   sectElement.density = occElement->value/256.0;
390                }
391
392                sectElement.confidence = occElement->totalWeight/256.0;
393                sectElement.valid = FALSE;
394                sectElement.realData = TRUE;
395                sectElement.set.ntris = 0;
396
397                rleSection->putNextElement(&sectElement);
398                xx += length;
399            }
400            else {
401                for (int i = 0; i < length; i++) {
402                    occElement = occGrid->getNextElement();
403                    totalWeight = occElement->totalWeight & 
404                        ~OccGridRLE::FALSE_DATA_BIT;
405                    if (UseValueWeightProduct) {
406                       if (totalWeight > 0) {
407                          sectElement.density = (occElement->value/256.0 - 128)
408                             *occElement->totalWeight + 128;
409                       } else {
410                          sectElement.density = occElement->value/256.0;
411                       }
412                    }
413                    else {
414                       sectElement.density = occElement->value/256.0;
415                    }
416                    sectElement.confidence = totalWeight/256.0;
417                    sectElement.valid = 
418                        totalWeight > OCC_CONF_THRESHOLD;
419                    sectElement.realData = 
420                        (occElement->totalWeight & 
421                          OccGridRLE::FALSE_DATA_BIT) == 0;
422
423                    sectElement.set.ntris = 0;
424                    rleSection->putNextElement(&sectElement);
425                    xx++;
426                }
427            }
428        }
429    }
430}
431
432
433static void
434copyElement(Point *sect, SectionElement *elem)
435{
436    sect->density = elem->density;
437    sect->nx = elem->nx;
438    sect->ny = elem->ny;
439    sect->nz = elem->nz;
440    sect->valid = elem->valid;
441    sect->realData = elem->realData;
442    sect->confidence = elem->confidence;
443}
444
445
446void
447computeGradients(SectionRLE **section, int level)
448{
449   int yy, xx, i;
450   int runType;
451   RunLength length;
452   SectionElement *element, *elem1, *elem2;   
453
454   for (yy = 1; yy < section[level]->ydim-1; yy++) {
455      section[level]->setScanline(yy);
456      xx = 0;
457      while (TRUE) {
458
459          length = section[level]->getNextLength();
460          runType = section[level]->getRunType(&length);
461
462          if (runType == SectionRLE::END_OF_RUN)
463              break;
464         
465          if (runType == SectionRLE::CONSTANT_DATA) {
466              element = section[level]->getNextElement();
467              xx += length;
468          }
469          else {
470              for (i=0; i<length; i++, xx++) {
471                  element = section[level]->getNextElement();
472                  if (xx == 0 || xx == section[level]->xdim-1)
473                      continue;
474                  element->nx = 
475                     -(section[level+1]->getElement(xx, yy)->density
476                       - section[level-1]->getElement(xx, yy)->density)/2.0/dx;
477                  element->ny = 
478                     (section[level]->getElement(xx, yy+1)->density
479                      - section[level]->getElement(xx, yy-1)->density)/2.0/dy;
480                  element->nz = 
481                     (section[level]->getElement(xx-1, yy)->density
482                      - section[level]->getElement(xx+1, yy)->density)/2.0/dz;
483              }
484          }
485      }
486  }
487}
488
489
490int
491loadCube(Cube cube, SectionRLE **section, int index, int xx, int yy)
492{
493    copyElement(&cube[0], section[index]->getElement(xx+1,yy));
494    if (!cube[0].valid) {
495        return 0;
496    }
497                 
498    copyElement(&cube[1], section[index]->getElement(xx+1,yy+1));
499    if (!cube[1].valid)
500        return 0;
501                 
502    copyElement(&cube[2], section[index]->getElement(xx,yy+1));
503    if (!cube[2].valid)
504        return 0;
505                 
506    copyElement(&cube[3], section[index]->getElement(xx,yy));
507    if (!cube[3].valid)
508        return 0;
509   
510    copyElement(&cube[4], section[index+1]->getElement(xx+1,yy));
511    if (!cube[4].valid)
512        return 0;
513   
514    copyElement(&cube[5], section[index+1]->getElement(xx+1,yy+1));
515    if (!cube[5].valid)
516        return 0;
517   
518    copyElement(&cube[6], section[index+1]->getElement(xx,yy+1));
519    if (!cube[6].valid)
520        return 0;
521   
522    copyElement(&cube[7], section[index+1]->getElement(xx,yy));
523    if (!cube[7].valid)
524        return 0;
525
526    return 1;
527}
528
529
530int
531loadCubeAllowPartiallyValid(Cube cube, SectionRLE **section, 
532                            int index, int xx, int yy)
533{
534   uchar valid_count = 0;
535
536   copyElement(&cube[0], section[index]->getElement(xx+1,yy));
537   valid_count += cube[0].valid;
538                 
539   copyElement(&cube[1], section[index]->getElement(xx+1,yy+1));
540   valid_count += cube[1].valid;
541                 
542   copyElement(&cube[2], section[index]->getElement(xx,yy+1));
543   valid_count += cube[2].valid;
544   
545   copyElement(&cube[3], section[index]->getElement(xx,yy));
546   valid_count += cube[3].valid;
547   
548   copyElement(&cube[4], section[index+1]->getElement(xx+1,yy));
549   valid_count += cube[4].valid;
550   
551   // First early exit
552   if (valid_count < 1)
553      return 0;
554     
555   copyElement(&cube[5], section[index+1]->getElement(xx+1,yy+1));
556   valid_count += cube[5].valid;
557   
558   if (valid_count < 2)
559      return 0;
560
561   copyElement(&cube[6], section[index+1]->getElement(xx,yy+1));
562   valid_count += cube[6].valid;
563   
564   if (valid_count < 3)
565      return 0;
566
567   copyElement(&cube[7], section[index+1]->getElement(xx,yy));
568   valid_count += cube[7].valid;
569
570   if (valid_count < 4)
571      return 0;
572   else
573      return 1;
574}
575                 
576
577static void
578writePlyFile(char *filename, float tx, float ty, float tz)
579{
580    int i, nvp;
581    float version;
582    PlyFile *ply;
583    PlyVertex aVert;
584    Vertex *vert;
585    Tri *currentTri;
586    PlyTri aTri;
587    Vec3f norm;
588
589    ply = ply_open_for_writing(filename, 2, elem_names, 
590                               PLY_BINARY_BE, &version);
591
592    nvp = 0;
593    vert_props[nvp].offset = offsetof(PlyVertex, x); nvp++;
594    vert_props[nvp].offset = offsetof(PlyVertex, y); nvp++;
595    vert_props[nvp].offset = offsetof(PlyVertex, z); nvp++;
596
597
598    vert_props[nvp].offset = offsetof(PlyVertex, confidence); nvp++;
599
600    vert_props[nvp].offset = offsetof(PlyVertex, realData); nvp++;
601
602    if (WriteNormals) {
603        vert_props[nvp].offset = offsetof(PlyVertex, nx); nvp++;
604        vert_props[nvp].offset = offsetof(PlyVertex, ny); nvp++;
605        vert_props[nvp].offset = offsetof(PlyVertex, nz); nvp++;
606    }
607
608
609    face_props[0].offset = offsetof(PlyTri, verts);
610    face_props[0].count_offset = offsetof(PlyTri, nverts);  /* count offset */
611   
612    ply_describe_element (ply, "vertex", numVerts, nvp, vert_props);
613
614    ply_describe_element (ply, "face", numTris, 1, face_props);
615
616    ply_header_complete (ply);
617   
618    /* set up and write the vertex elements */
619    ply_put_element_setup (ply, "vertex");
620
621    vertChunker->reset();
622    for (i = 0; i < numVerts; i++) {
623        vert = (Vertex *)vertChunker->nextElem(sizeof(Vertex));
624
625        aVert.x = -vert->z + tx;
626        aVert.y = vert->y + ty;
627        aVert.z = -vert->x + tz;
628
629        if (WriteNormals) {
630            norm.setValue(vert->nz, -vert->ny, vert->nx);
631            norm.normalize();
632            aVert.nx = norm.x;
633            aVert.ny = norm.y;
634            aVert.nz = norm.z;
635        }
636
637        aVert.confidence = vert->confidence;
638        aVert.realData = vert->realData;
639       
640        ply_put_element (ply, (void *) &aVert);
641    }
642
643    /* set up and write the face elements */
644    ply_put_element_setup (ply, "face");
645
646    int v[3];
647    aTri.nverts = 3;
648    aTri.verts = (int *)v;
649    triChunker->reset();
650    for (i = 0; i < numTris; i++) {
651        currentTri = (Tri *)triChunker->nextElem(sizeof(Tri));
652
653        aTri.verts[0] = currentTri->verts[1]->index;
654        aTri.verts[1] = currentTri->verts[0]->index;
655        aTri.verts[2] = currentTri->verts[2]->index;
656
657        ply_put_element (ply, (void *) &aTri);
658    }
659   
660    /* close the PLY file */
661    ply_close (ply);   
662}
Note: See TracBrowser for help on using the repository browser.