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);
41 //EED 2017-01-01 Migration VTK7
42 #if VTK_MAJOR_VERSION <= 5
44 // _imageSphere->SetScalarTypeToUnsignedShort();
45 _imageSphere->SetScalarType( _imageDataOriginal->GetScalarType() );
46 _imageSphere->AllocateScalars();
47 _imageSphere->Update();
49 _imageSphere->AllocateScalars(_imageDataOriginal->GetScalarType(),1);
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);
60 _transform = vtkTransform::New();
61 _transform1 = vtkTransform::New();
62 _transform2 = vtkTransform::New();
63 _transform ->Identity();
64 _transform1->Identity();
65 _transform2->Identity();
69 //-------------------------------------------------------------------
70 wxSphereView::~wxSphereView()
72 _transform -> Delete();
73 _transform1 -> Delete();
74 _transform2 -> Delete();
78 //----------------------------------------------------------------------------
79 void wxSphereView::SetImage()
82 _imageDataOriginal=_vtkmprbasedata->GetImageData();
85 //----------------------------------------------------------------------------
86 double wxSphereView::GetRadio()
91 //----------------------------------------------------------------------------
92 void wxSphereView::SetRadio(double radio)
101 //----------------------------------------------------------------------------
102 void wxSphereView::Configure()
104 wxVtk2DBaseView::Configure();
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 );
116 _vtkinteractorstylesphere = new vtkInteractorStyleSphere();
117 ((vtkInteractorStyleBaseView*)GetInteractorStyleBaseView())->AddInteractorStyleMaracas( _vtkinteractorstylesphere );
119 // EED purify 12/sep/2006
134 //----------------------------------------------------------------------------
135 void wxSphereView::RefreshPoint()
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);
153 //-------------------------------------------------------------------
154 void wxSphereView::Refresh( )
157 UpdateColorWindowLevel();
158 wxVtkBaseView::Refresh();
163 //----------------------------------------------------------------------------
164 void wxSphereView::RefreshView()
167 //EED 2017-01-01 Migration VTK7
168 #if VTK_MAJOR_VERSION <= 5
169 _imageViewer2XYZ->GetVtkImageViewer2()->SetInput( _imageSphere );
171 _imageViewer2XYZ->GetVtkImageViewer2()->SetInputData( _imageSphere );
174 // UpdateColorWindowLevel();
175 // wxVtk2DBaseView::Refresh();
179 //----------------------------------------------------------------------------
180 void wxSphereView::RotationEnd()
182 _transform1->RotateWXYZ(_ang,_vxb,_vyb,0);
183 _transform2->Identity();
188 //----------------------------------------------------------------------------
189 void wxSphereView::RotationStart(double vx, double vy, bool ok_v, bool ok_ang)
193 _ang = -sqrt( vx*vx + vy*vy ) / 1.0;
199 _transform2->Identity();
200 _transform2->RotateWXYZ(_ang,_vxb,_vyb,0);
205 //----------------------------------------------------------------------------
206 void wxSphereView::GetPointSphere(double p[3],double r1,double angA,double angB)
212 vtkTransform *transform = vtkTransform::New();
213 transform->Identity();
214 transform->RotateX(angB);
215 transform->RotateZ(angA);
216 transform->TransformPoint(in,out);
223 //----------------------------------------------------------------------------
224 void wxSphereView::RotatePointOverTheSphere( double pp[3], double p[3],double cc[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];
233 //----------------------------------------------------------------------------
234 void wxSphereView::ResetlstId()
236 int i,size=_lstId.size();
237 for (i=size-1;i>=0;i--)
244 //----------------------------------------------------------------------------
245 int wxSphereView::GetIdOfImage(double radio)
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++)
255 //idAlBeRa *tmp=_lstId[i]; // JPRx
256 if ((_lstId[i]->_radio==radio) && (_lstId[i]->_deltavoxel==_delta))
258 return _lstId[i]->_id;
261 if (size>sizeMaxList)
263 delete _lstId[size-1];
268 id = id % sizeMaxList;
272 FiltreImage(id,radio);
273 _lstId.insert(_lstId.begin(),1,new idAlBeRa(id,radio,_delta) );
277 //----------------------------------------------------------------------------
278 void wxSphereView::DefineImageSphere()
281 id=GetIdOfImage( _radio );
283 // GetVtkBaseData()->SetZ( id );
285 _imageViewer2XYZ->GetVtkImageViewer2()->SetSlice ( id );
288 //----------------------------------------------------------------------------
289 void wxSphereView::SetDeltaVoxel(int delta)
294 //----------------------------------------------------------------------------
295 void wxSphereView::SetVoxel(double i, double j, int delta,double id, double gris)
298 unsigned short *pRes;
300 _imageSphere->GetDimensions(dimRes);
303 for ( ii=(int)(i-delta2) ; ii<=(int)(i+delta2) ; ii++ )
305 for ( jj=(int)(j-delta2) ; jj<=(int)(j+delta2) ; jj++ )
307 if ( (ii>=0)&&(ii<dimRes[0]) &&
308 (jj>=0)&&(jj<dimRes[1]) )
311 // pRes = (unsigned short*)_imageSphere->GetScalarPointer( ii , jj , (int)id );
313 _imageSphere->SetScalarComponentFromDouble( ii , jj , (int)id , 0, gris );
319 //----------------------------------------------------------------------------
320 void wxSphereView::SetXYZtoParent(double i, double j)
323 double factor = 0.75;
324 double radio2 = _radio*_radio;
325 double pxx,pyy,d2x,d2y;
326 double cc[3],p[3],pp[3];
331 int dimRes[3],dimOrig[3];
332 _imageSphere->GetDimensions(dimRes);
335 _imageDataOriginal->GetDimensions(dimOrig);
336 p[0] = (i - d2x)*factor;
338 p[1] = (j - d2y)*factor;
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]) )
349 if (_vtkmprbasedata){
350 _vtkmprbasedata->SetX(pp[0]);
351 _vtkmprbasedata->SetY(pp[1]);
352 _vtkmprbasedata->SetZ(pp[2]);
358 //----------------------------------------------------------------------------
359 void wxSphereView::FiltreImageB(int id, double radio, bool ok,int deltaTMP)
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];
371 unsigned short *pOrig;
372 int dimRes[3],dimOrig[3];
375 _imageSphere->GetExtent(ext);
376 _imageSphere->GetDimensions(dimRes);
378 //_imageSphere->SetExtent(0,dimRes[0]-1,0,dimRes[1]-1,0,dimRes[2]-1);
379 _imageSphere->SetExtent(ext);
383 _imageDataOriginal->GetDimensions(dimOrig);
386 limitA = (int) ( (-radioB/factor)+d2x );
387 limitB = (int) ( (radioB/factor)+d2x );
396 for ( i=start ; i<end ; i=i+deltaTMP )
398 p[0] = (i - d2x)*factor;
400 for (j=start;j<end;j=j+deltaTMP)
402 p[1] = (j - d2y)*factor;
406 if (( ((ok==false) && (!((i>limitA) && (i<limitB) && (j>limitA) && (j<limitB)))) )
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]) )
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);
424 SetVoxel(i,j,deltaTMP,id,2000);
427 SetVoxel(i,j,deltaTMP,id,0);
432 _imageSphere->Modified();
433 //EED 2017-01-01 Migration VTK7
434 #if VTK_MAJOR_VERSION <= 5
435 _imageSphere->Update();
441 //----------------------------------------------------------------------------
442 void wxSphereView::FiltreImage(int id, double radio)
444 _transform -> Identity();
445 _transform -> Concatenate(_transform1);
446 _transform -> Concatenate(_transform2);
447 FiltreImageB(id,radio,false, _delta);
448 FiltreImageB(id,radio,true, 1);
451 //----------------------------------------------------------------------------
454 void wxSphereView::FiltreImage(int id, double radio)
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];
464 unsigned short *pOrig;
465 int dimRes[3],dimOrig[3];
467 _imageSphere->GetDimensions(dimRes);
468 _imageSphere->SetExtent(0,dimRes[0]-1,0,dimRes[1]-1,0,dimRes[2]-1);
471 double deltaTMP=_delta;
472 _imageDataOriginal->GetDimensions(dimOrig);
474 for ( i=0 ; i<dimRes[0] ; i=i+deltaTMP )
476 p[0] = (i - d2x)*0.75;
478 for (j=0;j<dimRes[1];j=j+deltaTMP)
480 p[1] = (j - d2y)*0.75;
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]) )
494 pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( pp[0] , pp[1] , pp[2] );
495 SetVoxel(i,j,deltaTMP,id,*pOrig);
497 SetVoxel(i,j,deltaTMP,id,2000);
500 SetVoxel(i,j,deltaTMP,id,0);
508 for ( i=0 ; i<dimRes[0] ; i=i+deltaTMP )
510 p[0] = (i - d2x)*0.75;
512 for (j=0;j<dimRes[1];j=j+deltaTMP)
514 p[1] = (j - d2y)*0.75;
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]) )
527 pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( pp[0] , pp[1] , pp[2] );
528 SetVoxel(i,j,deltaTMP,id,*pOrig);
530 SetVoxel(i,j,deltaTMP,id,2000);
533 SetVoxel(i,j,deltaTMP,id,0);
539 _imageSphere->Modified();
540 _imageSphere->Update();
544 void wxSphereView::FiltreImage(vtkImageData *imageSphere)
547 imageSphere->GetDimensions(dim);
548 for (i=0;i<dim[0];i++)
550 for (j=0;j<dim[1];j++)
552 for (k=0;k<dim[2];k++)
554 unsigned short *pRes=(unsigned short*)imageSphere->GetScalarPointer (i,j,k);
561 double cc[3],p1[3],p2[3],pp1[3],pp2[3];
565 double r1 = _sl_radio->GetValue() - _sl_thickness->GetValue()/2;
566 double r2 = _sl_radio->GetValue() + _sl_thickness->GetValue()/2;
571 double alpha= _sl_alpha->GetValue();
572 double beta = _sl_beta->GetValue();
575 for (angA=-deltaA;angA<deltaA;angA++)
577 for (angB=-deltaA;angB<deltaA;angB++)
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);
589 //----------------------------------------------------------------------------
591 void wxSphereView::TransferePoints(double pp1[3],double pp2[3],double AngX,double AngY,vtkImageData *image)
594 double difX = pp2[0]-pp1[0];
595 double difY = pp2[1]-pp1[1];
596 double difZ = pp2[2]-pp1[2];
601 _imageDataOriginal->GetDimensions(dimOrg);
602 image->GetDimensions(dimRes);
607 int xx=-1,yy=-1,zz=-1;
611 xx = (int) (x1+t*difX);
612 yy = (int) (y1+t*difY);
613 zz = (int) (z1+t*difZ);
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]) )
619 unsigned short *pOrg=(unsigned short*)_imageDataOriginal->GetScalarPointer (xx,yy,zz);
620 unsigned short *pRes=(unsigned short*)image->GetScalarPointer( (int)AngX , (int)AngY , z );
632 //----------------------------------------------------------------------------
634 void wxSphereView::InitSphere(double points[4][3])
637 double r = SphereFindCenter(points,cc); // 4-points , center
640 _centerX = (int)(cc[0]);
641 _centerY = (int)(cc[1]);
642 _centerZ = (int)(cc[2]);
645 _imageDataOriginal->GetDimensions(dim);
646 _centerX = (int)(dim[0]/2);
647 _centerY = (int)(dim[1]/2);
648 _centerZ = (int)(dim[2]/2);
652 //----------------------------------------------------------------------------
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])
660 double r, m11, m12, m13, m14, m15;
663 for (i = 0; i < 4; i++) // find minor 11
670 m11 = determinant( a, 4 );
672 for (i = 0; i < 4; i++) // find minor 12
674 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
679 m12 = determinant( a, 4 );
681 for (i = 0; i < 4; i++) // find minor 13
683 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
688 m13 = determinant( a, 4 );
690 for (i = 0; i < 4; i++) // find minor 14
692 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
697 m14 = determinant( a, 4 );
700 for (i = 0; i < 4; i++) // find minor 15
702 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
707 m15 = determinant( a, 4 );
716 cc[0] = 0.5*m12/m11; //cx
717 cc[1] = -0.5*m13/m11; //cy
718 cc[2] = 0.5*m14/m11; //cz
720 r = sqrt( cc[0]*cc[0] + cc[1]*cc[1] + cc[2]*cc[2] - m15/m11 );
723 return r; // the radius
725 //----------------------------------------------------------------------------
727 // Recursive definition of determinate using expansion by minors.
728 double wxSphereView::determinant(double a[4][4], int n)
742 if (n == 2) // terminate recursion
744 d = a[0][0]*a[1][1] - a[1][0]*a[0][1];
749 for (j1 = 0; j1 < n; j1++ ) // do each column
751 for (i = 1; i < n; i++) // create minor
754 for (j = 0; j < n; j++)
756 if (j == j1) continue;
757 m[i-1][j2] = a[i][j];
762 // sum (+/-)cofactor * minor
763 d = d + pow(-1.0, j1)*a[0][j1]*determinant( m, n-1 );