lib/maracasVisuLib/src/kernel/marContour.cpp

   1 /*=========================================================================
   2
   3  Program:   wxMaracas
   4  Module:    $RCSfile: marContour.cpp,v $
   5  Language:  C++
   6  Date:      $Date: 2008/10/31 16:32:54 $
   7  Version:   $Revision: 1.1 $
   8
   9   Copyright: (c) 2002, 2003
  10   License:
  11
  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.
  15
  16 =========================================================================*/
  17
  18
  19
  20 extern "C" {
  21 #include"idcnt.h"
  22 };
  23
  24 //#include <gsl/gsl_blas.h> //PS
  25 #include "marVector.h"
  26 #include <math.h> //PS
  27 #include "marMathConst.h" //PS
  28
  29 #include "marContour.h"
  30 #include <vtkPolyLine.h>
  31 #include <vtkKochanekSpline.h>
  32
  33 #include "load_snake.h"
  34
  35 // -------------------------------------------------------------------------
  36 marContour::marContour( double s, marParameters* p )
  37 : marObject( p ),
  38 _s( s )
  39 {
  40 }
  41
  42 // -------------------------------------------------------------------------
  43
  44 void marContour::addContourPoint(  double x, double y  )
  45 {
  46     double* np = new double[ 2 ];
  47
  48 //    memcpy( np, p, 2 * sizeof( double ) );
  49         np[0]=x;
  50         np[1]=y;
  51     _points.push_back( np );
  52 }
  53
  54 // -------------------------------------------------------------------------
  55 void marContour::resample( )
  56 { // TODO
  57 }
  58
  59 // -------------------------------------------------------------------------
  60 bool marContour::isCalculated( )
  61 {
  62     return( _points.size( ) > 0 );
  63 }
  64
  65 // -------------------------------------------------------------------------
  66 int marContour::getNumberOfPoints( )
  67 {
  68     return( _points.size( ) );
  69 }
  70
  71 // -------------------------------------------------------------------------
  72 void marContour::getPoint( double* p, int i )
  73 {
  74     memcpy( p, _points[ i ], 2 * sizeof( double ) );
  75 }
  76 // -------------------------------------------------------------------------
  77 void marContour::calculateVariables()
  78 {
  79
  80          _area                                  =       -1;
  81          _perimeter                             =       -1;
  82          _diameterFromArea              =       -1;
  83          _diameterFromPerimeter =       -1;
  84          _minimumDiameter               =       -1;
  85          _maximumDiameter               =       -1;
  86          _averageDiameter               =       -1;
  87
  88
  89         calculateArea();
  90         calculatePerimeter();
  91         calculateDiameterFromArea();
  92         calculateDiameterFromPerimeter();
  93         calculateMinimumMaximumDiameter();
  94         calculateAverageDiameter();
  95
  96 }
  97 /* METHOD DESCRIPTION ****************************************************
  98  *                                                                       *
  99  * marContour::do_spline (method)                                        *
 100  *                                                                       *
 101  * DESCRIPTION : Calculates a spline function in order to add points to  *
 102  *               actual points.                                          *
 103  *                                                                       *
 104  * SYNTAX : ctn.do_spline( )                                             *
 105  *                                                                       *
 106  ******************************************************* END DESCRIPTION */
 107 void marContour::do_spline( ) {
 108     unsigned int i, to;
 109     vtkKochanekSpline* CntSplineX = vtkKochanekSpline::New( );
 110     vtkKochanekSpline* CntSplineY = vtkKochanekSpline::New( );
 111
 112 //EED   temPoint< float, 2 >* tmp;
 113         double *tmp;
 114
 115 //EED    float t;
 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 ] );
 122     } //  rof
 123
 124     if( i > 0 ) {
 125         CntSplineX->AddPoint( i, ( _points[ 0 ] )[ 0 ] );
 126         CntSplineY->AddPoint( i, ( _points[ 0 ] )[ 1 ] );
 127     } // fi
 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( );
 132     } // fwhile
 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 );
 139     } // rof
 140
 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 );
 148 //EED    } // rof
 149 //EED    _updated_measures = false;
 150
 151         CntSplineX->Delete();
 152         CntSplineY->Delete();
 153
 154
 155 }
 156
 157 // -------------------------------------------------------------------------
 158 void marContour::calculateArea( )
 159 {
 160     double area;
 161     int i, j;
 162
 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++ ) {
 168
 169         j = ( i + 1 ) % _points.size( );
 170
 171         //  Area
 172         area +=
 173             ( _points[ i ][ 0 ] * _points[ j ][ 1 ] ) -
 174             ( _points[ j ][ 0 ] * _points[ i ][ 1 ] );
 175
 176
 177     } // rof
 178     area /= 2.0;
 179     area = fabs( area );
 180
 181     this->_area = area;
 182 }
 183
 184
 185
 186
 187
 188
 189
 190
 191 // -------------------------------------------------------------------------
 192 void marContour::calculatePerimeter( )
 193 {
 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);
 198     double L;
 199     unsigned int j;
 200
 201         int nbPoints=_points.size( );
 202
 203     for( j = 0, L = 0.0; j < nbPoints; j++ ) {
 204
 205                 // PS ->        //memcpy( pO->data, _points[ j ], 2 * sizeof( double ) ); //PS
 206                 (*pO)=_points[j];
 207                 // PS ->        //memcpy( pF->data, _points[ j + 1 ], 2 * sizeof( double ) ); //PS
 208                 (*pF)=_points[(j+1)%nbPoints];
 209
 210                 // PS ->        // gsl_vector_sub( pF, pO );//PS
 211                 (*pF)=(*pF)-(*pO);
 212                 // PS ->        //L += gsl_blas_dnrm2( pF );//PS
 213                 L+=pF->norm2();
 214
 215     } // rof
 216         /*
 217         // PS ->     //memcpy( pO->data, _points[ j - 1 ], 2 * sizeof( double ) ); //PS
 218      (*pO)=_points[j-1];
 219         // PS ->        //memcpy( pF->data, _points[ 0 ], 2 * sizeof( double ) ); //PS
 220         (*pF)=_points[0];
 221
 222
 223         // PS ->        // gsl_vector_sub( pF, pO );//PS
 224         (*pF)=(*pF)-(*pO);
 225         // PS ->        //L += gsl_blas_dnrm2( pF );//PS
 226         L+=pF->norm2();
 227
 228         // PS ->     //gsl_vector_free( pO );//PS
 229         // PS ->     //gsl_vector_free( pF );//PS
 230         */
 231         delete pO;
 232         delete pF;
 233
 234     _perimeter=L ;
 235 }
 236
 237 // -------------------------------------------------------------------------
 238 /* METHOD DESCRIPTION ****************************************************
 239  *                                                                       *
 240  * marContour::recalcule_measures (method)                               *
 241  *                                                                       *
 242  * DESCRIPTION : Does one step to calcule measures. Don't call this      *
 243  *               outside this file.                                      *
 244  *                                                                       *
 245  * SYNTAX : this->recalcule_measures( )                                  *
 246  *                                                                       *
 247  ******************************************************* END DESCRIPTION */
 248 void marContour::calculateMinimumMaximumDiameter( )
 249 {
 250
 251         double _minimum_diameter;
 252         double _maximum_diameter;
 253         double _p1_min[2];
 254         double _p2_min[2];
 255         double _p1_max[2];
 256         double _p2_max[2];
 257 //-----
 258
 259     unsigned int i, j;
 260     float  tmp1[2], tmp2[2];
 261     vtkPoints* points;
 262     vtkPolyLine* poly;
 263     vtkUnstructuredGrid* grid;
 264     PCONTOUR_FLOAT ctr, ctr2;
 265
 266     double mom_min;
 267     double mom_max;
 268
 269 //EED    if( !_updated_measures ) {
 270         //  Preparation
 271 //EED        _area = 0.0;
 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).
 283 //EED        //
 284 //EED        //  P = sum( |Pj_Pj+1| ); pO == pN
 285 //EED        //
 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( );
 291
 292         poly->GetPointIds( )->SetNumberOfIds( _points.size( ) );
 293         ctr2 = ( PCONTOUR_FLOAT ) IdCntAlloc( _points.size( ), CNT_FLOAT );
 294
 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( ) );
 308             //  Area & perimeter
 309 //EED            IdCntAddPointG( ctr2, tmp1[ 0 ], tmp1[ 1 ] );
 310 //EED            _area += ( ( tmp1[ 0 ] * tmp2[ 1 ] ) - ( tmp2[ 0 ] * tmp1[ 1 ] ) );
 311 //EED            _perimeter += tmp1 * tmp2;
 312         } // rof
 313 //EED        _area *= 0.5;
 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 );
 318         //  Libido contour
 319         ctr = ( PCONTOUR_FLOAT ) IdCntAlloc( _no_spline_points.size( ), CNT_FLOAT );
 320         for( i = 0; i < _no_spline_points.size( ); i++ ) {
 321             IdCntAddPointG(
 322                 ctr,
 323                 ( _no_spline_points[ i ] )[ 0 ],
 324                 ( _no_spline_points[ i ] )[ 1 ]
 325             );
 326         } // rof
 327         //  Diameters...
 328
 329         if( _points.size( ) > 0 ) {
 330             grid->Allocate( 1, 1 );
 331             grid->InsertNextCell( poly->GetCellType( ), poly->GetPointIds( ) );
 332             grid->SetPoints( points );
 333             axe_max_min(
 334                 ctr,
 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 ] ),
 342                 &mom_min, &mom_max
 343             );
 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 ];
 352
 353 //EED
 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;
 359
 360         } // fi
 361 //EED     _updated_measures = true;
 362 //EED   } // fi
 363 //}
 364 //-----
 365         _minimumDiameter = _minimum_diameter;
 366         _maximumDiameter = _maximum_diameter;
 367
 368         points->Delete();
 369         poly->Delete();
 370         grid->Delete();
 371
 372 }
 373
 374
 375 // -------------------------------------------------------------------------
 376 void marContour::getMaximumLine( double* pO, double* pF )
 377 {
 378     memcpy( pO, _max, 2 * sizeof( double ) );
 379     memcpy( pF, _max + 2, 2 * sizeof( double ) );
 380 }
 381
 382 // -------------------------------------------------------------------------
 383 void marContour::getMinimumLine( double* pO, double* pF )
 384 {
 385     memcpy( pO, _min, 2 * sizeof( double ) );
 386     memcpy( pF, _min + 2, 2 * sizeof( double ) );
 387 }
 388
 389 // -------------------------------------------------------------------------
 390 void marContour::reset( )
 391 {
 392     int i;
 393
 394     for( i = 0; i < _points.size( ); i++ )
 395                 delete _points[ i ];
 396     _points.clear( );
 397
 398     for( i = 0; i < _no_spline_points.size( ); i++ )
 399                 delete _no_spline_points[ i ];
 400     _no_spline_points.clear( );
 401
 402     _min[ 0 ] = _min[ 1 ] = _min[ 2 ] = _min[ 3 ] = 0.0;
 403     _max[ 0 ] = _max[ 1 ] = _max[ 2 ] = _max[ 3 ] = 0.0;
 404 }
 405
 406 // -------------------------------------------------------------------------
 407 void marContour::copyFrom( const marObject& from )
 408 {
 409     double* tmp;
 410     int i;
 411     marContour* f = &( ( marContour& )from );
 412
 413     reset( );
 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++ ) {
 417
 418                 tmp = new double[ 3 ];
 419                 memcpy( tmp, f->_points[ i ], 3 * sizeof( double ) );
 420                 _points.push_back( tmp );
 421
 422     } // rof
 423 }
 424
 425 // -------------------------------------------------------------------------
 426 bool marContour::save( std::ofstream& os )
 427 {
 428     int i;
 429
 430     i = _points.size( );
 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 ) );
 434
 435     return( true );
 436 }
 437
 438 // -------------------------------------------------------------------------
 439 bool marContour::load( std::ifstream& is )
 440 {
 441     int i, size;
 442     double* tmp;
 443
 444     reset( );
 445     is.read( ( char* )&size, sizeof( int ) );
 446     for( i = 0; i < size; i++ ) {
 447
 448                 tmp = new double[ 3 ];
 449                 is.read( ( char* )tmp, 3 * sizeof( double ) );
 450                 _points.push_back( tmp );
 451
 452     } // rof
 453
 454     return( true );
 455 }
 456
 457
 458 // ----------------------------------------------------------------------------
 459 vtkUnstructuredGrid* marContour::Draw( )
 460 {
 461         int taille;
 462         double p[2];
 463
 464         vtkPoints* PointsContour;
 465         vtkPolyLine* aPolyLine;
 466         int i;
 467
 468         taille = (int) (this->getNumberOfPoints( ) / 2);
 469
 470         PointsContour = vtkPoints::New( );
 471
 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 );
 478         }
 479
 480         _allData = vtkUnstructuredGrid::New( );
 481         _allData->Allocate( 1 , 1 );
 482         _allData->InsertNextCell( aPolyLine->GetCellType( ) ,
 483                 aPolyLine->GetPointIds( ) );
 484         _allData->SetPoints( PointsContour );
 485
 486         PointsContour->Delete( );
 487         aPolyLine->Delete( );
 488         return ( _allData );
 489 }
 490
 491 // ----------------------------------------------------------------------------
 492 void marContour::Delete( )
 493 {
 494         if(_allData) {
 495                 _allData->Delete();
 496                 _allData = NULL;
 497         }
 498
 499         reset();
 500 }
 501
 502
 503 // eof - contour.cxx