1 /*# ---------------------------------------------------------------------
3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
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
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.
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
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 # ------------------------------------------------------------------------ */
26 #include "wxSphereView.h"
29 wxSphereView::wxSphereView( wxWindow *parent, vtkMPRBaseData *vtkmprbasedata, vtkImageData *imageData )
30 : wxVtk2DBaseView(parent)
33 _vtkmprbasedata = vtkmprbasedata;
34 _imageDataOriginal = imageData;
36 _imageSphere = vtkImageData::New();
37 _imageSphere->SetDimensions (150,150,500);
38 _imageSphere->SetScalarTypeToUnsignedShort();
39 _imageSphere->AllocateScalars();
40 _imageSphere->Update();
43 vtkBaseData *vtkbasedata = new vtkBaseData();
44 vtkbasedata->SetMarImageData( new marImageData(_imageSphere) );
45 this->SetVtkBaseData(vtkbasedata);
47 _transform = vtkTransform::New();
48 _transform1 = vtkTransform::New();
49 _transform2 = vtkTransform::New();
50 _transform ->Identity();
51 _transform1->Identity();
52 _transform2->Identity();
57 //-------------------------------------------------------------------
59 wxSphereView::~wxSphereView()
61 _transform -> Delete();
62 _transform1 -> Delete();
63 _transform2 -> Delete();
67 //----------------------------------------------------------------------------
69 double wxSphereView::GetRadio()
74 //----------------------------------------------------------------------------
76 void wxSphereView::SetRadio(double radio)
85 //----------------------------------------------------------------------------
87 void wxSphereView::Configure()
89 wxVtk2DBaseView::Configure();
91 _vtkinteractorstylesphere = new vtkInteractorStyleSphere();
92 ((vtkInteractorStyleBaseView*)GetInteractorStyleBaseView())->AddInteractorStyleMaracas( _vtkinteractorstylesphere );
95 // EED purify 12/sep/2006
109 //----------------------------------------------------------------------------
111 void wxSphereView::RefreshPoint()
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) );
119 alpha = alpha*180/3.1416;
120 beta = beta*180/3.1416;
122 _transform1->Identity();
123 _transform1->RotateY(alpha);
124 _transform1->RotateX(-beta);
126 _radio= sqrt(x*x + y*y +z*z);
131 //----------------------------------------------------------------------------
133 void wxSphereView::RefreshView()
136 wxVtk2DBaseView::Refresh();
139 //----------------------------------------------------------------------------
141 void wxSphereView::RotationEnd()
143 _transform1->RotateWXYZ(_ang,_vxb,_vyb,0);
144 _transform2->Identity();
149 //----------------------------------------------------------------------------
151 void wxSphereView::RotationStart(double vx, double vy, bool ok_v, bool ok_ang)
155 _ang = -sqrt( vx*vx + vy*vy ) / 1.0;
163 _transform2->Identity();
164 _transform2->RotateWXYZ(_ang,_vxb,_vyb,0);
169 //----------------------------------------------------------------------------
171 void wxSphereView::GetPointSphere(double p[3],double r1,double angA,double angB)
177 vtkTransform *transform = vtkTransform::New();
178 transform->Identity();
179 transform->RotateX(angB);
180 transform->RotateZ(angA);
181 transform->TransformPoint(in,out);
188 //----------------------------------------------------------------------------
190 void wxSphereView::RotatePointOverTheSphere( double pp[3], double p[3],double cc[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];
201 //----------------------------------------------------------------------------
203 void wxSphereView::TransferePoints(double pp1[3],double pp2[3],double AngX,double AngY,vtkImageData *image)
206 double difX = pp2[0]-pp1[0];
207 double difY = pp2[1]-pp1[1];
208 double difZ = pp2[2]-pp1[2];
215 _imageDataOriginal->GetDimensions(dimOrg);
216 image->GetDimensions(dimRes);
222 int xx=-1,yy=-1,zz=-1;
227 xx = (int) (x1+t*difX);
228 yy = (int) (y1+t*difY);
229 zz = (int) (z1+t*difZ);
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]) )
235 unsigned short *pOrg=(unsigned short*)_imageDataOriginal->GetScalarPointer (xx,yy,zz);
236 unsigned short *pRes=(unsigned short*)image->GetScalarPointer( (int)AngX , (int)AngY , z );
242 //----------------------------------------------------------------------------
244 void wxSphereView::ResetlstId()
246 int i,size=_lstId.size();
247 for (i=size-1;i>=0;i--)
254 //----------------------------------------------------------------------------
256 int wxSphereView::GetIdOfImage(double radio)
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++)
266 //idAlBeRa *tmp=_lstId[i]; // JPRx
267 if ((_lstId[i]->_radio==radio) && (_lstId[i]->_deltavoxel==_delta))
269 return _lstId[i]->_id;
272 if (size>sizeMaxList)
274 delete _lstId[size-1];
279 id = id % sizeMaxList;
284 FiltreImage(id,radio);
285 _lstId.insert(_lstId.begin(),1,new idAlBeRa(id,radio,_delta) );
290 //----------------------------------------------------------------------------
292 void wxSphereView::DefineImageSphere()
295 id=GetIdOfImage( _radio );
296 GetVtkBaseData()->SetZ( id );
300 //----------------------------------------------------------------------------
301 void wxSphereView::SetDeltaVoxel(int delta)
306 //----------------------------------------------------------------------------
307 void wxSphereView::SetVoxel(double i, double j, int delta,double id, unsigned short gris)
310 unsigned short *pRes;
312 _imageSphere->GetDimensions(dimRes);
315 for ( ii=(int)(i-delta2) ; ii<=(int)(i+delta2) ; ii++ )
317 for ( jj=(int)(j-delta2) ; jj<=(int)(j+delta2) ; jj++ )
319 if ( (ii>=0)&&(ii<dimRes[0]) &&
320 (jj>=0)&&(jj<dimRes[1]) )
322 pRes = (unsigned short*)_imageSphere->GetScalarPointer( ii , jj , (int)id );
330 //----------------------------------------------------------------------------
332 void wxSphereView::SetXYZtoParent(double i, double j)
335 double factor = 0.75;
336 double radio2 = _radio*_radio;
337 double pxx,pyy,d2x,d2y;
338 double cc[3],p[3],pp[3];
343 int dimRes[3],dimOrig[3];
344 _imageSphere->GetDimensions(dimRes);
347 _imageDataOriginal->GetDimensions(dimOrig);
349 p[0] = (i - d2x)*factor;
351 p[1] = (j - d2y)*factor;
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]) )
362 if (_vtkmprbasedata){
363 _vtkmprbasedata->SetX(pp[0]);
364 _vtkmprbasedata->SetY(pp[1]);
365 _vtkmprbasedata->SetZ(pp[2]);
372 //----------------------------------------------------------------------------
374 void wxSphereView::FiltreImageB(int id, double radio, bool ok,int deltaTMP)
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];
385 unsigned short *pOrig;
386 int dimRes[3],dimOrig[3];
389 _imageSphere->GetExtent(ext);
390 _imageSphere->GetDimensions(dimRes);
392 //_imageSphere->SetExtent(0,dimRes[0]-1,0,dimRes[1]-1,0,dimRes[2]-1);
393 _imageSphere->SetExtent(ext);
397 // double deltaTMP=_delta;
398 _imageDataOriginal->GetDimensions(dimOrig);
402 limitA = (int) ( (-radioB/factor)+d2x );
403 limitB = (int) ( (radioB/factor)+d2x );
412 for ( i=start ; i<end ; i=i+deltaTMP )
414 p[0] = (i - d2x)*factor;
416 for (j=start;j<end;j=j+deltaTMP)
418 p[1] = (j - d2y)*factor;
422 if (( ((ok==false) && (!((i>limitA) && (i<limitB) && (j>limitA) && (j<limitB)))) )
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]) )
434 pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( (int)(pp[0]) , (int)(pp[1]) , (int)(pp[2]) );
435 SetVoxel(i,j,deltaTMP,id,*pOrig);
437 SetVoxel(i,j,deltaTMP,id,2000);
440 SetVoxel(i,j,deltaTMP,id,0);
446 _imageSphere->Modified();
447 _imageSphere->Update();
453 //----------------------------------------------------------------------------
455 void wxSphereView::FiltreImage(int id, double radio)
458 _transform -> Identity();
459 _transform -> Concatenate(_transform1);
460 _transform -> Concatenate(_transform2);
462 FiltreImageB(id,radio,false, _delta);
463 FiltreImageB(id,radio,true, 1);
467 //----------------------------------------------------------------------------
470 void wxSphereView::FiltreImage(int id, double radio)
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];
480 unsigned short *pOrig;
481 int dimRes[3],dimOrig[3];
483 _imageSphere->GetDimensions(dimRes);
484 _imageSphere->SetExtent(0,dimRes[0]-1,0,dimRes[1]-1,0,dimRes[2]-1);
487 double deltaTMP=_delta;
488 _imageDataOriginal->GetDimensions(dimOrig);
490 for ( i=0 ; i<dimRes[0] ; i=i+deltaTMP )
492 p[0] = (i - d2x)*0.75;
494 for (j=0;j<dimRes[1];j=j+deltaTMP)
496 p[1] = (j - d2y)*0.75;
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]) )
510 pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( pp[0] , pp[1] , pp[2] );
511 SetVoxel(i,j,deltaTMP,id,*pOrig);
513 SetVoxel(i,j,deltaTMP,id,2000);
516 SetVoxel(i,j,deltaTMP,id,0);
524 for ( i=0 ; i<dimRes[0] ; i=i+deltaTMP )
526 p[0] = (i - d2x)*0.75;
528 for (j=0;j<dimRes[1];j=j+deltaTMP)
530 p[1] = (j - d2y)*0.75;
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]) )
543 pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( pp[0] , pp[1] , pp[2] );
544 SetVoxel(i,j,deltaTMP,id,*pOrig);
546 SetVoxel(i,j,deltaTMP,id,2000);
549 SetVoxel(i,j,deltaTMP,id,0);
555 _imageSphere->Modified();
556 _imageSphere->Update();
560 void wxSphereView::FiltreImage(vtkImageData *imageSphere)
563 imageSphere->GetDimensions(dim);
564 for (i=0;i<dim[0];i++)
566 for (j=0;j<dim[1];j++)
568 for (k=0;k<dim[2];k++)
570 unsigned short *pRes=(unsigned short*)imageSphere->GetScalarPointer (i,j,k);
577 double cc[3],p1[3],p2[3],pp1[3],pp2[3];
581 double r1 = _sl_radio->GetValue() - _sl_thickness->GetValue()/2;
582 double r2 = _sl_radio->GetValue() + _sl_thickness->GetValue()/2;
587 double alpha= _sl_alpha->GetValue();
588 double beta = _sl_beta->GetValue();
591 for (angA=-deltaA;angA<deltaA;angA++)
593 for (angB=-deltaA;angB<deltaA;angB++)
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);
606 //----------------------------------------------------------------------------
608 void wxSphereView::InitSphere(double points[4][3])
611 double r = SphereFindCenter(points,cc); // 4-points , center
614 _centerX = (int)(cc[0]);
615 _centerY = (int)(cc[1]);
616 _centerZ = (int)(cc[2]);
619 _imageDataOriginal->GetDimensions(dim);
620 _centerX = (int)(dim[0]/2);
621 _centerY = (int)(dim[1]/2);
622 _centerZ = (int)(dim[2]/2);
626 //----------------------------------------------------------------------------
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])
634 double r, m11, m12, m13, m14, m15;
637 for (i = 0; i < 4; i++) // find minor 11
644 m11 = determinant( a, 4 );
646 for (i = 0; i < 4; i++) // find minor 12
648 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
653 m12 = determinant( a, 4 );
655 for (i = 0; i < 4; i++) // find minor 13
657 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
662 m13 = determinant( a, 4 );
664 for (i = 0; i < 4; i++) // find minor 14
666 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
671 m14 = determinant( a, 4 );
674 for (i = 0; i < 4; i++) // find minor 15
676 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
681 m15 = determinant( a, 4 );
690 cc[0] = 0.5*m12/m11; //cx
691 cc[1] = -0.5*m13/m11; //cy
692 cc[2] = 0.5*m14/m11; //cz
694 r = sqrt( cc[0]*cc[0] + cc[1]*cc[1] + cc[2]*cc[2] - m15/m11 );
697 return r; // the radius
699 //----------------------------------------------------------------------------
701 // Recursive definition of determinate using expansion by minors.
702 double wxSphereView::determinant(double a[4][4], int n)
716 if (n == 2) // terminate recursion
718 d = a[0][0]*a[1][1] - a[1][0]*a[0][1];
723 for (j1 = 0; j1 < n; j1++ ) // do each column
725 for (i = 1; i < n; i++) // create minor
728 for (j = 0; j < n; j++)
730 if (j == j1) continue;
731 m[i-1][j2] = a[i][j];
736 // sum (+/-)cofactor * minor
737 d = d + pow(-1.0, j1)*a[0][j1]*determinant( m, n-1 );