lib/maracasVisuLib/src/interface/wxWindows/widgets/wxSphereView.cxx

   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 #include "wxSphereView.h"
  27
  28 wxSphereView::wxSphereView( wxWindow *parent, vtkMPRBaseData *vtkmprbasedata/*, vtkImageData *imageData */)
  29 : wxVtk2DBaseView(parent)
  30 {
  31         _delta                          =       1;
  32         _vtkmprbasedata         =       vtkmprbasedata;
  33
  34 //EED 2016-08-31
  35 //      _imageDataOriginal      =       imageData;
  36
  37         SetImage();
  38         _imageSphere            =       vtkImageData::New();
  39         _imageSphere->SetDimensions (150,150,500);
  40
  41 //EED
  42 //      _imageSphere->SetScalarTypeToUnsignedShort();
  43         _imageSphere->SetScalarType( _imageDataOriginal->GetScalarType() );
  44
  45         _imageSphere->AllocateScalars();
  46         _imageSphere->Update();
  47
  48 //EED  ???? vtkBaseData  no esta compartido con los otros objetos .. PLOP
  49 //      vtkBaseData *vtkbasedata = new vtkBaseData();
  50 //      vtkbasedata->SetMarImageData( new marImageData(_imageSphere) );
  51 //      this->SetVtkBaseData(vtkbasedata);
  52         this->SetVtkBaseData(_vtkmprbasedata);
  53
  54     _transform                  =       vtkTransform::New();
  55     _transform1                 =       vtkTransform::New();
  56     _transform2                 =       vtkTransform::New();
  57         _transform ->Identity();
  58         _transform1->Identity();
  59         _transform2->Identity();
  60         _radio=25;
  61 }
  62
  63 //-------------------------------------------------------------------
  64 wxSphereView::~wxSphereView()
  65 {
  66         _transform  -> Delete();
  67         _transform1 -> Delete();
  68         _transform2 -> Delete();
  69         ResetlstId();
  70 }
  71
  72 //----------------------------------------------------------------------------
  73 void wxSphereView::SetImage()
  74 {
  75         ResetlstId();
  76         _imageDataOriginal=_vtkmprbasedata->GetImageData();
  77 }
  78
  79 //----------------------------------------------------------------------------
  80 double wxSphereView::GetRadio()
  81 {
  82         return _radio;
  83 }
  84
  85 //----------------------------------------------------------------------------
  86 void wxSphereView::SetRadio(double radio)
  87 {
  88         if (radio<0)
  89         {
  90                 radio=0;
  91         } // if
  92         _radio=radio;
  93 }
  94
  95 //----------------------------------------------------------------------------
  96 void wxSphereView::Configure()
  97 {
  98         wxVtk2DBaseView::Configure();
  99
 100
 101         vtkInteractorStyleBaseView2D *style2D = vtkInteractorStyleBaseView2D::New();
 102         manualInteractorWindowLevel *_manualinteractorwindowlevel= new manualInteractorWindowLevel();
 103         style2D->SetInteractorWindowLevel( _manualinteractorwindowlevel );
 104 //      vtkInteractorScrollZ *_vtkInteractorScrollZ = new vtkInteractorScrollZ();
 105 //      style2D->SetInteractorScrollZ(_vtkInteractorScrollZ);
 106         SetInteractorStyleImage( style2D );
 107
 108
 109
 110         _vtkinteractorstylesphere = new vtkInteractorStyleSphere();
 111         ((vtkInteractorStyleBaseView*)GetInteractorStyleBaseView())->AddInteractorStyleMaracas( _vtkinteractorstylesphere );
 112         double points[4][3];
 113 // EED purify 12/sep/2006
 114         int i,j;
 115         for (i=0;i<4;i++)
 116         {
 117                 for (j=0;j<3;j++)
 118                 {
 119                         points[i][j]=0;
 120                 } // for
 121         } // for
 122         InitSphere(points);
 123
 124
 125         DefineImageSphere();
 126 }
 127
 128 //----------------------------------------------------------------------------
 129 void wxSphereView::RefreshPoint()
 130 {
 131         double x        = _vtkmprbasedata->GetX() - _centerX;
 132         double y        = _vtkmprbasedata->GetY() - _centerY;
 133         double z        = _vtkmprbasedata->GetZ() - _centerZ;
 134         double alpha= atan2(x,z);
 135         double beta = atan2( y , sqrt(z*z+x*x) );
 136         alpha           = alpha*180/3.1416;
 137         beta            = beta*180/3.1416;
 138         _transform1->Identity();
 139         _transform1->RotateY(alpha);
 140         _transform1->RotateX(-beta);
 141         _radio= sqrt(x*x + y*y +z*z);
 142         RefreshView();
 143 }
 144
 145
 146
 147 //-------------------------------------------------------------------
 148 void wxSphereView::Refresh(  )
 149 {
 150 //      ExtractPlane();
 151         UpdateColorWindowLevel();
 152         wxVtkBaseView::Refresh();
 153 }
 154
 155
 156
 157 //----------------------------------------------------------------------------
 158 void wxSphereView::RefreshView()
 159 {
 160 //EED
 161         _imageViewer2XYZ->GetVtkImageViewer2()->SetInput ( _imageSphere );
 162         DefineImageSphere();
 163 //      UpdateColorWindowLevel();
 164 //      wxVtk2DBaseView::Refresh();
 165         Refresh();
 166 }
 167
 168 //----------------------------------------------------------------------------
 169 void wxSphereView::RotationEnd()
 170 {
 171         _transform1->RotateWXYZ(_ang,_vxb,_vyb,0);
 172         _transform2->Identity();
 173         SetDeltaVoxel(1);
 174         ResetlstId();
 175 }
 176
 177 //----------------------------------------------------------------------------
 178 void wxSphereView::RotationStart(double vx, double vy, bool ok_v, bool ok_ang)
 179 {
 180         if (ok_ang==false)
 181         {
 182                 _ang = -sqrt( vx*vx + vy*vy ) / 1.0;
 183         } // if ok_ang
 184         if (ok_v==false){
 185                 _vxb=-vy;
 186                 _vyb=vx;
 187         } // if ok_v
 188         _transform2->Identity();
 189         _transform2->RotateWXYZ(_ang,_vxb,_vyb,0);
 190         SetDeltaVoxel(3);
 191         ResetlstId();
 192 }
 193
 194 //----------------------------------------------------------------------------
 195 void wxSphereView::GetPointSphere(double p[3],double r1,double angA,double angB)
 196 {
 197         double in[3],out[3];
 198         in[0]=0;
 199         in[1]=r1;
 200         in[2]=0;
 201         vtkTransform *transform = vtkTransform::New();
 202         transform->Identity();
 203         transform->RotateX(angB);
 204         transform->RotateZ(angA);
 205         transform->TransformPoint(in,out);
 206         p[0]=out[0];
 207         p[1]=out[1];
 208         p[2]=out[2];
 209         transform->Delete();
 210 }
 211
 212 //----------------------------------------------------------------------------
 213 void wxSphereView::RotatePointOverTheSphere( double pp[3], double p[3],double cc[3])
 214 {
 215         double out[3];
 216         _transform->TransformPoint(p,out);
 217         pp[0] = out[0] + cc[0];
 218         pp[1] = out[1] + cc[1];
 219         pp[2] = out[2] + cc[2];
 220 }
 221
 222 //----------------------------------------------------------------------------
 223 void wxSphereView::ResetlstId()
 224 {
 225         int i,size=_lstId.size();
 226         for (i=size-1;i>=0;i--)
 227         {
 228                 delete _lstId[i];
 229         }
 230         _lstId.clear();
 231 }
 232
 233 //----------------------------------------------------------------------------
 234 int wxSphereView::GetIdOfImage(double radio)
 235 {
 236         int id=0;
 237         int dim[3];
 238         _imageSphere->GetDimensions(dim);
 239         int sizeMaxList = dim[2];
 240         // Search in list >> alpha beta radio
 241         int i,size=_lstId.size();
 242         for (i=0; i<size;i++)
 243         {
 244                 //idAlBeRa *tmp=_lstId[i]; // JPRx
 245                 if ((_lstId[i]->_radio==radio) && (_lstId[i]->_deltavoxel==_delta))
 246                 {
 247                         return _lstId[i]->_id;
 248                 }
 249         }
 250         if (size>sizeMaxList)
 251         {
 252                 delete _lstId[size-1];
 253                 _lstId.pop_back();
 254         }
 255         if (size!=0){
 256                 id=_lstId[0]->_id+1;
 257                 id = id % sizeMaxList;
 258         } else {
 259                 id = 0;
 260         }
 261         FiltreImage(id,radio);
 262         _lstId.insert(_lstId.begin(),1,new idAlBeRa(id,radio,_delta) );
 263         return id;
 264 }
 265
 266 //----------------------------------------------------------------------------
 267 void wxSphereView::DefineImageSphere()
 268 {
 269         int id;
 270         id=GetIdOfImage( _radio );
 271 //EED
 272 //      GetVtkBaseData()->SetZ( id );
 273
 274         _imageViewer2XYZ->GetVtkImageViewer2()->SetSlice ( id );
 275 }
 276
 277 //----------------------------------------------------------------------------
 278 void wxSphereView::SetDeltaVoxel(int delta)
 279 {
 280         _delta=delta;
 281 }
 282
 283 //----------------------------------------------------------------------------
 284 void wxSphereView::SetVoxel(double i, double j, int delta,double id,  double gris)
 285 {
 286         int ii,jj,delta2;
 287         unsigned short *pRes;
 288         int dimRes[3];
 289         _imageSphere->GetDimensions(dimRes);
 290
 291         delta2=delta-1;
 292         for ( ii=(int)(i-delta2) ; ii<=(int)(i+delta2) ; ii++ )
 293         {
 294                 for ( jj=(int)(j-delta2) ; jj<=(int)(j+delta2) ; jj++ )
 295                 {
 296                         if ( (ii>=0)&&(ii<dimRes[0]) &&
 297                                  (jj>=0)&&(jj<dimRes[1]) )
 298                         {
 299 //EED
 300 //                              pRes = (unsigned short*)_imageSphere->GetScalarPointer( ii , jj , (int)id );
 301 //                              *pRes=gris;
 302                                 _imageSphere->SetScalarComponentFromDouble( ii , jj , (int)id , 0, gris );
 303                         }
 304                 } // for jj
 305         } // for ii
 306 }
 307
 308 //----------------------------------------------------------------------------
 309 void wxSphereView::SetXYZtoParent(double i, double j)
 310 {
 311
 312         double factor = 0.75;
 313         double radio2   = _radio*_radio;
 314         double pxx,pyy,d2x,d2y;
 315         double cc[3],p[3],pp[3];
 316         cc[0]=_centerX;
 317         cc[1]=_centerY;
 318         cc[2]=_centerZ;
 319         double aa;
 320         int dimRes[3],dimOrig[3];
 321         _imageSphere->GetDimensions(dimRes);
 322         d2x=dimRes[0]/2;
 323         d2y=dimRes[1]/2;
 324         _imageDataOriginal->GetDimensions(dimOrig);
 325         p[0]  = (i - d2x)*factor;
 326         pxx=p[0]*p[0];
 327         p[1]  = (j - d2y)*factor;
 328         pyy=p[1]*p[1];
 329         aa = pxx + pyy;
 330         if (radio2>aa){
 331                 aa=radio2-aa;
 332                 p[2]  = sqrt(aa);
 333                 RotatePointOverTheSphere( pp, p,cc);
 334                 if ( (pp[0]>=0) && (pp[0]<dimOrig[0]) &&
 335                          (pp[1]>=0) && (pp[1]<dimOrig[1]) &&
 336                          (pp[2]>=0) && (pp[2]<dimOrig[2]) )
 337                 {
 338                         if (_vtkmprbasedata){
 339                                 _vtkmprbasedata->SetX(pp[0]);
 340                                 _vtkmprbasedata->SetY(pp[1]);
 341                                 _vtkmprbasedata->SetZ(pp[2]);
 342                         }
 343                 } // if pp
 344         } // if radio
 345 }
 346
 347 //----------------------------------------------------------------------------
 348 void wxSphereView::FiltreImageB(int id, double radio, bool ok,int deltaTMP)
 349 {
 350         double value;
 351         double factor = 0.75;
 352         double radioB   = radio/3;
 353         double radio2   = radio*radio;
 354         double pxx,pyy,d2x,d2y;
 355         double cc[3],p[3],pp[3];
 356         cc[0]=_centerX;
 357         cc[1]=_centerY;
 358         cc[2]=_centerZ;
 359         double aa;
 360         unsigned short *pOrig;
 361         int dimRes[3],dimOrig[3];
 362         double i,j;
 363         int ext[6];
 364         _imageSphere->GetExtent(ext);
 365         _imageSphere->GetDimensions(dimRes);
 366         //JCP 24 - 04 -09
 367         //_imageSphere->SetExtent(0,dimRes[0]-1,0,dimRes[1]-1,0,dimRes[2]-1);
 368         _imageSphere->SetExtent(ext);
 369         //JCP 24 - 04 -09
 370         d2x=dimRes[0]/2;
 371         d2y=dimRes[1]/2;
 372         _imageDataOriginal->GetDimensions(dimOrig);
 373         int start,end;
 374         int limitA,limitB;
 375         limitA  = (int) ( (-radioB/factor)+d2x );
 376         limitB  = (int) ( (radioB/factor)+d2x );
 377         if (ok==true){
 378                 start   = limitA;
 379                 end             = limitB;
 380         } else {
 381                 start=0;
 382                 end=dimRes[0];
 383         }
 384
 385         for ( i=start ; i<end ; i=i+deltaTMP )
 386         {
 387                 p[0]  = (i - d2x)*factor;
 388                 pxx=p[0]*p[0];
 389                 for (j=start;j<end;j=j+deltaTMP)
 390                 {
 391                         p[1]  = (j - d2y)*factor;
 392                         pyy=p[1]*p[1];
 393                         aa = pxx + pyy;
 394
 395                         if  (( ((ok==false) && (!((i>limitA) && (i<limitB) && (j>limitA) && (j<limitB)))) )
 396                                     ||
 397                                         (ok==true))
 398                         {
 399                                 if (radio2>aa){
 400                                         aa=radio2-aa;
 401                                         p[2]  = sqrt(aa);
 402                                         RotatePointOverTheSphere( pp, p,cc);
 403                                         if ( (pp[0]>=0) && (pp[0]<dimOrig[0]) &&
 404                                                  (pp[1]>=0) && (pp[1]<dimOrig[1]) &&
 405                                                  (pp[2]>=0) && (pp[2]<dimOrig[2]) )
 406                                         {
 407 //EED
 408 //                                              pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( (int)(pp[0]) , (int)(pp[1]) , (int)(pp[2]) );
 409 //                                              SetVoxel(i,j,deltaTMP,id,*pOrig);
 410                                                 value=_imageDataOriginal->GetScalarComponentAsDouble( (int)(pp[0]) , (int)(pp[1]) , (int)(pp[2]),  0);
 411                                                 SetVoxel(i,j,deltaTMP,id,value);
 412                                         } else {
 413                                                 SetVoxel(i,j,deltaTMP,id,2000);
 414                                         }
 415                                 } else {
 416                                         SetVoxel(i,j,deltaTMP,id,0);
 417                                 }
 418                         }
 419                 }
 420         }
 421         _imageSphere->Modified();
 422         _imageSphere->Update();
 423 }
 424
 425 //----------------------------------------------------------------------------
 426 void wxSphereView::FiltreImage(int id, double radio)
 427 {
 428         _transform -> Identity();
 429         _transform -> Concatenate(_transform1);
 430         _transform -> Concatenate(_transform2);
 431         FiltreImageB(id,radio,false, _delta);
 432         FiltreImageB(id,radio,true, 1);
 433 }
 434
 435 //----------------------------------------------------------------------------
 436
 437 /*
 438 void wxSphereView::FiltreImage(int id, double radio)
 439 {
 440         double radio2   = radio*radio;
 441         double radio2TMP= (radio/2)*(radio/2);
 442         double pxx,pyy,d2x,d2y;
 443         double cc[3],p[3],pp[3];
 444         cc[0]=_centerX;
 445         cc[1]=_centerY;
 446         cc[2]=_centerZ;
 447         double aa;
 448         unsigned short *pOrig;
 449         int dimRes[3],dimOrig[3];
 450         double i,j;
 451         _imageSphere->GetDimensions(dimRes);
 452         _imageSphere->SetExtent(0,dimRes[0]-1,0,dimRes[1]-1,0,dimRes[2]-1);
 453         d2x=dimRes[0]/2;
 454         d2y=dimRes[1]/2;
 455         double deltaTMP=_delta;
 456         _imageDataOriginal->GetDimensions(dimOrig);
 457
 458         for ( i=0 ; i<dimRes[0] ; i=i+deltaTMP )
 459         {
 460                 p[0]  = (i - d2x)*0.75;
 461                 pxx=p[0]*p[0];
 462                 for (j=0;j<dimRes[1];j=j+deltaTMP)
 463                 {
 464                         p[1]  = (j - d2y)*0.75;
 465                         pyy=p[1]*p[1];
 466                         aa = pxx + pyy;
 467
 468                         if (aa>radio2TMP)
 469                         {
 470                                 if (radio2>aa){
 471                                         aa=radio2-aa;
 472                                         p[2]  = sqrt(aa);
 473                                         RotatePointOverTheSphere( pp, p,cc);
 474                                         if ( (pp[0]>=0) && (pp[0]<dimOrig[0]) &&
 475                                                  (pp[1]>=0) && (pp[1]<dimOrig[1]) &&
 476                                                  (pp[2]>=0) && (pp[2]<dimOrig[2]) )
 477                                         {
 478                                                 pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( pp[0] , pp[1] , pp[2] );
 479                                                 SetVoxel(i,j,deltaTMP,id,*pOrig);
 480                                         } else {
 481                                                 SetVoxel(i,j,deltaTMP,id,2000);
 482                                         }
 483                                 } else {
 484                                         SetVoxel(i,j,deltaTMP,id,0);
 485                                 }
 486                         }
 487                 }
 488         }
 489
 490
 491         deltaTMP=1;
 492         for ( i=0 ; i<dimRes[0] ; i=i+deltaTMP )
 493         {
 494                 p[0]  = (i - d2x)*0.75;
 495                 pxx=p[0]*p[0];
 496                 for (j=0;j<dimRes[1];j=j+deltaTMP)
 497                 {
 498                         p[1]  = (j - d2y)*0.75;
 499                         pyy=p[1]*p[1];
 500                         aa = pxx + pyy;
 501                         if (aa<=radio2TMP)
 502                         {
 503                                 if (radio2>aa){
 504                                         aa=radio2-aa;
 505                                         p[2]  = sqrt(aa);
 506                                         RotatePointOverTheSphere( pp, p,cc);
 507                                         if ( (pp[0]>=0) && (pp[0]<dimOrig[0]) &&
 508                                                  (pp[1]>=0) && (pp[1]<dimOrig[1]) &&
 509                                                  (pp[2]>=0) && (pp[2]<dimOrig[2]) )
 510                                         {
 511                                                 pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( pp[0] , pp[1] , pp[2] );
 512                                                 SetVoxel(i,j,deltaTMP,id,*pOrig);
 513                                         } else {
 514                                                 SetVoxel(i,j,deltaTMP,id,2000);
 515                                         }
 516                                 } else {
 517                                         SetVoxel(i,j,deltaTMP,id,0);
 518                                 }
 519                         }
 520                 }
 521         }
 522
 523         _imageSphere->Modified();
 524         _imageSphere->Update();
 525 }
 526 */
 527 /*
 528 void wxSphereView::FiltreImage(vtkImageData *imageSphere)
 529 {
 530         int dim[3],i,j,k;
 531         imageSphere->GetDimensions(dim);
 532         for (i=0;i<dim[0];i++)
 533         {
 534                 for (j=0;j<dim[1];j++)
 535                 {
 536                         for (k=0;k<dim[2];k++)
 537                         {
 538                                 unsigned short *pRes=(unsigned short*)imageSphere->GetScalarPointer (i,j,k);
 539                                 *pRes=0;
 540                         }
 541                 }
 542         }
 543
 544         double deltaA=90;
 545         double cc[3],p1[3],p2[3],pp1[3],pp2[3];
 546         cc[0]=_centerX;
 547         cc[1]=_centerY;
 548         cc[2]=_centerZ;
 549         double r1       = _sl_radio->GetValue() - _sl_thickness->GetValue()/2;
 550         double r2       = _sl_radio->GetValue() + _sl_thickness->GetValue()/2;
 551         if (r1<10)
 552         {
 553                 r1=10;
 554         }
 555         double alpha= _sl_alpha->GetValue();
 556         double beta     = _sl_beta->GetValue();
 557
 558         double angA,angB;
 559         for (angA=-deltaA;angA<deltaA;angA++)
 560         {
 561                 for (angB=-deltaA;angB<deltaA;angB++)
 562                 {
 563                         GetPointSphere(p1,r1,angA,angB);
 564                         GetPointSphere(p2,r2,angA,angB);
 565                         RotatePointOverTheSphere( pp1, alpha, beta, p1 ,cc );
 566                         RotatePointOverTheSphere( pp2, alpha, beta, p2 ,cc );
 567                         TransferePoints(pp1,pp2,angA+alpha+180,angB+beta+90,imageSphere);
 568                 }
 569         }
 570 }
 571
 572
 573 //----------------------------------------------------------------------------
 574
 575 void wxSphereView::TransferePoints(double pp1[3],double pp2[3],double AngX,double AngY,vtkImageData *image)
 576 {
 577         double t;
 578         double difX = pp2[0]-pp1[0];
 579         double difY = pp2[1]-pp1[1];
 580         double difZ = pp2[2]-pp1[2];
 581         double  max             = 200;
 582         int dimOrg[3];
 583         int dimRes[3];
 584         int z;
 585         _imageDataOriginal->GetDimensions(dimOrg);
 586         image->GetDimensions(dimRes);
 587         int i;
 588         double x1=pp1[0];
 589         double y1=pp1[1];
 590         double z1=pp1[2];
 591         int xx=-1,yy=-1,zz=-1;
 592         for (i=0;i<max;i++)
 593         {
 594                 t  = i/max;
 595                 xx = (int) (x1+t*difX);
 596                 yy = (int) (y1+t*difY);
 597                 zz = (int) (z1+t*difZ);
 598
 599                 z=i;
 600                 if ((xx>=0) && (xx<dimOrg[0]) && (yy>=0) && (yy<dimOrg[1]) && (zz>=0) && (zz<dimOrg[2]) &&
 601                         (AngX>=0) && (AngX<dimRes[0]) && (AngY>=0) && (AngY<dimRes[1]) && (z>=0) && (z<dimRes[2]) )
 602                 {
 603                         unsigned short *pOrg=(unsigned short*)_imageDataOriginal->GetScalarPointer (xx,yy,zz);
 604                         unsigned short *pRes=(unsigned short*)image->GetScalarPointer( (int)AngX , (int)AngY , z );
 605                         *pRes=*pOrg;
 606                 }
 607         }
 608 }
 609
 610
 611
 612
 613 */
 614
 615
 616 //----------------------------------------------------------------------------
 617
 618 void wxSphereView::InitSphere(double points[4][3])
 619 {
 620         double cc[3];
 621     double r = SphereFindCenter(points,cc); // 4-points , center
 622     if (r > 0)
 623     {
 624         _centerX        = (int)(cc[0]);
 625         _centerY        = (int)(cc[1]);
 626         _centerZ        = (int)(cc[2]);
 627     } else {
 628                 int dim[3];
 629                 _imageDataOriginal->GetDimensions(dim);
 630         _centerX        = (int)(dim[0]/2);
 631         _centerY        = (int)(dim[1]/2);
 632         _centerZ        = (int)(dim[2]/2);
 633         }
 634 }
 635
 636 //----------------------------------------------------------------------------
 637
 638 // Calculate center and radius of sphere given four points
 639 // http://home.att.net/~srschmitt/script_sphere_solver.html
 640 // source code HTML <script language=JavaScript>
 641 double wxSphereView::SphereFindCenter(double P[4][3], double cc[3])
 642 {
 643     int i;
 644     double r, m11, m12, m13, m14, m15;
 645         double a[4][4];
 646
 647     for (i = 0; i < 4; i++)                    // find minor 11
 648     {
 649         a[i][0] = P[i][0];
 650         a[i][1] = P[i][1];
 651         a[i][2] = P[i][2];
 652         a[i][3] = 1;
 653     }
 654     m11 = determinant( a, 4 );
 655
 656     for (i = 0; i < 4; i++)                    // find minor 12
 657     {
 658         a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
 659         a[i][1] = P[i][1];
 660         a[i][2] = P[i][2];
 661         a[i][3] = 1;
 662     }
 663     m12 = determinant( a, 4 );
 664
 665     for (i = 0; i < 4; i++)                    // find minor 13
 666     {
 667         a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
 668         a[i][1] = P[i][0];
 669         a[i][2] = P[i][2];
 670         a[i][3] = 1;
 671     }
 672     m13 = determinant( a, 4 );
 673
 674     for (i = 0; i < 4; i++)                    // find minor 14
 675     {
 676         a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
 677         a[i][1] = P[i][0];
 678         a[i][2] = P[i][1];
 679         a[i][3] = 1;
 680     }
 681     m14 = determinant( a, 4 );
 682
 683
 684     for (i = 0; i < 4; i++)                    // find minor 15
 685     {
 686         a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
 687         a[i][1] = P[i][0];
 688         a[i][2] = P[i][1];
 689         a[i][3] = P[i][2];
 690     }
 691     m15 = determinant( a, 4 );
 692
 693     if (m11 == 0)
 694     {
 695         r = 0;
 696     }
 697     else
 698     {
 699                 // center of sphere
 700         cc[0] =  0.5*m12/m11;   //cx
 701         cc[1] = -0.5*m13/m11;   //cy
 702         cc[2] =  0.5*m14/m11;   //cz
 703                 // Sphere radio
 704         r  = sqrt( cc[0]*cc[0] + cc[1]*cc[1] + cc[2]*cc[2] - m15/m11 );
 705     }
 706
 707     return r;                                  // the radius
 708 }
 709 //----------------------------------------------------------------------------
 710
 711 //  Recursive definition of determinate using expansion by minors.
 712 double wxSphereView::determinant(double a[4][4], int n)
 713 {
 714     int i, j, j1, j2;
 715     double d;
 716         double m[4][4];
 717
 718         for (i=0;i<4;i++)
 719         {
 720                 for (j=0;j<4;j++)
 721                 {
 722                         m[i][j]=0;
 723                 }
 724         }
 725
 726     if (n == 2)                                // terminate recursion
 727     {
 728         d = a[0][0]*a[1][1] - a[1][0]*a[0][1];
 729     }
 730     else
 731     {
 732         d = 0;
 733         for (j1 = 0; j1 < n; j1++ )            // do each column
 734         {
 735             for (i = 1; i < n; i++)            // create minor
 736             {
 737                 j2 = 0;
 738                 for (j = 0; j < n; j++)
 739                 {
 740                     if (j == j1) continue;
 741                     m[i-1][j2] = a[i][j];
 742                     j2++;
 743                 }
 744             }
 745
 746             // sum (+/-)cofactor * minor
 747             d = d + pow(-1.0, j1)*a[0][j1]*determinant( m, n-1 );
 748         }
 749     }
 750
 751     return d;
 752 }
 753