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