[37] | 1 | /* |
---|
| 2 | |
---|
| 3 | Name: cube.c |
---|
| 4 | |
---|
| 5 | Coded: Paul Ning |
---|
| 6 | |
---|
| 7 | Modified by: Brian Curless |
---|
| 8 | Computer Graphics Laboratory |
---|
| 9 | Stanford University |
---|
| 10 | |
---|
| 11 | Comment: Processes a single cube. Is passed absolute i,j,k position. |
---|
| 12 | |
---|
| 13 | Copyright (1997) The Board of Trustees of the Leland Stanford Junior |
---|
| 14 | University. Except for commercial resale, lease, license or other |
---|
| 15 | commercial transactions, permission is hereby given to use, copy, |
---|
| 16 | modify this software for academic purposes only. No part of this |
---|
| 17 | software or any derivatives thereof may be used in the production of |
---|
| 18 | computer models for resale or for use in a commercial |
---|
| 19 | product. STANFORD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND |
---|
| 20 | CONCERNING THIS SOFTWARE. No support is implied or provided. |
---|
| 21 | |
---|
| 22 | */ |
---|
| 23 | |
---|
| 24 | |
---|
| 25 | #include <stdio.h> |
---|
| 26 | #include <iostream> |
---|
| 27 | #include <math.h> |
---|
| 28 | #include <assert.h> |
---|
| 29 | #include "mc.h" |
---|
| 30 | |
---|
| 31 | static TriangleVertex VertexOnEdge[12]; |
---|
| 32 | |
---|
| 33 | |
---|
| 34 | static void Interpolate(int n, Cube cube, int i, int j, int k); |
---|
| 35 | |
---|
| 36 | |
---|
| 37 | TriangleVertex* |
---|
| 38 | NewDoCube(Cube cube, int i, int j, int k, unsigned char *pEdgeTableIndex) |
---|
| 39 | { |
---|
| 40 | int n; |
---|
| 41 | unsigned char EdgeTableIndex = 0x00; |
---|
| 42 | |
---|
| 43 | /* form 8 bit index into edge table */ |
---|
| 44 | for (n=0;n<8;n++) |
---|
| 45 | if (cube[n].density > threshold) |
---|
| 46 | EdgeTableIndex = EdgeTableIndex | (1 << n); |
---|
| 47 | |
---|
| 48 | /* if edge table entry indicates triangles, */ |
---|
| 49 | if (TheEdgeTable[EdgeTableIndex].Ntriangles != 0) { |
---|
| 50 | |
---|
| 51 | /* interpolate the active edges, */ |
---|
| 52 | for (n=0;n<12;n++) |
---|
| 53 | if ((TheEdgeTable[EdgeTableIndex].edge)[n]) |
---|
| 54 | Interpolate(n,cube,i,j,k); |
---|
| 55 | |
---|
| 56 | TotalTriangles += TheEdgeTable[EdgeTableIndex].Ntriangles; |
---|
| 57 | |
---|
| 58 | } /* if (Ntriangles != 0) */ |
---|
| 59 | |
---|
| 60 | *pEdgeTableIndex = EdgeTableIndex; |
---|
| 61 | return VertexOnEdge; |
---|
| 62 | |
---|
| 63 | } /* DoCube */ |
---|
| 64 | |
---|
| 65 | |
---|
| 66 | /* |
---|
| 67 | * Interpolate along one edge of cube, setting VertexOnEdge[n] |
---|
| 68 | */ |
---|
| 69 | static void |
---|
| 70 | Interpolate(int n, Cube cube, int i, int j, int k) |
---|
| 71 | { |
---|
| 72 | float alpha=0,beta; |
---|
| 73 | float x=0,y=0,z=0,nx,ny,nz,mag,confidence; |
---|
| 74 | float NrmScale; /* scales nx,ny,nz to just within -32768 to 32767 */ |
---|
| 75 | Index a=0,b=0; |
---|
| 76 | unsigned char realData, valid=0; |
---|
| 77 | |
---|
| 78 | realData = TRUE; |
---|
| 79 | |
---|
| 80 | /* set the spatial components */ |
---|
| 81 | switch (n) { |
---|
| 82 | |
---|
| 83 | /* |
---|
| 84 | * four edges in k direction |
---|
| 85 | */ |
---|
| 86 | |
---|
| 87 | case 0: |
---|
| 88 | if (cube[1].density == cube[0].density) |
---|
| 89 | alpha = 0.5; |
---|
| 90 | else |
---|
| 91 | alpha=(threshold - cube[0].density)/ |
---|
| 92 | (cube[1].density - cube[0].density); |
---|
| 93 | x = xstart - i*dx; |
---|
| 94 | y = ystart + (k+alpha)*dy; |
---|
| 95 | z = zstart - (j+1)*dz; |
---|
| 96 | a = 0; |
---|
| 97 | b = 1; |
---|
| 98 | realData = cube[0].realData && cube[1].realData; |
---|
| 99 | valid = cube[0].valid && cube[1].valid; |
---|
| 100 | break; |
---|
| 101 | |
---|
| 102 | case 2: |
---|
| 103 | if (cube[2].density == cube[3].density) |
---|
| 104 | alpha = 0.5; |
---|
| 105 | else |
---|
| 106 | alpha=(threshold - cube[3].density)/ |
---|
| 107 | (cube[2].density - cube[3].density); |
---|
| 108 | x = xstart - i*dx; |
---|
| 109 | y = ystart + (k+alpha)*dy; |
---|
| 110 | z = zstart - j*dz; |
---|
| 111 | a = 3; |
---|
| 112 | b = 2; |
---|
| 113 | realData = cube[2].realData && cube[3].realData; |
---|
| 114 | valid = cube[2].valid && cube[3].valid; |
---|
| 115 | break; |
---|
| 116 | |
---|
| 117 | case 4: |
---|
| 118 | if (cube[5].density == cube[4].density) |
---|
| 119 | alpha = 0.5; |
---|
| 120 | else |
---|
| 121 | alpha=(threshold - cube[4].density)/ |
---|
| 122 | (cube[5].density - cube[4].density); |
---|
| 123 | x = xstart - (i+1)*dx; |
---|
| 124 | y = ystart + (k+alpha)*dy; |
---|
| 125 | z = zstart - (j+1)*dz; |
---|
| 126 | a = 4; |
---|
| 127 | b = 5; |
---|
| 128 | realData = cube[4].realData && cube[5].realData; |
---|
| 129 | valid = cube[4].valid && cube[5].valid; |
---|
| 130 | break; |
---|
| 131 | |
---|
| 132 | case 6: |
---|
| 133 | if (cube[6].density == cube[7].density) |
---|
| 134 | alpha = 0.5; |
---|
| 135 | else |
---|
| 136 | alpha=(threshold - cube[7].density)/ |
---|
| 137 | (cube[6].density - cube[7].density); |
---|
| 138 | x = xstart - (i+1)*dx; |
---|
| 139 | y = ystart + (k+alpha)*dy; |
---|
| 140 | z = zstart - j*dz; |
---|
| 141 | a = 7; |
---|
| 142 | b = 6; |
---|
| 143 | realData = cube[6].realData && cube[7].realData; |
---|
| 144 | valid = cube[6].valid && cube[7].valid; |
---|
| 145 | break; |
---|
| 146 | |
---|
| 147 | /* |
---|
| 148 | * four edges in j direction |
---|
| 149 | */ |
---|
| 150 | |
---|
| 151 | case 1: |
---|
| 152 | if (cube[1].density == cube[2].density) |
---|
| 153 | alpha = 0.5; |
---|
| 154 | else |
---|
| 155 | alpha=(threshold - cube[2].density)/ |
---|
| 156 | (cube[1].density - cube[2].density); |
---|
| 157 | x = xstart - i*dx; |
---|
| 158 | y = ystart + (k+1)*dy; |
---|
| 159 | z = zstart - (j+alpha)*dz; |
---|
| 160 | a = 2; |
---|
| 161 | b = 1; |
---|
| 162 | realData = cube[1].realData && cube[2].realData; |
---|
| 163 | valid = cube[1].valid && cube[2].valid; |
---|
| 164 | break; |
---|
| 165 | |
---|
| 166 | case 3: |
---|
| 167 | if (cube[0].density == cube[3].density) |
---|
| 168 | alpha = 0.5; |
---|
| 169 | else |
---|
| 170 | alpha=(threshold - cube[3].density)/ |
---|
| 171 | (cube[0].density - cube[3].density); |
---|
| 172 | x = xstart - i*dx; |
---|
| 173 | y = ystart + k*dy; |
---|
| 174 | z = zstart - (j+alpha)*dz; |
---|
| 175 | a = 3; |
---|
| 176 | b = 0; |
---|
| 177 | realData = cube[0].realData && cube[3].realData; |
---|
| 178 | valid = cube[0].valid && cube[3].valid; |
---|
| 179 | break; |
---|
| 180 | |
---|
| 181 | case 5: |
---|
| 182 | if (cube[5].density == cube[6].density) |
---|
| 183 | alpha = 0.5; |
---|
| 184 | else |
---|
| 185 | alpha=(threshold - cube[6].density)/ |
---|
| 186 | (cube[5].density - cube[6].density); |
---|
| 187 | x = xstart - (i+1)*dx; |
---|
| 188 | y = ystart + (k+1)*dy; |
---|
| 189 | z = zstart - (j+alpha)*dz; |
---|
| 190 | a = 6; |
---|
| 191 | b = 5; |
---|
| 192 | realData = cube[5].realData && cube[6].realData; |
---|
| 193 | valid = cube[5].valid && cube[6].valid; |
---|
| 194 | break; |
---|
| 195 | |
---|
| 196 | case 7: |
---|
| 197 | if (cube[4].density == cube[7].density) |
---|
| 198 | alpha = 0.5; |
---|
| 199 | else |
---|
| 200 | alpha=(threshold - cube[7].density)/ |
---|
| 201 | (cube[4].density - cube[7].density); |
---|
| 202 | x = xstart - (i+1)*dx; |
---|
| 203 | y = ystart + k*dy; |
---|
| 204 | z = zstart - (j+alpha)*dz; |
---|
| 205 | a = 7; |
---|
| 206 | b = 4; |
---|
| 207 | realData = cube[4].realData && cube[7].realData; |
---|
| 208 | valid = cube[4].valid && cube[7].valid; |
---|
| 209 | break; |
---|
| 210 | |
---|
| 211 | /* |
---|
| 212 | * four edges in i direction |
---|
| 213 | */ |
---|
| 214 | |
---|
| 215 | case 8: |
---|
| 216 | if (cube[4].density == cube[0].density) |
---|
| 217 | alpha = 0.5; |
---|
| 218 | else |
---|
| 219 | alpha=(threshold - cube[0].density)/ |
---|
| 220 | (cube[4].density - cube[0].density); |
---|
| 221 | x = xstart - (i+alpha)*dx; |
---|
| 222 | y = ystart + k*dy; |
---|
| 223 | z = zstart - (j+1)*dz; |
---|
| 224 | a = 0; |
---|
| 225 | b = 4; |
---|
| 226 | realData = cube[0].realData && cube[4].realData; |
---|
| 227 | valid = cube[0].valid && cube[4].valid; |
---|
| 228 | break; |
---|
| 229 | |
---|
| 230 | case 9: |
---|
| 231 | if (cube[5].density == cube[1].density) |
---|
| 232 | alpha = 0.5; |
---|
| 233 | else |
---|
| 234 | alpha=(threshold - cube[1].density)/ |
---|
| 235 | (cube[5].density - cube[1].density); |
---|
| 236 | x = xstart - (i+alpha)*dx; |
---|
| 237 | y = ystart + (k+1)*dy; |
---|
| 238 | z = zstart - (j+1)*dz; |
---|
| 239 | a = 1; |
---|
| 240 | b = 5; |
---|
| 241 | realData = cube[1].realData && cube[5].realData; |
---|
| 242 | valid = cube[1].valid && cube[5].valid; |
---|
| 243 | break; |
---|
| 244 | |
---|
| 245 | case 10: |
---|
| 246 | if (cube[7].density == cube[3].density) |
---|
| 247 | alpha = 0.5; |
---|
| 248 | else |
---|
| 249 | alpha=(threshold - cube[3].density)/ |
---|
| 250 | (cube[7].density - cube[3].density); |
---|
| 251 | x = xstart - (i+alpha)*dx; |
---|
| 252 | y = ystart + k*dy; |
---|
| 253 | z = zstart - j*dz; |
---|
| 254 | a = 3; |
---|
| 255 | b = 7; |
---|
| 256 | realData = cube[3].realData && cube[7].realData; |
---|
| 257 | valid = cube[3].valid && cube[7].valid; |
---|
| 258 | break; |
---|
| 259 | |
---|
| 260 | case 11: |
---|
| 261 | if (cube[6].density == cube[2].density) |
---|
| 262 | alpha = 0.5; |
---|
| 263 | else |
---|
| 264 | alpha=(threshold - cube[2].density)/ |
---|
| 265 | (cube[6].density - cube[2].density); |
---|
| 266 | x = xstart - (i+alpha)*dx; |
---|
| 267 | y = ystart + (k+1)*dy; |
---|
| 268 | z = zstart - j*dz; |
---|
| 269 | a = 2; |
---|
| 270 | b = 6; |
---|
| 271 | realData = cube[2].realData && cube[6].realData; |
---|
| 272 | valid = cube[2].valid && cube[6].valid; |
---|
| 273 | break; |
---|
| 274 | } |
---|
| 275 | |
---|
| 276 | /* compute the normal components and magnitude */ |
---|
| 277 | beta = 1 - alpha; |
---|
| 278 | nx = beta*cube[a].nx + alpha*cube[b].nx; |
---|
| 279 | ny = beta*cube[a].ny + alpha*cube[b].ny; |
---|
| 280 | nz = beta*cube[a].nz + alpha*cube[b].nz; |
---|
| 281 | confidence = beta*cube[a].confidence + alpha*cube[b].confidence; |
---|
| 282 | mag = sqrt(nx*nx + ny*ny + nz*nz); |
---|
| 283 | |
---|
| 284 | if (mag == 0) { |
---|
| 285 | nx = 0; |
---|
| 286 | ny = 0; |
---|
| 287 | nz = 0; |
---|
| 288 | } |
---|
| 289 | else { |
---|
| 290 | nx = nx/mag; |
---|
| 291 | ny = ny/mag; |
---|
| 292 | nz = nz/mag; |
---|
| 293 | } |
---|
| 294 | |
---|
| 295 | VertexOnEdge[n].x = x; |
---|
| 296 | VertexOnEdge[n].y = y; |
---|
| 297 | VertexOnEdge[n].z = z; |
---|
| 298 | VertexOnEdge[n].nx = nx; |
---|
| 299 | VertexOnEdge[n].ny = ny; |
---|
| 300 | VertexOnEdge[n].nz = nz; |
---|
| 301 | VertexOnEdge[n].realData = realData; |
---|
| 302 | VertexOnEdge[n].valid = valid; |
---|
| 303 | |
---|
| 304 | if (SaveGradientAsConfidence) |
---|
| 305 | VertexOnEdge[n].confidence = mag; |
---|
| 306 | else |
---|
| 307 | VertexOnEdge[n].confidence = confidence; |
---|
| 308 | |
---|
| 309 | VertexOnEdge[n].tex1 = (unsigned char) round(mag); |
---|
| 310 | VertexOnEdge[n].tex2 = 0; |
---|
| 311 | |
---|
| 312 | } /* Interpolate */ |
---|
| 313 | |
---|
| 314 | |
---|