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"
28 wxSphereView::wxSphereView( wxWindow *parent, vtkMPRBaseData *vtkmprbasedata/*, vtkImageData *imageData */)
29 : wxVtk2DBaseView(parent)
32 _vtkmprbasedata = vtkmprbasedata;
35 // _imageDataOriginal = imageData;
38 _imageSphere = vtkImageData::New();
39 _imageSphere->SetDimensions (150,150,500);
42 // _imageSphere->SetScalarTypeToUnsignedShort();
43 _imageSphere->SetScalarType( _imageDataOriginal->GetScalarType() );
45 _imageSphere->AllocateScalars();
46 _imageSphere->Update();
48 //EED ???? vtkBaseData no esta compartido con los otros objetos .. PLOP
49 // vtkBaseData *vtkbasedata = new vtkBaseData();
50 // vtkbasedata->SetMarImageData( new marImageData(_imageSphere) );
51 // this->SetVtkBaseData(vtkbasedata);
52 this->SetVtkBaseData(_vtkmprbasedata);
54 _transform = vtkTransform::New();
55 _transform1 = vtkTransform::New();
56 _transform2 = vtkTransform::New();
57 _transform ->Identity();
58 _transform1->Identity();
59 _transform2->Identity();
63 //-------------------------------------------------------------------
64 wxSphereView::~wxSphereView()
66 _transform -> Delete();
67 _transform1 -> Delete();
68 _transform2 -> Delete();
72 //----------------------------------------------------------------------------
73 void wxSphereView::SetImage()
76 _imageDataOriginal=_vtkmprbasedata->GetImageData();
79 //----------------------------------------------------------------------------
80 double wxSphereView::GetRadio()
85 //----------------------------------------------------------------------------
86 void wxSphereView::SetRadio(double radio)
95 //----------------------------------------------------------------------------
96 void wxSphereView::Configure()
98 wxVtk2DBaseView::Configure();
101 vtkInteractorStyleBaseView2D *style2D = vtkInteractorStyleBaseView2D::New();
102 manualInteractorWindowLevel *_manualinteractorwindowlevel= new manualInteractorWindowLevel();
103 style2D->SetInteractorWindowLevel( _manualinteractorwindowlevel );
104 // vtkInteractorScrollZ *_vtkInteractorScrollZ = new vtkInteractorScrollZ();
105 // style2D->SetInteractorScrollZ(_vtkInteractorScrollZ);
106 SetInteractorStyleImage( style2D );
110 _vtkinteractorstylesphere = new vtkInteractorStyleSphere();
111 ((vtkInteractorStyleBaseView*)GetInteractorStyleBaseView())->AddInteractorStyleMaracas( _vtkinteractorstylesphere );
113 // EED purify 12/sep/2006
128 //----------------------------------------------------------------------------
129 void wxSphereView::RefreshPoint()
131 double x = _vtkmprbasedata->GetX() - _centerX;
132 double y = _vtkmprbasedata->GetY() - _centerY;
133 double z = _vtkmprbasedata->GetZ() - _centerZ;
134 double alpha= atan2(x,z);
135 double beta = atan2( y , sqrt(z*z+x*x) );
136 alpha = alpha*180/3.1416;
137 beta = beta*180/3.1416;
138 _transform1->Identity();
139 _transform1->RotateY(alpha);
140 _transform1->RotateX(-beta);
141 _radio= sqrt(x*x + y*y +z*z);
147 //-------------------------------------------------------------------
148 void wxSphereView::Refresh( )
151 UpdateColorWindowLevel();
152 wxVtkBaseView::Refresh();
157 //----------------------------------------------------------------------------
158 void wxSphereView::RefreshView()
161 _imageViewer2XYZ->GetVtkImageViewer2()->SetInput ( _imageSphere );
163 // UpdateColorWindowLevel();
164 // wxVtk2DBaseView::Refresh();
168 //----------------------------------------------------------------------------
169 void wxSphereView::RotationEnd()
171 _transform1->RotateWXYZ(_ang,_vxb,_vyb,0);
172 _transform2->Identity();
177 //----------------------------------------------------------------------------
178 void wxSphereView::RotationStart(double vx, double vy, bool ok_v, bool ok_ang)
182 _ang = -sqrt( vx*vx + vy*vy ) / 1.0;
188 _transform2->Identity();
189 _transform2->RotateWXYZ(_ang,_vxb,_vyb,0);
194 //----------------------------------------------------------------------------
195 void wxSphereView::GetPointSphere(double p[3],double r1,double angA,double angB)
201 vtkTransform *transform = vtkTransform::New();
202 transform->Identity();
203 transform->RotateX(angB);
204 transform->RotateZ(angA);
205 transform->TransformPoint(in,out);
212 //----------------------------------------------------------------------------
213 void wxSphereView::RotatePointOverTheSphere( double pp[3], double p[3],double cc[3])
216 _transform->TransformPoint(p,out);
217 pp[0] = out[0] + cc[0];
218 pp[1] = out[1] + cc[1];
219 pp[2] = out[2] + cc[2];
222 //----------------------------------------------------------------------------
223 void wxSphereView::ResetlstId()
225 int i,size=_lstId.size();
226 for (i=size-1;i>=0;i--)
233 //----------------------------------------------------------------------------
234 int wxSphereView::GetIdOfImage(double radio)
238 _imageSphere->GetDimensions(dim);
239 int sizeMaxList = dim[2];
240 // Search in list >> alpha beta radio
241 int i,size=_lstId.size();
242 for (i=0; i<size;i++)
244 //idAlBeRa *tmp=_lstId[i]; // JPRx
245 if ((_lstId[i]->_radio==radio) && (_lstId[i]->_deltavoxel==_delta))
247 return _lstId[i]->_id;
250 if (size>sizeMaxList)
252 delete _lstId[size-1];
257 id = id % sizeMaxList;
261 FiltreImage(id,radio);
262 _lstId.insert(_lstId.begin(),1,new idAlBeRa(id,radio,_delta) );
266 //----------------------------------------------------------------------------
267 void wxSphereView::DefineImageSphere()
270 id=GetIdOfImage( _radio );
272 // GetVtkBaseData()->SetZ( id );
274 _imageViewer2XYZ->GetVtkImageViewer2()->SetSlice ( id );
277 //----------------------------------------------------------------------------
278 void wxSphereView::SetDeltaVoxel(int delta)
283 //----------------------------------------------------------------------------
284 void wxSphereView::SetVoxel(double i, double j, int delta,double id, double gris)
287 unsigned short *pRes;
289 _imageSphere->GetDimensions(dimRes);
292 for ( ii=(int)(i-delta2) ; ii<=(int)(i+delta2) ; ii++ )
294 for ( jj=(int)(j-delta2) ; jj<=(int)(j+delta2) ; jj++ )
296 if ( (ii>=0)&&(ii<dimRes[0]) &&
297 (jj>=0)&&(jj<dimRes[1]) )
300 // pRes = (unsigned short*)_imageSphere->GetScalarPointer( ii , jj , (int)id );
302 _imageSphere->SetScalarComponentFromDouble( ii , jj , (int)id , 0, gris );
308 //----------------------------------------------------------------------------
309 void wxSphereView::SetXYZtoParent(double i, double j)
312 double factor = 0.75;
313 double radio2 = _radio*_radio;
314 double pxx,pyy,d2x,d2y;
315 double cc[3],p[3],pp[3];
320 int dimRes[3],dimOrig[3];
321 _imageSphere->GetDimensions(dimRes);
324 _imageDataOriginal->GetDimensions(dimOrig);
325 p[0] = (i - d2x)*factor;
327 p[1] = (j - d2y)*factor;
333 RotatePointOverTheSphere( pp, p,cc);
334 if ( (pp[0]>=0) && (pp[0]<dimOrig[0]) &&
335 (pp[1]>=0) && (pp[1]<dimOrig[1]) &&
336 (pp[2]>=0) && (pp[2]<dimOrig[2]) )
338 if (_vtkmprbasedata){
339 _vtkmprbasedata->SetX(pp[0]);
340 _vtkmprbasedata->SetY(pp[1]);
341 _vtkmprbasedata->SetZ(pp[2]);
347 //----------------------------------------------------------------------------
348 void wxSphereView::FiltreImageB(int id, double radio, bool ok,int deltaTMP)
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];
360 unsigned short *pOrig;
361 int dimRes[3],dimOrig[3];
364 _imageSphere->GetExtent(ext);
365 _imageSphere->GetDimensions(dimRes);
367 //_imageSphere->SetExtent(0,dimRes[0]-1,0,dimRes[1]-1,0,dimRes[2]-1);
368 _imageSphere->SetExtent(ext);
372 _imageDataOriginal->GetDimensions(dimOrig);
375 limitA = (int) ( (-radioB/factor)+d2x );
376 limitB = (int) ( (radioB/factor)+d2x );
385 for ( i=start ; i<end ; i=i+deltaTMP )
387 p[0] = (i - d2x)*factor;
389 for (j=start;j<end;j=j+deltaTMP)
391 p[1] = (j - d2y)*factor;
395 if (( ((ok==false) && (!((i>limitA) && (i<limitB) && (j>limitA) && (j<limitB)))) )
402 RotatePointOverTheSphere( pp, p,cc);
403 if ( (pp[0]>=0) && (pp[0]<dimOrig[0]) &&
404 (pp[1]>=0) && (pp[1]<dimOrig[1]) &&
405 (pp[2]>=0) && (pp[2]<dimOrig[2]) )
408 // pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( (int)(pp[0]) , (int)(pp[1]) , (int)(pp[2]) );
409 // SetVoxel(i,j,deltaTMP,id,*pOrig);
410 value=_imageDataOriginal->GetScalarComponentAsDouble( (int)(pp[0]) , (int)(pp[1]) , (int)(pp[2]), 0);
411 SetVoxel(i,j,deltaTMP,id,value);
413 SetVoxel(i,j,deltaTMP,id,2000);
416 SetVoxel(i,j,deltaTMP,id,0);
421 _imageSphere->Modified();
422 _imageSphere->Update();
425 //----------------------------------------------------------------------------
426 void wxSphereView::FiltreImage(int id, double radio)
428 _transform -> Identity();
429 _transform -> Concatenate(_transform1);
430 _transform -> Concatenate(_transform2);
431 FiltreImageB(id,radio,false, _delta);
432 FiltreImageB(id,radio,true, 1);
435 //----------------------------------------------------------------------------
438 void wxSphereView::FiltreImage(int id, double radio)
440 double radio2 = radio*radio;
441 double radio2TMP= (radio/2)*(radio/2);
442 double pxx,pyy,d2x,d2y;
443 double cc[3],p[3],pp[3];
448 unsigned short *pOrig;
449 int dimRes[3],dimOrig[3];
451 _imageSphere->GetDimensions(dimRes);
452 _imageSphere->SetExtent(0,dimRes[0]-1,0,dimRes[1]-1,0,dimRes[2]-1);
455 double deltaTMP=_delta;
456 _imageDataOriginal->GetDimensions(dimOrig);
458 for ( i=0 ; i<dimRes[0] ; i=i+deltaTMP )
460 p[0] = (i - d2x)*0.75;
462 for (j=0;j<dimRes[1];j=j+deltaTMP)
464 p[1] = (j - d2y)*0.75;
473 RotatePointOverTheSphere( pp, p,cc);
474 if ( (pp[0]>=0) && (pp[0]<dimOrig[0]) &&
475 (pp[1]>=0) && (pp[1]<dimOrig[1]) &&
476 (pp[2]>=0) && (pp[2]<dimOrig[2]) )
478 pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( pp[0] , pp[1] , pp[2] );
479 SetVoxel(i,j,deltaTMP,id,*pOrig);
481 SetVoxel(i,j,deltaTMP,id,2000);
484 SetVoxel(i,j,deltaTMP,id,0);
492 for ( i=0 ; i<dimRes[0] ; i=i+deltaTMP )
494 p[0] = (i - d2x)*0.75;
496 for (j=0;j<dimRes[1];j=j+deltaTMP)
498 p[1] = (j - d2y)*0.75;
506 RotatePointOverTheSphere( pp, p,cc);
507 if ( (pp[0]>=0) && (pp[0]<dimOrig[0]) &&
508 (pp[1]>=0) && (pp[1]<dimOrig[1]) &&
509 (pp[2]>=0) && (pp[2]<dimOrig[2]) )
511 pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( pp[0] , pp[1] , pp[2] );
512 SetVoxel(i,j,deltaTMP,id,*pOrig);
514 SetVoxel(i,j,deltaTMP,id,2000);
517 SetVoxel(i,j,deltaTMP,id,0);
523 _imageSphere->Modified();
524 _imageSphere->Update();
528 void wxSphereView::FiltreImage(vtkImageData *imageSphere)
531 imageSphere->GetDimensions(dim);
532 for (i=0;i<dim[0];i++)
534 for (j=0;j<dim[1];j++)
536 for (k=0;k<dim[2];k++)
538 unsigned short *pRes=(unsigned short*)imageSphere->GetScalarPointer (i,j,k);
545 double cc[3],p1[3],p2[3],pp1[3],pp2[3];
549 double r1 = _sl_radio->GetValue() - _sl_thickness->GetValue()/2;
550 double r2 = _sl_radio->GetValue() + _sl_thickness->GetValue()/2;
555 double alpha= _sl_alpha->GetValue();
556 double beta = _sl_beta->GetValue();
559 for (angA=-deltaA;angA<deltaA;angA++)
561 for (angB=-deltaA;angB<deltaA;angB++)
563 GetPointSphere(p1,r1,angA,angB);
564 GetPointSphere(p2,r2,angA,angB);
565 RotatePointOverTheSphere( pp1, alpha, beta, p1 ,cc );
566 RotatePointOverTheSphere( pp2, alpha, beta, p2 ,cc );
567 TransferePoints(pp1,pp2,angA+alpha+180,angB+beta+90,imageSphere);
573 //----------------------------------------------------------------------------
575 void wxSphereView::TransferePoints(double pp1[3],double pp2[3],double AngX,double AngY,vtkImageData *image)
578 double difX = pp2[0]-pp1[0];
579 double difY = pp2[1]-pp1[1];
580 double difZ = pp2[2]-pp1[2];
585 _imageDataOriginal->GetDimensions(dimOrg);
586 image->GetDimensions(dimRes);
591 int xx=-1,yy=-1,zz=-1;
595 xx = (int) (x1+t*difX);
596 yy = (int) (y1+t*difY);
597 zz = (int) (z1+t*difZ);
600 if ((xx>=0) && (xx<dimOrg[0]) && (yy>=0) && (yy<dimOrg[1]) && (zz>=0) && (zz<dimOrg[2]) &&
601 (AngX>=0) && (AngX<dimRes[0]) && (AngY>=0) && (AngY<dimRes[1]) && (z>=0) && (z<dimRes[2]) )
603 unsigned short *pOrg=(unsigned short*)_imageDataOriginal->GetScalarPointer (xx,yy,zz);
604 unsigned short *pRes=(unsigned short*)image->GetScalarPointer( (int)AngX , (int)AngY , z );
616 //----------------------------------------------------------------------------
618 void wxSphereView::InitSphere(double points[4][3])
621 double r = SphereFindCenter(points,cc); // 4-points , center
624 _centerX = (int)(cc[0]);
625 _centerY = (int)(cc[1]);
626 _centerZ = (int)(cc[2]);
629 _imageDataOriginal->GetDimensions(dim);
630 _centerX = (int)(dim[0]/2);
631 _centerY = (int)(dim[1]/2);
632 _centerZ = (int)(dim[2]/2);
636 //----------------------------------------------------------------------------
638 // Calculate center and radius of sphere given four points
639 // http://home.att.net/~srschmitt/script_sphere_solver.html
640 // source code HTML <script language=JavaScript>
641 double wxSphereView::SphereFindCenter(double P[4][3], double cc[3])
644 double r, m11, m12, m13, m14, m15;
647 for (i = 0; i < 4; i++) // find minor 11
654 m11 = determinant( a, 4 );
656 for (i = 0; i < 4; i++) // find minor 12
658 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
663 m12 = determinant( a, 4 );
665 for (i = 0; i < 4; i++) // find minor 13
667 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
672 m13 = determinant( a, 4 );
674 for (i = 0; i < 4; i++) // find minor 14
676 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
681 m14 = determinant( a, 4 );
684 for (i = 0; i < 4; i++) // find minor 15
686 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
691 m15 = determinant( a, 4 );
700 cc[0] = 0.5*m12/m11; //cx
701 cc[1] = -0.5*m13/m11; //cy
702 cc[2] = 0.5*m14/m11; //cz
704 r = sqrt( cc[0]*cc[0] + cc[1]*cc[1] + cc[2]*cc[2] - m15/m11 );
707 return r; // the radius
709 //----------------------------------------------------------------------------
711 // Recursive definition of determinate using expansion by minors.
712 double wxSphereView::determinant(double a[4][4], int n)
726 if (n == 2) // terminate recursion
728 d = a[0][0]*a[1][1] - a[1][0]*a[0][1];
733 for (j1 = 0; j1 < n; j1++ ) // do each column
735 for (i = 1; i < n; i++) // create minor
738 for (j = 0; j < n; j++)
740 if (j == j1) continue;
741 m[i-1][j2] = a[i][j];
746 // sum (+/-)cofactor * minor
747 d = d + pow(-1.0, j1)*a[0][j1]*determinant( m, n-1 );