Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members  

guge_marchcube.cpp

Go to the documentation of this file.
00001 /* LIBGUL - Geometry Utility Library
00002  * Copyright (C) 1998-1999 Norbert Irmer
00003  *
00004  * This library is free software; you can redistribute it and/or
00005  * modify it under the terms of the GNU Library General Public
00006  * License as published by the Free Software Foundation; either
00007  * version 2 of the License, or (at your option) any later version.
00008  *
00009  * This library is distributed in the hope that it will be useful,
00010  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012  * Library General Public License for more details.
00013  *
00014  * You should have received a copy of the GNU Library General Public
00015  * License along with this library; if not, write to the
00016  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
00017  * Boston, MA 02111-1307, USA.
00018  */
00019 
00020 #include "stdafx.h"
00021 
00022 #include <iostream>
00023  
00024 #include "gul_types.h"
00025 #include "gul_vector.h"
00026 #include "guge_mcpattern.h"
00027 #include "guge_marchcube.h"
00028 
00029 namespace guge {
00030 
00031 using gul::normalize;
00032 using gul::cross_product;
00033 
00034 #define VEC(v,a,b,c) {(v).x=(a); (v).y=(b); (v).z=(c);}
00035 #define DER(v,k,i,j) {\
00036             (v).x = (F[k][(i)][(j+1)] - F[k][(i)][(j-1)]) / dx2;\
00037             (v).y = (F[k][(i+1)][(j)] - F[k][(i-1)][(j)]) / dy2;\
00038             (v).z = (F[k+1][(i)][(j)] - F[k-1][(i)][(j)]) / dz2;}
00039 
00040 #define MPT(a,b) med[ EdgePerms[perm[(a)]][perm[(b)]] ]
00041 #define NPT(a,b) norm[EdgePerms[perm[(a)]][perm[(b)]] ]
00042 
00043 #define FUNC(a,b,c,d,e,f) RenumberVerts((a),(b),(c),(d),(e),(f),\
00044                                         inv,func,fdat)
00045 #define CON_X1 0x01
00046 #define CON_X2 0x02
00047 #define CON_Y1 0x04
00048 #define CON_Y2 0x08
00049 #define CON_Z1 0x10
00050 #define CON_Z2 0x20
00051 #define CON_SPECIAL 0x40
00052 #define DRAW_ALT 0x80
00053 
00054 /*----------------------------------------------------------------------------
00055   Renumber the vertices of a triangle, so that its normal has the same
00056   direction as the middled gradient
00057 ---------------------------------------------------------------------------*/
00058 template< class T >
00059 void RenumberVerts(
00060               const point<T>& P1, const point<T>& P2, const point<T>& P3,
00061               const point<T>& T1, const point<T>& T2, const point<T>& T3,
00062               bool inv,
00063               void (*func) (
00064                 bool, 
00065                 const point<T> *, const point<T> *, const point<T> *,
00066                 const point<T> *, const point<T> *, const point<T> *, void * ),
00067               void *fdat )
00068 {
00069   point<T> a,b,N;
00070   T s;
00071  
00072   if( !inv )
00073   {
00074     a = P2 - P1;
00075     b = P3 - P1;
00076     N = cross_product( a, b );
00077   
00078     s=(T1.x+T2.x+T3.x)*N.x + (T1.y+T2.y+T3.y)*N.y + (T1.z+T2.z+T3.z)*N.z;
00079     
00080     if( s > 0.0 )
00081       func(0,&P1,&P2,&P3,&T1,&T2,&T3,fdat);
00082     else
00083       func(1,&P1,&P3,&P2,&T1,&T3,&T2,fdat); /* front and back faces swapped */
00084   }
00085   else
00086   {
00087     a = P3 - P1;
00088     b = P2 - P1;
00089     N = cross_product( a, b );
00090   
00091     s=(T1.x+T2.x+T3.x)*N.x + (T1.y+T2.y+T3.y)*N.y + (T1.z+T2.z+T3.z)*N.z;
00092 
00093     if( s > 0.0 )
00094       func(0,&P1,&P3,&P2,&T1,&T3,&T2,fdat);
00095     else
00096       func(1,&P1,&P2,&P3,&T1,&T2,&T3,fdat);
00097   }
00098 }        
00099 
00100 
00101 /* -------------------------------------------------------------------------
00102   Functions for drawing the 14 patterns
00103 --------------------------------------------------------------------------- */
00104               
00105 /* pattern 1
00106     6 7    0 0
00107    4 5    0 0
00108     2 3    0 0
00109    0 1    x 0
00110 */
00111 template< class T > void Draw1( 
00112               unsigned char pat, 
00113               point<T> *med, 
00114               point<T> *norm,
00115               void (*func) (
00116                 bool, 
00117                 const point<T> *, const point<T> *, const point<T> *,
00118                 const point<T> *, const point<T> *, const point<T> *, void * ),
00119               void *fdat )
00120 {
00121   int *perm;
00122   bool inv = PatternTable[pat].nbits > 4;
00123   perm = VertexPerms[PatternTable[pat].perm];
00124   FUNC( MPT(0,2), MPT(0,4), MPT(0,1), NPT(0,2), NPT(0,4), NPT(0,1) );   
00125 }
00126 
00127 /* pattern 2
00128     6 7    0 0
00129    4 5    0 0
00130     2 3    0 0
00131    0 1    x x
00132 */
00133 template< class T > void Draw2( 
00134               unsigned char pat, 
00135               point<T> *med, 
00136               point<T> *norm,
00137               void (*func) (
00138                 bool, 
00139                 const point<T> *, const point<T> *, const point<T> *,
00140                 const point<T> *, const point<T> *, const point<T> *, void * ),
00141               void *fdat )
00142 {
00143   int *perm;
00144   bool inv = PatternTable[pat].nbits > 4;
00145   perm = VertexPerms[PatternTable[pat].perm];
00146   FUNC( MPT(0,2), MPT(1,5), MPT(1,3), NPT(0,2), NPT(1,5), NPT(1,3) );   
00147   FUNC( MPT(0,2), MPT(0,4), MPT(1,5), NPT(0,2), NPT(0,4), NPT(1,5) );   
00148 }
00149 
00150 /* pattern 3
00151     6 7    0 0
00152    4 5    0 x
00153     2 3    0 0
00154    0 1    x 0
00155 */
00156 template< class T > void Draw3( 
00157               unsigned char pat, 
00158               point<T> *med, 
00159               point<T> *norm,
00160               void (*func) (
00161                 bool, 
00162                 const point<T> *, const point<T> *, const point<T> *,
00163                 const point<T> *, const point<T> *, const point<T> *, void * ),
00164               void *fdat )
00165 {
00166   int *perm;
00167   bool inv = PatternTable[pat].nbits > 4;
00168   perm = VertexPerms[PatternTable[pat].perm];
00169   FUNC( MPT(0,2), MPT(0,4), MPT(0,1), NPT(0,2), NPT(0,4), NPT(0,1) );   
00170   FUNC( MPT(4,5), MPT(5,7), MPT(1,5), NPT(4,5), NPT(5,7), NPT(1,5) );   
00171 }
00172 template< class T > void Draw3a( 
00173               unsigned char pat, 
00174               point<T> *med, 
00175               point<T> *norm,
00176               void (*func) (
00177                 bool, 
00178                 const point<T> *, const point<T> *, const point<T> *,
00179                 const point<T> *, const point<T> *, const point<T> *, void * ),
00180               void *fdat )
00181 {
00182   int *perm;
00183   bool inv = PatternTable[pat].nbits > 4;
00184   perm = VertexPerms[PatternTable[pat].perm];
00185   FUNC( MPT(0,4), MPT(4,5), MPT(0,2), NPT(0,4), NPT(4,5), NPT(0,2) );   
00186   FUNC( MPT(4,5), MPT(5,7), MPT(0,2), NPT(4,5), NPT(5,7), NPT(0,2) );   
00187   FUNC( MPT(0,2), MPT(5,7), MPT(1,5), NPT(0,2), NPT(5,7), NPT(1,5) );
00188   FUNC( MPT(1,5), MPT(0,1), MPT(0,2), NPT(1,5), NPT(0,1), NPT(0,2) );
00189 }
00190 
00191 /* pattern 4
00192     6 7    0 x
00193    4 5    0 0
00194     2 3    0 0
00195    0 1    x 0
00196 */
00197 template< class T > void Draw4( 
00198               unsigned char pat, 
00199               point<T> *med, 
00200               point<T> *norm,
00201               void (*func) (
00202                 bool, 
00203                 const point<T> *, const point<T> *, const point<T> *,
00204                 const point<T> *, const point<T> *, const point<T> *, void * ),
00205               void *fdat )
00206 {
00207   int *perm;
00208   bool inv = PatternTable[pat].nbits > 4;
00209   perm = VertexPerms[PatternTable[pat].perm];
00210   FUNC( MPT(0,2), MPT(0,4), MPT(0,1), NPT(0,2), NPT(0,4), NPT(0,1) );   
00211   FUNC( MPT(6,7), MPT(3,7), MPT(5,7), NPT(6,7), NPT(3,7), NPT(5,7) );   
00212 }
00213 
00214 /* pattern 5
00215     6 7    0 0
00216    4 5    0 0
00217     2 3    0 x
00218    0 1    x x
00219 */
00220 template< class T > void Draw5( 
00221               unsigned char pat, 
00222               point<T> *med, 
00223               point<T> *norm,
00224               void (*func) (
00225                 bool, 
00226                 const point<T> *, const point<T> *, const point<T> *,
00227                 const point<T> *, const point<T> *, const point<T> *, void * ),
00228               void *fdat )
00229 {
00230   int *perm;
00231   bool inv = PatternTable[pat].nbits > 4;
00232   perm = VertexPerms[PatternTable[pat].perm];
00233   FUNC( MPT(0,2), MPT(0,4), MPT(2,3), NPT(0,2), NPT(0,4), NPT(2,3)  );   
00234   FUNC( MPT(0,4), MPT(3,7), MPT(2,3), NPT(0,4), NPT(3,7), NPT(2,3) );   
00235   FUNC( MPT(0,4), MPT(1,5), MPT(3,7), NPT(0,4), NPT(1,5), NPT(3,7) );   
00236 }
00237 
00238 /* pattern 6
00239     6 7    0 x
00240    4 5    0 0
00241     2 3    0 0
00242    0 1    x x
00243 */
00244 template< class T > void Draw6( 
00245               unsigned char pat, 
00246               point<T> *med, 
00247               point<T> *norm,
00248               void (*func) (
00249                 bool, 
00250                 const point<T> *, const point<T> *, const point<T> *,
00251                 const point<T> *, const point<T> *, const point<T> *, void * ),
00252               void *fdat )
00253 {
00254   int *perm;
00255   bool inv = PatternTable[pat].nbits > 4;
00256   perm = VertexPerms[PatternTable[pat].perm];
00257   FUNC( MPT(0,2), MPT(1,5), MPT(1,3), NPT(0,2), NPT(1,5), NPT(1,3) );   
00258   FUNC( MPT(0,2), MPT(0,4), MPT(1,5), NPT(0,2), NPT(0,4), NPT(1,5) );   
00259   FUNC( MPT(6,7), MPT(3,7), MPT(5,7), NPT(6,7), NPT(3,7), NPT(5,7) );   
00260 }
00261 template< class T > void Draw6a( 
00262               unsigned char pat, 
00263               point<T> *med, 
00264               point<T> *norm,
00265               void (*func) (
00266                 bool, 
00267                 const point<T> *, const point<T> *, const point<T> *,
00268                 const point<T> *, const point<T> *, const point<T> *, void * ),
00269               void *fdat )
00270 {
00271   int *perm;
00272   bool inv = PatternTable[pat].nbits > 4;
00273   perm = VertexPerms[PatternTable[pat].perm];
00274   FUNC( MPT(0,2), MPT(3,7), MPT(1,3), NPT(0,2), NPT(3,7), NPT(1,3) );   
00275   FUNC( MPT(0,2), MPT(6,7), MPT(3,7), NPT(0,2), NPT(6,7), NPT(3,7) );   
00276   FUNC( MPT(0,2), MPT(0,4), MPT(6,7), NPT(0,2), NPT(0,4), NPT(6,7) );   
00277   FUNC( MPT(0,4), MPT(1,5), MPT(6,7), NPT(0,4), NPT(1,5), NPT(6,7) );   
00278   FUNC( MPT(1,5), MPT(5,7), MPT(6,7), NPT(1,5), NPT(5,7), NPT(6,7) );   
00279 }
00280 
00281 
00282 /* pattern 7
00283     6 7    0 0
00284    4 5    0 x
00285     2 3    0 x
00286    0 1    x 0
00287 */
00288 template< class T > void Draw7( 
00289               unsigned char pat, 
00290               point<T> *med, 
00291               point<T> *norm,
00292               void (*func) (
00293                 bool, 
00294                 const point<T> *, const point<T> *, const point<T> *,
00295                 const point<T> *, const point<T> *, const point<T> *, void * ),
00296               void *fdat )
00297 {
00298   int *perm;
00299   bool inv = PatternTable[pat].nbits > 4;
00300   perm = VertexPerms[PatternTable[pat].perm];
00301   FUNC( MPT(0,2), MPT(0,4), MPT(0,1), NPT(0,2), NPT(0,4), NPT(0,1)  );   
00302   FUNC( MPT(2,3), MPT(1,3), MPT(3,7), NPT(2,3), NPT(1,3), NPT(3,7)  );   
00303   FUNC( MPT(1,5), MPT(4,5), MPT(5,7), NPT(1,5), NPT(4,5), NPT(5,7)  );
00304 }
00305 template< class T > void Draw7a( 
00306               unsigned char pat, 
00307               point<T> *med, 
00308               point<T> *norm,
00309               void (*func) (
00310                 bool, 
00311                 const point<T> *, const point<T> *, const point<T> *,
00312                 const point<T> *, const point<T> *, const point<T> *, void * ),
00313               void *fdat )
00314 {
00315   int *perm;
00316   bool inv = PatternTable[pat].nbits > 4;
00317   perm = VertexPerms[PatternTable[pat].perm];
00318   FUNC( MPT(2,3), MPT(5,7), MPT(3,7), NPT(2,3), NPT(5,7), NPT(3,7)  );   
00319   FUNC( MPT(0,2), MPT(5,7), MPT(2,3), NPT(0,2), NPT(5,7), NPT(2,3)  );   
00320   FUNC( MPT(5,7), MPT(0,2), MPT(4,5), NPT(5,7), NPT(0,2), NPT(4,5)  );   
00321   FUNC( MPT(4,5), MPT(0,2), MPT(0,4), NPT(4,5), NPT(0,2), NPT(0,4)  );   
00322   FUNC( MPT(0,1), MPT(1,3), MPT(1,5), NPT(0,1), NPT(1,3), NPT(1,5)  );   
00323 }
00324 
00325 
00326 
00327 /* pattern 8
00328     6 7    0 0
00329    4 5    0 0
00330     2 3    x x
00331    0 1    x x
00332 */
00333 template< class T > void Draw8( 
00334               unsigned char pat, 
00335               point<T> *med, 
00336               point<T> *norm,
00337               void (*func) (
00338                 bool, 
00339                 const point<T> *, const point<T> *, const point<T> *,
00340                 const point<T> *, const point<T> *, const point<T> *, void * ),
00341               void *fdat )
00342 {
00343   int *perm;
00344   bool inv = PatternTable[pat].nbits > 4;
00345   perm = VertexPerms[PatternTable[pat].perm];
00346   FUNC( MPT(0,4), MPT(1,5), MPT(3,7), NPT(0,4), NPT(1,5), NPT(3,7) );   
00347   FUNC( MPT(0,4), MPT(3,7), MPT(2,6), NPT(0,4), NPT(3,7), NPT(2,6) );   
00348 }
00349 
00350 /* pattern 9
00351     6 7    0 0
00352    4 5    x 0
00353     2 3    x 0
00354    0 1    x x
00355 */
00356 template< class T > void Draw9( 
00357               unsigned char pat, 
00358               point<T> *med, 
00359               point<T> *norm,
00360               void (*func) (
00361                 bool, 
00362                 const point<T> *, const point<T> *, const point<T> *,
00363                 const point<T> *, const point<T> *, const point<T> *, void * ),
00364               void *fdat )
00365 {
00366   int *perm;
00367   bool inv = PatternTable[pat].nbits > 4;
00368   perm = VertexPerms[PatternTable[pat].perm];
00369   FUNC( MPT(4,5), MPT(1,5), MPT(1,3), NPT(4,5), NPT(1,5), NPT(1,3) );   
00370   FUNC( MPT(4,5), MPT(1,3), MPT(4,6), NPT(4,5), NPT(1,3), NPT(4,6) );   
00371   FUNC( MPT(4,6), MPT(1,3), MPT(2,3), NPT(4,6), NPT(1,3), NPT(2,3) );   
00372   FUNC( MPT(4,6), MPT(2,3), MPT(2,6), NPT(4,6), NPT(2,3), NPT(2,6) );   
00373 }
00374 
00375 /* pattern 10
00376     6 7    x x
00377    4 5    0 0
00378     2 3    0 0
00379    0 1    x x
00380 */
00381 template< class T > void Draw10( 
00382               unsigned char pat, 
00383               point<T> *med, 
00384               point<T> *norm,
00385               void (*func) (
00386                 bool, 
00387                 const point<T> *, const point<T> *, const point<T> *,
00388                 const point<T> *, const point<T> *, const point<T> *, void * ),
00389               void *fdat )
00390 {
00391   int *perm;
00392   bool inv = PatternTable[pat].nbits > 4;
00393   perm = VertexPerms[PatternTable[pat].perm];
00394   FUNC( MPT(0,4), MPT(1,3), MPT(0,2), NPT(0,4), NPT(1,3), NPT(0,2)  );   
00395   FUNC( MPT(0,4), MPT(1,5), MPT(1,3), NPT(0,4), NPT(1,5), NPT(1,3) );   
00396   FUNC( MPT(4,6), MPT(2,6), MPT(5,7), NPT(4,6), NPT(2,6), NPT(5,7)  );   
00397   FUNC( MPT(5,7), MPT(2,6), MPT(3,7), NPT(5,7), NPT(2,6), NPT(3,7)  );   
00398 }
00399 template< class T > void Draw10a( 
00400               unsigned char pat, 
00401               point<T> *med, 
00402               point<T> *norm,
00403               void (*func) (
00404                 bool, 
00405                 const point<T> *, const point<T> *, const point<T> *,
00406                 const point<T> *, const point<T> *, const point<T> *, void * ),
00407               void *fdat )
00408 {
00409   int *perm;
00410   bool inv = PatternTable[pat].nbits > 4;
00411   perm = VertexPerms[PatternTable[pat].perm];
00412   FUNC( MPT(0,2), MPT(2,6), MPT(3,7), NPT(0,2), NPT(2,6), NPT(3,7) );   
00413   FUNC( MPT(0,2), MPT(3,7), MPT(1,3), NPT(0,2), NPT(3,7), NPT(1,3) );   
00414   FUNC( MPT(4,6), MPT(0,4), MPT(5,7), NPT(4,6), NPT(0,4), NPT(5,7) );   
00415   FUNC( MPT(0,4), MPT(1,5), MPT(5,7), NPT(0,4), NPT(1,5), NPT(5,7) );   
00416 
00417 }
00418 
00419 
00420 /* pattern 11
00421     6 7    x 0
00422    4 5    0 0
00423     2 3    x 0
00424    0 1    x x
00425 */
00426 template< class T > void Draw11( 
00427               unsigned char pat, 
00428               point<T> *med, 
00429               point<T> *norm,
00430               void (*func) (
00431                 bool, 
00432                 const point<T> *, const point<T> *, const point<T> *,
00433                 const point<T> *, const point<T> *, const point<T> *, void * ),
00434               void *fdat )
00435 {
00436   int *perm;
00437   bool inv = PatternTable[pat].nbits > 4;
00438   perm = VertexPerms[PatternTable[pat].perm];
00439   FUNC( MPT(4,6), MPT(0,4), MPT(1,5), NPT(4,6), NPT(0,4), NPT(1,5) );
00440   FUNC( MPT(4,6), MPT(1,5), MPT(6,7), NPT(4,6), NPT(1,5), NPT(6,7) );   
00441   FUNC( MPT(6,7), MPT(1,5), MPT(1,3), NPT(6,7), NPT(1,5), NPT(1,3) );   
00442   FUNC( MPT(1,3), MPT(2,3), MPT(6,7), NPT(1,3), NPT(2,3), NPT(6,7) );   
00443 }
00444 
00445 
00446 /* pattern 12
00447     6 7    0 x
00448    4 5    0 0
00449     2 3    x 0
00450    0 1    x x
00451 */
00452 template< class T > void Draw12( 
00453               unsigned char pat, 
00454               point<T> *med, 
00455               point<T> *norm,
00456               void (*func) (
00457                 bool, 
00458                 const point<T> *, const point<T> *, const point<T> *,
00459                 const point<T> *, const point<T> *, const point<T> *, void * ),
00460               void *fdat )
00461 {
00462   int *perm;
00463   bool inv = PatternTable[pat].nbits > 4;
00464   perm = VertexPerms[PatternTable[pat].perm];
00465   FUNC( MPT(0,4), MPT(1,5), MPT(2,6), NPT(0,4), NPT(1,5), NPT(2,6) );   
00466   FUNC( MPT(2,6), MPT(1,3), MPT(2,3), NPT(2,6), NPT(1,3), NPT(2,3) );   
00467   FUNC( MPT(1,3), MPT(2,6), MPT(1,5), NPT(1,3), NPT(2,6), NPT(1,5) );   
00468   FUNC( MPT(6,7), MPT(3,7), MPT(5,7), NPT(6,7), NPT(3,7), NPT(5,7) );
00469 }
00470 template< class T > void Draw12a( 
00471               unsigned char pat, 
00472               point<T> *med, 
00473               point<T> *norm,
00474               void (*func) (
00475                 bool, 
00476                 const point<T> *, const point<T> *, const point<T> *,
00477                 const point<T> *, const point<T> *, const point<T> *, void * ),
00478               void *fdat )
00479 {
00480   int *perm;
00481   bool inv = PatternTable[pat].nbits > 4;
00482   perm = VertexPerms[PatternTable[pat].perm];
00483   FUNC( MPT(0,4), MPT(1,5), MPT(2,6), NPT(0,4), NPT(1,5), NPT(2,6) );   
00484   FUNC( MPT(1,3), MPT(2,3), MPT(3,7), NPT(1,3), NPT(2,3), NPT(3,7) );   
00485   FUNC( MPT(1,5), MPT(5,7), MPT(6,7), NPT(1,5), NPT(5,7), NPT(6,7) );   
00486   FUNC( MPT(6,7), MPT(2,6), MPT(1,5), NPT(6,7), NPT(2,6), NPT(1,5) );   
00487 }
00488 
00489 /* pattern 13
00490     6 7    x 0
00491    4 5    0 x
00492     2 3    0 x
00493    0 1    x 0
00494 */
00495 template< class T > void Draw13( 
00496               unsigned char pat, 
00497               point<T> *med, 
00498               point<T> *norm,
00499               void (*func) (
00500                 bool, 
00501                 const point<T> *, const point<T> *, const point<T> *,
00502                 const point<T> *, const point<T> *, const point<T> *, void * ),
00503               void *fdat )
00504 {
00505   int *perm;
00506   bool inv = PatternTable[pat].nbits > 4;
00507   perm = VertexPerms[PatternTable[pat].perm];
00508   FUNC( MPT(0,2), MPT(0,4), MPT(0,1), NPT(0,2), NPT(0,4), NPT(0,1) );   
00509   FUNC( MPT(2,3), MPT(1,3), MPT(3,7), NPT(2,3), NPT(1,3), NPT(3,7)  );   
00510   FUNC( MPT(4,5), MPT(5,7), MPT(1,5), NPT(4,5), NPT(5,7), NPT(1,5) );   
00511   FUNC( MPT(4,6), MPT(2,6), MPT(6,7), NPT(4,6), NPT(2,6), NPT(6,7) );   
00512 }
00513 template< class T > void Draw13a( 
00514               unsigned char pat, 
00515               point<T> *med, 
00516               point<T> *norm,
00517               void (*func) (
00518                 bool, 
00519                 const point<T> *, const point<T> *, const point<T> *,
00520                 const point<T> *, const point<T> *, const point<T> *, void * ),
00521               void *fdat )
00522 {
00523   int *perm;
00524   bool inv = PatternTable[pat].nbits > 4;
00525   perm = VertexPerms[PatternTable[pat].perm];
00526   FUNC( MPT(0,1), MPT(1,3), MPT(1,5), NPT(0,1), NPT(1,3), NPT(1,5) );   
00527   FUNC( MPT(0,2), MPT(2,6), MPT(2,3), NPT(0,2), NPT(2,6), NPT(2,3) );   
00528   FUNC( MPT(4,6), MPT(0,4), MPT(4,5), NPT(4,6), NPT(0,4), NPT(4,5) );   
00529   FUNC( MPT(6,7), MPT(5,7), MPT(3,7), NPT(6,7), NPT(5,7), NPT(3,7) );   
00530 }
00531 
00532 /* pattern 14
00533     6 7    0 0
00534    4 5    0 x
00535     2 3    x 0
00536    0 1    x x
00537 */
00538 template< class T > void Draw14( 
00539               unsigned char pat, 
00540               point<T> *med, 
00541               point<T> *norm,
00542               void (*func) (
00543                 bool, 
00544                 const point<T> *, const point<T> *, const point<T> *,
00545                 const point<T> *, const point<T> *, const point<T> *, void * ),
00546               void *fdat )
00547 {
00548   int *perm;
00549   bool inv = PatternTable[pat].nbits > 4;
00550   perm = VertexPerms[PatternTable[pat].perm];
00551   FUNC( MPT(5,7), MPT(2,6), MPT(0,4), NPT(5,7), NPT(2,6), NPT(0,4) );   
00552   FUNC( MPT(5,7), MPT(2,3), MPT(2,6), NPT(5,7), NPT(2,3), NPT(2,6) );   
00553   FUNC( MPT(5,7), MPT(1,3), MPT(2,3), NPT(5,7), NPT(1,3), NPT(2,3) );   
00554   FUNC( MPT(0,4), MPT(4,5), MPT(5,7), NPT(0,4), NPT(4,5), NPT(5,7) );   
00555 }
00556 
00557 
00558 template< class T >
00559 class MethodTable
00560 {
00561 public:
00562   typedef 
00563     void (*func_t)( 
00564       unsigned char, 
00565       point<T> *, 
00566       point<T> *,
00567       void (*) ( bool, 
00568         const point<T> *, const point<T> *, const point<T> *,
00569         const point<T> *, const point<T> *, const point<T> *, void * ),
00570       void *);
00571 
00572   static func_t m_table[15];
00573 
00574   static func_t get(size_t i) { return m_table[i]; } 
00575 };
00576 template class MethodTable<float>;
00577 template class MethodTable<double>;
00578 
00579 template<>
00580 MethodTable<float>::func_t MethodTable<float>::m_table[15] =
00581 {
00582   0,
00583   Draw1<float>,  
00584   Draw2<float>,  
00585   Draw3<float>,  
00586   Draw4<float>,  
00587   Draw5<float>,  
00588   Draw6<float>,  
00589   Draw7<float>,  
00590   Draw8<float>,  
00591   Draw9<float>,  
00592   Draw10<float>,  
00593   Draw11<float>,  
00594   Draw12<float>,  
00595   Draw13<float>,  
00596   Draw14<float>
00597 }; 
00598 template<>
00599 MethodTable<double>::func_t MethodTable<double>::m_table[15] =
00600 {
00601   0,
00602   Draw1<double>,  
00603   Draw2<double>,  
00604   Draw3<double>,  
00605   Draw4<double>,  
00606   Draw5<double>,  
00607   Draw6<double>,  
00608   Draw7<double>,  
00609   Draw8<double>,  
00610   Draw9<double>,  
00611   Draw10<double>,  
00612   Draw11<double>,  
00613   Draw12<double>,  
00614   Draw13<double>,  
00615   Draw14<double>
00616 }; 
00617 
00618 template< class T >
00619 class AltMethodTable
00620 {
00621 public:
00622   typedef 
00623     void (*func_t)( 
00624       unsigned char, 
00625       point<T> *, 
00626       point<T> *,
00627       void (*) ( bool, 
00628         const point<T> *, const point<T> *, const point<T> *,
00629         const point<T> *, const point<T> *, const point<T> *, void * ),
00630       void *);
00631 
00632   static func_t m_table[15];
00633 
00634   static func_t get(size_t i) { return m_table[i]; } 
00635 };
00636 template class AltMethodTable<float>;
00637 template class AltMethodTable<double>;
00638 
00639 template<>
00640 AltMethodTable<float>::func_t AltMethodTable<float>::m_table[15] =
00641 {
00642   0,
00643   Draw1<float>,  
00644   Draw2<float>,  
00645   Draw3a<float>,  
00646   Draw4<float>,  
00647   Draw5<float>,  
00648   Draw6a<float>,  
00649   Draw7a<float>,  
00650   Draw8<float>,  
00651   Draw9<float>,  
00652   Draw10a<float>,  
00653   Draw11<float>,  
00654   Draw12a<float>,  
00655   Draw13a<float>,  
00656   Draw14<float>  
00657 };
00658 template<>
00659 AltMethodTable<double>::func_t AltMethodTable<double>::m_table[15] =
00660 {
00661   0,
00662   Draw1<double>,  
00663   Draw2<double>,  
00664   Draw3a<double>,  
00665   Draw4<double>,  
00666   Draw5<double>,  
00667   Draw6a<double>,  
00668   Draw7a<double>,  
00669   Draw8<double>,  
00670   Draw9<double>,  
00671   Draw10a<double>,  
00672   Draw11<double>,  
00673   Draw12a<double>,  
00674   Draw13a<double>,  
00675   Draw14<double>  
00676 };
00677 
00678 /*-----------------------------------------------------------------------------
00679   Process the first plane of function values for Marching Cube algorithm
00680 -----------------------------------------------------------------------------*/
00681 template< class T >
00682 void FirstPlane( 
00683          int                        nx,
00684          int                        ny,
00685          Ptr< Ptr<T> >              F0,       /* for getting function values  */
00686          Ptr< Ptr<unsigned char> >  Pat1      /* for storing patterns         */
00687        )
00688 {
00689   T f;
00690   int ix,iy;
00691 
00692   /* --- initial code for iy=0,ix=0 --- */
00693 
00694   Pat1[0][0] = 0x00;   /*initialize new pattern */
00695 
00696   f = F0[0][0]; /* get function value */
00697   if( f < 0.0 )
00698   {
00699     Pat1[0][0] |= 0x01;     /* influenced patterns */
00700   }
00701   /* --- initial code for iy=0 and ix = 1,..,nx ---- */
00702 
00703   for( ix=1; ix<=nx; ix++ ) 
00704   {
00705     Pat1[0][ix] = 0x00;   /*initialize new pattern */
00706 
00707     f = F0[0][ix]; /* get function value */
00708     if( f < 0.0 )
00709     {
00710       Pat1[0][ix] |= 0x01;     /* influenced patterns */
00711       Pat1[0][ix-1] |= 0x02;
00712     }
00713   }
00714   for( iy=1; iy<=ny; iy++ ) 
00715   {
00716     /* --- initial code for ix = 0 --- */
00717     
00718     Pat1[iy][0] = 0x00;   /*initialize new pattern */
00719  
00720     f = F0[iy][0]; /* get function value */
00721     if( f < 0.0 )
00722     { 
00723       Pat1[iy][0] |= 0x01;     /* influenced patterns */
00724       Pat1[iy-1][0] |= 0x04;
00725     }  
00726     /* ---- code for ix=1,..,nx ---- */
00727     
00728     for( ix=1; ix<=nx; ix++ ) 
00729     {
00730       Pat1[iy][ix] = 0x00;   /*initialize new pattern */
00731       
00732       f = F0[iy][ix]; /* get function value */            
00733       if( f < 0.0 )
00734       {
00735         Pat1[iy][ix] |= 0x01;     /* influenced patterns */
00736         Pat1[iy][ix-1] |= 0x02;
00737         Pat1[iy-1][ix] |= 0x04;
00738         Pat1[iy-1][ix-1] |= 0x08;
00739       }
00740     }      
00741   }
00742 }
00743 
00744 /*-----------------------------------------------------------------------------
00745   Process the next plane of function values for Marching Cube algorithm
00746 -----------------------------------------------------------------------------*/
00747 template< class T >
00748 void NextPlane(
00749          int       nx,
00750          int       ny,
00751          Ptr< Ptr<T> >              Fi,       /* for getting function values  */
00752          Ptr< Ptr<unsigned char> >  Pat0,     /* for storing patterns         */
00753          Ptr< Ptr<unsigned char> >  Pat1      /* for storing patterns         */
00754        )
00755 {
00756   T f;
00757   int ix,iy;
00758 
00759   /* --- initial code for iy=0,ix=0 --- */
00760 
00761   Pat1[0][0] = 0x00;   /*initialize new pattern */
00762 
00763   f = Fi[0][0]; /* get function value */
00764   if( f < 0.0 )
00765   {
00766     Pat1[0][0] |= 0x01;     /* influenced patterns */
00767   }
00768   /* --- initial code for iy=0 and ix = 1,..,nx ---- */
00769 
00770   for( ix=1; ix<=nx; ix++ ) 
00771   {
00772     Pat1[0][ix] = 0x00;   /*initialize new pattern */
00773 
00774     f = Fi[0][ix]; /* get function value */
00775     if( f < 0.0 )
00776     {
00777       Pat1[0][ix] |= 0x01;     /* influenced patterns */
00778       Pat1[0][ix-1] |= 0x02;
00779       Pat0[0][ix] |= 0x10;
00780       Pat0[0][ix-1] |= 0x20;
00781     }
00782   }
00783   for( iy=1; iy<=ny; iy++ ) 
00784   {
00785     /* --- initial code for ix = 0 --- */
00786     
00787     Pat1[iy][0] = 0x00;   /*initialize new pattern */
00788  
00789     f = Fi[iy][0]; /* get function value */
00790     if( f < 0.0 )
00791     { 
00792       Pat1[iy][0] |= 0x01;     /* influenced patterns */
00793       Pat1[iy-1][0] |= 0x04;
00794       Pat0[iy][0] |= 0x10;
00795       Pat0[iy-1][0] |= 0x40;
00796     }
00797     /* ---- code for ix=1,..,nx ---- */
00798     
00799     for( ix=1; ix<=nx; ix++ ) 
00800     {
00801       Pat1[iy][ix] = 0x00;   /*initialize new pattern */
00802       
00803       f = Fi[iy][ix]; /* get function value */
00804       if( f < 0.0 )
00805       {
00806         Pat1[iy][ix] |= 0x01;     /* influenced patterns */
00807         Pat1[iy][ix-1] |= 0x02;
00808         Pat1[iy-1][ix] |= 0x04;
00809         Pat1[iy-1][ix-1] |= 0x08;
00810         Pat0[iy][ix] |= 0x10;
00811         Pat0[iy][ix-1] |= 0x20;
00812         Pat0[iy-1][ix] |= 0x40;
00813         Pat0[iy-1][ix-1] |= 0x80;
00814       }
00815     }      
00816   }
00817 }
00818 
00819 /*------------------------------------------------------------------------
00820   Rasterize the function F(x,y,z), and draw the surface F(x,y,z)=0,
00821   using Marching Cube algorithm
00822 ------------------------------------------------------------------------*/
00823 
00824 void Harmonize( 
00825       int k, int j, int i, int nx, int ny, int nz,
00826       Ptr< Ptr< Ptr<unsigned char> > > flags )
00827 {
00828   if( (flags[i][j][k]&(DRAW_ALT|CON_SPECIAL)) != 0 )
00829     return;
00830 
00831   flags[i][j][k] |= DRAW_ALT;
00832   
00833   if( ((flags[i][j][k]&CON_X1)!=0) && (k>0) )
00834     Harmonize( k-1, j, i, nx, ny, nz, flags );
00835 
00836   if( ((flags[i][j][k]&CON_X2)!=0) && (k<nx-1) )
00837     Harmonize( k+1, j, i, nx, ny, nz, flags );
00838 
00839   if( ((flags[i][j][k]&CON_Y1)!=0) && (j>0) )
00840     Harmonize( k, j-1, i, nx, ny, nz, flags );
00841       
00842   if( ((flags[i][j][k]&CON_Y2)!=0) && (j<ny-1) )
00843     Harmonize( k, j+1, i, nx, ny, nz, flags );
00844 
00845   if( ((flags[i][j][k]&CON_Z1)!=0) && (i>0) )
00846     Harmonize( k, j, i-1, nx, ny, nz, flags );
00847 
00848   if( ((flags[i][j][k]&CON_Z2)!=0) && (i<nz-1) )
00849     Harmonize( k, j, i+1, nx, ny, nz, flags );
00850 }  
00851 
00852 template< class T >
00853 GULAPI void MarchingCube( 
00854         T    x0,
00855         T    y0,
00856         T    z0,
00857         T    dx,
00858         T    dy,
00859         T    dz,
00860         int  nx,
00861         int  ny,
00862         int  nz,
00863         Ptr< Ptr< Ptr<T> > > F, 
00864         void (*trifunc)(
00865                  const bool,
00866                  const point<T> *, const point<T> *, const point<T> *, 
00867                  const point<T> *, const point<T> *, const point<T> *, void *),
00868         void *tridata
00869      )
00870 {
00871   Ptr< Ptr< Ptr<unsigned char> > > Pat,Flags;
00872   unsigned char pat,pat0; 
00873   int i,j,k,iz;
00874   T dx2,dy2,dz2,bx0,by0,bz0,bx1,by1,bz1,f,a;
00875   point<T> M[12],G[12],D[8];
00876   unsigned char ders;
00877 
00878   dx2 = dx+dx;
00879   dy2 = dy+dy;
00880   dz2 = dz+dz;
00881 
00882   Pat.reserve_pool(nz+1);
00883   Flags.reserve_pool(nz+1);
00884 
00885   for( i = 0; i <= nz; i++ )
00886   {
00887     Pat[i].reserve_pool(ny+1);
00888     Flags[i].reserve_pool(ny+1);
00889 
00890     for( j = 0; j <= ny; j++ )
00891     {
00892       Pat[i][j].reserve_pool(nx+1);
00893       Flags[i][j].reserve_pool(nx+1);
00894     }
00895   }
00896 
00897   FirstPlane( nx, ny, F[0], Pat[0] );
00898 
00899   for( i=1; i <= nz; i++ )
00900   {
00901     NextPlane( nx, ny, F[i], Pat[i-1], Pat[i] );
00902   }  
00903 
00904   for( k = 0; k < nz; k++ )
00905     for( i = 0; i < ny; i++ )
00906       for( j = 0; j < nx; j++ )
00907         Flags[k][i][j] = 0x00;
00908 
00909   for( k = 1; k < nz; k++ )
00910   {
00911     for( i = 1; i < ny; i++ )
00912     {
00913       for( j = 1; j < nx; j++ )
00914       {
00915         pat = Pat[k][i][j];
00916 
00917         if( ((pat & 0x33)==0x21)||((pat & 0x33)==0x12) )
00918         {
00919           Flags[k][i][j] |= CON_Y1;
00920           Flags[k][i-1][j] |= CON_Y2;
00921 
00922           pat0 = Pat[k][i-1][j];
00923           if(
00924       ((PatternTable[pat].nbits>4)&&(PatternTable[pat0].nbits<=4))||
00925       ((PatternTable[pat].nbits<=4)&&(PatternTable[pat0].nbits>4))
00926           )
00927           {
00928             Flags[k][i][j] |= CON_SPECIAL;
00929             Flags[k][i-1][j] |= CON_SPECIAL;
00930 
00931             if( (k&1) )
00932             {
00933               if( ((!(i&1)) && (!(j&1))) ||
00934                   (((i&1)) && ((j&1))) )
00935                 Flags[k][i][j] |= DRAW_ALT;
00936               else
00937                 Flags[k][i-1][j] |= DRAW_ALT;
00938             }
00939             else
00940             {
00941               if( ((!(i&1)) && (!(j&1))) ||
00942                   ((i&1) && (j&1)) )
00943                 Flags[k][i-1][j] |= DRAW_ALT;
00944               else
00945                 Flags[k][i][j] |= DRAW_ALT;
00946             }                                    
00947           }
00948         }
00949 
00950         if( ((pat & 0x55)==0x41)||((pat & 0x55)==0x14) )
00951         {
00952           Flags[k][i][j] |= CON_X1;
00953           Flags[k][i][j-1] |= CON_X2;
00954 
00955           pat0 = Pat[k][i][j-1];
00956           if(
00957       ((PatternTable[pat].nbits>4)&&(PatternTable[pat0].nbits<=4))||
00958       ((PatternTable[pat].nbits<=4)&&(PatternTable[pat0].nbits>4))
00959           )
00960           {
00961             Flags[k][i][j] |= CON_SPECIAL;
00962             Flags[k][i][j-1] |= CON_SPECIAL;
00963 
00964             if( (k&1) )
00965             {
00966               if( ((!(i&1)) && (!(j&1))) ||
00967                   ((i&1) && (j&1)) )
00968                 Flags[k][i][j] |= DRAW_ALT;
00969               else
00970                 Flags[k][i][j-1] |= DRAW_ALT;
00971             }
00972             else
00973             {
00974               if( ((!(i&1)) && (!(j&1))) ||
00975                   ((i&1) && (j&1)) )
00976                 Flags[k][i][j-1] |= DRAW_ALT;
00977               else
00978                 Flags[k][i][j] |= DRAW_ALT;
00979             }
00980           }
00981         }
00982 
00983         if( ((pat & 0x0f)==0x09)||((pat & 0x0f)==0x06) )
00984         {
00985           Flags[k][i][j] |= CON_Z1;
00986           Flags[k-1][i][j] |= CON_Z2;
00987 
00988           pat0 = Pat[k-1][i][j];
00989           if(
00990       ((PatternTable[pat].nbits>4)&&(PatternTable[pat0].nbits<=4))||
00991       ((PatternTable[pat].nbits<=4)&&(PatternTable[pat0].nbits>4))
00992           )
00993           {
00994             Flags[k][i][j] |= CON_SPECIAL;
00995             Flags[k-1][i][j] |= CON_SPECIAL;
00996 
00997             if( (k&1) )
00998             {
00999               if( ((!(i&1)) && (!(j&1))) ||
01000                   ((i&1) && (j&1)) )
01001                 Flags[k][i][j] |= DRAW_ALT;
01002               else
01003                 Flags[k-1][i][j] |= DRAW_ALT;
01004             }
01005             else
01006             {
01007               if( ((!(i&1)) && (!(j&1))) ||
01008                   ((i&1) && (j&1)) )
01009                 Flags[k-1][i][j] |= DRAW_ALT;
01010               else
01011                 Flags[k][i][j] |= DRAW_ALT;
01012             }
01013           }
01014         }
01015 
01016 
01017       }
01018     }
01019   }
01020 
01021   for( i = 0; i < nz; i++ )
01022     for( j = 0; j < ny; j++ )
01023       for( k = 0; k < nx; k++ )
01024       {
01025         if( ((Flags[i][j][k]&CON_X1)!=0) && (k>0) )
01026           Harmonize( k-1, j, i, nx, ny, nz, Flags );
01027 
01028         if( ((Flags[i][j][k]&CON_X2)!=0) && (k<nx-1) )
01029           Harmonize( k+1, j, i, nx, ny, nz, Flags );
01030  
01031         if( ((Flags[i][j][k]&CON_Y1)!=0) && (j>0) )
01032           Harmonize( k, j-1, i, nx, ny, nz, Flags );
01033       
01034         if( ((Flags[i][j][k]&CON_Y2)!=0) && (j<ny-1) )
01035           Harmonize( k, j+1, i, nx, ny, nz, Flags );
01036 
01037         if( ((Flags[i][j][k]&CON_Z1)!=0) && (i>0) )
01038           Harmonize( k, j, i-1, nx, ny, nz, Flags );
01039 
01040         if( ((Flags[i][j][k]&CON_Z2)!=0) && (i<nz-1) )
01041           Harmonize( k, j, i+1, nx, ny, nz, Flags );
01042       }
01043 
01044 
01045   for( iz=1,bz0=z0+dz; iz < nz-1; iz++,bz0+=dz )
01046   {
01047     for( i=1,by0=y0+dy; i < ny-1; i++,by0+=dy )
01048     {
01049       for( j=1,bx0=x0+dx; j < nx-1; j++,bx0+=dx )
01050       {
01051         pat = Pat[iz][i][j];
01052         
01053         if( (pat != 0x00) && (pat != 0xff) )
01054         {
01055           bx1 = bx0+dx;
01056           by1 = by0+dy;
01057           bz1 = bz0+dz; 
01058           ders = 0;
01059 
01060           if( ((pat & 0x03) == 0x01) || ((pat & 0x03) == 0x02) )
01061           {
01062             DER( D[0], iz, i, j );
01063             DER( D[1], iz, i, j+1 );
01064             ders |= 0x03;
01065 
01066             f=F[iz][i][j+1];      /* intersection */  
01067             a=F[iz][i][j];        
01068             VEC( M[0], bx0+a*dx/(a-f), by0, bz0 );
01069 
01070             G[0].x = D[0].x + (M[0].x-bx0)*(D[1].x-D[0].x)/dx;
01071             G[0].y = D[0].y + (M[0].x-bx0)*(D[1].y-D[0].y)/dx;
01072             G[0].z = D[0].z + (M[0].x-bx0)*(D[1].z-D[0].z)/dx;
01073             normalize(&G[0],G[0]);
01074           }
01075           if( ((pat & 0x05) == 0x01) || ((pat & 0x05) == 0x04) )
01076           {
01077             if( (ders & 0x01) == 0 )
01078               DER( D[0], iz, i, j );
01079             if( (ders & 0x04) == 0 )
01080               DER( D[2], iz, i+1, j );
01081             ders |= 0x05;
01082 
01083             f=F[iz][i+1][j];      /* intersection */  
01084             a=F[iz][i][j];        
01085             VEC( M[2], bx0, by0+a*dy/(a-f), bz0 );
01086             
01087             G[2].x = D[0].x + (M[2].y-by0)*(D[2].x-D[0].x)/dy;
01088             G[2].y = D[0].y + (M[2].y-by0)*(D[2].y-D[0].y)/dy;
01089             G[2].z = D[0].z + (M[2].y-by0)*(D[2].z-D[0].z)/dy;
01090             normalize(&G[2],G[2]);
01091           }
01092           if( ((pat & 0x0a) == 0x02) || ((pat & 0x0a) == 0x08) )
01093           {
01094             if( (ders & 0x02) == 0 )
01095               DER( D[1], iz, i, j+1 );
01096             if( (ders & 0x08) == 0 )
01097               DER( D[3], iz, i+1, j+1 );
01098             ders |= 0x0a;
01099               
01100             f=F[iz][i+1][j+1];      /* intersection */  
01101             a=F[iz][i][j+1];        
01102             VEC( M[3], bx1, by0+a*dy/(a-f), bz0 );
01103 
01104             G[3].x = D[1].x + (M[3].y-by0)*(D[3].x-D[1].x)/dy;
01105             G[3].y = D[1].y + (M[3].y-by0)*(D[3].y-D[1].y)/dy;
01106             G[3].z = D[1].z + (M[3].y-by0)*(D[3].z-D[1].z)/dy;
01107             normalize(&G[3],G[3]);
01108           }
01109           if( ((pat & 0x0c) == 0x04) || ((pat & 0x0c) == 0x08) )
01110           {
01111             if( (ders & 0x04) == 0 )
01112               DER( D[2], iz, i+1, j );
01113             if( (ders & 0x08) == 0 )
01114               DER( D[3], iz, i+1, j+1 );
01115             ders |= 0x0c;
01116 
01117             f=F[iz][i+1][j+1];      /* intersection */  
01118             a=F[iz][i+1][j];        
01119             VEC( M[1], bx0+a*dx/(a-f), by1, bz0 );
01120               
01121             G[1].x = D[2].x + (M[1].x-bx0)*(D[3].x-D[2].x)/dx;
01122             G[1].y = D[2].y + (M[1].x-bx0)*(D[3].y-D[2].y)/dx;
01123             G[1].z = D[2].z + (M[1].x-bx0)*(D[3].z-D[2].z)/dx;
01124             normalize(&G[1],G[1]);
01125           }
01126 
01127 
01128           if( ((pat & 0x30) == 0x10) || ((pat & 0x30) == 0x20) )
01129           {
01130             if( (ders & 0x10) == 0 )
01131               DER( D[4], iz+1, i, j );
01132             if( (ders & 0x20) == 0 )
01133               DER( D[5], iz+1, i, j+1 );
01134             ders |= 0x30;
01135              
01136             f=F[iz+1][i][j+1];      /* intersection */  
01137             a=F[iz+1][i][j];        
01138             VEC( M[4], bx0+a*dx/(a-f), by0, bz1 );
01139 
01140             G[4].x = D[4].x + (M[4].x-bx0)*(D[5].x-D[4].x)/dx;
01141             G[4].y = D[4].y + (M[4].x-bx0)*(D[5].y-D[4].y)/dx;
01142             G[4].z = D[4].z + (M[4].x-bx0)*(D[5].z-D[4].z)/dx;
01143             normalize(&G[4],G[4]);
01144           }
01145           if( ((pat & 0x50) == 0x10) || ((pat & 0x50) == 0x40) )
01146           {
01147             if( (ders & 0x10) == 0 )
01148               DER( D[4], iz+1, i, j );
01149             if( (ders & 0x40) == 0 )
01150               DER( D[6], iz+1, i+1, j );
01151             ders |= 0x50;
01152               
01153             f=F[iz+1][i+1][j];      /* intersection */  
01154             a=F[iz+1][i][j];        
01155             VEC( M[6], bx0, by0+a*dy/(a-f), bz1 );
01156 
01157             G[6].x = D[4].x + (M[6].y-by0)*(D[6].x-D[4].x)/dy;
01158             G[6].y = D[4].y + (M[6].y-by0)*(D[6].y-D[4].y)/dy; 
01159             G[6].z = D[4].z + (M[6].y-by0)*(D[6].z-D[4].z)/dy; 
01160             normalize(&G[6],G[6]);
01161           }
01162           if( ((pat & 0xa0) == 0x20) || ((pat & 0xa0) == 0x80) )
01163           {
01164             if( (ders & 0x20) == 0 )
01165               DER( D[5], iz+1, i, j+1 );
01166             if( (ders & 0x80) == 0 )
01167               DER( D[7], iz+1, i+1, j+1 );
01168             ders |= 0xa0;
01169               
01170             f=F[iz+1][i+1][j+1];     /* intersection */  
01171             a=F[iz+1][i][j+1];        
01172             VEC( M[7], bx1, by0+a*dy/(a-f), bz1 );
01173 
01174             G[7].x = D[5].x + (M[7].y-by0)*(D[7].x-D[5].x)/dy;  
01175             G[7].y = D[5].y + (M[7].y-by0)*(D[7].y-D[5].y)/dy;  
01176             G[7].z = D[5].z + (M[7].y-by0)*(D[7].z-D[5].z)/dy;  
01177             normalize(&G[7],G[7]);
01178           }
01179           if( ((pat & 0xc0) == 0x40) || ((pat & 0xc0) == 0x80) )
01180           {
01181             if( (ders & 0x40) == 0 )
01182               DER( D[6], iz+1, i+1, j );
01183             if( (ders & 0x80) == 0 )
01184               DER( D[7], iz+1, i+1, j+1 );
01185             ders |= 0xc0;
01186               
01187             f=F[iz+1][i+1][j+1];      /* intersection */  
01188             a=F[iz+1][i+1][j];        
01189             VEC( M[5], bx0+a*dx/(a-f), by1, bz1 );
01190             
01191             G[5].x = D[6].x + (M[5].x-bx0)*(D[7].x-D[6].x)/dx;  
01192             G[5].y = D[6].y + (M[5].x-bx0)*(D[7].y-D[6].y)/dx;  
01193             G[5].z = D[6].z + (M[5].x-bx0)*(D[7].z-D[6].z)/dx;  
01194             normalize(&G[5],G[5]);
01195           }
01196 
01197 
01198           if( ((pat & 0x11) == 0x01) || ((pat & 0x11) == 0x10) )
01199           {
01200             if( (ders & 0x01) == 0 )
01201               DER( D[0], iz, i, j );
01202             if( (ders & 0x10) == 0 )
01203               DER( D[4], iz+1, i, j );
01204             ders |= 0x11;
01205      
01206             f=F[iz+1][i][j];      /* intersection */  
01207             a=F[iz][i][j];        
01208             VEC( M[8], bx0, by0, bz0+a*dz/(a-f) );
01209 
01210             G[8].x = D[0].x + (M[8].z-bz0)*(D[4].x-D[0].x)/dz;   
01211             G[8].y = D[0].y + (M[8].z-bz0)*(D[4].y-D[0].y)/dz;   
01212             G[8].z = D[0].z + (M[8].z-bz0)*(D[4].z-D[0].z)/dz;   
01213             normalize(&G[8],G[8]);
01214           }
01215           if( ((pat & 0x22) == 0x02) || ((pat & 0x22) == 0x20) )
01216           {
01217             if( (ders & 0x02) == 0 )
01218               DER( D[1], iz, i, j+1 );
01219             if( (ders & 0x20) == 0 )
01220               DER( D[5], iz+1, i, j+1 );
01221             ders |= 0x22;
01222 
01223             f=F[iz+1][i][j+1];      /* intersection */  
01224             a=F[iz][i][j+1];        
01225             VEC( M[9], bx1, by0, bz0+a*dz/(a-f) );
01226               
01227             G[9].x = D[1].x + (M[9].z-bz0)*(D[5].x-D[1].x)/dz;   
01228             G[9].y = D[1].y + (M[9].z-bz0)*(D[5].y-D[1].y)/dz;   
01229             G[9].z = D[1].z + (M[9].z-bz0)*(D[5].z-D[1].z)/dz;   
01230             normalize(&G[9],G[9]);
01231           }
01232           if( ((pat & 0x44) == 0x04) || ((pat & 0x44) == 0x40) )
01233           {
01234             if( (ders & 0x04) == 0 )
01235               DER( D[2], iz, i+1, j );
01236             if( (ders & 0x40) == 0 )
01237               DER( D[6], iz+1, i+1, j );
01238             ders |= 0x44;
01239               
01240             f=F[iz+1][i+1][j];      /* intersection */  
01241             a=F[iz][i+1][j];        
01242             VEC( M[10], bx0, by1, bz0+a*dz/(a-f) );
01243 
01244             G[10].x = D[2].x + (M[10].z-bz0)*(D[6].x-D[2].x)/dz;   
01245             G[10].y = D[2].y + (M[10].z-bz0)*(D[6].y-D[2].y)/dz;   
01246             G[10].z = D[2].z + (M[10].z-bz0)*(D[6].z-D[2].z)/dz;   
01247             normalize(&G[10],G[10]);
01248           }
01249           if( ((pat & 0x88) == 0x08) || ((pat & 0x88) == 0x80) )
01250           {
01251             if( (ders & 0x08) == 0 )
01252               DER( D[3], iz, i+1, j+1 );
01253             if( (ders & 0x80) == 0 )
01254               DER( D[7], iz+1, i+1, j+1 );
01255             ders |= 0x88;
01256 
01257             f=F[iz+1][i+1][j+1];      /* intersection */  
01258             a=F[iz][i+1][j+1];        
01259             VEC( M[11], bx1, by1, bz0+a*dz/(a-f) );
01260               
01261             G[11].x = D[3].x + (M[11].z-bz0)*(D[7].x-D[3].x)/dz;   
01262             G[11].y = D[3].y + (M[11].z-bz0)*(D[7].y-D[3].y)/dz;   
01263             G[11].z = D[3].z + (M[11].z-bz0)*(D[7].z-D[3].z)/dz;   
01264             normalize(&G[11],G[11]);
01265           }
01266 
01267           if( (Flags[iz][i][j]&DRAW_ALT)!=0 )
01268           {
01269             if( AltMethodTable<T>::get(PatternTable[pat].method) != 0 )       
01270               AltMethodTable<T>::get(PatternTable[pat].method)
01271                                            ( pat, M, G, trifunc, tridata );
01272           }
01273           else
01274           {
01275             if( MethodTable<T>::get(PatternTable[pat].method) != NULL )       
01276               MethodTable<T>::get(PatternTable[pat].method)
01277                                            ( pat, M, G, trifunc, tridata );
01278           }
01279         }
01280       }
01281     }
01282   }
01283 } 
01284 template GULAPI void MarchingCube( 
01285         float    x0,
01286         float    y0,
01287         float    z0,
01288         float    dx,
01289         float    dy,
01290         float    dz,
01291         int  nx,
01292         int  ny,
01293         int  nz,
01294         Ptr< Ptr< Ptr<float> > > F, 
01295         void (*trifunc)(
01296      const bool,
01297      const point<float> *, const point<float> *, const point<float> *, 
01298      const point<float> *, const point<float> *, const point<float> *, void *),
01299         void *tridata
01300      );
01301 template GULAPI void MarchingCube( 
01302         double    x0,
01303         double    y0,
01304         double    z0,
01305         double    dx,
01306         double    dy,
01307         double    dz,
01308         int  nx,
01309         int  ny,
01310         int  nz,
01311         Ptr< Ptr< Ptr<double> > > F, 
01312         void (*trifunc)(
01313   const bool,
01314   const point<double> *, const point<double> *, const point<double> *, 
01315   const point<double> *, const point<double> *, const point<double> *, void *),
01316         void *tridata
01317      );
01318 }
Generated on Mon Jan 21 04:17:27 2002 for GUL 0.6 - Geometry Utility Library by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001