svn merge -r 12716:12856 https://svn.blender.org/svnroot/bf-blender/trunk/blender
[blender.git] / intern / elbeem / intern / isosurface.cpp
index 9925565b85d67f55fd6785d172db95ad5b33db9b..14aa22d06b59d4f4c1640004bc73c1a10baec028 100644 (file)
 #include <algorithm>
 #include <stdio.h>
 
 #include <algorithm>
 #include <stdio.h>
 
+#if !defined(linux) && (defined (__sparc) || defined (__sparc__))
+#include <ieeefp.h>
+#endif
+
+
 // just use default rounding for platforms where its not available
 #ifndef round
 #define round(x) (x)
 #endif
 
 // just use default rounding for platforms where its not available
 #ifndef round
 #define round(x) (x)
 #endif
 
+#if PARALLEL==1        
+#include <omp.h>
+#endif
+
 /******************************************************************************
  * Constructor
  *****************************************************************************/
 /******************************************************************************
  * Constructor
  *****************************************************************************/
@@ -156,13 +165,6 @@ void IsoSurface::triangulate( void )
                mpEdgeVerticesZ[i] = -1;
        }
 
                mpEdgeVerticesZ[i] = -1;
        }
 
-       ntlVec3Gfx pos[8];
-       float value[8];
-       int cubeIndex;      // index entry of the cube 
-       int triIndices[12]; // vertex indices 
-       int *eVert[12];
-       IsoLevelVertex ilv;
-
        // edges between which points?
        const int mcEdges[24] = { 
                0,1,  1,2,  2,3,  3,0,
        // edges between which points?
        const int mcEdges[24] = { 
                0,1,  1,2,  2,3,  3,0,
@@ -189,7 +191,12 @@ void IsoSurface::triangulate( void )
                                px = mStart[0]-gsx*0.5;
                                for(int i=1;i<(mSizex-2);i++) {
                                        px += gsx;
                                px = mStart[0]-gsx*0.5;
                                for(int i=1;i<(mSizex-2);i++) {
                                        px += gsx;
-
+                                       int cubeIndex;      // index entry of the cube 
+                                       float value[8];
+                                       int triIndices[12]; // vertex indices 
+                                       int *eVert[12];
+                                       IsoLevelVertex ilv;
+                                       
                                        value[0] = *getData(i  ,j  ,k  );
                                        value[1] = *getData(i+1,j  ,k  );
                                        value[2] = *getData(i+1,j+1,k  );
                                        value[0] = *getData(i  ,j  ,k  );
                                        value[1] = *getData(i+1,j  ,k  );
                                        value[2] = *getData(i+1,j+1,k  );
@@ -235,6 +242,7 @@ void IsoSurface::triangulate( void )
                                        eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
 
                                        // grid positions
                                        eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
 
                                        // grid positions
+                                       ntlVec3Gfx pos[8];
                                        pos[0] = ntlVec3Gfx(px    ,py    ,pz);
                                        pos[1] = ntlVec3Gfx(px+gsx,py    ,pz);
                                        pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz);
                                        pos[0] = ntlVec3Gfx(px    ,py    ,pz);
                                        pos[1] = ntlVec3Gfx(px+gsx,py    ,pz);
                                        pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz);
@@ -340,10 +348,7 @@ void IsoSurface::triangulate( void )
                if(mUseFullEdgeArrays) {
                        errMsg("IsoSurface::triangulate","Disabling mUseFullEdgeArrays!");
                }
                if(mUseFullEdgeArrays) {
                        errMsg("IsoSurface::triangulate","Disabling mUseFullEdgeArrays!");
                }
-
-               // subdiv local arrays
-               gfxReal orgval[8];
-               gfxReal subdAr[2][11][11]; // max 10 subdivs!
+               
                ParticleObject* *arppnt = new ParticleObject*[mSizez*mSizey*mSizex];
 
                // construct pointers
                ParticleObject* *arppnt = new ParticleObject*[mSizez*mSizey*mSizex];
 
                // construct pointers
@@ -404,13 +409,25 @@ void IsoSurface::triangulate( void )
 
                debMsgStd("IsoSurface::triangulate",DM_MSG,"Starting. Parts in use:"<<pInUse<<", Subdivs:"<<mSubdivs, 9);
                pz = mStart[2]-(double)(0.*gsz)-0.5*orgGsz;
 
                debMsgStd("IsoSurface::triangulate",DM_MSG,"Starting. Parts in use:"<<pInUse<<", Subdivs:"<<mSubdivs, 9);
                pz = mStart[2]-(double)(0.*gsz)-0.5*orgGsz;
+       
                for(int ok=1;ok<(mSizez-2)*mSubdivs;ok++) {
                        pz += gsz;
                        const int k = ok/mSubdivs;
                        if(k<=0) continue; // skip zero plane
                for(int ok=1;ok<(mSizez-2)*mSubdivs;ok++) {
                        pz += gsz;
                        const int k = ok/mSubdivs;
                        if(k<=0) continue; // skip zero plane
+#if PARALLEL==1        
+#pragma omp parallel for
+#endif 
                        for(int j=1;j<(mSizey-2);j++) {
                                for(int i=1;i<(mSizex-2);i++) {
                        for(int j=1;j<(mSizey-2);j++) {
                                for(int i=1;i<(mSizex-2);i++) {
-
+                                       float value[8];
+                                       ntlVec3Gfx pos[8];
+                                       int cubeIndex;      // index entry of the cube 
+                                       int triIndices[12]; // vertex indices 
+                                       int *eVert[12];
+                                       IsoLevelVertex ilv;
+                                       gfxReal orgval[8];
+                                       gfxReal subdAr[2][11][11]; // max 10 subdivs!
+                                       
                                        orgval[0] = *getData(i  ,j  ,k  );
                                        orgval[1] = *getData(i+1,j  ,k  );
                                        orgval[2] = *getData(i+1,j+1,k  ); // with subdivs
                                        orgval[0] = *getData(i  ,j  ,k  );
                                        orgval[1] = *getData(i+1,j  ,k  );
                                        orgval[2] = *getData(i+1,j+1,k  ); // with subdivs
@@ -422,6 +439,7 @@ void IsoSurface::triangulate( void )
 
                                        // prebuild subsampled array slice
                                        const int sdkOffset = ok-k*mSubdivs; 
 
                                        // prebuild subsampled array slice
                                        const int sdkOffset = ok-k*mSubdivs; 
+
                                        for(int sdk=0; sdk<2; sdk++) 
                                                for(int sdj=0; sdj<mSubdivs+1; sdj++) 
                                                        for(int sdi=0; sdi<mSubdivs+1; sdi++) {
                                        for(int sdk=0; sdk<2; sdk++) 
                                                for(int sdj=0; sdj<mSubdivs+1; sdj++) 
                                                        for(int sdi=0; sdi<mSubdivs+1; sdi++) {
@@ -576,8 +594,13 @@ void IsoSurface::triangulate( void )
 
                                                                                // init isolevel vertex
                                                                                ilv.v = p1 + (p2-p1)*mu; // with subdivs
 
                                                                                // init isolevel vertex
                                                                                ilv.v = p1 + (p2-p1)*mu; // with subdivs
+#if PARALLEL==1        
+#pragma omp critical
+#endif
+                                                                               {
                                                                                mPoints.push_back( ilv );
                                                                                triIndices[e] = (mPoints.size()-1);
                                                                                mPoints.push_back( ilv );
                                                                                triIndices[e] = (mPoints.size()-1);
+                                                                               }
                                                                                // store vertex 
                                                                                *eVert[ e ] = triIndices[e]; 
                                                                        }       else {
                                                                                // store vertex 
                                                                                *eVert[ e ] = triIndices[e]; 
                                                                        }       else {
@@ -587,23 +610,27 @@ void IsoSurface::triangulate( void )
                                                                } // along all edges 
                                                        }
                                                        // removed cutoff treatment...
                                                                } // along all edges 
                                                        }
                                                        // removed cutoff treatment...
-
+                                                       
                                                        // Create the triangles... 
                                                        // Create the triangles... 
+#if PARALLEL==1        
+#pragma omp critical
+#endif
+                                                       {
                                                        for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
                                                                mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] );
                                                                mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] ); // with subdivs
                                                                mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
                                                                //errMsg("TTT"," i1"<<mIndices[mIndices.size()-3]<<" "<< " i2"<<mIndices[mIndices.size()-2]<<" "<< " i3"<<mIndices[mIndices.size()-1]<<" "<< mIndices.size() );
                                                        for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
                                                                mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] );
                                                                mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] ); // with subdivs
                                                                mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
                                                                //errMsg("TTT"," i1"<<mIndices[mIndices.size()-3]<<" "<< " i2"<<mIndices[mIndices.size()-2]<<" "<< " i3"<<mIndices[mIndices.size()-1]<<" "<< mIndices.size() );
-                                                       }
-
                                                        } // triangles in edge table?
                                                        } // triangles in edge table?
+                                                       }
+                                                       }
                                                        
                                                }//si
                                        }// sj
 
                                }//i
                        }// j
                                                        
                                                }//si
                                        }// sj
 
                                }//i
                        }// j
-
+                       
                        // copy edge arrays
                        for(int j=0;j<(mSizey-0)*mSubdivs;j++) 
                                for(int i=0;i<(mSizex-0)*mSubdivs;i++) {
                        // copy edge arrays
                        for(int j=0;j<(mSizey-0)*mSubdivs;j++) 
                                for(int i=0;i<(mSizex-0)*mSubdivs;i++) {