lib/maracasVisuLib/src/interface/wxWindows/widgets/UtilVtk3DGeometriSelection.cxx

   1 /*# ---------------------------------------------------------------------
   2 #
   3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
   4 #                        pour la Sant�)
   5 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
   6 # Previous Authors : Laurent Guigues, Jean-Pierre Roux
   7 # CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
   8 #
   9 #  This software is governed by the CeCILL-B license under French law and
  10 #  abiding by the rules of distribution of free software. You can  use,
  11 #  modify and/ or redistribute the software under the terms of the CeCILL-B
  12 #  license as circulated by CEA, CNRS and INRIA at the following URL
  13 #  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  14 #  or in the file LICENSE.txt.
  15 #
  16 #  As a counterpart to the access to the source code and  rights to copy,
  17 #  modify and redistribute granted by the license, users are provided only
  18 #  with a limited warranty  and the software's author,  the holder of the
  19 #  economic rights,  and the successive licensors  have only  limited
  20 #  liability.
  21 #
  22 #  The fact that you are presently reading this means that you have had
  23 #  knowledge of the CeCILL-B license and that you accept its terms.
  24 # ------------------------------------------------------------------------ */
  25
  26 /*=========================================================================
  27
  28   Program:   wxMaracas
  29   Module:    $RCSfile: UtilVtk3DGeometriSelection.cxx,v $
  30   Language:  C++
  31   Date:      $Date: 2012/11/15 14:14:35 $
  32   Version:   $Revision: 1.5 $
  33
  34   Copyright: (c) 2002, 2003
  35   License:
  36
  37      This software is distributed WITHOUT ANY WARRANTY; without even
  38      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  39      PURPOSE.  See the above copyright notice for more information.
  40
  41 =========================================================================*/
  42
  43
  44 #include "UtilVtk3DGeometriSelection.h"
  45
  46 #include <vtkPolyData.h>
  47 #include <vtkCell.h>
  48 #include <vtkCellLocator.h>
  49 #include <vtkPointData.h>
  50
  51 #include <vector>
  52 #include "matrix.h"
  53
  54
  55 //----------------------------------------------------------------------------
  56 //----------------------------------------------------------------------------
  57 //----------------------------------------------------------------------------
  58 UtilVtk3DGeometriSelection::UtilVtk3DGeometriSelection()
  59 {
  60         _width  =       0;
  61         _height =       0;
  62         _depth  =       0;
  63         _mCubes =       NULL;
  64 }
  65 //----------------------------------------------------------------------------
  66 UtilVtk3DGeometriSelection::~UtilVtk3DGeometriSelection()
  67 {
  68 }
  69 //----------------------------------------------------------------------------
  70 void UtilVtk3DGeometriSelection::SetDimentions(int w,int h,int d)
  71 {
  72         _width  =       w;
  73         _height =       h;
  74         _depth  =       d;
  75 }
  76 //----------------------------------------------------------------------------
  77 void UtilVtk3DGeometriSelection::SetMarchingCube(vtkMarchingCubes *mCubes)
  78 {
  79         _mCubes = mCubes;
  80 }
  81 //----------------------------------------------------------------------------
  82 bool UtilVtk3DGeometriSelection::FindCubePointsFromPoints( double* pO, double* pF, double* pickPoint, double* cameraPos )
  83 {
  84     gtm::TVector< double > p1( 3 ), p2( 3 ), p3( 3 ), p4( 3 ), p5( 3 ), p6( 3 );
  85     gtm::TVector< double > c( 3 );
  86     gtm::TVector< double >* swp;
  87     gtm::TVector< double > tpO( pO, 3, false ), tpF( pF, 3, false );
  88     std::vector< gtm::TVector< double >* > points;
  89     double x1[ 3 ], x2[ 3 ], x3[ 3 ], d1, d2;
  90     int i, j;
  91
  92     // 1st plane intersection
  93     x1[ 0 ] = 0; x1[ 1 ] = _height; x1[ 2 ] =      0;
  94     x2[ 0 ] = 0; x2[ 1 ] =       0; x2[ 2 ] =      0;
  95     x3[ 0 ] = 0; x3[ 1 ] =       0; x3[ 2 ] = _depth;
  96     if( this->IntersectPlaneWithLine( p1.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
  97         points.push_back( &p1 );
  98
  99     // 2nd plane intersection
 100     x1[ 0 ] =      0; x1[ 1 ] = _height; x1[ 2 ] = _depth;
 101     x2[ 0 ] =      0; x2[ 1 ] =       0; x2[ 2 ] = _depth;
 102     x3[ 0 ] = _width; x3[ 1 ] =       0; x3[ 2 ] = _depth;
 103     if( this->IntersectPlaneWithLine( p2.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
 104         points.push_back( &p2 );
 105
 106     // 3rd plane intersection
 107     x1[ 0 ] =      0; x1[ 1 ] = 0; x1[ 2 ] = _depth;
 108     x2[ 0 ] =      0; x2[ 1 ] = 0; x2[ 2 ] =      0;
 109     x3[ 0 ] = _width; x3[ 1 ] = 0; x3[ 2 ] =      0;
 110     if( this->IntersectPlaneWithLine( p3.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
 111         points.push_back( &p3 );
 112
 113     // 4th plane intersection
 114     x1[ 0 ] = _width; x1[ 1 ] = _height; x1[ 2 ] = _depth;
 115     x2[ 0 ] = _width; x2[ 1 ] =       0; x2[ 2 ] = _depth;
 116     x3[ 0 ] = _width; x3[ 1 ] =       0; x3[ 2 ] =      0;
 117     if( this->IntersectPlaneWithLine( p4.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
 118         points.push_back( &p4 );
 119
 120     // 5th plane intersection
 121     x1[ 0 ] = _width; x1[ 1 ] =       0; x1[ 2 ] = 0;
 122     x2[ 0 ] =      0; x2[ 1 ] =       0; x2[ 2 ] = 0;
 123     x3[ 0 ] =      0; x3[ 1 ] = _height; x3[ 2 ] = 0;
 124     if( this->IntersectPlaneWithLine( p5.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
 125         points.push_back( &p5 );
 126
 127     // 6th plane intersection
 128     x1[ 0 ] =      0; x1[ 1 ] = _height; x1[ 2 ] =      0;
 129     x2[ 0 ] =      0; x2[ 1 ] = _height; x2[ 2 ] = _depth;
 130     x3[ 0 ] = _width; x3[ 1 ] = _height; x3[ 2 ] = _depth;
 131     if( this->IntersectPlaneWithLine( p6.GetAnsiRef( ), x1, x2, x3, cameraPos, pickPoint ) )
 132         points.push_back( &p6 );
 133
 134     if( points.size( ) >= 2 ) { // Did I find at least 2 points?
 135
 136         c( 0 ) = ( double )_width  / 2.0;
 137         c( 1 ) = ( double )_height / 2.0;
 138         c( 2 ) = ( double )_depth  / 2.0;
 139
 140         // Sort with bubble sort. Only 30 iterations!
 141         for( i = 0; i <(int) (points.size( )); i++ ) {
 142             for( j = 0; j < (int)(points.size( )) - 1; j++ ) {
 143
 144                 d1 = ( c - *points[ j ] ).GetNorm( );
 145                 d2 = ( c - *points[ j + 1 ] ).GetNorm( );
 146                 if( d2 < d1 ) {
 147
 148                     swp = points[ j ];
 149                     points[ j ] = points[ j + 1 ];
 150                     points[ j + 1 ] = swp;
 151
 152                 } // fi
 153
 154     } // rof
 155
 156         } // rof
 157
 158         // Order the two points according to distance to camera.
 159         c = cameraPos;
 160         d1 = ( c - *points[ 0 ] ).GetNorm( );
 161         d2 = ( c - *points[ 1 ] ).GetNorm( );
 162         tpO = ( d1 < d2 )? *points[ 0 ]: *points[ 1 ];
 163         tpF = ( d1 > d2 )? *points[ 0 ]: *points[ 1 ];
 164         return( true );
 165
 166     } else return( false );
 167
 168 }
 169 //----------------------------------------------------------------------------
 170 bool UtilVtk3DGeometriSelection::GetPointAndNormalIntersection( double* p, double* n, double* pO, double* pF )
 171 {
 172     gtm::TVector< double > n1( 3 ), n2( 3 ), n3( 3 );
 173
 174         //EED 01/03/2010
 175         //int     cellId,
 176         vtkIdType cellId;
 177
 178         int subId,  returnVal;
 179
 180     double t, pcoords[ 3 ], x[ 3 ];
 181     double fpO[ 3 ], fpF[ 3 ];
 182     double p1[ 3 ], p2[ 3 ], p3[ 3 ];
 183     vtkPolyData* data = _mCubes->GetOutput( );
 184     vtkCellLocator* locator = vtkCellLocator::New( );
 185
 186     locator->SetDataSet( data );
 187     locator->Initialize( );
 188     locator->Update( );
 189
 190     fpO[ 0 ] = pO[ 0 ]; fpO[ 1 ] = pO[ 1 ]; fpO[ 2 ] = pO[ 2 ];
 191     fpF[ 0 ] = pF[ 0 ]; fpF[ 1 ] = pF[ 1 ]; fpF[ 2 ] = pF[ 2 ];
 192     returnVal = locator->IntersectWithLine( fpO, fpF, 0.1, t, x, pcoords, subId, cellId );
 193     locator->Delete( );
 194
 195     if( returnVal )
 196     {
 197         vtkCell* cell = data->GetCell( cellId );
 198         vtkPoints* points = cell->GetPoints( );
 199
 200         data->GetPointData( )->GetNormals( )->GetTuple( cell->GetPointIds( )->GetId( 0 ), n1.GetAnsiRef( ) );
 201         data->GetPointData( )->GetNormals( )->GetTuple( cell->GetPointIds( )->GetId( 1 ), n2.GetAnsiRef( ) );
 202         data->GetPointData( )->GetNormals( )->GetTuple( cell->GetPointIds( )->GetId( 2 ), n3.GetAnsiRef( ) );
 203
 204         n1 += n2 + n3;
 205         n1 *= ( 1.0 / 3.0 );
 206         n1.Normalize( );
 207         n[ 0 ] = n1( 0 ); n[ 1 ] = n1( 1 ); n[ 2 ] = n1( 2 );
 208
 209         points->GetPoint( 0, p1 );
 210         points->GetPoint( 1, p2 );
 211         points->GetPoint( 2, p3 );
 212         this->IntersectPlaneWithLine( p, p1, p2, p3, pO, pF );
 213         return( true );
 214     } else return( false );
 215
 216 }
 217
 218 //----------------------------------------------------------------------------
 219 bool UtilVtk3DGeometriSelection::IntersectPlaneWithLine( double* p, double* x1, double* x2, double* x3, double* x4, double* x5 )
 220 {
 221     gtm::TVector< double > vx1( 3 ), vx2( 3 ), vx3( 3 );
 222     gtm::TVector< double > vx4( 3 ), vx5( 3 ), vx6( p, 3, false );
 223     gtm::TMatrix< double > mU( 4, 4 ), mD( 4, 4 );
 224     double t;
 225     bool ret;
 226
 227     vx1 = x1;
 228     vx2 = x2;
 229     vx3 = x3;
 230     vx4 = x4;
 231     vx5 = x5;
 232     vx6 = vx5 - vx4;
 233
 234     mD( 0, 0 ) = mU( 0, 0 ) = 1;
 235     mD( 0, 1 ) = mU( 0, 1 ) = vx1( 0 );
 236     mD( 0, 2 ) = mU( 0, 2 ) = vx1( 1 );
 237     mD( 0, 3 ) = mU( 0, 3 ) = vx1( 2 );
 238     mD( 1, 0 ) = mU( 1, 0 ) = 1;
 239     mD( 1, 1 ) = mU( 1, 1 ) = vx2( 0 );
 240     mD( 1, 2 ) = mU( 1, 2 ) = vx2( 1 );
 241     mD( 1, 3 ) = mU( 1, 3 ) = vx2( 2 );
 242     mD( 2, 0 ) = mU( 2, 0 ) = 1;
 243     mD( 2, 1 ) = mU( 2, 1 ) = vx3( 0 );
 244     mD( 2, 2 ) = mU( 2, 2 ) = vx3( 1 );
 245     mD( 2, 3 ) = mU( 2, 3 ) = vx3( 2 );
 246     mU( 3, 0 ) = 1;
 247     mU( 3, 1 ) = vx4( 0 );
 248     mU( 3, 2 ) = vx4( 1 );
 249     mU( 3, 3 ) = vx4( 2 );
 250     mD( 3, 0 ) = 0;
 251     mD( 3, 1 ) = vx6( 0 );
 252     mD( 3, 2 ) = vx6( 1 );
 253     mD( 3, 3 ) = vx6( 2 );
 254
 255     ret = ( mD.Det( ) != 0 );
 256     if( ret ) {
 257
 258         t = mU.Det( ) / mD.Det( );
 259         vx6 = ( ( vx4 - vx5 ) * t ) + vx4;
 260
 261     } // fi
 262     return( ret );
 263
 264 }
 265
 266
 267
 268 //----------------------------------------------------------------------------
 269 // pResult  Result
 270 // n        Normal to plane
 271 // p        Point of the plane
 272 // pA       point1 of the line
 273 // pB       point2 of the line
 274 void UtilVtk3DGeometriSelection::IntersectionPlaneAndLine(double *pResult,double *n,double *p,double *pA,double *pB)
 275 {
 276         double u,A,B,D;
 277         D = - (n[0]*p[0] + n[1]*p[1] + n[2]*p[2]) ;
 278         A= n[0]*pA[0] +  n[1]*pA[1] + n[2]*pA[2] + D;
 279         B= n[0]*(pA[0]-pB[0]) +  n[1]*(pA[1]-pB[1]) + n[2]*(pA[2]-pB[2]);
 280         if (B!=0)
 281         {
 282                 u = A / B ;
 283         } else {
 284                 u=-1;
 285         }
 286         pResult[0] = pA[0] + u*(pB[0] - pA[0]);
 287         pResult[1] = pA[1] + u*(pB[1] - pA[1]);
 288         pResult[2] = pA[2] + u*(pB[2] - pA[2]);
 289 }
 290
 291 //----------------------------------------------------------------------------
 292 // p   point
 293 // pA  point1 of the line
 294 // pB  point2 of the line
 295 double UtilVtk3DGeometriSelection::DistanceMinPointToLine(double *p,double *pA, double *pB)
 296 {
 297         double pResult[3];
 298         double normalToPlane[3];
 299         normalToPlane[0]=pA[0]-pB[0];
 300         normalToPlane[1]=pA[1]-pB[1];
 301         normalToPlane[2]=pA[2]-pB[2];
 302     IntersectionPlaneAndLine(pResult,normalToPlane,p,pA,pB);
 303
 304         double dx = p[0]-pResult[0];
 305         double dy = p[1]-pResult[1];
 306         double dz = p[2]-pResult[2];
 307
 308     double result=sqrt( dx*dx + dy*dy + dz*dz );
 309
 310         return result;
 311 }
 312
 313
 314
 315