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

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