]> Creatis software - creaMaracasVisu.git/blob - lib/maracasVisuLib/src/interface/wxWindows/widgets/wxSphereView.cxx
*** empty log message ***
[creaMaracasVisu.git] / 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