lib/maracasVisuLib/src/kernel/include/marVector.cpp

   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 // marVector.cpp: implementation of the marVector class.
  27 //
  28 //////////////////////////////////////////////////////////////////////
  29 #include "marVector.h"
  30 #include <assert.h>
  31 #include <math.h>
  32
  33
  34 //////////////////////////////////////////////////////////////////////
  35 // Constructeurs
  36 //////////////////////////////////////////////////////////////////////
  37 marVector::marVector(size_t s)
  38         :shallowCopy(false)
  39 {
  40         if (s!=0)
  41         {
  42                 _data=new double[s];
  43                 _size=s;
  44         }
  45         else
  46         {
  47                 _data=NULL;
  48                 _size=0;
  49                 assert(false);
  50         }
  51 }
  52
  53 marVector::marVector(const marVector &v)
  54         :shallowCopy(false)
  55 {
  56         _size=v._size;
  57         _data=new double[_size];
  58         int i;
  59         for (i=0;i<_size;i++)
  60         {
  61                 _data[i]=v._data[i];
  62         }
  63 }
  64
  65 marVector::marVector(double*v,size_t s)
  66         :shallowCopy(true)
  67 {
  68         if (s!=0)
  69         {
  70                 _size=s;
  71                 _data=v;
  72 // PS ->                _data=new double[_size];
  73 // PS ->                for (int i=0;i<_size;i++)
  74 // PS ->                {
  75 // PS ->                        _data[i]=v[i];
  76 // PS ->                }
  77
  78         }
  79         else
  80         {
  81                 _data=NULL;
  82                 _size=0;
  83                 assert(false);
  84         }
  85
  86 }
  87
  88 //////////////////////////////////////////////////////////////////////
  89 // Destructeur
  90 //////////////////////////////////////////////////////////////////////
  91 marVector::~marVector()
  92 {
  93         if (!shallowCopy)
  94         {
  95                 delete _data;
  96         }
  97         else
  98         {
  99                 _data=NULL;
 100         }
 101 }
 102
 103 //////////////////////////////////////////////////////////////////////
 104 // Op�rateurs
 105 //////////////////////////////////////////////////////////////////////
 106 // ---------------------------------------------------------------------
 107 /*Affichage*/
 108 std::ostream& operator<<(std::ostream& os,const marVector& v)
 109 {
 110         int i;
 111     for(i=0;i < v.size();i++)
 112         {
 113                 os << " " << v(i);
 114         }
 115         os <<"\n";
 116     return(os);
 117 }
 118
 119 // ---------------------------------------------------------------------
 120 /*Casting*/
 121 marVector::operator double*() const
 122 {
 123         return(_data);
 124 }
 125
 126 // ---------------------------------------------------------------------
 127 /*Indexation*/
 128 double& marVector::operator()(size_t i)
 129 {
 130         return(_data[i]);
 131 }
 132
 133 const double& marVector::operator()(size_t i) const
 134 {
 135         return(_data[i]);
 136 }
 137
 138 // ---------------------------------------------------------------------
 139 /*Affectation*/
 140 marVector& marVector::operator=(const marVector& o)
 141 {
 142         _size=o._size;
 143         if (!(_data==NULL)&&(!shallowCopy))
 144         {
 145                 delete _data;
 146         }
 147         shallowCopy=false;
 148         _data=new double[_size];
 149         int i;
 150         for (i=0;i<_size;i++)
 151         {
 152                 _data[i]=o._data[i];
 153         }
 154         return (*this);
 155 }
 156
 157 marVector& marVector::operator=(double o)
 158 {
 159         if (_size!=0)
 160         {
 161                 if (!(_data==NULL)&&(!shallowCopy))
 162                 {
 163                         delete _data;
 164                 }
 165                 shallowCopy=false;
 166                 _data=new double[_size];
 167                 int i;
 168                 for (i=0;i<_size;i++)
 169                 {
 170                         _data[i]=o;
 171                 }
 172         }
 173         else
 174         {
 175                 assert(false);
 176         }
 177         return (*this);
 178 }
 179
 180 marVector& marVector::operator=(double*o)
 181 {
 182         if (_size!=0)
 183         {
 184                 if (!(_data==NULL)&&(!shallowCopy))
 185                 {
 186                         delete _data;
 187                 }
 188                 shallowCopy=false;
 189                 _data=new double[_size];
 190                 int i;
 191                 for (i=0;i<_size;i++)
 192                 {
 193                         _data[i]=o[i];
 194                 }
 195         }
 196         else
 197         {
 198                 assert(false);
 199         }
 200         return (*this);
 201 }
 202
 203 // ---------------------------------------------------------------------
 204 /*Comparaison*/
 205 bool marVector::operator==(const marVector& o) const
 206 {
 207     bool equal=true;
 208
 209     if(_size!=o._size)
 210         {
 211                 return(false);
 212         }
 213
 214         int i;
 215         for(i=0;i<_size && equal;i++)
 216         {
 217                 equal=equal && (_data[i]==o._data[i]);
 218         }
 219     return(equal);
 220 }
 221
 222 bool marVector::operator!=(const marVector& o) const
 223 {
 224         return(!((*this)==o));
 225 }
 226
 227 // ---------------------------------------------------------------------
 228 /*Addition*/
 229 marVector marVector::operator+(const marVector& o)
 230 {
 231         marVector result(*this);
 232         size_t s=result._size;
 233
 234         if ((o._size!=s) || (s==0))
 235         {
 236                 assert(false);
 237         }
 238         else
 239         {
 240                 for (int i=0;i < s;i++)
 241                 {
 242                         result._data[i] +=o._data[i];
 243                 }
 244         }
 245         return result;
 246 }
 247
 248 marVector marVector::operator+(double o)
 249 {
 250         marVector result(*this);
 251         size_t s=result._size;
 252
 253         if (s==0)
 254         {
 255                 assert(false);
 256         }
 257         else
 258         {
 259                 for (int i=0;i < s;i++)
 260                 {
 261                         result._data[i] +=o;
 262                 }
 263         }
 264         return result;
 265 }
 266
 267 marVector marVector::operator+(double*o)
 268 {
 269         marVector result(*this);
 270         size_t s=result._size;
 271
 272         if (s==0)
 273         {
 274                 assert(false);
 275         }
 276         else
 277         {
 278                 for (int i=0;i < s;i++)
 279                 {
 280                         result._data[i] +=o[i];
 281                 }
 282         }
 283         return result;
 284 }
 285
 286 // ---------------------------------------------------------------------
 287 /*Addition + Affectation*/
 288 marVector& marVector::operator+=(const marVector& o)
 289 {
 290         if ((o._size!=_size) || (_size==0))
 291         {
 292                 assert(false);
 293         }
 294         else
 295         {
 296                 for (int i=0;i < _size;i++)
 297                 {
 298                         _data[i] +=o._data[i];
 299                 }
 300         }
 301         return (*this);
 302 }
 303
 304 marVector& marVector::operator+=(double o)
 305 {
 306         if (_size==0)
 307         {
 308                 assert(false);
 309         }
 310         else
 311         {
 312                 for (int i=0;i < _size;i++)
 313                 {
 314                         _data[i] +=o;
 315                 }
 316         }
 317     return(*this);
 318 }
 319
 320 marVector& marVector::operator+=(double*o)
 321 {
 322         if (_size==0)
 323         {
 324                 assert(false);
 325         }
 326         else
 327         {
 328                 for (int i=0;i < _size;i++)
 329                 {
 330                         _data[i] +=o[i];
 331                 }
 332         }
 333     return(*this);
 334 }
 335
 336 // ---------------------------------------------------------------------
 337 /*Soustraction*/
 338 marVector marVector::operator-(const marVector& o)
 339 {
 340         marVector result(*this);
 341         size_t s=result._size;
 342
 343         if ((o._size !=s) || (s==0))
 344         {
 345                 assert(false);
 346         }
 347         else
 348         {
 349                 for (int i=0;i < s;i++)
 350                 {
 351                         result._data[i] -=o._data[i];
 352                 }
 353         }
 354         return result;
 355 }
 356
 357 marVector marVector::operator-(double o)
 358 {
 359         marVector result(*this);
 360         size_t s=result._size;
 361
 362         if (s==0)
 363         {
 364                 assert(false);
 365         }
 366         else
 367         {
 368                 for (int i=0;i < s;i++)
 369                 {
 370                         result._data[i] -=o;
 371                 }
 372         }
 373         return result;
 374 }
 375
 376 marVector marVector::operator-(double*o)
 377 {
 378         marVector result(*this);
 379         size_t s=result._size;
 380
 381         if (s==0)
 382         {
 383                 assert(false);
 384         }
 385         else
 386         {
 387                 for (int i=0;i < s;i++)
 388                 {
 389                         result._data[i] -=o[i];
 390                 }
 391         }
 392         return result;
 393 }
 394
 395 // ---------------------------------------------------------------------
 396 /*Sooustraction + Affection*/
 397 marVector& marVector::operator-=(const marVector& o)
 398 {
 399         if ((o._size!=_size) || (_size==0))
 400         {
 401                 assert(false);
 402         }
 403         else
 404         {
 405                 for (int i=0;i < _size;i++)
 406                 {
 407                         _data[i] -=o._data[i];
 408                 }
 409         }
 410         return (*this);
 411 }
 412
 413 marVector& marVector::operator-=(double o)
 414 {
 415         if (_size==0)
 416         {
 417                 assert(false);
 418         }
 419         else
 420         {
 421                 for (int i=0;i < _size;i++)
 422                 {
 423                         _data[i] -=o;
 424                 }
 425         }
 426     return(*this);
 427 }
 428
 429 marVector& marVector::operator-=(double*o)
 430 {
 431         if (_size==0)
 432         {
 433                 assert(false);
 434         }
 435         else
 436         {
 437                 for (int i=0;i < _size;i++)
 438                 {
 439                         _data[i] -=o[i];
 440                 }
 441         }
 442     return(*this);
 443 }
 444
 445 // ---------------------------------------------------------------------
 446 /*Multiplication (produit scalaire)*/
 447 marVector marVector::operator*(double o)
 448 {
 449         marVector result(*this);
 450         size_t s=result._size;
 451
 452         if (s==0)
 453         {
 454                 assert(false);
 455         }
 456         else
 457         {
 458                 for (int i=0;i < s;i++)
 459                 {
 460                         result._data[i]*=o;
 461                 }
 462         }
 463         return result;
 464 }
 465
 466 // ---------------------------------------------------------------------
 467 /*Multiplication (produit scalaire) + Affectation*/
 468 marVector& marVector::operator*=(double o)
 469 {
 470         if (_size==0)
 471         {
 472                 assert(false);
 473         }
 474         else
 475         {
 476                 for (int i=0;i < _size;i++)
 477                 {
 478                         _data[i]*=o;
 479                 }
 480         }
 481     return(*this);
 482 }
 483
 484 // ---------------------------------------------------------------------
 485 /*Division (division scalaire)*/
 486 marVector marVector::operator/(double o)
 487 {
 488         marVector result(*this);
 489         size_t s=result._size;
 490
 491         if (s==0)
 492         {
 493                 assert(false);
 494         }
 495         else
 496         {
 497                 for (int i=0;i < s;i++)
 498                 {
 499                         result._data[i] /=o;
 500                 }
 501         }
 502         return result;
 503 }
 504
 505 // ---------------------------------------------------------------------
 506 /*Division (division scalaire) + Affectation*/
 507 marVector& marVector::operator/=(double o)
 508 {
 509         if (_size==0)
 510         {
 511                 assert(false);
 512         }
 513         else
 514         {
 515                 for (int i=0;i < _size;i++)
 516                 {
 517                         _data[i] /=o;
 518                 }
 519         }
 520     return(*this);
 521 }
 522
 523 //////////////////////////////////////////////////////////////////////
 524 // Op�rations
 525 //////////////////////////////////////////////////////////////////////
 526 // ---------------------------------------------------------------------
 527 /*Produit scalaire*/
 528 double marVector::dot(const marVector& o)
 529 {
 530     double result;
 531         if ((_size!=o._size) || (_size==0))
 532         {
 533                 assert(false);
 534         }
 535         else
 536         {
 537                 result=0.0;
 538                 for (int i=0;i<_size;i++)
 539                 {
 540                         result+=(_data[i]*o._data[i]);
 541                 }
 542         }
 543         return result;
 544 }
 545
 546 double marVector::dot(double*o)
 547 {
 548     double result;
 549         if (_size==0)
 550         {
 551                 assert(false);
 552         }
 553         else
 554         {
 555                 result=0.0;
 556                 for (int i=0;i<_size;i++)
 557                 {
 558                         result+=(_data[i]*o[i]);
 559                 }
 560         }
 561         return result;
 562 }
 563
 564 // ---------------------------------------------------------------------
 565 /*Produit vectoriel*/
 566 marVector marVector::cross(const marVector& o)
 567 {
 568     marVector result(*this);
 569         if ((o._size!=_size) || (_size!=3))
 570         {
 571                 assert(false);
 572         }
 573         else
 574         {
 575                 double S,s;
 576
 577                 // 1st element
 578                 S=_data[1]*o._data[2];
 579                 s=_data[2]*o._data[1];
 580                 result._data[0]=S - s;
 581
 582                 // 2nd element
 583                 S=_data[2]*o._data[0];
 584                 s=_data[0]*o._data[2];
 585                 result._data[1]=S - s;
 586
 587                 // 3rd element
 588                 S=_data[0]*o._data[1];
 589                 s=_data[1]*o._data[0];
 590                 result._data[2]=S - s;
 591         }
 592     return(result);
 593 }
 594
 595 marVector marVector::cross(double*o)
 596 {
 597     marVector result(*this);
 598         if (_size!=3)
 599         {
 600                 assert(false);
 601         }
 602         else
 603         {
 604                 double S,s;
 605
 606                 // 1st element
 607                 S=_data[1]*o[2];
 608                 s=_data[2]*o[1];
 609                 result._data[0]=S - s;
 610
 611                 // 2nd element
 612                 S=_data[2]*o[0];
 613                 s=_data[0]*o[2];
 614                 result._data[1]=S - s;
 615
 616                 // 3rd element
 617                 S=_data[0]*o[1];
 618                 s=_data[1]*o[0];
 619                 result._data[2]=S - s;
 620         }
 621     return(result);
 622 }
 623
 624 int marVector::scross(const marVector& o)
 625 {
 626         if ((o._size!=_size) && (_size!=3))
 627         {
 628                 return (1);
 629         }
 630         else
 631         {
 632                 *this=cross(o);
 633                 return(0);
 634         }
 635 }
 636
 637 int marVector::scross(double*o)
 638 {
 639         if (_size!=3)
 640         {
 641                 return (1);
 642         }
 643         else
 644         {
 645                 *this=cross(o);
 646                 return(0);
 647         }
 648 }
 649
 650 // ---------------------------------------------------------------------
 651 /*Norme euclidienne*/
 652 double marVector::norm2()
 653 // Calcule la norme euclidienne du vecteur
 654 // en utilisant l'agorithme de Blue
 655 {
 656
 657         double norme2=0.0;
 658         if (_size==0)
 659         {
 660                 assert(false);
 661         }
 662         else
 663         {
 664                 const double seuil1=pow( 2., 154. ),seuil2=1/seuil1;
 665                 double n1=0; /* les sommes partielles */
 666                 double n2=0;
 667                 double n3=0;
 668                 double n4;
 669                 int i;
 670                 for ( i=0; i<_size; i++ )
 671                 {
 672                         double x=_data[i];      /* x=abs(vecteur[i]) */
 673
 674                         if (x<0)
 675                         {
 676                                 x=-x;
 677                         }
 678                         /* sommation par classe */
 679                         if ( x > seuil1 )
 680                         {
 681                                 x *=seuil2;
 682                                 n1 +=x*x;
 683                         }
 684                         else if ( x < seuil2 )
 685                         {
 686                                 x  *=seuil1;
 687                                 n2 +=x*x;
 688                         }
 689                         else
 690                         {
 691                                 n3 +=x*x;
 692                         }
 693                 }
 694
 695                 n3=sqrt(n3);
 696                 if (n1>0)
 697                 {
 698                         n4=seuil1 * sqrt(n1);
 699                 }
 700                 else
 701                 {
 702                         n4=seuil2 * sqrt(n2);
 703                 }
 704
 705 // EED
 706 /*
 707                 if (n3 < n4)
 708                 {
 709                         n3 /=n4;
 710                         norme2=n4 * sqrt( 1 + n3*n3 );
 711                 }
 712                 else
 713                 {
 714                         n4 /=n3;
 715                         norme2=n3 * sqrt( 1 + n4*n4 );
 716                 }
 717 */
 718                 norme2=n3;
 719         }
 720         return norme2;
 721 }
 722
 723 // ---------------------------------------------------------------------
 724 /*Normalisation*/
 725 marVector marVector::normalize( )
 726 {
 727         marVector result=marVector(*this);
 728         result *= (1/result.norm2());
 729         return result;
 730 }
 731
 732 int marVector::snormalize( )
 733 {
 734         (*this)*=(1/this->norm2());
 735         return(0);
 736 }// renvoie 0 si OK, 1 sinon
 737
 738 //////////////////////////////////////////////////////////////////////
 739 // M�thodes
 740 //////////////////////////////////////////////////////////////////////
 741
 742 size_t marVector::size() const
 743 {
 744         return _size;
 745 }