1 /*=========================================================================
4 Module: $RCSfile: marContour.cpp,v $
6 Date: $Date: 2008/10/31 16:32:54 $
7 Version: $Revision: 1.1 $
9 Copyright: (c) 2002, 2003
12 This software is distributed WITHOUT ANY WARRANTY; without even
13 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14 PURPOSE. See the above copyright notice for more information.
16 =========================================================================*/
24 //#include <gsl/gsl_blas.h> //PS
25 #include "marVector.h"
26 #include <math.h> //PS
27 #include "marMathConst.h" //PS
29 #include "marContour.h"
30 #include <vtkPolyLine.h>
31 #include <vtkKochanekSpline.h>
33 #include "load_snake.h"
35 // -------------------------------------------------------------------------
36 marContour::marContour( double s, marParameters* p )
42 // -------------------------------------------------------------------------
44 void marContour::addContourPoint( double x, double y )
46 double* np = new double[ 2 ];
48 // memcpy( np, p, 2 * sizeof( double ) );
51 _points.push_back( np );
54 // -------------------------------------------------------------------------
55 void marContour::resample( )
59 // -------------------------------------------------------------------------
60 bool marContour::isCalculated( )
62 return( _points.size( ) > 0 );
65 // -------------------------------------------------------------------------
66 int marContour::getNumberOfPoints( )
68 return( _points.size( ) );
71 // -------------------------------------------------------------------------
72 void marContour::getPoint( double* p, int i )
74 memcpy( p, _points[ i ], 2 * sizeof( double ) );
76 // -------------------------------------------------------------------------
77 void marContour::calculateVariables()
82 _diameterFromArea = -1;
83 _diameterFromPerimeter = -1;
84 _minimumDiameter = -1;
85 _maximumDiameter = -1;
86 _averageDiameter = -1;
91 calculateDiameterFromArea();
92 calculateDiameterFromPerimeter();
93 calculateMinimumMaximumDiameter();
94 calculateAverageDiameter();
97 /* METHOD DESCRIPTION ****************************************************
99 * marContour::do_spline (method) *
101 * DESCRIPTION : Calculates a spline function in order to add points to *
104 * SYNTAX : ctn.do_spline( ) *
106 ******************************************************* END DESCRIPTION */
107 void marContour::do_spline( ) {
109 vtkKochanekSpline* CntSplineX = vtkKochanekSpline::New( );
110 vtkKochanekSpline* CntSplineY = vtkKochanekSpline::New( );
112 //EED temPoint< float, 2 >* tmp;
116 unsigned int nps = 1000; // TODO: 1000?????
117 CntSplineX->SetDefaultTension( 0 ); CntSplineX->SetDefaultBias( 0 ); CntSplineX->SetDefaultContinuity( 0 );
118 CntSplineY->SetDefaultTension( 0 ); CntSplineY->SetDefaultBias( 0 ); CntSplineY->SetDefaultContinuity( 0 );
119 for( i = 0; i < _points.size( ); i++ ) {
120 CntSplineX->AddPoint( i, ( _points[ i ] )[ 0 ] );
121 CntSplineY->AddPoint( i, ( _points[ i ] )[ 1 ] );
125 CntSplineX->AddPoint( i, ( _points[ 0 ] )[ 0 ] );
126 CntSplineY->AddPoint( i, ( _points[ 0 ] )[ 1 ] );
128 to = _points.size( );
129 while( _no_spline_points.size( ) > 0 ) {
130 delete _no_spline_points[ _no_spline_points.size( ) - 1 ];
131 _no_spline_points.pop_back( );
133 for( i = 0; i < _points.size( ); i++ ) {
134 //EED tmp = new temPoint< float, 2 >;
135 tmp = (double*)malloc ( 2*sizeof(double) );
136 ( tmp )[ 0 ] = ( _points[ i ] )[ 0 ];
137 ( tmp )[ 1 ] = ( _points[ i ] )[ 1 ];
138 _no_spline_points.push_back( tmp );
141 //EED this->clean( );
142 //EED for( i = 0; i < nps; i++ ) {
143 //EED tmp = new temPoint< float, 2 >;
144 //EED t = ( float ) ( to - 1 ) / ( float ) ( nps - 1 ) * ( float ) i;
145 //EED ( *tmp )[ 0 ] = CntSplineX->Evaluate( t );
146 //EED ( *tmp )[ 1 ] = CntSplineY->Evaluate( t );
147 //EED _points.push_back( tmp );
149 //EED _updated_measures = false;
151 CntSplineX->Delete();
152 CntSplineY->Delete();
157 // -------------------------------------------------------------------------
158 void marContour::calculateArea( )
163 // This uses Green's theorem:
164 // A = 1/2 * sum( xiyi+1 - xi+1yi); pO == pN
165 // A < 0 -> A = |A| (a negative value could raise because points are
166 // given in clockwise order).
167 for( i = 0, area = 0.0; i < _points.size( ); i++ ) {
169 j = ( i + 1 ) % _points.size( );
173 ( _points[ i ][ 0 ] * _points[ j ][ 1 ] ) -
174 ( _points[ j ][ 0 ] * _points[ i ][ 1 ] );
191 // -------------------------------------------------------------------------
192 void marContour::calculatePerimeter( )
194 // PS -> //gsl_vector* pO = gsl_vector_alloc( 2 ); //PS
195 marVector* pO = new marVector(2);
196 // PS -> //gsl_vector* pF = gsl_vector_alloc( 2 ); //PS
197 marVector* pF = new marVector(2);
201 int nbPoints=_points.size( );
203 for( j = 0, L = 0.0; j < nbPoints; j++ ) {
205 // PS -> //memcpy( pO->data, _points[ j ], 2 * sizeof( double ) ); //PS
207 // PS -> //memcpy( pF->data, _points[ j + 1 ], 2 * sizeof( double ) ); //PS
208 (*pF)=_points[(j+1)%nbPoints];
210 // PS -> // gsl_vector_sub( pF, pO );//PS
212 // PS -> //L += gsl_blas_dnrm2( pF );//PS
217 // PS -> //memcpy( pO->data, _points[ j - 1 ], 2 * sizeof( double ) ); //PS
219 // PS -> //memcpy( pF->data, _points[ 0 ], 2 * sizeof( double ) ); //PS
223 // PS -> // gsl_vector_sub( pF, pO );//PS
225 // PS -> //L += gsl_blas_dnrm2( pF );//PS
228 // PS -> //gsl_vector_free( pO );//PS
229 // PS -> //gsl_vector_free( pF );//PS
237 // -------------------------------------------------------------------------
238 /* METHOD DESCRIPTION ****************************************************
240 * marContour::recalcule_measures (method) *
242 * DESCRIPTION : Does one step to calcule measures. Don't call this *
243 * outside this file. *
245 * SYNTAX : this->recalcule_measures( ) *
247 ******************************************************* END DESCRIPTION */
248 void marContour::calculateMinimumMaximumDiameter( )
251 double _minimum_diameter;
252 double _maximum_diameter;
260 float tmp1[2], tmp2[2];
263 vtkUnstructuredGrid* grid;
264 PCONTOUR_FLOAT ctr, ctr2;
269 //EED if( !_updated_measures ) {
272 //EED _perimeter = 0.0;
273 _minimum_diameter = 0.0;
274 _maximum_diameter = 0.0;
275 _p1_min[ 0 ] = _p1_min[ 1 ] = 0.0;
276 _p2_min[ 0 ] = _p2_min[ 1 ] = 0.0;
277 _p1_max[ 0 ] = _p1_max[ 1 ] = 0.0;
278 _p2_max[ 0 ] = _p2_max[ 1 ] = 0.0;
279 //EED // This uses Green's theorem:
280 //EED // A = 1/2 * sum( xiyi+1 - xi+1yi); pO == pN
281 //EED // A < 0 -> A = |A| (a negative value is because points are given in
282 //EED // clockwise order).
284 //EED // P = sum( |Pj_Pj+1| ); pO == pN
286 //EED // Prepares data for maximum & minimum process.
287 //EED // The scale factor must be applied.
288 points = vtkPoints::New( );
289 poly = vtkPolyLine::New( );
290 grid = vtkUnstructuredGrid::New( );
292 poly->GetPointIds( )->SetNumberOfIds( _points.size( ) );
293 ctr2 = ( PCONTOUR_FLOAT ) IdCntAlloc( _points.size( ), CNT_FLOAT );
295 for( i = 0, j = 1; i < _points.size( ); i++, j = ( j + 1 ) % _points.size( ) ) {
296 tmp1[0] = _points[ i ][0];
297 tmp1[1] = _points[ i ][1];
298 //EED tmp1 = *_points[ i ];
299 tmp2[0] = _points[ j ][0];
300 tmp2[1] = _points[ j ][1];
301 //EED tmp2 = *_points[ j ];
302 // VTK stuff for maximum & minimum
303 poly->GetPointIds( )->SetId( i, i );
304 points->InsertNextPoint( tmp1[ 0 ], tmp1[ 1 ], 1 );
305 // Point scale factor
306 //EED tmp1 *= ( _parameters->get_dim_ima( ) / ( float ) _parameters->get_size_ima( ) );
307 //EED tmp2 *= ( _parameters->get_dim_ima( ) / ( float ) _parameters->get_size_ima( ) );
309 //EED IdCntAddPointG( ctr2, tmp1[ 0 ], tmp1[ 1 ] );
310 //EED _area += ( ( tmp1[ 0 ] * tmp2[ 1 ] ) - ( tmp2[ 0 ] * tmp1[ 1 ] ) );
311 //EED _perimeter += tmp1 * tmp2;
314 //EED _area = ( float ) fabs( ( double ) _area );
315 // TODO : Just testing: LibIDO is better than green theorem's?
316 // which is the difference?
317 //EED _area = ( float ) IdCntCalculAireG( ( PCONTOUR ) ctr2 );
319 ctr = ( PCONTOUR_FLOAT ) IdCntAlloc( _no_spline_points.size( ), CNT_FLOAT );
320 for( i = 0; i < _no_spline_points.size( ); i++ ) {
323 ( _no_spline_points[ i ] )[ 0 ],
324 ( _no_spline_points[ i ] )[ 1 ]
329 if( _points.size( ) > 0 ) {
330 grid->Allocate( 1, 1 );
331 grid->InsertNextCell( poly->GetCellType( ), poly->GetPointIds( ) );
332 grid->SetPoints( points );
335 _no_spline_points.size( ),
336 ( int ) grid->GetCell( 0 ),
337 &_minimum_diameter, &_maximum_diameter,
338 &( _p1_min[ 0 ] ), &( _p1_min[ 1 ] ),
339 &( _p2_min[ 0 ] ), &( _p2_min[ 1 ] ),
340 &( _p1_max[ 0 ] ), &( _p1_max[ 1 ] ),
341 &( _p2_max[ 0 ] ), &( _p2_max[ 1 ] ),
344 _min[0]=_p1_min[ 0 ];
345 _min[1]=_p1_min[ 1 ];
346 _min[2]=_p2_min[ 0 ];
347 _min[3]=_p2_min[ 1 ];
348 _max[0]=_p1_max[ 0 ];
349 _max[1]=_p1_max[ 1 ];
350 _max[2]=_p2_max[ 0 ];
351 _max[3]=_p2_max[ 1 ];
354 // float param_dim =(float) _parameters->get_dim_ima( );
355 // float param_size =(float) _parameters->get_size_ima( );
356 // float relation=param_dim / param_size
357 // _minimum_diameter *= relation;
358 // _maximum_diameter *= relation;
361 //EED _updated_measures = true;
365 _minimumDiameter = _minimum_diameter;
366 _maximumDiameter = _maximum_diameter;
375 // -------------------------------------------------------------------------
376 void marContour::getMaximumLine( double* pO, double* pF )
378 memcpy( pO, _max, 2 * sizeof( double ) );
379 memcpy( pF, _max + 2, 2 * sizeof( double ) );
382 // -------------------------------------------------------------------------
383 void marContour::getMinimumLine( double* pO, double* pF )
385 memcpy( pO, _min, 2 * sizeof( double ) );
386 memcpy( pF, _min + 2, 2 * sizeof( double ) );
389 // -------------------------------------------------------------------------
390 void marContour::reset( )
394 for( i = 0; i < _points.size( ); i++ )
398 for( i = 0; i < _no_spline_points.size( ); i++ )
399 delete _no_spline_points[ i ];
400 _no_spline_points.clear( );
402 _min[ 0 ] = _min[ 1 ] = _min[ 2 ] = _min[ 3 ] = 0.0;
403 _max[ 0 ] = _max[ 1 ] = _max[ 2 ] = _max[ 3 ] = 0.0;
406 // -------------------------------------------------------------------------
407 void marContour::copyFrom( const marObject& from )
411 marContour* f = &( ( marContour& )from );
414 memcpy( _min, f->_min, 4 * sizeof( double ) );
415 memcpy( _max, f->_max, 4 * sizeof( double ) );
416 for( i = 0; i < f->_points.size( ); i++ ) {
418 tmp = new double[ 3 ];
419 memcpy( tmp, f->_points[ i ], 3 * sizeof( double ) );
420 _points.push_back( tmp );
425 // -------------------------------------------------------------------------
426 bool marContour::save( std::ofstream& os )
431 os.write( ( const char* )&i, sizeof( int ) );
432 for( i = 0; i < _points.size( ); i++ )
433 os.write( ( const char* )_points[ i ], 3 * sizeof( double ) );
438 // -------------------------------------------------------------------------
439 bool marContour::load( std::ifstream& is )
445 is.read( ( char* )&size, sizeof( int ) );
446 for( i = 0; i < size; i++ ) {
448 tmp = new double[ 3 ];
449 is.read( ( char* )tmp, 3 * sizeof( double ) );
450 _points.push_back( tmp );
458 // ----------------------------------------------------------------------------
459 vtkUnstructuredGrid* marContour::Draw( )
464 vtkPoints* PointsContour;
465 vtkPolyLine* aPolyLine;
468 taille = (int) (this->getNumberOfPoints( ) / 2);
470 PointsContour = vtkPoints::New( );
472 aPolyLine = vtkPolyLine::New( );
473 ( aPolyLine->GetPointIds( ) )->SetNumberOfIds( taille );
474 for ( i = 0; i < taille; i++) {
475 ( aPolyLine->GetPointIds( ) )->SetId( i , i);
476 getPoint((double *) p, i*2);
477 PointsContour->InsertNextPoint( p[ 0 ], p[ 1 ], 0 );
480 _allData = vtkUnstructuredGrid::New( );
481 _allData->Allocate( 1 , 1 );
482 _allData->InsertNextCell( aPolyLine->GetCellType( ) ,
483 aPolyLine->GetPointIds( ) );
484 _allData->SetPoints( PointsContour );
486 PointsContour->Delete( );
487 aPolyLine->Delete( );
491 // ----------------------------------------------------------------------------
492 void marContour::Delete( )