/*========================================================================= Program: wxMaracas Module: $RCSfile: marAxis.cpp,v $ Language: C++ Date: $Date: 2008/10/31 16:32:54 $ Version: $Revision: 1.1 $ Copyright: (c) 2002, 2003 License: This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notice for more information. =========================================================================*/ #include #include #include #include #include #include #include #include #include "marAxis.h" #include "marVector.h" #include "ExtractionAxe.h" #include "FonctionsGrales.h" #include #include #include #include #include #include #include #include #include #include #include #include #include // ---------------------------------------------------------------------------- marAxis::marAxis( marParameters* p ) : marObject( p ), kCurve( 3, INDX_count - 3 ), _healthySlice( -1 ), _startQuant( -1 ), _finishQuant( -1 ), _actualQuant( -1 ) { _totalAxisLenght=-1; _subAxisLenght=-1; _referenceArea=-1; _referenceAverDiam=-1; _allData = NULL; _healthySliceStart=-1; _healthySliceEnd=-1; reset( ); calibration = false; } // ---------------------------------------------------------------------------- void marAxis::addAxisPoint( double* p ) { double pt[ INDX_count ]; memcpy( pt, p, sizeof( POINTAXE ) ); pt[ INDX_SIGNALVALUE ] = 0; addControlPoint( pt, pt + 3 ); } // ---------------------------------------------------------------------------- double* marAxis::getNormal( unsigned int slice ) { double* ret = new double[ 3 ]; if ( slice == 0 ) { ret[0] = _points[ slice + 1 ][0] - _points[ slice ][0]; ret[1] = _points[ slice + 1 ][1] - _points[ slice ][1]; ret[2] = _points[ slice + 1 ][2] - _points[ slice ][2]; } else if ( slice == _points.size( ) - 1 ) { ret[0] = _points[ slice - 1 ][0] - _points[ slice ][0]; ret[1] = _points[ slice - 1 ][1] - _points[ slice ][1]; ret[2] = _points[ slice - 1 ][2] - _points[ slice ][2]; } else if ( slice > 0 && slice < _points.size( ) - 1 ) { ret[0] = _points[ slice - 1 ][0] - _points[ slice + 1 ][0]; ret[1] = _points[ slice - 1 ][1] - _points[ slice + 1 ][1]; ret[2] = _points[ slice - 1 ][2] - _points[ slice + 1 ][2]; } return( ret ); } // ---------------------------------------------------------------------------- void marAxis::changeAxisResolution( ) {/* int i; int nOut = ( int )ceil( getParameters( )-> getDoubleParam( marParameters::e_axis_discret_step ) * ( double )( GetNumberOfPoints( ) ) ); i = GetNumberOfPoints( ); double d = getParameters( )-> getDoubleParam( marParameters::e_axis_discret_step ); SetResolution( nOut ); for( i = 0; i < _contours.size( ); i++ ) delete _contours[ i ]; _contours.clear( ); for( i = 0; i < nOut; i++ ) _contours.push_back( new marContour( getParameters( ) ) ); */ } // ---------------------------------------------------------------------------- void marAxis::calculateSignal( kVolume* vol ) { /** * PSIGNAL_USHORT <-> unsigned short * */ /* PSIGNAL_FLOAT PutSignalValueFloat(PSIGNAL_FLOAT signal, int pos, float value) { signal[pos] = value; return signal; }*/ ushort maxlocal; double *p; double rayon; unsigned int i, mask_size = getParameters( )-> getIntParam( marParameters::e_mask_size ); while( _signal.size( ) > 0 ) _signal.pop_back( ); // Intensity signal PSIGNAL_USHORT signal = ( PSIGNAL_USHORT ) IdSigAlloc( _points.size( ), SIG_USHORT ); for( i = 0; i < _points.size( ); i++ ) { p = _points[i]; rayon = RayonLocal( ( int ) ( p[ 0 ] ), ( int ) ( p[ 1 ] ), ( int ) ( p[ 2 ] ), _points_disc, getNumberOfControlPoints( )); maxlocal = vol->GetMaxIntSphere2( p, rayon ); PutSignalValue( signal, i, maxlocal ); } // rof PSIGNAL_FLOAT mask_signal = ( PSIGNAL_FLOAT )IdSigAlloc( mask_size, SIG_FLOAT ); for( i = 0; i < mask_size; i++ ) PutSignalValueFloat( mask_signal, i, 1 ); PSIGNAL_USHORT signalfilt = ( PSIGNAL_USHORT )IdSigAlloc( _points.size( ), SIG_USHORT ); signalfilt = MonSigConvolve( signal, mask_signal, signalfilt, 1.0 / ( double )mask_size, 0 ); for( i = 0; i < _points.size( ); i++ ){ _signal.push_back( GetSignalValue( signalfilt, i ) ); } IdSigFree( signal ); IdSigFree( mask_signal ); IdSigFree( signalfilt ); } // ---------------------------------------------------------------------------- void marAxis::start( ) { /* TODO int swap, n = GetNumberOfPoints( ); swap = GSL_MIN( _startQuant, _finishQuant ); _finishQuant = GSL_MAX( _startQuant, _finishQuant ); _startQuant = swap; _startQuant = ( _startQuant >= 0 )? _startQuant: 0; _finishQuant = ( _finishQuant >= n )? _finishQuant: n - 1; _actualQuant = _startQuant; */ } // ---------------------------------------------------------------------------- void marAxis::cut( int slice, bool up ) { // TODO } // ---------------------------------------------------------------------------- void marAxis::doSpline( ) { unsigned int i; unsigned int nop, nip; float length=0, t; float p1[ 3 ], p2[ 3 ]; double p[marAxis::INDX_count]; vtkKochanekSpline* spX = vtkKochanekSpline::New( ); vtkKochanekSpline* spY = vtkKochanekSpline::New( ); vtkKochanekSpline* spZ = vtkKochanekSpline::New( ); for( i = 0, length = 0.0; i < _controlPoints.size( ); i++ ) { getControlPoint(i, p, p+3); p2[ 0 ] = (float) p[marAxis::INDX_X]; p2[ 1 ] = (float) p[marAxis::INDX_Y]; p2[ 2 ] = (float) p[marAxis::INDX_Z]; spX->AddPoint( (float) i, p2[ 0 ] ); spY->AddPoint( (float) i, p2[ 1 ] ); spZ->AddPoint( (float) i, p2[ 2 ] ); if( i > 0 ) length += (float) (sqrt( (p2[0]-p1[0])*(p2[0]-p1[0]) + (p2[1]-p1[1])*(p2[1]-p1[1]) + (p2[2]-p1[2])*(p2[2]-p1[2]) ) ); p1[ 0 ] = p2[ 0 ]; p1[ 1 ] = p2[ 1 ]; p1[ 2 ] = p2[ 2 ]; } // rof nop = ( unsigned int ) ( getParameters( )->getDoubleParam( marParameters::e_axis_discret_step ) * length ); nip = _controlPoints.size( ); _healthySlice = -1; while( _points.size( ) > 0 ) { delete _points[ _points.size( ) - 1 ]; _points.pop_back( ); } // fwhile while( _contours.size( ) > 0 ) { delete _contours[ _contours.size( ) - 1 ]; _contours.pop_back( ); } // fwhile for( i = 0; i < nop; i++ ) { t = (float) (( ( double ) nip - 1.0 ) / ( ( double ) nop - 1.0 ) * i); double* np = new double[ 3 ]; np[ 0 ] = 0; np[ 1 ] = 0; np[ 2 ] = 0; np[ 0 ] = (double) (spX->Evaluate( t )); np[ 1 ] = (double) (spY->Evaluate( t )); np[ 2 ] = (double) (spZ->Evaluate( t )); _points.push_back( np ); _contours.push_back( new marContour( i, getParameters( ) ) ); } // rof spX->Delete( ); spY->Delete( ); spZ->Delete( ); } // ---------------------------------------------------------------------------- void marAxis::sliceVolumeAxis( kVolume* vol, bool forceCnt ) { /* int nCnts = ( int ) getNumberOfSplinePoints( ); int sizeIma = getParameters( )->getSizeIma( ); double dimIma = getParameters( )->getDimIma( ); double voxSize = getParameters( )->getVoxelSize( ); if( forceCnt ) for( unsigned int i = 0; i < _contours.size( ); i++ ) delete _contours[ i ]; _contours.clear( ); for( unsigned int i = 0; i < _slices.size( ); i++ ) { delete _slices[ i ]; _3Dslices[ i ]->Delete( ); _quantificationImages[ i ]->Delete( ); } _slices.clear( ); _3Dslices.clear( ); _quantificationImages.clear( ); // "Cutter" initialization vtkStructuredPoints* stPoints = vtkStructuredPoints::New( ); // Perform "cut" double *p, *n; for( int k = 0; k < nCnts; k++ ) { //s = ( double )k / ( double )( nCnts - 1 ); p = _points[k]; n = getNormal( k ); vtkPlaneSource* pSource = vtkPlaneSource::New( ); pSource->SetOrigin( p[ 0 ], p[ 1 ], p[ 2 ] ); pSource->SetPoint1( p[ 0 ] + dimIma - 1.0, p[ 1 ], p[ 2 ] ); pSource->SetPoint2( p[ 0 ], p[ 1 ], p[ 2 ] + dimIma - 1.0 ); pSource->SetResolution( sizeIma - 1, sizeIma - 1 ); pSource->Update( ); pSource->SetCenter( p[ 0 ], p[ 1 ], p[ 2 ] ); pSource->SetNormal( n[ 0 ], n[ 1 ], n[ 2 ] ); pSource->Update( ); _3Dslices.push_back( vtkProbeFilter::New( ) ); _3Dslices[ k ]->SetInput( ( vtkDataSet* )pSource->GetOutput( ) ); _3Dslices[ k ]->SetSource( vol->castVtk( ) ); _3Dslices[ k ]->Update( ); pSource->Delete( ); stPoints->GetPointData( )-> SetScalars( _3Dslices[ k ]->GetOutput( )-> GetPointData( )->GetScalars( ) ); stPoints->SetDimensions( sizeIma, sizeIma, 1 ); stPoints->SetScalarType( vol->castVtk( )->GetDataObjectType( ) ); stPoints->Update( ); if( forceCnt ) { marContour* tmp_marContour = new marContour( k, getParameters( ) ) ; tmp_marContour->calculateVariables(); _contours.push_back( tmp_marContour ); } _slices.push_back( new kVolume( kVolume::UCHAR, 1, 1, 1 ) ); *( _slices[ k ] ) = ( vtkImageData* )stPoints; //copy _quantificationImages.push_back( (_slices[ k ])->castVtk( ) ); } // rof stPoints->Delete( ); */ } // ---------------------------------------------------------------------------- int marAxis::getNumberOfContours( ) { return( _contours.size( ) ); } // ---------------------------------------------------------------------------- int marAxis::getNumberOfSplinePoints( ) { return( _points.size( ) ); } /* // EED Esto toca limpiarlo o reimplementarlo en varias partes o algo // ---------------------------------------------------------------------------- vtkPolyData* marAxis::setContour( int i, int x, int y, std::vector< double* >* points, marContour* c ) { if( c != NULL ) { c->setS( c->getS( ) ); _contours[ i ]->copyFrom( *c ); } else { // PS -> // int numberOfPoints, numberOfCells; // PS // PS -> // int id; // PS // PS -> // gslobj_vector p( 3 ); // PS // PS -> // double new_p[3]; // PS // PS -> // float pf[ 3 ]; // PS vtkImageData* imagedata; // Contour parameters // PS -> int typ = getParameters( )->getIntParam( marParameters::e_algorithm_type ); double thr = getParameters( )->getDoubleParam( marParameters::e_threshold_isocontour ); // PS -> double thr = ( typ == marParameters::ISOCONTOURS )? // PS -> getParameters( )-> // PS -> getDoubleParam( marParameters::e_threshold_isocontour ): // PS -> ( typ == marParameters::SNAKE ) ? // PS -> getParameters( )-> // PS -> getDoubleParam( marParameters::e_threshold_snake_isocontour ): // PS -> 30; int sizeIma = getParameters( )->getSizeIma( ); double localint = getSignal( i ); double sigma = getParameters( )->getDoubleParam( marParameters::e_sigma ); double threshold = localint * thr / 100.0; int dimIma = ( int ) ( ( double ) sizeIma /getParameters( )->getDoubleParam( marParameters::e_scale ) ); int vmin = ( int )threshold; int vmax = ( int )threshold; x = ( x != -1 )? x: sizeIma / 2; y = ( y != -1 )? y: sizeIma / 2; // Isocontour imagedata = (vtkImageData*) ((getSlice( i , NULL ))->castVtk( )); float *range = imagedata->GetScalarRange(); //vtkKitwareContourFilter* cntVTK = vtkKitwareContourFilter::New( ); vtkContourFilter* cntVTK = vtkContourFilter::New( ); cntVTK->SetInput( imagedata ); cntVTK->SetNumberOfContours( 1 ); //cntVTK->SetValue( 0, vmin ); cntVTK->SetValue( 0, (range[1]*thr/100) ); //cntVTK->SetValue( 1, vmax ); cntVTK->Update( ); vtkCleanPolyData* cpd = vtkCleanPolyData::New( ); cpd->SetInput( cntVTK->GetOutput( ) ); cpd->ConvertLinesToPointsOff( ); cpd->Update( ); vtkPolyDataConnectivityFilter* conn = vtkPolyDataConnectivityFilter::New( ); //conn->SetMaxRecursionDepth( 3000 ); // PS -> //conn->SetInput( cntVTK->GetOutput( ) ); PS conn->SetInput( cpd->GetOutput( ) ); conn->SetClosestPoint( x, y, 0 ); conn->SetExtractionModeToClosestPointRegion( ); conn->Update( ); vtkCleanPolyData* cpd2 = vtkCleanPolyData::New( ); cpd2->SetInput( conn->GetOutput( ) ); cpd2->Update(); vtkPolyData* polyDataResult = cpd2->GetOutput() ; double ta = _contours[ i ]->getS( ); _contours[ i ]->reset( ); // PS -> gslobj_vector o( 3 ), n( 3 ), y( 3 ), vr( 3 ); int id; int numberOfCells; int numberOfPoints,numberOfPoints2; vtkCell* cell; vtkIdList* pids; marVector p( 3 ); // PS -> double new_p[3]; float pf[ 3 ]; numberOfCells = polyDataResult->GetNumberOfCells( ); cell = polyDataResult->GetCell( 0 ); pids = cell->GetPointIds( ); numberOfPoints = pids->GetNumberOfIds( ); numberOfPoints2=polyDataResult->GetNumberOfPoints(); marVector o(3),n(3),y(3),vr(3); double cost, sint; y( 0 ) = y( 2 ) = 0.0; y( 1 ) = 1.0; getPoint( ( double* )o, ta ); getTangent( ( double* )n, ta ); vr = y.cross( n ).normalize( ); sint = y.cross( n ).norm2( ); cost = y.dot( n ); for( int j = 0; j < numberOfPoints; j++ ) { id = pids->GetId( j ); polyDataResult->GetPoint( id, pf ); p( 0 ) = pf[ 0 ]; p( 1 ) = 0.0; p( 2 ) = pf[ 1 ]; // 3D transformation p = ( p * cost ) + ( vr * ( vr.dot( p ) * ( 1.0 - cost ) ) ) - ( p.cross( vr ) * sint ) + o; // PS -> double p_cross_vr[3]; // PS -> vtkMath::Cross( p, vr, p_cross_vr); // PS -> for(int k=0; k<3; k++) // PS -> { // PS -> new_p[i] = (new_p[k] * cost) + ( vr(k) * ( vtkMath::Dot(vr, p) * ( 1.0 - cost ) ) ) - // PS -> ( p_cross_vr[k] * sint ) + o(i); // PS -> } //TODO: check if p and new_p are the same... _contours[ i ]->addContourPoint( p ); } // rof // PS -> int nbPoints=polyDataResult->GetNumberOfPoints(); // PS -> // PS -> float* vtkContourPoint; // PS -> vtkContourPoint= new float[3]; // PS -> double* marContourPoint; // PS -> marContourPoint= new double[3]; // PS -> // PS -> for (int j=0;j< nbPoints;j++) // PS -> { // PS -> vtkContourPoint=polyDataResult->GetPoint(j); // PS -> marContourPoint[0]=(double)vtkContourPoint[0]; // PS -> marContourPoint[1]=(double)vtkContourPoint[1]; // PS -> marContourPoint[2]=(double)vtkContourPoint[2]; // PS -> _contours[ i ]->addContourPoint( marContourPoint ); // PS -> } marContour * tutu=_contours[ i ]; _contours[ i ]->getArea(); _contours[ i ]->getPerimeter(); return( polyDataResult ); //return( cpd->GetOutput( ) ); //return( cntVTK->GetOutput( ) ); //EED1 vtkStripper* isoStrips = vtkStripper::New( ); //EED1 isoStrips->SetInput( cpd2->GetOutput( ) ); //EED1 isoStrips->Update( ); //EED1 //EED1 vtkPolyData* data = isoStrips->GetOutput( ); //EED1 data->Update( ); //EED1 vtkCell* cell; //EED1 vtkIdList* pids; //EED1 //EED1 numberOfCells = data->GetNumberOfCells( ); //EED1 cell = data->GetCell( 0 ); //EED1 pids = cell->GetPointIds( ); //EED1 numberOfPoints = pids->GetNumberOfIds( ); //EED2 _contours[ i ]->reset( ); //EED2 if( typ == marParameters::ISOCONTOURS ) { //EED2 gslobj_vector o( 3 ), n( 3 ), y( 3 ), vr( 3 ); //EED2 double cost, sint; //EED2 y( 0 ) = y( 2 ) = 0.0; //EED2 y( 1 ) = 1.0; //EED2 getPoint( ( double* )o, ta ); //EED2 getTangent( ( double* )n, ta ); //EED2 vr = y.cross( n ).normalize( ); //EED2 sint = y.cross( n ).norm2( ); //EED2 cost = y.dot( n ); //EED2 //EED2 for( int j = 0; j < numberOfPoints; j++ ) { //EED2 id = pids->GetId( j ); //EED2 data->GetPoint( id, pf ); //EED2 p( 0 ) = pf[ 0 ]; //EED2 p( 1 ) = 0.0; //EED2 p( 2 ) = pf[ 1 ]; //EED2 //EED2 // 3D transformation //EED2 p = ( p * cost ) + ( vr * ( vr.dot( p ) * ( 1.0 - cost ) ) ) - //EED2 ( p.cross( vr ) * sint ) + o; //EED2 //EED2 double p_cross_vr[3]; //EED2 vtkMath::Cross( p, vr, p_cross_vr); //EED2 for(int ii=0; ii<3; ii++) { //EED2 new_p[i] = (new_p[ii] * cost) + ( vr(ii) * ( vtkMath::Dot(vr, p) * ( 1.0 - cost ) ) ) - //EED2 ( p_cross_vr[ii] * sint ) + o(i); //EED2 } //EED2 //EED2 //TODO: check if p and new_p are the same... //EED2 _contours[ i ]->addContourPoint( p ); //EED2 } // rof //EED2 } else if( typ == marParameters::SNAKE || typ == marParameters::DERICHE ) { //EED2 // ONLY Isocontours at the moment. //EED2 } // fi //EED2 //EED2 cntVTK->Delete( ); //EED2 cpd->Delete( ); //EED2 conn->Delete( ); //EED2 cpd2->Delete( ); //EED2 isoStrips->Delete( ); } // fi } */ // ---------------------------------------------------------------------------- void marAxis::createContour(int i, kVolume* vol) { int j,id; int numberOfPoints,numberOfCells; vtkCell* cell; vtkIdList* pids; double p[ 3 ]; vtkPolyData* vtkPolyData_2DContour; double x,y; int sizeIma = getParameters( )->getSizeIma( ); double dimIma = getParameters( )->getDimIma( ); // EED 24 Dec 2007 dimIma=sizeIma/3; _contours[i] = new marContour( i, getParameters( ) ) ; vtkPolyData_2DContour = get2Dcontour(i,vol); numberOfCells = vtkPolyData_2DContour->GetNumberOfCells( ); cell = vtkPolyData_2DContour->GetCell( 0 ); pids = cell->GetPointIds( ); numberOfPoints = pids->GetNumberOfIds( ); for( j = 0; j < numberOfPoints; j++ ) { id = pids->GetId( j ); vtkPolyData_2DContour->GetPoint( id, p); x=p[0]-64.0; y=p[1]-64.0; x=x * dimIma / ( float ) sizeIma; y=y * dimIma / ( float ) sizeIma; _contours[ i ]->addContourPoint( x , y ); } _contours[i]->do_spline(); _contours[i]->calculateVariables(); } // ---------------------------------------------------------------------------- void marAxis::EraseContour(int i){ if (_3Dcontour[i]!=NULL){ _3Dcontour[i]->Delete(); _3Dcontour[i]=NULL; } if (_2DDiameterMax[i]!=NULL){ _2DDiameterMax[i]->Delete(); _2DDiameterMax[i]=NULL; } if (_2DDiameterMin[i]!=NULL){ _2DDiameterMin[i]->Delete(); _2DDiameterMin[i]=NULL; } if (_2Dcontours[i]!=NULL){ _2Dcontours[i]->Delete(); _2Dcontours[i]=NULL; } if (_contours[i]!=NULL){ delete _contours[i]; _contours[i]=NULL; } } // ---------------------------------------------------------------------------- void marAxis::replaceContour2D(int i,int size,double *vx,double *vy){ EraseContour(i); vtkPoints *_pts = vtkPoints::New(); _pts->SetNumberOfPoints(size); int j; for (j=0 ; jSetPoint(j, vx[j] , vy[j] , 0 ); } // _pts->SetPoint(0, vx[0] , vy[0] , 0 ); vtkCellArray *lines = vtkCellArray::New(); lines->InsertNextCell( size ); for ( j=0 ; jInsertCellPoint(j % size ); } vtkPolyData *_pd = vtkPolyData::New(); _pd->SetPoints( _pts ); _pd->SetLines( lines ); lines->Delete(); //do not delete lines ?? _pts->Delete(); _2Dcontours[i]=_pd; createContour(i,NULL); } // ---------------------------------------------------------------------------- void marAxis::createSlice(int i , kVolume* vol){ int k=i; int sizeIma = getParameters( )->getSizeIma( ); vtkStructuredPoints* stPoints = vtkStructuredPoints::New( ); double *p, *n; p = _points[k]; n = getNormal( k ); stPoints->GetPointData( )->SetScalars( get3DSlice(k,vol)->GetOutput()->GetPointData()->GetScalars() ); stPoints->SetDimensions( sizeIma, sizeIma, 1 ); stPoints->SetScalarType( vol->castVtk( )->GetDataObjectType( ) ); stPoints->SetScalarTypeToShort (); stPoints->Update(); vtkImageChangeInformation* change = vtkImageChangeInformation::New(); change->SetInput( stPoints ); change->Update(); //important if (_slices[k]!=NULL) { delete _slices[k]; } _slices[k] = new kVolume( change->GetOutput() ) ; stPoints->Delete( ); change->Delete(); } // ---------------------------------------------------------------------------- void marAxis::create3DSlice(int i , kVolume* vol ) { int k=i; int sizeIma = getParameters( )->getSizeIma( ); double dimIma = getParameters( )->getDimIma( ); dimIma=sizeIma/3; // double voxSize = getParameters( )->getVoxelSize( ); // if( forceCnt ) // for( unsigned int i = 0; i < _contours.size( ); i++ ) // delete _contours[ i ]; // _contours.clear( ); // for( unsigned int i = 0; i < _slices.size( ); i++ ) { // delete _slices[ i ]; // _3Dslices[ i ]->Delete( ); // _quantificationImages[ i ]->Delete( ); // } // _slices.clear( ); // _3Dslices.clear( ); // _quantificationImages.clear( ); // "Cutter" initialization // vtkStructuredPoints* stPoints = vtkStructuredPoints::New( ); // Perform "cut" double *p, *n; // for( int k = 0; k < nCnts; k++ ) { //s = ( double )k / ( double )( nCnts - 1 ); p = _points[k]; n = getNormal( k ); vtkPlaneSource* pSource = vtkPlaneSource::New( ); pSource->SetOrigin( p[ 0 ], p[ 1 ], p[ 2 ] ); pSource->SetPoint1( p[ 0 ] + dimIma - 1.0, p[ 1 ], p[ 2 ] ); pSource->SetPoint2( p[ 0 ], p[ 1 ], p[ 2 ] + dimIma - 1.0 ); pSource->SetResolution( sizeIma - 1, sizeIma - 1 ); pSource->Update( ); pSource->SetCenter( p[ 0 ], p[ 1 ], p[ 2 ] ); pSource->SetNormal( n[ 0 ], n[ 1 ], n[ 2 ] ); pSource->Update( ); if (_3Dslices[ k ]!=NULL){ _3Dslices[ k ]->Delete(); } _3Dslices[ k ] = vtkProbeFilter::New( ) ; _3Dslices[ k ]->SetInput( ( vtkDataSet* )pSource->GetOutput( ) ); _3Dslices[ k ]->SetSource( vol->castVtk( ) ); _3Dslices[ k ]->Update( ); pSource->Delete( ); // stPoints->GetPointData( )-> // SetScalars( _3Dslices[ k ]->GetOutput( )-> // GetPointData( )->GetScalars( ) ); // stPoints->SetDimensions( sizeIma, sizeIma, 1 ); // stPoints->SetScalarType( vol->castVtk( )->GetDataObjectType( ) ); // stPoints->Update( ); // // if( forceCnt ) // { // marContour* tmp_marContour = new marContour( k, getParameters( ) ) ; // tmp_marContour->calculateVariables(); // _contours.push_back( tmp_marContour ); // } // _slices.push_back( new kVolume( kVolume::UCHAR, 1, 1, 1 ) ); // *( _slices[ k ] ) = ( vtkImageData* )stPoints; //copy // _quantificationImages.push_back( (_slices[ k ])->castVtk( ) ); // } // rof // stPoints->Delete( ); // _3Dslices[i]=NULL; } // ---------------------------------------------------------------------------- void marAxis::create3Dcontour(int i, kVolume* vol) { vtkPoints *_vtkPoints; double *o; double *n; vtkMatrix4x4* mat; vtkTransform* trans; double nn[3]; double c[3],p1[3],p2[3]; double nA,nB,nC; _vtkPoints = vtkPoints::New(); o = _points[i]; n = getNormal( i ); mat = vtkMatrix4x4::New(); trans = vtkTransform::New(); nn[0]=n[0]; nn[1]=n[1]; nn[2]=n[2]; nC=sqrt( nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2] ); nn[0]=nn[0]/nC; nn[1]=nn[1]/nC; nn[2]=nn[2]/nC; vtkPlaneSource* pSource = vtkPlaneSource::New( ); pSource->SetOrigin( 0, 0 , 0 ); pSource->SetPoint1( 1, 0 , 0 ); pSource->SetPoint2( 0, 0 , 1.0 ); pSource->SetCenter( 0,0,0 ); pSource->SetNormal( nn[ 0 ], nn[ 1 ], nn[ 2 ] ); pSource->Update( ); // pSource->Update( ); pSource->GetOrigin(c); pSource->GetPoint1(p1); pSource->GetPoint2(p2); pSource->Delete(); p1[0]=p1[0]-c[0]; p1[1]=p1[1]-c[1]; p1[2]=p1[2]-c[2]; p2[0]=p2[0]-c[0]; p2[1]=p2[1]-c[1]; p2[2]=p2[2]-c[2]; nA=sqrt( p1[0]*p1[0] + p1[1]*p1[1] + p1[2]*p1[2] ); nB=sqrt( p2[0]*p2[0] + p2[1]*p2[1] + p2[2]*p2[2] ); p1[0]=p1[0]/nA; p1[1]=p1[1]/nA; p1[2]=p1[2]/nA; p2[0]=p2[0]/nB; p2[1]=p2[1]/nB; p2[2]=p2[2]/nB; mat->SetElement (0,0, nn[0]); mat->SetElement (1,0, nn[1]); mat->SetElement (2,0, nn[2]); mat->SetElement (3,0, 0); mat->SetElement (0,1, p2[0]); mat->SetElement (1,1, p2[1]); mat->SetElement (2,1, p2[2]); mat->SetElement (3,1, 0); mat->SetElement (0,2, p1[0]); mat->SetElement (1,2, p1[1]); mat->SetElement (2,2, p1[2]); mat->SetElement (3,2, 0); mat->SetElement (0,3, 0); mat->SetElement (1,3, 0); mat->SetElement (2,3, 0); mat->SetElement (3,3, 1); // double deter=mat->Determinant(mat); trans->Identity(); trans->SetMatrix(mat); float ang; ang=-90; trans->RotateWXYZ ( ang,0,1,0); // EED Borrame // double scale = 1; // trans->Scale ( scale , scale, scale); trans->Update(); vtkMatrix4x4 *m=vtkMatrix4x4::New(); trans->GetMatrix ( m ); int j,numberOfPoints; marContour* contour2D = getContour(i,vol); numberOfPoints = contour2D->getNumberOfPoints( ); double fpA[3]; double pA[3],ppp[3]; for( j = 0; j <= numberOfPoints; j++ ) { contour2D->getPoint( fpA,j % numberOfPoints); pA[0] = fpA[0] + 0.2; // correction EED pA[1] = fpA[1] + 0.2; // correction EED pA[2] = 0; // Why this does not works...(release version)???? // trans->TransformPoint(pA,ppp); // or // trans->TransformVector(pA,ppp); // so.. ppp[0] = m->GetElement(0,0)*pA[0] + m->GetElement(0,1)*pA[1] + m->GetElement(0,2)*pA[2] + m->GetElement(0,3); ppp[1] = m->GetElement(1,0)*pA[0] + m->GetElement(1,1)*pA[1] + m->GetElement(1,2)*pA[2] + m->GetElement(1,3); ppp[2] = m->GetElement(2,0)*pA[0] + m->GetElement(2,1)*pA[1] + m->GetElement(2,2)*pA[2] + m->GetElement(2,3); // ppp[0]=ppp[0]+o[0]-c[0]; ppp[1]=ppp[1]+o[1]-c[1]; ppp[2]=ppp[2]+o[2]-c[2]; ppp[0]=ppp[0]+o[0]; ppp[1]=ppp[1]+o[1]; ppp[2]=ppp[2]+o[2]; // ppp[0]=ppp[0]+o[0]-c[0]*0.781; ppp[1]=ppp[1]+o[1]-c[1]*0.781;; ppp[2]=ppp[2]+o[2]-c[2]*0.781;; _vtkPoints->InsertNextPoint( (float)ppp[0],(float)ppp[1],(float)ppp[2] ); } m->Delete(); mat->Delete(); trans->Delete(); if (_3Dcontour[i]!=NULL){ _3Dcontour[i]->Delete(); } _3Dcontour[i]=_vtkPoints; } // ---------------------------------------------------------------------------- void marAxis::createSliceImage(int i , kVolume* vol ) { _quantificationImages[i] = getSlice(i,vol)->castVtk( ) ; } // ---------------------------------------------------------------------------- void marAxis::create2Dcontour(int i , kVolume* vol ) { vtkImageData* imagedata; double thr = getParameters( )->getDoubleParam( marParameters::e_threshold_isocontour ); int sizeIma = getParameters( )->getSizeIma( ); double localint = getSignal( i ); double sigma = getParameters( )->getDoubleParam( marParameters::e_sigma ); double threshold = localint * thr / 100.0; int dimIma = ( int ) ( ( double ) sizeIma /getParameters( )->getDoubleParam( marParameters::e_scale ) ); int vmin = ( int )threshold; int vmax = ( int )threshold; int x = -1; int y = -1; x = ( x != -1 )? x: sizeIma / 2; y = ( y != -1 )? y: sizeIma / 2; // Isocontour imagedata = (vtkImageData*) ((getSlice( i , vol ))->castVtk( )); double *range = imagedata->GetScalarRange(); //vtkKitwareContourFilter* cntVTK = vtkKitwareContourFilter::New( ); vtkContourFilter* cntVTK = vtkContourFilter::New( ); cntVTK->SetInput( imagedata ); cntVTK->SetNumberOfContours( 1 ); //cntVTK->SetValue( 0, vmin ); cntVTK->SetValue( 0, (range[1]*thr/100) ); //cntVTK->SetValue( 1, vmax ); cntVTK->Update( ); vtkCleanPolyData* cpd = vtkCleanPolyData::New( ); cpd->SetInput( cntVTK->GetOutput( ) ); cpd->ConvertLinesToPointsOff( ); cpd->Update( ); vtkPolyDataConnectivityFilter* conn = vtkPolyDataConnectivityFilter::New( ); //conn->SetMaxRecursionDepth( 3000 ); // PS -> //conn->SetInput( cntVTK->GetOutput( ) ); PS conn->SetInput( cpd->GetOutput( ) ); conn->SetClosestPoint( x, y, 0 ); conn->SetExtractionModeToClosestPointRegion( ); conn->Update( ); vtkCleanPolyData* cpd2 = vtkCleanPolyData::New( ); cpd2->SetInput( conn->GetOutput( ) ); cpd2->Update(); vtkStripper* vtkstripper = vtkStripper::New( ); vtkstripper->SetInput( cpd2->GetOutput() ); vtkstripper->Update(); vtkPolyData* polyDataResult = vtkstripper->GetOutput(); polyDataResult->Update( ); cntVTK -> Delete(); cpd2 -> Delete(); cpd -> Delete(); conn -> Delete(); /* // EED 27 Oct 2007 double *p; int ii,size=polyDataResult->GetNumberOfPoints(); printf("EED marAxis::create2Dcontour size %d\n",size); for (ii=0;iiGetPoint(ii); printf("EED marAxis::create2Dcontour %d >> %f %f %f\n",ii,p[0], p[1], p[2]); p[2]=-100; } for (ii=0;iiGetPoint(ii); printf("EED marAxis::create2Dcontour xxxxx %d >> %f %f %f\n",ii,p[0], p[1], p[2]); } */ if ( _2Dcontours[i]!=NULL ) { _2Dcontours[i]->Delete(); } _2Dcontours[i]=polyDataResult; } // ---------------------------------------------------------------------------- void marAxis::create2DDiameterMin( int i , kVolume* vol ){ double p1[2],p2[2]; marContour *marcontour=getContour(i,vol); marcontour->getMinimumLine(p1,p2); vtkPoints *_vtkPoints = vtkPoints::New(); int sizeIma = getParameters( )->getSizeIma( ); double dimIma = getParameters( )->getDimIma( ); // EED 24 Dec 2007 dimIma=sizeIma/3; double factor=( float ) sizeIma / dimIma; p1[0]=p1[0]*factor+sizeIma/2; p1[1]=p1[1]*factor+sizeIma/2; p2[0]=p2[0]*factor+sizeIma/2; p2[1]=p2[1]*factor+sizeIma/2; _vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 ); _vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 ); if ( _2DDiameterMin[i]!=NULL ) { _2DDiameterMin[i]->Delete(); } _2DDiameterMin[i] = _vtkPoints; } // ---------------------------------------------------------------------------- void marAxis::create2DDiameterMax( int i , kVolume* vol ){ double p1[2],p2[2]; marContour *marcontour=getContour(i,vol); marcontour->getMaximumLine(p1,p2); vtkPoints *_vtkPoints = vtkPoints::New(); int sizeIma = getParameters( )->getSizeIma( ); double dimIma = getParameters( )->getDimIma( ); // EED 24 Dec 2007 dimIma=sizeIma/3; double factor=( float ) sizeIma / dimIma; p1[0]=p1[0]*factor+sizeIma/2; p1[1]=p1[1]*factor+sizeIma/2; p2[0]=p2[0]*factor+sizeIma/2; p2[1]=p2[1]*factor+sizeIma/2; _vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 ); _vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 ); if ( _2DDiameterMax[i]!=NULL ) { _2DDiameterMax[i]->Delete(); } _2DDiameterMax[i] = _vtkPoints; } // ---------------------------------------------------------------------------- void marAxis::createEmptyVectors() { int nCnts = ( int ) getNumberOfSplinePoints( ); clearAllVectors(); int i; // int xxxx=getNumberOfContours(); for (i=0;iDelete(); _3Dslices[i]=NULL; } } _3Dslices.clear(); for (i=0;i<_3Dcontour.size();i++){ // VISUALISATION_VTK 3D if (_3Dcontour[i]!=NULL) { _3Dcontour[i]->Delete(); _3Dcontour[i]=NULL; } } _3Dcontour.clear(); for (i=0;i<_quantificationImages.size();i++){ // VISUALISATION_VTK 2D if (_quantificationImages[i]!=NULL) { // _quantificationImages[i]->Delete(); // _quantificationImages[i]=NULL; } } _quantificationImages.clear(); for (i=0;i<_2Dcontours.size();i++){ // VISUALISATION_VTK 2D if (_2Dcontours[i]!=NULL) { // EED ??? This object was allready erased bye the VTK pipeline // _2Dcontours[i]->Delete(); _2Dcontours[i]=NULL; } } _2Dcontours.clear(); for (i=0;i<_2DDiameterMin.size();i++){ // VISUALISATION_VTK 2D if (_2DDiameterMin[i]!=NULL) { _2DDiameterMin[i]->Delete() ; _2DDiameterMin[i]=NULL; } } _2DDiameterMin.clear(); for (i=0;i<_2DDiameterMax.size();i++){ // VISUALISATION_VTK 2D if (_2DDiameterMax[i]!=NULL) { _2DDiameterMax[i]->Delete(); _2DDiameterMax[i]=NULL; } } _2DDiameterMax.clear(); } // ---------------------------------------------------------------------------- void marAxis::eraseContourVectorsContent() { unsigned int i; for (i=0;i<_contours.size();i++){ if (_contours[i]!=NULL) { // DATA-MODEL-2D delete _contours[i]; _contours[i]=NULL; } } for (i=0;i<_2Dcontours.size();i++){ // VISUALISATION_VTK 2D if (_2Dcontours[i]!=NULL) { _2Dcontours[i]->Delete(); _2Dcontours[i]=NULL; } } for (i=0;i<_2DDiameterMin.size();i++){ // VISUALISATION_VTK 2D if (_2DDiameterMin[i]!=NULL) { _2DDiameterMin[i]->Delete(); _2DDiameterMin[i]=NULL; } } for (i=0;i<_2DDiameterMax.size();i++){ // VISUALISATION_VTK 2D if (_2DDiameterMax[i]!=NULL) { _2DDiameterMax[i]->Delete(); _2DDiameterMax[i]=NULL; } } for (i=0;i<_3Dcontour.size();i++){ // VISUALISATION_VTK 3D if (_3Dcontour[i]!=NULL) { _3Dcontour[i]->Delete(); _3Dcontour[i]=NULL; } } } // ---------------------------------------------------------------------------- marContour* marAxis::getContour(int i, kVolume* vol) // DATA-MODEL 2D { if (_contours[i]==NULL){ createContour(i, vol); } return _contours[i]; } // ---------------------------------------------------------------------------- kVolume* marAxis::getSlice(int i, kVolume* vol) // DATA-MODEL-Voxel XxYx1 { if (_slices[i]==NULL){ createSlice( i , vol ); } return _slices[i]; } // ---------------------------------------------------------------------------- bool marAxis::if3DcontourExist(int i) { bool result=true; if (_3Dcontour[i]==NULL) { result=false; } return result; } // ---------------------------------------------------------------------------- void marAxis::Save3Dcontour(FILE *ff,int id) { int i,size; double point[3]; if (_3Dcontour[id]!=NULL) { //Ramiro fprintf(ff,"contour_id: %d \n", id); size = _3Dcontour[id]->GetNumberOfPoints(); //Ramiro fprintf(ff,"numberOfPoints: %d \n", size); for ( i=0 ; iGetPoint( i , point ); fprintf(ff,"%f %f %f %d\n", point[0], point[1],point[2],id); } } } // ---------------------------------------------------------------------------- void marAxis::SaveExisting3DContours(FILE *ff) { if (ff!=NULL) { // counting existing 3Dcontours int acum=0; int i,size=_3Dcontour.size(); for (i=0; ipEnd){ int tmp=pIni; pIni=pEnd; pEnd=tmp; } marVector* pO = new marVector(2); marVector* pF = new marVector(2); double L; unsigned int j; for( j = pIni, L = 0.0; j < pEnd-1; j++ ) { (*pO)=_points[j]; (*pF)=_points[j+1]; (*pF)=(*pF)-(*pO); L+=pF->norm2(); } // rof delete pO; delete pF; return L; } // ---------------------------------------------------------------------------- void marAxis::calculateTotalAxisLenght(){ _totalAxisLenght = getAxisLenght( 0 , _points.size( ) ); } // ---------------------------------------------------------------------------- double marAxis::getTotalLength(){ if (_totalAxisLenght==-1){ calculateTotalAxisLenght(); } return _totalAxisLenght; } // ---------------------------------------------------------------------------- void marAxis::calculateSubAxisLength(){ _subAxisLenght = getAxisLenght( _startQuant , _finishQuant ); } // ---------------------------------------------------------------------------- double marAxis::getSubAxisLength(){ if (_subAxisLenght==-1){ calculateSubAxisLength(); } return _subAxisLenght; } // ---------------------------------------------------------------------------- double marAxis::getAverageArea(int pIni, int pEnd, kVolume* vol){ marContour *marcontourHealthy; int ihealthySlice,itemp; double acumArea = 0; for (ihealthySlice = pIni ; ihealthySlice <= pEnd ; ihealthySlice++ ){ itemp=ihealthySlice; if (itemp<0) { itemp=0; } if (itemp>=_points.size( )) { itemp=_points.size( )-1; } marcontourHealthy = getContour( itemp , vol ); acumArea = acumArea + marcontourHealthy->getArea(); } return acumArea / (pEnd - pIni + 1); } // ---------------------------------------------------------------------------- void marAxis::calculateReferenceArea(kVolume* vol){ _referenceArea = getAverageArea(_healthySliceStart,_healthySliceEnd, vol); } // ---------------------------------------------------------------------------- double marAxis::getReferenceArea(kVolume* vol){ if (_referenceArea==-1){ calculateReferenceArea(vol); } return _referenceArea; } // ---------------------------------------------------------------------------- void marAxis::calculateReferenceAverDiam(kVolume* vol){ marContour *marcontourHealthy; int ihealthySlice; double acumDiam=0; for (ihealthySlice=_healthySliceStart; ihealthySlice<=_healthySliceEnd; ihealthySlice++){ marcontourHealthy = getContour( ihealthySlice , vol ); acumDiam = acumDiam + marcontourHealthy->getAverageDiameter(); } _referenceAverDiam = acumDiam / (_healthySliceEnd-_healthySliceStart+1); } // ---------------------------------------------------------------------------- double marAxis::getReferenceAverDiam(kVolume* vol){ if (_referenceAverDiam==-1){ calculateReferenceAverDiam(vol); } return _referenceAverDiam; } // ---------------------------------------------------------------------------- void marAxis::reset( ) { //EED Borrame // unsigned int i; kCurve::reset( ); _description = ""; _healthySlice = -1; _startQuant = -1; _finishQuant = -1; _actualQuant = -1; _subAxisLenght = -1; _referenceArea = -1; _referenceAverDiam = -1; Delete( ); // EED Borrame /* for( i = 0; i < _contours.size( ); i++ ) delete _contours[ i ]; _contours.clear( ); for( i = 0; i < _slices.size( ); i++ ) delete _slices[ i ]; _slices.clear( ); */ } // ---------------------------------------------------------------------------- void marAxis::copyFrom( const marObject& from ) { // TODO } // ---------------------------------------------------------------------------- bool marAxis::save( std::ofstream& os ) { double p[ INDX_count ]; unsigned int i; i = _description.length( ); os.write( ( const char* )&i, sizeof( int ) ); os.write( ( char* )_description.c_str( ), i * sizeof( char ) ); i = getNumberOfControlPoints( ); os.write( ( const char* )&i, sizeof( int ) ); for( i = 0; i < getNumberOfControlPoints( ); i++ ) { memcpy( p, _controlPoints[ i ], 3 * sizeof( double ) ); memcpy( p + 3, _controlState[ i ], ( INDX_count - 3 ) * sizeof( double ) ); os.write( ( const char* )p, INDX_count * sizeof( double ) ); } // rof i = getNumberOfContours( ); os.write( ( const char* )&i, sizeof( int ) ); for( i = 0; i < getNumberOfContours( ); i++ ) _contours[ i ]->save( os ); return( true ); } // ---------------------------------------------------------------------------- bool marAxis::load( std::ifstream& is ) { double p[ INDX_count ]; int i, n; reset( ); is.read( ( char* )&n, sizeof( int ) ); _description.resize( n ); is.read( ( char* )_description.c_str( ), n * sizeof( char ) ); is.read( ( char* )&n, sizeof( int ) ); for( i = 0; i < n; i++ ) { is.read( ( char* )p, INDX_count * sizeof( double ) ); addControlPoint( p, p + 3 ); } // rof is.read( ( char* )&n, sizeof( int ) ); for( i = 0; i < n; i++ ) { _contours.push_back( new marContour( 0, getParameters( ) ) ); _contours[ i ]->load( is ); } // rof return( true ); } // ---------------------------------------------------------------------------- vtkPolyData* marAxis::Draw( ) { unsigned int i, j; double *p; vtkPoints* allPoints = vtkPoints::New( ); vtkCellArray* allTopology = vtkCellArray::New( ); allTopology->InsertNextCell( _points.size( ) ); for( i = 0, j=0; i < _points.size( ); i++, j++ ) { p = _points[i]; allPoints->InsertNextPoint( p[ 0 ], p[ 1 ], p[ 2 ] ); allTopology->InsertCellPoint( j ); } // rof _allData = vtkPolyData::New( ); _allData->SetPoints( allPoints ); _allData->SetLines( allTopology ); allPoints->Delete(); allTopology->Delete(); return ( _allData ); } // ---------------------------------------------------------------------------- vtkPolyData *marAxis::GetAxisData( ) { unsigned int i; double point[ 3 ]; double radio; double p[marAxis::INDX_count]; vtkPoints *allPoints = vtkPoints::New( ); vtkCellArray *allTopology = vtkCellArray::New( ); vtkDoubleArray *allRadios = vtkDoubleArray::New( ); allRadios->SetName("radio"); for( i = 0; i < _controlPoints.size( ); i++ ) { getControlPoint(i, p, p+3); point[ 0 ] = p[marAxis::INDX_X]; point[ 1 ] = p[marAxis::INDX_Y]; point[ 2 ] = p[marAxis::INDX_Z]; radio = p[marAxis::INDX_RAYON]*2.0; allPoints -> InsertPoint( i, point[ 0 ] , point[ 1 ] , point[ 2 ] ); //para saber exactamante el indice que se le asigno allRadios -> InsertValue(i,radio); if(i>0){ allTopology->InsertNextCell(2); allTopology->InsertCellPoint(i); allTopology->InsertCellPoint(i-1); } } // rof _allData = vtkPolyData::New( ); _allData -> SetPoints( allPoints ); _allData -> SetLines( allTopology ); _allData -> GetPointData()->SetScalars(allRadios); return ( _allData ); } // ---------------------------------------------------------------------------- void marAxis::Delete( ){ clearAllVectors(); int i; for (i=0;i<_points.size();i++){ if (_points[i]!=NULL) { delete _points[i]; _points[i]=NULL; } } _points.clear(); _signal.clear(); if(_allData) { _allData->Delete(); _allData = NULL; } // reset(); } // ---------------------------------------------------------------------------- void marAxis::setHealthySlice( int hsS, int hs, int hsE ){ _healthySliceStart = hsS; _healthySlice = hs; _healthySliceEnd = hsE; _referenceArea=-1; _referenceAverDiam=-1; }; // ---------------------------------------------------------------------------- void marAxis::setStartQuant( int sq ){ _startQuant = sq; _subAxisLenght = -1; }; // ---------------------------------------------------------------------------- void marAxis::setFinishQuant( int fq ){ _finishQuant = fq; _subAxisLenght = -1; }; // ---------------------------------------------------------------------------- void marAxis::AddPointToList(double x, double y, double z, int signal) { double *p; p=(double *)malloc(sizeof(double)*3); p[0]=x; p[1]=y; p[2]=z; _points.push_back( p ); _signal.push_back( signal ); } // eof - axis.cxx