File List
amroc/converters/HDFToV3/V3VisualGridData.h

00001 // -*- C++ -*-
00002 
00003 #ifndef AMROC_V3VISUALGRIDDATA_H
00004 #define AMROC_V3VISUALGRIDDATA_H
00005 
00013 #include "VisualGridData.h" 
00014 
00015 float tol=FLT_EPSILON*100.0;
00016 
00017 template <class VectorType> class V3VisualGrid;
00018 
00025 template <class VectorType>
00026 class V3VisualGridData : public VisualGridData<VectorType> {
00027   typedef VisualGridData<VectorType> base;
00028   typedef V3VisualGridData<VectorType> self;
00029   typedef typename VectorType::InternalDataType DataType;
00030 
00031 public:
00032   typedef V3VisualGrid<VectorType> grid_type;
00033   typedef typename base::grid_type base_grid_type;
00034   typedef typename base::evaluator_type evaluator_type;
00035   typedef typename base::data_block_type data_block_type;
00036   typedef typename base::vector_block_type vector_block_type;
00037 
00038   V3VisualGridData(grid_type& gr, data_block_type& db, BBox &bb, int comp, evaluator_type& ev,
00039                    int lev, int proc, float* pdx, float* pgeom) : 
00040     base((base_grid_type &) gr, db, bb, comp, ev, lev, proc, pdx, pgeom) { 
00041     CrossPlane = (float *)0; 
00042   }
00043 
00044   V3VisualGridData(grid_type& gr, vector_block_type& db, BBox &bb, evaluator_type& ev,
00045                    int lev, int proc, float* pdx, float* pgeom) : 
00046     base((base_grid_type &) gr, db, bb, ev, lev, proc, pdx, pgeom) {
00047     CrossPlane = (float *)0; 
00048   }
00049 
00050   void Draw3DGrids()
00051     {
00052       float orii,orij,orik;
00053       float exti,extj,extk;
00054       exti=base::DB().stepsize(0)*base::DB().extents(0)*base::dx[0];
00055       extj=base::DB().stepsize(1)*base::DB().extents(1)*base::dx[1];
00056       extk=base::DB().stepsize(2)*base::DB().extents(2)*base::dx[2];
00057       
00058       orii=base::DB().lower(0)*base::dx[0]+base::geom[0];
00059       orij=base::DB().lower(1)*base::dx[1]+base::geom[2];
00060       orik=base::DB().lower(2)*base::dx[2]+base::geom[4];     
00061 
00062       float radii = 0.;
00063       int intzero = 0;
00064       int itype = 4;
00065       int icolor = 1;
00066       int np = 5;
00067       float color;
00068       switch (base::_Level){
00069         case 0:
00070           color=.54;
00071           break;
00072         case 1:
00073           color=.44;
00074           break;
00075         case 2:
00076          color=.34;
00077           break;
00078         case 3:
00079           color=.24;
00080           break;
00081         case 4:
00082           color=.14;
00083           break;
00084         default:
00085           color=.04;
00086           break;
00087         }
00088       
00089       /* change first entry of colormap for our drawing and reset afterwards*/
00090       long iv[1]={1};
00091       float rv[3]={color,color,color};
00092       float saveentry[3];
00093       int opt=18;
00094       char blank=' ';
00095       V3_GETSTATE(&opt,iv,saveentry,&blank,1);
00096       V3_SETSTATE(&opt,iv,rv,&blank,1);
00097      
00098       itype=4;
00099       np=5;
00100       float col[8]={0.,0.,0.,0.,0.,0.,0.,0.};
00101       float line1[5*3] = { orii, orij, orik, orii, orij+extj, orik, 
00102                            orii+exti, orij+extj, orik, orii+exti, orij, orik,
00103                            orii, orij, orik};   
00104       
00105       V3_OBJECT3D(&itype, &icolor, line1, &radii, col, &intzero, &np); 
00106       
00107       float line2[5*3] = { orii, orij, orik+extk,
00108                            orii, orij+extj, orik+extk, 
00109                            orii+exti, orij+extj, orik+extk,
00110                            orii+exti, orij, orik+extk, 
00111                            orii, orij, orik+extk};    
00112        
00113       V3_OBJECT3D(&itype, &icolor, line2, &radii, col, &intzero, &np); 
00114       
00115       itype = 3; 
00116       np = 2*4; 
00117       float lines[2*4*3] = { orii, orij, orik,
00118                              orii, orij, orik+extk,  
00119                              orii, orij+extj, orik,
00120                              orii, orij+extj, orik+extk, 
00121                              
00122                              orii+exti, orij+extj, orik,
00123                              orii+exti, orij+extj, orik+extk,  
00124                              orii+exti, orij, orik,
00125                              orii+exti, orij, orik+extk};   
00126       
00127       V3_OBJECT3D(&itype, &icolor, lines, &radii, col, &intzero, &np);
00128       V3_SETSTATE(&opt,iv,saveentry,&blank,1);
00129     }
00130  
00131     void Draw3DCrossSections(){
00132       if (CrossPlane == (float *)0)
00133         return;
00134       
00135       float ext[3]={base::DB().stepsize(0)*base::DB().extents(0)*base::dx[0],
00136                     base::DB().stepsize(1)*base::DB().extents(1)*base::dx[1],
00137                     base::DB().stepsize(2)*base::DB().extents(2)*base::dx[2]};
00138       float ori[3]={base::DB().lower(0)*base::dx[0]+base::geom[0],
00139                     base::DB().lower(1)*base::dx[1]+base::geom[2],
00140                     base::DB().lower(2)*base::dx[2]+base::geom[4]};
00141       float vector1[3]={CrossPlane[0]-CrossPlane[3],
00142                         CrossPlane[1]-CrossPlane[4],
00143                         CrossPlane[2]-CrossPlane[5]};
00144       float vector2[3]={CrossPlane[0]-CrossPlane[6],
00145                         CrossPlane[1]-CrossPlane[7],
00146                         CrossPlane[2]-CrossPlane[8]};
00147       float ortho[3]={vector1[1]*vector2[2]-vector1[2]*vector2[1],
00148                       vector1[2]*vector2[0]-vector1[0]*vector2[2],
00149                       vector1[0]*vector2[1]-vector1[1]*vector2[0]};
00150       float linevector[12*3]={0,0,ext[2],
00151                               0,0,ext[2],
00152                               0,0,ext[2],
00153                               0,0,ext[2],
00154                               0,ext[1],0,
00155                               0,ext[1],0,
00156                               0,ext[1],0,
00157                               0,ext[1],0,
00158                               ext[0],0,0,
00159                               ext[0],0,0,
00160                               ext[0],0,0,
00161                               ext[0],0,0};
00162       float lineori[12*3]={ori[0],ori[1],ori[2],
00163                            ori[0]+ext[0],ori[1],ori[2],
00164                            ori[0]+ext[0],ori[1]+ext[1],ori[2],
00165                            ori[0],ori[1]+ext[1],ori[2],
00166 
00167                            ori[0],ori[1],ori[2],
00168                            ori[0]+ext[0],ori[1],ori[2],
00169                            ori[0]+ext[0],ori[1],ori[2]+ext[2],
00170                            ori[0],ori[1],ori[2]+ext[2],
00171 
00172                            ori[0],ori[1],ori[2],
00173                            ori[0],ori[1],ori[2]+ext[2],
00174                            ori[0],ori[1]+ext[1],ori[2]+ext[2],
00175                            ori[0],ori[1]+ext[1],ori[2]};
00176       float polygon[4*3];
00177       int index=0;
00178 
00179       for(int i=0;i<12;i++){
00180 
00181         float A[9]={vector1[0],vector2[0],-linevector[3*i],
00182                     vector1[1],vector2[1],-linevector[3*i+1],
00183                     vector1[2],vector2[2],-linevector[3*i+2]};
00184         float b[3]={lineori[3*i]-CrossPlane[0],
00185                     lineori[3*i+1]-CrossPlane[1],
00186                     lineori[3*i+2]-CrossPlane[2]};
00187         float parallel=linevector[3*i]*ortho[0]
00188                       +linevector[3*i+1]*ortho[1]
00189                       +linevector[3*i+2]*ortho[2];
00190     
00191         if(std::fabs(parallel)>tol){
00192           gauss(A,b);   
00193           int k; 
00194           for(int j=0;j<3;j++)
00195             polygon[index*3+j]=b[2]*linevector[3*i+j]+lineori[3*i+j];
00196 
00197           if(tol<b[2] && b[2]<1.0-tol)
00198             index++;
00199           if((std::fabs(b[2])<tol || std::fabs(1.0-b[2])<tol)){
00200             if(index>0){
00201               /*check if we've alread got one of the corners*/
00202               for(k=0;k<index;k++){
00203                 if(std::fabs(polygon[index*3]-polygon[k*3])<tol &&
00204                    std::fabs(polygon[index*3+1]-polygon[k*3+1])<tol &&
00205                    std::fabs(polygon[index*3+2]-polygon[k*3+2])<tol)
00206                   break;
00207               }
00208               if(k==index)
00209                 index++;
00210             }
00211             else
00212               index++;
00213 
00214           }
00215           
00216           
00217         }
00218       }
00219       
00220       if(index>2){
00221         float radii = 0.;
00222         int intzero = 0;
00223         int itype = 1;
00224         int icolor = 1;
00225         int np=3*3;
00226         if(index==3)
00227           np=3;
00228         float color;
00229         switch (base::_Level){
00230         case 0:
00231           color=.6;
00232           break;
00233         case 1:
00234           color=.5;
00235           break;
00236         case 2:
00237           color=.4;
00238           break;
00239         case 3:
00240           color=.3;
00241           break;
00242         case 4:
00243           color=.2;
00244           break;
00245         default:
00246           color=.1;
00247           break;
00248         }
00249         /* change first entry of colormap for our drawing and reset afterwards*/
00250         long iv[1]={1};
00251         float rv[3]={color,color,color};
00252         float saveentry[3];
00253         int opt=18;
00254         char blank=' ';
00255         V3_GETSTATE(&opt,iv,saveentry,&blank,1);
00256         V3_SETSTATE(&opt,iv,rv,&blank,1);
00257         color=0;
00258         float triang[3*3*3]={polygon[0],polygon[1],polygon[2],
00259                              polygon[3],polygon[4],polygon[5],
00260                              polygon[6],polygon[7],polygon[8],
00261                              
00262                              polygon[3],polygon[4],polygon[5],
00263                              polygon[6],polygon[7],polygon[8],
00264                              polygon[9],polygon[10],polygon[11],
00265                              
00266                              polygon[6],polygon[7],polygon[8],
00267                              polygon[9],polygon[10],polygon[11],
00268                              polygon[0],polygon[1],polygon[2]};
00269         float triangcol[3*3]={color,color,color,
00270                               color,color,color,
00271                               color,color,color}; 
00272         V3_OBJECT3D(&itype, &icolor, triang, &radii, triangcol, &intzero, &np);
00273         V3_SETSTATE(&opt,iv,saveentry,&blank,1);
00274       }
00275     }
00276 
00277 
00278     int FindMatchingNode(int ijk[]){
00279     
00280     if (ijk[0]<base::DB().lower(0) || ijk[0]>base::DB().upper(0)+base::DB().stepsize(0) ||
00281         ijk[1]<base::DB().lower(1) || ijk[1]>base::DB().upper(1)+base::DB().stepsize(1) ||
00282         ijk[2]<base::DB().lower(2) || ijk[2]>base::DB().upper(2)+base::DB().stepsize(2))
00283         return -1;
00284     
00285     if (ijk[0]%base::DB().stepsize(0) != 0 ||
00286         ijk[1]%base::DB().stepsize(1) != 0 ||
00287         ijk[2]%base::DB().stepsize(2) != 0)
00288       return -1;
00289     
00290     int exti=base::DB().extents(0)+1;
00291     int extj=base::DB().extents(1)+1;
00292     return ( (ijk[2]-base::DB().lower(2))/base::DB().stepsize(2)*exti*extj +
00293              (ijk[1]-base::DB().lower(1))/base::DB().stepsize(1)*exti +
00294              (ijk[0]-base::DB().lower(0))/base::DB().stepsize(0) + 1 + base::_nodeOffset );
00295   }
00296   
00297   inline const int& NodeOffset() const { return base::_nodeOffset; }
00298   inline int NodeOffset() { return base::_nodeOffset; }  
00299   
00300   int ExternalSurface(Coords& slc, Coords& suc, BBoxList& surflist) {
00301     for (BBox *bb = surflist.first(); bb ; bb = surflist.next()) 
00302       if (bb->inside(slc) && bb->inside(suc))
00303         return 0;
00304     return 1;
00305   }
00306   
00307   void SetSurfaces(int dir, BBoxList& surflist, int& ksurf, int scel[][4],
00308                    int set) {
00309     
00310     int exti=base::DB().extents(0)+1;
00311     int extj=base::DB().extents(1)+1;
00312     int extk=base::DB().extents(2)+1;
00313     int i, j, k, kn;
00314     const Coords& base = base::DB().lower();
00315     const Coords& step = base::DB().stepsize();
00316 
00317     switch (dir) {
00318     case 0:
00319       for (k = 0; k < extk-1; k++){
00320         for (j = 0; j < extj-1; j++){
00321           kn = k*exti*extj + j*exti + 1+base::_nodeOffset;
00322 
00323           Coords slc(3,base(0),base(1)+ j   *step(1),base(2)+ k   *step(2));
00324           Coords suc(3,base(0),base(1)+(j+1)*step(1),base(2)+(k+1)*step(2));
00325           if (ExternalSurface(slc,suc,surflist)) {
00326             if (set) {
00327               scel[ksurf][0] = kn;
00328               scel[ksurf][1] = kn + exti;
00329               scel[ksurf][2] = kn + exti + exti*extj;
00330               scel[ksurf][3] = kn + exti*extj;
00331             }
00332             ksurf++;
00333           }       
00334         }
00335       }
00336       break;
00337     case 1:
00338       for (k = 0; k < extk-1; k++){
00339         for (j = 0; j < extj-1; j++){
00340           kn = k*exti*extj + j*exti + exti+base::_nodeOffset;
00341 
00342           Coords slc(3,base(0)+(exti-1)*step(0),base(1)+ j   *step(1),
00343                      base(2)+ k   *step(2));
00344           Coords suc(3,base(0)+(exti-1)*step(0),base(1)+(j+1)*step(1),
00345                      base(2)+(k+1)*step(2));
00346           if (ExternalSurface(slc,suc,surflist)) {
00347             if (set) {
00348               scel[ksurf][0] = kn;
00349               scel[ksurf][1] = kn + exti;
00350               scel[ksurf][2] = kn + exti + exti*extj;
00351               scel[ksurf][3] = kn + exti*extj;  
00352             }
00353             ksurf++;
00354           }  
00355         }
00356       }
00357       break;
00358     case 2:
00359       for (k = 0; k < extk-1; k++){
00360         for (i = 0; i < exti-1; i++){     
00361           kn = k*exti*extj + i + 1+base::_nodeOffset;
00362 
00363           Coords slc(3,base(0)+ i   *step(0),base(1),base(2)+ k   *step(2));
00364           Coords suc(3,base(0)+(i+1)*step(0),base(1),base(2)+(k+1)*step(2));
00365           if (ExternalSurface(slc,suc,surflist)) {
00366             if (set) {
00367               scel[ksurf][0] = kn;
00368               scel[ksurf][1] = kn + 1;
00369               scel[ksurf][2] = kn + 1 + exti*extj;
00370               scel[ksurf][3] = kn + exti*extj;
00371             }
00372             ksurf++;
00373           }
00374         }
00375       }     
00376       break;
00377     case 3:
00378       for (k = 0; k < extk-1; k++){
00379         for (i = 0; i < exti-1; i++){     
00380           kn = k*exti*extj + (extj-1)*exti + i + 1+base::_nodeOffset;
00381 
00382           Coords slc(3,base(0)+ i   *step(0),base(1)+(extj-1)*step(1),
00383                      base(2)+ k   *step(2));
00384           Coords suc(3,base(0)+(i+1)*step(0),base(1)+(extj-1)*step(1),
00385                      base(2)+(k+1)*step(2));
00386           if (ExternalSurface(slc,suc,surflist)) {
00387             if (set) {
00388               scel[ksurf][0] = kn;
00389               scel[ksurf][1] = kn + 1;
00390               scel[ksurf][2] = kn + 1 + exti*extj;
00391               scel[ksurf][3] = kn + exti*extj;
00392             }
00393             ksurf++;
00394           }
00395         }
00396       }
00397       break;
00398     case 4:
00399       for (j = 0; j < extj-1; j++){
00400         for (i = 0; i < exti-1; i++){     
00401           kn = j*exti + i + 1+base::_nodeOffset;
00402 
00403           Coords slc(3,base(0)+ i   *step(0),base(1)+ j   *step(1),base(2));
00404           Coords suc(3,base(0)+(i+1)*step(0),base(1)+(j+1)*step(1),base(2));
00405           if (ExternalSurface(slc,suc,surflist)) {
00406             if (set) {
00407               scel[ksurf][0] = kn;
00408               scel[ksurf][1] = kn + 1;
00409               scel[ksurf][2] = kn + 1 + exti;
00410               scel[ksurf][3] = kn + exti;
00411             }
00412             ksurf++;
00413           }
00414         }
00415       }
00416       break;
00417     case 5:
00418       for (j = 0; j < extj-1; j++){
00419         for (i = 0; i < exti-1; i++){     
00420           kn = (extk-1)*extj*exti + j*exti + i + 1+base::_nodeOffset;
00421 
00422           Coords slc(3,base(0)+ i   *step(0),base(1)+ j   *step(1),
00423                      base(2)+(extk-1)*step(2));
00424           Coords suc(3,base(0)+(i+1)*step(0),base(1)+(j+1)*step(1),
00425                      base(2)+(extk-1)*step(2));
00426           if (ExternalSurface(slc,suc,surflist)) {
00427             if (set) {
00428               scel[ksurf][0] = kn;
00429               scel[ksurf][1] = kn + 1;
00430               scel[ksurf][2] = kn + 1 + exti;
00431               scel[ksurf][3] = kn + exti;
00432             }
00433             ksurf++;
00434           }
00435         }
00436       }
00437       break;
00438 
00439     default:
00440       break;
00441     }
00442   }
00443 
00444   void SetCrossPlane(float* cp) { CrossPlane = cp; }
00445 
00446 private:  
00447   void gauss(float *a,float *b) {
00448     int i,j,k;  
00449     float dummy;
00450         
00451     for(k=0;k<2;k++){
00452       for(i=k;i<2;i++){ 
00453         if(std::fabs(a[i*3])>tol)
00454           break;
00455       }
00456       for(j=0;j<3;j++){
00457         dummy=a[i*3+j];
00458         a[i*3+j]=a[k*3+j];
00459         a[k*3+j]=dummy;
00460       } 
00461       dummy=b[i];
00462       b[i]=b[k];
00463       b[k]=dummy;
00464       
00465     }
00466     
00467     if(std::fabs(a[3])>tol){
00468       a[4]=a[3]*a[1]-a[0]*a[4];
00469       a[5]=a[3]*a[2]-a[0]*a[5];
00470       b[1]=a[3]*b[0]-a[0]*b[1];
00471     }
00472     if(std::fabs(a[6])>tol){
00473       a[7]=a[6]*a[1]-a[0]*a[7];
00474       a[8]=a[6]*a[2]-a[0]*a[8];
00475       b[2]=a[6]*b[0]-a[0]*b[2];
00476     }
00477     
00478     b[2]= (a[4]*b[2]-a[7]*b[1])/(a[4]*a[8]-a[7]*a[5]);
00479   } 
00480 
00481 private:
00482   float* CrossPlane;
00483 };
00484 
00485 #endif
Generated on Thu Jun 30 02:19:02 2016 for AMROC's HDF Tools by 1.6.3