1 /*# ---------------------------------------------------------------------
3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
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
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.
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
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 # ------------------------------------------------------------------------ */
26 /*=========================================================================
29 Module: $RCSfile: UtilVtk3DGeometriSelection.cxx,v $
31 Date: $Date: 2012/11/15 14:14:35 $
32 Version: $Revision: 1.5 $
34 Copyright: (c) 2002, 2003
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.
41 =========================================================================*/
44 #include "UtilVtk3DGeometriSelection.h"
46 #include <vtkPolyData.h>
48 #include <vtkCellLocator.h>
49 #include <vtkPointData.h>
55 //----------------------------------------------------------------------------
56 //----------------------------------------------------------------------------
57 //----------------------------------------------------------------------------
58 UtilVtk3DGeometriSelection::UtilVtk3DGeometriSelection()
65 //----------------------------------------------------------------------------
66 UtilVtk3DGeometriSelection::~UtilVtk3DGeometriSelection()
69 //----------------------------------------------------------------------------
70 void UtilVtk3DGeometriSelection::SetDimentions(int w,int h,int d)
76 //----------------------------------------------------------------------------
77 void UtilVtk3DGeometriSelection::SetMarchingCube(vtkMarchingCubes *mCubes)
81 //----------------------------------------------------------------------------
82 bool UtilVtk3DGeometriSelection::FindCubePointsFromPoints( double* pO, double* pF, double* pickPoint, double* cameraPos )
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;
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 );
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 );
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 );
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 );
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 );
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 );
134 if( points.size( ) >= 2 ) { // Did I find at least 2 points?
136 c( 0 ) = ( double )_width / 2.0;
137 c( 1 ) = ( double )_height / 2.0;
138 c( 2 ) = ( double )_depth / 2.0;
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++ ) {
144 d1 = ( c - *points[ j ] ).GetNorm( );
145 d2 = ( c - *points[ j + 1 ] ).GetNorm( );
149 points[ j ] = points[ j + 1 ];
150 points[ j + 1 ] = swp;
158 // Order the two points according to distance to camera.
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 ];
166 } else return( false );
169 //----------------------------------------------------------------------------
170 bool UtilVtk3DGeometriSelection::GetPointAndNormalIntersection( double* p, double* n, double* pO, double* pF )
172 gtm::TVector< double > n1( 3 ), n2( 3 ), n3( 3 );
178 int subId, returnVal;
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( );
186 locator->SetDataSet( data );
187 locator->Initialize( );
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 );
197 vtkCell* cell = data->GetCell( cellId );
198 vtkPoints* points = cell->GetPoints( );
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( ) );
207 n[ 0 ] = n1( 0 ); n[ 1 ] = n1( 1 ); n[ 2 ] = n1( 2 );
209 points->GetPoint( 0, p1 );
210 points->GetPoint( 1, p2 );
211 points->GetPoint( 2, p3 );
212 this->IntersectPlaneWithLine( p, p1, p2, p3, pO, pF );
214 } else return( false );
218 //----------------------------------------------------------------------------
219 bool UtilVtk3DGeometriSelection::IntersectPlaneWithLine( double* p, double* x1, double* x2, double* x3, double* x4, double* x5 )
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 );
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 );
247 mU( 3, 1 ) = vx4( 0 );
248 mU( 3, 2 ) = vx4( 1 );
249 mU( 3, 3 ) = vx4( 2 );
251 mD( 3, 1 ) = vx6( 0 );
252 mD( 3, 2 ) = vx6( 1 );
253 mD( 3, 3 ) = vx6( 2 );
255 ret = ( mD.Det( ) != 0 );
258 t = mU.Det( ) / mD.Det( );
259 vx6 = ( ( vx4 - vx5 ) * t ) + vx4;
268 //----------------------------------------------------------------------------
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)
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]);
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]);
291 //----------------------------------------------------------------------------
293 // pA point1 of the line
294 // pB point2 of the line
295 double UtilVtk3DGeometriSelection::DistanceMinPointToLine(double *p,double *pA, double *pB)
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);
304 double dx = p[0]-pResult[0];
305 double dy = p[1]-pResult[1];
306 double dz = p[2]-pResult[2];
308 double result=sqrt( dx*dx + dy*dy + dz*dz );