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