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 /*=========================================================================
29 =========================================================================*/
30 #include "axisExtractor02.h"
31 #include "vtkPolyData.h"
32 #include "vtkObjectFactory.h"
33 #include "vtkImageThreshold.h"
34 #include "vtkImageCast.h"
35 #include "vtkImageSeedConnectivity.h"
36 #include "vtkImageData.h"
37 #include "vtkMarchingCubes.h"
38 #include "vtkDoubleArray.h"
39 #include "vtkPointData.h"
40 #include "vtkExtractVOI.h"
41 #include "vtkImageConstantPad.h"
42 #include "vtkImageTranslateExtent.h"
44 #include "vtkTransformPolyDataFilter.h"
45 #include "vtkTransform.h"
52 vtkStandardNewMacro(axisExtractor02);
55 //----------------------------------------------------------------------------
56 void axisExtractor02::SetParam(double value)
63 //----------------------------------------------------------------------------
64 double axisExtractor02::GetParam()
72 //----------------------------------------------------------------------------
73 void axisExtractor02::SetParam2(double value)
80 //----------------------------------------------------------------------------
81 double axisExtractor02::GetParam2()
89 //----------------------------------------------------------------------------
90 void axisExtractor02::SetParam3(double value)
97 //----------------------------------------------------------------------------
98 double axisExtractor02::GetParam3()
105 //----------------------------------------------------------------------------
106 void axisExtractor02::SetMaxant(int value)
113 //----------------------------------------------------------------------------
114 int axisExtractor02::GetMaxant()
121 //----------------------------------------------------------------------------
122 void axisExtractor02::SetMinant(int value)
129 //----------------------------------------------------------------------------
130 int axisExtractor02::GetMinant()
136 //----------------------------------------------------------------------------
137 void axisExtractor02::SetPoint(double value[3])
143 this->GetInput()->GetSpacing (spa);
145 value[0]=value[0]+maxant*spa[0];
146 value[1]=value[1]+maxant*spa[1];
147 value[2]=value[2]+maxant*spa[2];
149 this->m_Stack0.push_front(value[0]);
150 this->m_Stack1.push_front(value[1]);
151 this->m_Stack2.push_front(value[2]);
152 this->m_Stack3.push_front(value[0]);
153 this->m_Stack4.push_front(value[1]);
154 this->m_Stack5.push_front(value[2]);
155 this->m_Stack6.push_front(value[0]);
156 this->m_Stack7.push_front(value[1]);
157 this->m_Stack8.push_front(value[2]);
158 this->m_Stack.push_front(0);
159 this->m_Stackr.push_front(0);
160 this->m_Stackra.push_front(0);
161 this->m_Stackrp.push_front(0);
166 //----------------------------------------------------------------------------
167 vtkImageData *axisExtractor02::GetVolumen()
170 vtkImageTranslateExtent *trans;
172 trans = vtkImageTranslateExtent::New();
174 trans->SetInput(this->data4);
176 trans->SetTranslation (maxant, maxant, maxant);
180 return trans->GetOutput();
187 //----------------------------------------------------------------------------
188 vtkPolyData *axisExtractor02::GetOutput ()
194 vtkTransform *transL1;
195 vtkTransformPolyDataFilter *trans;
197 transL1 = vtkTransform::New();
198 trans = vtkTransformPolyDataFilter::New();
202 this->GetInput()->GetSpacing (spa);
203 transL1->Translate(-maxant*spa[0], -maxant*spa[1], -maxant*spa[2]);
204 // transL1->Translate(-maxant, -maxant, -maxant);
207 trans->SetInput((vtkPolyData *)(this->Outputs[0]));
209 trans->SetTransform (transL1);
213 return trans->GetOutput();
222 //----------------------------------------------------------------------------
223 axisExtractor02::axisExtractor02() {
226 this->NumberOfRequiredInputs = 1;
232 /*this->resample= vtkImageResample::New();
233 this->resample->SetDimensionality (3);*/
235 this->data4=vtkImageData::New();
237 this->data1=vtkImageData::New();
238 this->data2=vtkImageData::New();
239 this->data3=vtkImageData::New();
240 this->data6=vtkImageData::New();
242 this->extrac = vtkExtractVOI::New();
243 // this->extrac->SetInput(resample->GetOutput());
244 this->extrac->SetSampleRate(1, 1, 1);
246 this->connect = vtkImageSeedConnectivity::New();
247 this->connect->SetInput(this->data1);
248 this->connect->SetInputConnectValue(255);
249 this->connect->SetOutputConnectedValue(255);
250 this->connect->SetOutputUnconnectedValue(0);
252 this->distance= vtkImageEuclideanDistance::New();
253 this->distance->SetInput(this->connect->GetOutput());
254 this->distance->InitializeOn();
255 this->distance->ConsiderAnisotropyOff();
262 //----------------------------------------------------------------------------
263 void axisExtractor02::SetInput(vtkImageData *input)
266 vtkImageConstantPad *pad;
267 vtkImageTranslateExtent *trans;
269 pad = vtkImageConstantPad::New();
270 trans = vtkImageTranslateExtent::New();
272 pad->SetInput(input);
273 trans->SetInput(pad->GetOutput());
278 pad->SetInput(input);
282 input->GetExtent(ext);
284 ext[0]=ext[0]-maxant;
285 ext[1]=ext[1]+maxant;
286 ext[2]=ext[2]-maxant;
287 ext[3]=ext[3]+maxant;
288 ext[4]=ext[4]-maxant;
289 ext[5]=ext[5]+maxant;
291 pad->SetOutputWholeExtent(ext);
293 trans->SetTranslation (maxant, maxant, maxant);
298 //this->vtkProcessObject::SetNthInput(0, input);
299 this->vtkProcessObject::SetNthInput(0, trans->GetOutput());
301 //this->GetInput()->SetOrigin(0,0,0);
302 //this->resample->SetInput(this->GetInput());
303 this->extrac->SetInput(this->GetInput());
309 //----------------------------------------------------------------------------
310 void axisExtractor02::distanciaejes(vtkPolyData *eje1, vtkPolyData *eje2)
312 vtkIdType totalpuntos1=eje1->GetNumberOfPoints();
313 vtkIdType totalpuntos2=eje2->GetNumberOfPoints();
314 vtkIdType totallineas=eje2->GetNumberOfLines();
320 double point[3], point2[3], point3[3], point4[3];
324 double x, dist, min, y;
327 ff=fopen("comparacion.txt","w");
331 for(i=0;i<totalpuntos1;i++){
332 eje1->GetPoint(i, point);
334 for(j=0;j<totallineas;j++){
335 lineas=eje2->GetCell(j);
337 lineas->GetPoints()->GetPoint(0, point3);
338 lineas->GetPoints()->GetPoint(1, point2);
339 v[0]=point3[0]-point2[0];
340 v[1]=point3[1]-point2[1];
341 v[2]=point3[2]-point2[2];
343 x=((v[0]*point[0])+(v[1]*point[1])+(v[2]*point[2])-(v[0]*point2[0])-(v[1]*point2[1])-(v[2]*point2[2]))/((v[0]*v[0])+(v[1]*v[1])+(v[2]*v[2]));
346 point4[0]=point2[0]+(x*v[0]);
347 point4[1]=point2[1]+(x*v[1]);
348 point4[2]=point2[2]+(x*v[2]);
350 v2[0]=point4[0]-point[0];
351 v2[1]=point4[1]-point[1];
352 v2[2]=point4[2]-point[2];
354 y=(v[0]*v2[0])+(v[1]*v2[1])+(v[2]*v2[2]);
357 dist=sqrt(((point4[0]-point[0])*(point4[0]-point[0]))+((point4[1]-point[1])*(point4[1]-point[1]))+((point4[2]-point[2])*(point4[2]-point[2])));
362 /*printf("punto: %f %f %f\n", point[0], point[1], point[2]);
363 printf("punto2: %f %f %f\n", point2[0], point2[1], point2[2]);
364 printf("punto3: %f %f %f\n", point3[0], point3[1], point3[2]);
365 printf("punto4: %f %f %f\n", point4[0], point4[1], point4[2]);
366 printf("v: %f %f %f\n", v[0], v[1], v[2]);
367 printf("v2: %f %f %f\n", v2[0], v2[1], v2[2]);
380 fprintf(ff,"%f\n", min);
396 //----------------------------------------------------------------------------
397 void axisExtractor02::Execute()
413 this->GetInput()->GetExtent(extprin0);
414 this->GetInput()->GetExtent(extprin);
415 this->GetInput()->GetSpacing(espprin);
418 // stream = fopen( "resultado.txt", "w" );
419 // fprintf(stream, "accion\tindexp\tbuenos\tradiotrabajo\tradioactual\tradiopred\tradioanterior\tmejordst\tmejorrad\tmejorcant\tflag2\tflag3\tflag4\tflag5\tflag6\tflag7\tflag8\tflag9\tflag10\tflag11\tflag12\tflag13\tflag14\tflag15\tflag16\tflag17\tflag18\tflag19\tflag20\tflagg\tflagg2\tcantidad\tcantidadb\tburned\tpuntoactual\tpuntoanterior\tpuntopre\tmejor\n");
421 /*minspac=min(espprin[0],min(espprin[1],espprin[2]));
423 this->resample->SetAxisOutputSpacing( 0, minspac);
424 this->resample->SetAxisOutputSpacing( 1, minspac);
425 this->resample->SetAxisOutputSpacing( 2, minspac);
426 resample->Update();*/
429 this->data4->SetScalarTypeToUnsignedChar();
430 this->data4->SetExtent(this->GetInput()->GetExtent());
431 this->data4->SetSpacing(this->GetInput()->GetSpacing());
433 this->blanquear(this->data4);
435 this->points = vtkPoints::New();
437 this->lineas = vtkCellArray::New();
446 ((vtkPolyData *)this->Outputs[0])->SetPoints (this->points);
447 ((vtkPolyData *)this->Outputs[0])->SetLines(this->lineas);
449 dif = difftime (end,start);
451 // ff=fopen("time.txt","w");
452 // fprintf(ff,"%d \n",this->buenos);
453 // fprintf(ff,"%.2lf \n",dif);
454 // fprintf(ff,"%.2lf \n",(double)this->buenos/dif);
464 //----------------------------------------------------------------------------
465 vtkImageData *axisExtractor02::GetInput()
467 if (this->NumberOfInputs < 1){
471 return (vtkImageData *)(this->Inputs[0]);
480 //----------------------------------------------------------------------------
481 void axisExtractor02::PrintSelf(ostream& os, vtkIndent indent)
483 this->Superclass::PrintSelf(os,indent);
489 //----------------------------------------------------------------------------
490 void axisExtractor02::realtoreal(double a[3], double b[3] )
493 b[0]=a[0]+(espprin[0]*extprin[0]);
494 b[1]=a[1]+(espprin[1]*extprin[2]);
503 //----------------------------------------------------------------------------
504 void axisExtractor02::realtoreal2(double a[3], double b[3] )
507 b[0]=a[0]-(espprin[0]*extprin[0]);
508 b[1]=a[1]-(espprin[1]*extprin[2]);
514 //----------------------------------------------------------------------------
515 void axisExtractor02::realtoindex(double a[3], int b[3] )
521 // EED 15 Mai 2007 .NET
522 // minspac=min(espprin[0],min(espprin[1],espprin[2]));
525 if (espprin[1]<minspac) { minspac=espprin[1]; }
526 if (espprin[2]<minspac) { minspac=espprin[2]; }
533 b[0]=(int)(a[0]/minspac);
534 b[1]=(int)(a[1]/minspac);
535 b[2]=(int)(a[2]/minspac);
556 //----------------------------------------------------------------------------
557 void axisExtractor02::indextoreal(int a[3], double b[3] )
561 // EED 15 Mai 2007 .NET
562 // minspac=min(espprin[0],min(espprin[1],espprin[2]));
564 if (espprin[1]<minspac) { minspac=espprin[1]; }
565 if (espprin[2]<minspac) { minspac=espprin[2]; }
574 //----------------------------------------------------------------------------
575 void axisExtractor02::indextoreal(double a[3], double b[3] )
579 // EED 15 Mai 2007 .NET
580 // minspac=min(espprin[0],min(espprin[1],espprin[2]));
583 if (espprin[1]<minspac) { minspac=espprin[1]; }
584 if (espprin[2]<minspac) { minspac=espprin[2]; }
593 //----------------------------------------------------------------------------
594 double axisExtractor02::distanciaejepunto(double point[3], double point2[3], double point3[3])
600 double v[3], v2[3], v3[3];
605 v[0]=point3[0]-point2[0];
606 v[1]=point3[1]-point2[1];
607 v[2]=point3[2]-point2[2];
609 v3[0]=point[0]-point2[0];
610 v3[1]=point[1]-point2[1];
611 v3[2]=point[2]-point2[2];
613 x=((v[0]*v3[0])+(v[1]*v3[1])+(v[2]*v3[2]))/((v[0]*v[0])+(v[1]*v[1])+(v[2]*v[2]));
616 point4[0]=point2[0]+(x*v[0]);
617 point4[1]=point2[1]+(x*v[1]);
618 point4[2]=point2[2]+(x*v[2]);
620 v2[0]=point4[0]-point[0];
621 v2[1]=point4[1]-point[1];
622 v2[2]=point4[2]-point[2];
624 y=(v[0]*v2[0])+(v[1]*v2[1])+(v[2]*v2[2]);
626 dist=sqrt((v2[0]*v2[0])+(v2[1]*v2[1])+(v2[2]*v2[2]));
634 //----------------------------------------------------------------------------
635 double axisExtractor02::proporcioejepunto(double point[3], double point2[3], double point3[3])
643 v[0]=point3[0]-point2[0];
644 v[1]=point3[1]-point2[1];
645 v[2]=point3[2]-point2[2];
647 v3[0]=point[0]-point2[0];
648 v3[1]=point[1]-point2[1];
649 v3[2]=point[2]-point2[2];
651 x=((v[0]*v3[0])+(v[1]*v3[1])+(v[2]*v3[2]))/((v[0]*v[0])+(v[1]*v[1])+(v[2]*v[2]));
667 //----------------------------------------------------------------------------
668 void axisExtractor02::searc(int i, int j, int k, vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
679 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
680 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
684 data->GetExtent(ext);
686 radio=(ext[1]-ext[0])/2;
693 vector[label-1][0]+=1;
694 vector[label-1][1]+=i;
695 vector[label-1][2]+=j;
696 vector[label-1][3]+=k;
698 for(i3=-1;i3<=1;i3++){
699 for(j3=-1;j3<=1;j3++){
700 for(k3=-1;k3<=1;k3++){
701 if(i+i3>=ext[0] && i+i3<=ext[1] && j+j3>=ext[2] && j+j3<=ext[3] && k+k3>=ext[4] && k+k3<=ext[5] ){
702 ptr=(unsigned char *) data->GetScalarPointer(i+i3,j+j3,k+k3);
703 ptr2=(unsigned char *) data2->GetScalarPointer(i+i3,j+j3,k+k3);
704 if(*ptr==255 && *ptr2==0 && (radio*radio)>=(((i2+i3)*(i2+i3))+((j2+j3)*(j2+j3))+((k2+k3)*(k2+k3))) && ((radio-1)*(radio-1))<=(((i2+i3)*(i2+i3))+((j2+j3)*(j2+j3))+((k2+k3)*(k2+k3)))){
705 searc(i+i3, j+j3, k+k3, data, data2, label, vector);
717 //----------------------------------------------------------------------------
718 void axisExtractor02::searcb(int i, int j, int k, vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
729 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
730 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
734 data->GetExtent(ext);
736 radio=(ext[1]-ext[0])/2;
743 vector[label-1][0]+=1;
744 vector[label-1][1]+=i;
745 vector[label-1][2]+=j;
746 vector[label-1][3]+=k;
748 for(i3=-1;i3<=1;i3++){
749 for(j3=-1;j3<=1;j3++){
750 for(k3=-1;k3<=1;k3++){
751 if(i+i3>=ext[0] && i+i3<=ext[1] && j+j3>=ext[2] && j+j3<=ext[3] && k+k3>=ext[4] && k+k3<=ext[5] ){
752 ptr=(unsigned char *) data->GetScalarPointer(i+i3,j+j3,k+k3);
753 ptr2=(unsigned char *) data2->GetScalarPointer(i+i3,j+j3,k+k3);
754 if(*ptr==0 && *ptr2==0 && (radio*radio)>=(((i2+i3)*(i2+i3))+((j2+j3)*(j2+j3))+((k2+k3)*(k2+k3))) && ((radio-1)*(radio-1))<=(((i2+i3)*(i2+i3))+((j2+j3)*(j2+j3))+((k2+k3)*(k2+k3)))){
755 searcb(i+i3, j+j3, k+k3, data, data2, label, vector);
767 //----------------------------------------------------------------------------
768 unsigned char axisExtractor02::find_components(vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
774 data->GetExtent(ext);
775 double *ff= data->GetOrigin();
781 radio=(ext[1]-ext[0])/2;
783 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
784 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
785 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
786 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
787 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
789 if(*ptr==255 && *ptr2==0 && (radio*radio)>=((i2*i2)+(j2*j2)+(k2*k2)) && ((radio-1)*(radio-1))<=((i2*i2)+(j2*j2)+(k2*k2)) ){
791 vector[label-1][0]=0;
792 vector[label-1][1]=0;
793 vector[label-1][2]=0;
794 vector[label-1][3]=0;
795 searc(i, j, k, data, data2, label, vector);
808 //----------------------------------------------------------------------------
809 unsigned char axisExtractor02::find_componentsb(vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
815 data->GetExtent(ext);
820 radio=(ext[1]-ext[0])/2;
822 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
823 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
824 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
825 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
826 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
828 if(*ptr==0 && *ptr2==0 && (radio*radio)>=((i2*i2)+(j2*j2)+(k2*k2)) && ((radio-1)*(radio-1))<=((i2*i2)+(j2*j2)+(k2*k2)) ){
830 vector[label-1][0]=0;
831 vector[label-1][1]=0;
832 vector[label-1][2]=0;
833 vector[label-1][3]=0;
834 searcb(i, j, k, data, data2, label, vector);
848 //----------------------------------------------------------------------------
849 int axisExtractor02::proporcion(vtkImageData *data )
856 data->GetExtent(ext);
860 for(i=ext[0];i<=ext[1];i++){
861 for(j=ext[2];j<=ext[3];j++){
862 for(k=ext[4];k<=ext[5];k++){
863 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
872 return (max*100)/total;
882 //----------------------------------------------------------------------------
883 bool axisExtractor02::border(vtkImageData *data, int p1[3] )
894 data->GetExtent(ext);
896 centro[0]=(ext[1]+ext[0])/2;
897 centro[1]=(ext[3]+ext[2])/2;
898 centro[2]=(ext[5]+ext[4])/2;
902 ptr=(unsigned char *) data->GetScalarPointer(p1);
916 if(ext[0]<=p2[0] && ext[1]>=p2[0] && ext[2]<=p2[1] && ext[3]>=p2[1] && ext[4]<=p2[2] && ext[5]>=p2[2] ){
917 ptr=(unsigned char *) data->GetScalarPointer(p2);
932 /*if(p1[0]==centro[0] && p1[1]==centro[1] && p1[2]==centro[2]){
944 //----------------------------------------------------------------------------
945 void axisExtractor02::optim(vtkImageData *data, vtkImageData *data2 )
955 double inerciaxx, inerciayy, inerciazz, inerciaxy, inerciayz, inerciazx;
956 double sumx, sumy, sumz;
957 double sumxx, sumyy, sumzz, sumxy, sumyz, sumzx;
959 double sumix, sumiy, sumiz;
960 double sumixx, sumiyy, sumizz, sumixy, sumiyz, sumizx;
962 double inerciaixx, inerciaiyy, inerciaizz, inerciaixy, inerciaiyz, inerciaizx;
964 double inerciaxx2, inerciayy2, inerciazz2, inerciaxy2, inerciayz2, inerciazx2;
965 double sumx2, sumy2, sumz2;
966 double sumxx2, sumyy2, sumzz2, sumxy2, sumyz2, sumzx2;
968 double sumix2, sumiy2, sumiz2;
969 double sumixx2, sumiyy2, sumizz2, sumixy2, sumiyz2, sumizx2;
971 //double sumixp, sumiyp, sumizp;
972 //double sumix2p, sumiy2p, sumiz2p;
974 double inerciaxxp, inerciayyp, inerciazzp, inerciaxyp, inerciayzp, inerciazxp;
975 double sumxp, sumyp, sumzp;
976 double sumxxp, sumyyp, sumzzp, sumxyp, sumyzp, sumzxp;
978 double inerciaxx2p, inerciayy2p, inerciazz2p, inerciaxy2p, inerciayz2p, inerciazx2p;
979 double sumx2p, sumy2p, sumz2p;
980 double sumxx2p, sumyy2p, sumzz2p, sumxy2p, sumyz2p, sumzx2p;
982 //double inerciaixx2p, inerciaiyy2p, inerciaizz2p, inerciaixy2p, inerciaiyz2p, inerciaizx2p;
990 unsigned long cant2p;
1002 double im, jm, km, lm;
1073 data->GetExtent(ext);
1075 unsigned short *ptr;
1076 unsigned char *ptr2;
1078 radio=(ext[1]-ext[0])/2;
1080 centro[0]=(ext[1]+ext[0])/2;
1081 centro[1]=(ext[3]+ext[2])/2;
1082 centro[2]=(ext[5]+ext[4])/2;
1086 for(i=ext[0];i<=ext[1];i++){
1087 for(j=ext[2];j<=ext[3];j++){
1088 for(k=ext[4];k<=ext[5];k++){
1089 tmpsqrt=((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
1090 if( radio>=sqrt( tmpsqrt) ){
1091 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
1105 for(i=ext[0];i<=ext[1];i++){
1106 for(j=ext[2];j<=ext[3];j++){
1107 for(k=ext[4];k<=ext[5];k++){
1108 tmpsqrt=((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
1109 if( radio>=sqrt(tmpsqrt) ){
1110 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
1111 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
1116 im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
1117 jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
1118 km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
1119 lm=2*(((double)*ptr-(double)min)/((double)max-(double)min));
1180 centi=(double)sumi/(double)cant;
1181 inerciaii= ((double)sumii/(double)cant)-(centi*centi);
1183 centi2=(double)sumi2/(double)cant2;
1184 inerciaii2= ((double)sumii2/(double)cant2)-(centi2*centi2);
1186 centx=(double)sumx/(double)cant;
1187 centy=(double)sumy/(double)cant;
1188 centz=(double)sumz/(double)cant;
1190 centx2=(double)sumx2/(double)cant2;
1191 centy2=(double)sumy2/(double)cant2;
1192 centz2=(double)sumz2/(double)cant2;
1194 inerciaxx= ((double)sumxx/(double)cant)-centx*centx;
1195 inerciayy= ((double)sumyy/(double)cant)-centy*centy;
1196 inerciazz= ((double)sumzz/(double)cant)-centz*centz;
1198 inerciaxx2= ((double)sumxx2/(double)cant2)-centx2*centx2;
1199 inerciayy2= ((double)sumyy2/(double)cant2)-centy2*centy2;
1200 inerciazz2= ((double)sumzz2/(double)cant2)-centz2*centz2;
1202 inerciaxy= ((double)sumxy/(double)cant)-centx*centy;
1203 inerciayz= ((double)sumyz/(double)cant)-centy*centz;
1204 inerciazx= ((double)sumzx/(double)cant)-centz*centx;
1206 inerciaxy2= ((double)sumxy2/(double)cant2)-centx2*centy2;
1207 inerciayz2= ((double)sumyz2/(double)cant2)-centy2*centz2;
1208 inerciazx2= ((double)sumzx2/(double)cant2)-centz2*centx2;
1211 centix=((double)sumix+(double)sumix2)/((double)sumi+(double)sumi2);
1212 centiy=((double)sumiy+(double)sumiy2)/((double)sumi+(double)sumi2);
1213 centiz=((double)sumiz+(double)sumiz2)/((double)sumi+(double)sumi2);
1215 inerciaixx= ( ((double)sumixx+(double)sumixx2) / ((double)sumi+(double)sumi2) )-(centix*centix);
1216 inerciaiyy= ( ((double)sumiyy+(double)sumiyy2) / ((double)sumi+(double)sumi2) )-(centiy*centiy);
1217 inerciaizz= ( ((double)sumizz+(double)sumizz2) / ((double)sumi+(double)sumi2) )-(centiz*centiz);
1220 inerciaixy= ( ((double)sumixy+(double)sumixy2) / ((double)sumi+(double)sumi2) )-(centix*centiy);
1221 inerciaiyz= ( ((double)sumiyz+(double)sumiyz2) / ((double)sumi+(double)sumi2) )-(centiy*centiz);
1222 inerciaizx= ( ((double)sumizx+(double)sumizx2) / ((double)sumi+(double)sumi2) )-(centiz*centix);
1228 A[0][1]=A[1][0]=inerciaxy;
1229 A[0][2]=A[2][0]=inerciazx;
1230 A[1][2]=A[2][1]=inerciayz;
1232 A2[0][0]=inerciaxx2;
1233 A2[1][1]=inerciayy2;
1234 A2[2][2]=inerciazz2;
1235 A2[0][1]=A2[1][0]=inerciaxy2;
1236 A2[0][2]=A2[2][0]=inerciazx2;
1237 A2[1][2]=A2[2][1]=inerciayz2;
1239 Ai[0][0]=inerciaixx;
1240 Ai[1][1]=inerciaiyy;
1241 Ai[2][2]=inerciaizz;
1242 Ai[0][1]=Ai[1][0]=inerciaixy;
1243 Ai[0][2]=Ai[2][0]=inerciaizx;
1244 Ai[1][2]=Ai[2][1]=inerciaiyz;
1246 vtkMath::Diagonalize3x3 (A, w, V);
1247 vtkMath::Diagonalize3x3 (A2, w2, V2);
1248 vtkMath::Diagonalize3x3 (Ai, wi, Vi);
1256 inerciary=w[0]+w[2];
1257 inerciarz=w[0]+w[1];
1260 inerciary=w[0]+w[1];
1261 inerciarz=w[0]+w[2];
1268 inerciary=w[0]+w[2];
1269 inerciarz=w[1]+w[2];
1272 inerciary=w[1]+w[2];
1273 inerciarz=w[0]+w[2];
1282 inerciary=w[1]+w[0];
1283 inerciarz=w[1]+w[2];
1286 inerciary=w[1]+w[2];
1287 inerciarz=w[1]+w[0];
1294 inerciary=w[0]+w[2];
1295 inerciarz=w[1]+w[2];
1298 inerciary=w[1]+w[2];
1299 inerciarz=w[0]+w[2];
1308 inercia2r=w2[1]+w2[2];
1312 inercia2r=w2[1]+w2[0];
1318 inercia2r=w2[2]+w2[0];
1322 inercia2r=w2[1]+w2[0];
1330 inerciari=wi[1]+wi[2];
1332 inerciariy=wi[0]+wi[2];
1333 inerciariz=wi[0]+wi[1];
1336 inerciariy=wi[0]+wi[1];
1337 inerciariz=wi[0]+wi[2];
1342 inerciari=wi[1]+wi[0];
1344 inerciariy=wi[0]+wi[2];
1345 inerciariz=wi[1]+wi[2];
1348 inerciariy=wi[1]+wi[2];
1349 inerciariz=wi[0]+wi[2];
1357 inerciari=wi[2]+wi[0];
1359 inerciariy=wi[1]+wi[0];
1360 inerciariz=wi[1]+wi[2];
1363 inerciariy=wi[1]+wi[2];
1364 inerciariz=wi[1]+wi[0];
1369 inerciari=wi[1]+wi[0];
1371 inerciariy=wi[0]+wi[2];
1372 inerciariz=wi[1]+wi[2];
1375 inerciariy=wi[1]+wi[2];
1376 inerciariz=wi[0]+wi[2];
1383 for(i=ext[0];i<=ext[1];i++){
1384 for(j=ext[2];j<=ext[3];j++){
1385 for(k=ext[4];k<=ext[5];k++){
1386 if( radio>=sqrt(((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]))) ){
1387 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
1388 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
1391 im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
1392 jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
1393 km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
1399 point2[0]=centx+V[0][ejemin];
1400 point2[1]=centy+V[1][ejemin];
1401 point2[2]=centz+V[2][ejemin];
1403 point3[0]=centx-V[0][ejemin];
1404 point3[1]=centy-V[1][ejemin];
1405 point3[2]=centz-V[2][ejemin];
1407 //distcent=1+sqrt(centx*centx+centy*centy+centz*centz);
1412 sumitotal+=(distcent-fabs(distanciaejepunto(point1, point2, point3)))*(distcent-fabs(distanciaejepunto(point1, point2, point3)));
1420 inercia2r=sumitotal/cant2;
1426 /*printf("centx: %f\n", centx);
1427 printf("centy: %f\n", centy);
1428 printf("centz: %f\n", centz);
1430 printf("w[0]: %f\n", w[0]);
1431 printf("w[1]: %f\n", w[1]);
1432 printf("w[2]: %f\n", w[2]);
1433 printf("inerciar: %f\n", inerciar);
1435 printf("w[ejemin]: %f\n", w[ejemin]);
1437 printf("w[ejemin]/inerciar: %f\n", w[ejemin]/inerciar);
1439 printf("V[0][0]: %f\n", V[0][0]);
1440 printf("V[1][0]: %f\n", V[1][0]);
1441 printf("V[2][0]: %f\n", V[2][0]);
1443 printf("V[0][1]: %f\n", V[0][1]);
1444 printf("V[1][1]: %f\n", V[1][1]);
1445 printf("V[2][1]: %f\n", V[2][1]);
1447 printf("V[0][2]: %f\n", V[0][2]);
1448 printf("V[1][2]: %f\n", V[1][2]);
1449 printf("V[2][2]: %f\n", V[2][2]);*/
1451 w0=(double)cant/((double)cant+(double)cant2);
1452 w1=(double)cant2/((double)cant+(double)cant2);
1456 ut=w0*centi+w1*centi2;
1457 intervar=w0*inerciaii+w1*inerciaii2;
1465 //costo=w0*inerciaii+w1*inerciaii2+param*inerciar;
1466 //costo=param*(w0*param3*inerciaii+w1*(1-param3)*inerciaii2)+(1-param)*(w0*param2*inerciar+w1*(1-param2)*inercia2r);
1467 costo=param*w0*inerciaii+w1*param2*inerciaii2+w0*param3*inerciar+w1*param4*inercia2r;
1473 while(changes>3 && total<500){
1476 for(i=ext[0];i<=ext[1];i++){
1477 for(j=ext[2];j<=ext[3];j++){
1478 for(k=ext[4];k<=ext[5];k++){
1479 tmpsqrt = ((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
1480 if( radio>=sqrt(tmpsqrt) ){
1481 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
1482 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
1487 im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
1488 jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
1489 km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
1490 lm=2*(((double)*ptr-(double)min)/((double)max-(double)min));
1492 if(border(data2, point)){
1501 sumii2p=sumii2+lm*lm;
1518 sumxx2p=sumxx2+im*im;
1519 sumyy2p=sumyy2+jm*jm;
1520 sumzz2p=sumzz2+km*km;
1521 sumxy2p=sumxy2+im*jm;
1522 sumyz2p=sumyz2+jm*km;
1523 sumzx2p=sumzx2+km*im;
1525 centip=(double)sumip/(double)cantp;
1526 inerciaiip= ((double)sumiip/(double)cantp)-(centip*centip);
1528 centi2p=(double)sumi2p/(double)cant2p;
1529 inerciaii2p= ((double)sumii2p/(double)cant2p)-(centi2p*centi2p);
1531 centxp=(double)sumxp/(double)cantp;
1532 centyp=(double)sumyp/(double)cantp;
1533 centzp=(double)sumzp/(double)cantp;
1535 centx2p=(double)sumx2p/(double)cant2p;
1536 centy2p=(double)sumy2p/(double)cant2p;
1537 centz2p=(double)sumz2p/(double)cant2p;
1539 inerciaxxp= ((double)sumxxp/(double)cantp)-centxp*centxp;
1540 inerciayyp= ((double)sumyyp/(double)cantp)-centyp*centyp;
1541 inerciazzp= ((double)sumzzp/(double)cantp)-centzp*centzp;
1543 inerciaxx2p= ((double)sumxx2p/(double)cant2p)-centx2p*centx2p;
1544 inerciayy2p= ((double)sumyy2p/(double)cant2p)-centy2p*centy2p;
1545 inerciazz2p= ((double)sumzz2p/(double)cant2p)-centz2p*centz2p;
1547 inerciaxyp= ((double)sumxyp/(double)cantp)-centxp*centyp;
1548 inerciayzp= ((double)sumyzp/(double)cantp)-centyp*centzp;
1549 inerciazxp= ((double)sumzxp/(double)cantp)-centzp*centxp;
1551 inerciaxy2p= ((double)sumxy2p/(double)cant2p)-centx2p*centy2p;
1552 inerciayz2p= ((double)sumyz2p/(double)cant2p)-centy2p*centz2p;
1553 inerciazx2p= ((double)sumzx2p/(double)cant2p)-centz2p*centx2p;
1555 Ap[0][0]=inerciaxxp;
1556 Ap[1][1]=inerciayyp;
1557 Ap[2][2]=inerciazzp;
1558 Ap[0][1]=Ap[1][0]=inerciaxyp;
1559 Ap[0][2]=Ap[2][0]=inerciazxp;
1560 Ap[1][2]=Ap[2][1]=inerciayzp;
1563 A2p[0][0]=inerciaxx2p;
1564 A2p[1][1]=inerciayy2p;
1565 A2p[2][2]=inerciazz2p;
1566 A2p[0][1]=A2p[1][0]=inerciaxy2p;
1567 A2p[0][2]=A2p[2][0]=inerciazx2p;
1568 A2p[1][2]=A2p[2][1]=inerciayz2p;
1572 w0p=(double)cantp/((double)cantp+(double)cant2p);
1573 w1p=(double)cant2p/((double)cantp+(double)cant2p);
1578 vtkMath::Diagonalize3x3 (Ap, wp, Vp);
1579 vtkMath::Diagonalize3x3 (A2p, w2p, V2p);
1585 inerciarp=wp[2]+wp[1];
1587 inerciarpy=wp[0]+wp[2];
1588 inerciarpz=wp[0]+wp[1];
1591 inerciarpy=wp[0]+wp[1];
1592 inerciarpz=wp[0]+wp[2];
1597 inerciarp=wp[1]+wp[0];
1599 inerciarpy=wp[0]+wp[2];
1600 inerciarpz=wp[1]+wp[2];
1603 inerciarpy=wp[1]+wp[2];
1604 inerciarpz=wp[0]+wp[2];
1611 inerciarp=wp[2]+wp[0];
1613 inerciarpy=wp[1]+wp[0];
1614 inerciarpz=wp[1]+wp[2];
1617 inerciarpy=wp[1]+wp[2];
1618 inerciarpz=wp[1]+wp[0];
1623 inerciarp=wp[1]+wp[0];
1625 inerciarpy=wp[0]+wp[2];
1626 inerciarpz=wp[1]+wp[2];
1629 inerciarpy=wp[1]+wp[2];
1630 inerciarpz=wp[0]+wp[2];
1639 inercia2rp=w2p[2]+w2p[1];
1643 inercia2rp=w2p[1]+w2p[0];
1649 inercia2rp=w2p[2]+w2p[0];
1653 inercia2rp=w2p[1]+w2p[0];
1661 point2[0]=centxp+Vp[0][ejeminp];
1662 point2[1]=centyp+Vp[1][ejeminp];
1663 point2[2]=centzp+Vp[2][ejeminp];
1665 point3[0]=centxp-Vp[0][ejeminp];
1666 point3[1]=centyp-Vp[1][ejeminp];
1667 point3[2]=centzp-Vp[2][ejeminp];*/
1669 //distcent=1+sqrt(centx*centx+centy*centy+centz*centz);
1672 sumitotalp=sumitotal+(distcent-fabs(distanciaejepunto(point1, point2, point3)))*(distcent-fabs(distanciaejepunto(point1, point2, point3)));
1674 inercia2rp=sumitotalp/cant2p;*/
1676 //costop=w0p*inerciaiip+w1p*inerciaii2p+param*inerciarp;
1678 //costop=param*(w0p*param3*inerciaiip+w1p*(1-param3)*inerciaii2p)+(1-param)*(w0p*param2*inerciarp+w1p*(1-param2)*inercia2rp);
1681 costop=param*w0p*inerciaiip+w1p*param2*inerciaii2p+w0p*param3*inerciarp+w1p*param4*inercia2rp;
1683 if(costop<costo && centip>centi2p){
1717 inerciaii= inerciaiip;
1720 inerciaii2p= inerciaii2p;
1730 inerciaxx= inerciaxxp;
1731 inerciayy= inerciayyp;
1732 inerciazz= inerciazzp;
1734 inerciaxx2= inerciaxx2p;
1735 inerciayy2= inerciayy2p;
1736 inerciazz2= inerciazz2p;
1738 inerciaxy= inerciaxyp;
1739 inerciayz= inerciayzp;
1740 inerciazx= inerciazxp;
1742 inerciaxy2= inerciaxy2p;
1743 inerciayz2= inerciayz2p;
1744 inerciazx2= inerciazx2p;
1746 // sumitotal=sumitotalp;
1759 sumii2p=sumii2-lm*lm;
1776 sumxx2p=sumxx2-im*im;
1777 sumyy2p=sumyy2-jm*jm;
1778 sumzz2p=sumzz2-km*km;
1779 sumxy2p=sumxy2-im*jm;
1780 sumyz2p=sumyz2-jm*km;
1781 sumzx2p=sumzx2-km*im;
1783 centip=(double)sumip/(double)cantp;
1784 inerciaiip= ((double)sumiip/(double)cantp)-(centip*centip);
1786 centi2p=(double)sumi2p/(double)cant2p;
1787 inerciaii2p= ((double)sumii2p/(double)cant2p)-(centi2p*centi2p);
1789 centxp=(double)sumxp/(double)cantp;
1790 centyp=(double)sumyp/(double)cantp;
1791 centzp=(double)sumzp/(double)cantp;
1793 centx2p=(double)sumx2p/(double)cant2p;
1794 centy2p=(double)sumy2p/(double)cant2p;
1795 centz2p=(double)sumz2p/(double)cant2p;
1798 inerciaxxp= ((double)sumxxp/(double)cantp)-centxp*centxp;
1799 inerciayyp= ((double)sumyyp/(double)cantp)-centyp*centyp;
1800 inerciazzp= ((double)sumzzp/(double)cantp)-centzp*centzp;
1802 inerciaxx2p= ((double)sumxx2p/(double)cant2p)-centx2p*centx2p;
1803 inerciayy2p= ((double)sumyy2p/(double)cant2p)-centy2p*centy2p;
1804 inerciazz2p= ((double)sumzz2p/(double)cant2p)-centz2p*centz2p;
1806 inerciaxyp= ((double)sumxyp/(double)cantp)-centxp*centyp;
1807 inerciayzp= ((double)sumyzp/(double)cantp)-centyp*centzp;
1808 inerciazxp= ((double)sumzxp/(double)cantp)-centzp*centxp;
1810 inerciaxy2p= ((double)sumxy2p/(double)cant2p)-centx2p*centy2p;
1811 inerciayz2p= ((double)sumyz2p/(double)cant2p)-centy2p*centz2p;
1812 inerciazx2p= ((double)sumzx2p/(double)cant2p)-centz2p*centx2p;
1814 Ap[0][0]=inerciaxxp;
1815 Ap[1][1]=inerciayyp;
1816 Ap[2][2]=inerciazzp;
1817 Ap[0][1]=Ap[1][0]=inerciaxyp;
1818 Ap[0][2]=Ap[2][0]=inerciazxp;
1819 Ap[1][2]=Ap[2][1]=inerciayzp;
1821 A2p[0][0]=inerciaxx2p;
1822 A2p[1][1]=inerciayy2p;
1823 A2p[2][2]=inerciazz2p;
1824 A2p[0][1]=A2p[1][0]=inerciaxy2p;
1825 A2p[0][2]=A2p[2][0]=inerciazx2p;
1826 A2p[1][2]=A2p[2][1]=inerciayz2p;
1828 w0p=(double)cantp/((double)cantp+(double)cant2p);
1829 w1p=(double)cant2p/((double)cantp+(double)cant2p);
1834 vtkMath::Diagonalize3x3 (Ap, wp, Vp);
1835 vtkMath::Diagonalize3x3 (A2p, w2p, V2p);
1840 inerciarp=wp[2]+wp[1];
1842 inerciarpy=wp[0]+wp[2];
1843 inerciarpz=wp[0]+wp[1];
1846 inerciarpy=wp[0]+wp[1];
1847 inerciarpz=wp[0]+wp[2];
1852 inerciarp=wp[1]+wp[0];
1854 inerciarpy=wp[0]+wp[2];
1855 inerciarpz=wp[1]+wp[2];
1858 inerciarpy=wp[1]+wp[2];
1859 inerciarpz=wp[0]+wp[2];
1866 inerciarp=wp[2]+wp[0];
1868 inerciarpy=wp[1]+wp[0];
1869 inerciarpz=wp[1]+wp[2];
1872 inerciarpy=wp[1]+wp[2];
1873 inerciarpz=wp[1]+wp[0];
1878 inerciarp=wp[1]+wp[0];
1880 inerciarpy=wp[0]+wp[2];
1881 inerciarpz=wp[1]+wp[2];
1884 inerciarpy=wp[1]+wp[2];
1885 inerciarpz=wp[0]+wp[2];
1894 inercia2rp=w2p[2]+w2p[1];
1898 inercia2rp=w2p[1]+w2p[0];
1904 inercia2rp=w2p[2]+w2p[0];
1908 inercia2rp=w2p[1]+w2p[0];
1917 point2[0]=centxp+Vp[0][ejeminp];
1918 point2[1]=centyp+Vp[1][ejeminp];
1919 point2[2]=centzp+Vp[2][ejeminp];
1921 point3[0]=centxp-Vp[0][ejeminp];
1922 point3[1]=centyp-Vp[1][ejeminp];
1923 point3[2]=centzp-Vp[2][ejeminp];*/
1925 //distcent=1+sqrt(centx*centx+centy*centy+centz*centz);
1928 sumitotalp=sumitotal-(distcent-fabs(distanciaejepunto(point1, point2, point3)))*(distcent-fabs(distanciaejepunto(point1, point2, point3)));
1930 inercia2rp=sumitotalp/cant2p;*/
1933 //costop=w0p*inerciaiip+w1p*inerciaii2p+param*inerciarp;
1934 //costop=param*(w0p*param3*inerciaiip+w1p*(1-param3)*inerciaii2p)+(1-param)*(w0p*param2*inerciarp+w1p*(1-param2)*inercia2rp);
1937 costop=param*w0p*inerciaiip+w1p*param2*inerciaii2p+w0p*param3*inerciarp+w1p*param4*inercia2rp;
1939 if(costop<costo && centip>centi2p){
1973 inerciaii= inerciaiip;
1976 inerciaii2p= inerciaii2p;
1986 inerciaxx= inerciaxxp;
1987 inerciayy= inerciayyp;
1988 inerciazz= inerciazzp;
1990 inerciaxx2= inerciaxx2p;
1991 inerciayy2= inerciayy2p;
1992 inerciazz2= inerciazz2p;
1994 inerciaxy= inerciaxyp;
1995 inerciayz= inerciayzp;
1996 inerciazx= inerciazxp;
1998 inerciaxy2= inerciaxy2p;
1999 inerciayz2= inerciayz2p;
2000 inerciazx2= inerciazx2p;
2002 // sumitotal=sumitotalp;
2015 // printf("changes: %d\n", changes);
2020 //printf("total: %d\n", total);
2029 //----------------------------------------------------------------------------
2030 void axisExtractor02::costominimo(vtkImageData *data, vtkImageData *data2 )
2037 data->GetExtent(ext);
2040 unsigned char *ptr2;
2044 radio=(ext[1]-ext[0])/2;
2046 centro[0]=(ext[1]+ext[0])/2;
2047 centro[1]=(ext[3]+ext[2])/2;
2048 centro[2]=(ext[5]+ext[4])/2;
2055 for(i=ext[0];i<=ext[1];i++){
2056 for(j=ext[2];j<=ext[3];j++){
2057 for(k=ext[4];k<=ext[5];k++){
2058 ptr=(double *)data->GetScalarPointer(i,j,k);
2059 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
2064 if(minis[ind-1]<=*ptr && *ptr!=0 && ind<11){
2080 //----------------------------------------------------------------------------
2081 void axisExtractor02::costominimo2(vtkImageData *data, vtkImageData *data3, int p1[3], int p2[3], int p3[3])
2094 unsigned char *ptr3;
2095 unsigned char *ptr4;
2097 unsigned char *ptr6;
2101 vtkImageData *data4;
2103 data4=vtkImageData::New();
2104 data4->SetScalarTypeToUnsignedChar();
2105 data4->SetExtent(data->GetExtent());
2106 data4->SetSpacing(data->GetSpacing());
2108 vtkImageData *data6;
2110 data6=vtkImageData::New();
2111 data6->SetScalarTypeToUnsignedChar();
2112 data6->SetExtent(data->GetExtent());
2113 data6->SetSpacing(data->GetSpacing());
2116 vtkImageData *data2;
2118 data2=vtkImageData::New();
2119 data2->SetScalarTypeToDouble();
2120 data2->SetExtent(data->GetExtent());
2121 data2->SetSpacing(data->GetSpacing());
2123 vtkImageData *data8;
2125 data8=vtkImageData::New();
2126 data8->SetScalarTypeToDouble();
2127 data8->SetExtent(data->GetExtent());
2128 data8->SetSpacing(data->GetSpacing());
2131 double mincst, mincstb;
2137 data->GetExtent(ext);
2139 radio=(ext[1]-ext[0])/2;
2141 for(i=ext[0];i<=ext[1];i++){
2142 for(j=ext[2];j<=ext[3];j++){
2143 for(k=ext[4];k<=ext[5];k++){
2144 ptr2=(double *)data2->GetScalarPointer(i,j,k);
2145 ptr4=(unsigned char *)data4->GetScalarPointer(i,j,k);
2146 ptr6=(unsigned char *)data6->GetScalarPointer(i,j,k);
2147 ptr8=(double *)data8->GetScalarPointer(i,j,k);
2167 ptr2=(double *)data2->GetScalarPointer(p1[0],p1[1],p1[2]);
2170 ptr4=(unsigned char *)data4->GetScalarPointer(p1[0],p1[1],p1[2]);
2173 ptr8=(double *)data8->GetScalarPointer(p2[0],p2[1],p2[2]);
2176 ptr6=(unsigned char *)data6->GetScalarPointer(p2[0],p2[1],p2[2]);
2186 ptr5=(double *)data2->GetScalarPointer(punto[0],punto[1],punto[2]);
2187 ptr7=(double *)data8->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
2192 if(i!=0 || j!=0 || k!=0){
2193 punto2[0]=punto[0]+i;
2194 punto2[1]=punto[1]+j;
2195 punto2[2]=punto[2]+k;
2197 puntob2[0]=puntob[0]+i;
2198 puntob2[1]=puntob[1]+j;
2199 puntob2[2]=puntob[2]+k;
2201 if(punto2[0]>=ext[0] && punto2[0]<=ext[1] && punto2[1]>=ext[2] && punto2[1]<=ext[3] && punto2[2]>=ext[4] && punto2[2]<=ext[5] ){
2203 ptr=(double *)data->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
2204 ptr2=(double *)data2->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
2205 ptr3=(unsigned char *)data3->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
2206 ptr4=(unsigned char *)data4->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
2209 if(*ptr3!=0 && *ptr4!=2){
2210 costoactual=*ptr5+(1/(*ptr));
2211 if(*ptr2>costoactual){
2218 if(puntob2[0]>=ext[0] && puntob2[0]<=ext[1] && puntob2[1]>=ext[2] && puntob2[1]<=ext[3] && puntob2[2]>=ext[4] && puntob2[2]<=ext[5] ){
2220 ptr=(double *)data->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
2221 ptr8=(double *)data8->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
2222 ptr3=(unsigned char *)data3->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
2223 ptr6=(unsigned char *)data6->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
2225 if(*ptr3!=0 && *ptr6!=2){
2226 costoactual=*ptr7+(1/(*ptr));
2227 if(*ptr8>costoactual){
2238 ptr4=(unsigned char *)data4->GetScalarPointer(punto[0],punto[1],punto[2]);
2241 ptr6=(unsigned char *)data6->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
2247 for(i=ext[0];i<=ext[1];i++){
2248 for(j=ext[2];j<=ext[3];j++){
2249 for(k=ext[4];k<=ext[5];k++){
2250 ptr2=(double *)data2->GetScalarPointer(i,j,k);
2251 ptr8=(double *)data8->GetScalarPointer(i,j,k);
2252 ptr4=(unsigned char *)data4->GetScalarPointer(i,j,k);
2253 ptr6=(unsigned char *)data6->GetScalarPointer(i,j,k);
2276 punto[0]=puntomin[0];
2277 punto[1]=puntomin[1];
2278 punto[2]=puntomin[2];
2280 puntob[0]=puntobmin[0];
2281 puntob[1]=puntobmin[1];
2282 puntob[2]=puntobmin[2];
2284 ptr4=(unsigned char *)data4->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
2285 ptr6=(unsigned char *)data6->GetScalarPointer(punto[0],punto[1],punto[2]);
2313 //----------------------------------------------------------------------------
2314 void axisExtractor02::invertir(vtkImageData *data )
2320 data->GetExtent(ext);
2324 radio=(ext[1]-ext[0])/2;
2328 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2329 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2330 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2331 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2332 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
2333 if( radio<sqrt(tmpsqrt) ){
2348 //----------------------------------------------------------------------------
2349 void axisExtractor02::redondear(vtkImageData *data )
2355 data->GetExtent(ext);
2359 radio=(ext[1]-ext[0])/2;
2363 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2364 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2365 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2366 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2367 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
2368 if( radio<sqrt(tmpsqrt) ){
2371 /*if( i2==0 && j2==0 && k2==0 ){
2383 //----------------------------------------------------------------------------
2384 void axisExtractor02::redondear2(vtkImageData *data )
2390 data->GetExtent(ext);
2394 radio=(ext[1]-ext[0])/2;
2398 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2399 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2400 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2401 ptr=(double *) data->GetScalarPointer(i,j,k);
2402 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
2403 if( radio<sqrt(tmpsqrt) ){
2416 //----------------------------------------------------------------------------
2417 void axisExtractor02::redondear3(vtkImageData *data )
2423 data->GetExtent(ext);
2427 radio=(ext[1]-ext[0])/2;
2431 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2432 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2433 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2434 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2435 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
2436 if( radio<sqrt(tmpsqrt) ){
2449 //----------------------------------------------------------------------------
2450 double axisExtractor02::distancia(double a[3], double b[3] )
2452 double tmpsqrt=((a[0]-b[0])*(a[0]-b[0]))+((a[1]-b[1])*(a[1]-b[1]))+((a[2]-b[2])*(a[2]-b[2]));
2453 return sqrt(tmpsqrt);
2460 //----------------------------------------------------------------------------
2461 double axisExtractor02::distancia(int a[3], int b[3] )
2463 double tmpsqrt=((a[0]-b[0])*(a[0]-b[0]))+((a[1]-b[1])*(a[1]-b[1]))+((a[2]-b[2])*(a[2]-b[2]));
2464 return sqrt(tmpsqrt);
2473 //----------------------------------------------------------------------------
2474 void axisExtractor02::blanquear(vtkImageData *data )
2479 data->GetExtent(ext);
2483 for(i=ext[0];i<=ext[1];i++){
2484 for(j=ext[2];j<=ext[3];j++){
2485 for(k=ext[4];k<=ext[5];k++){
2486 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2495 //----------------------------------------------------------------------------
2496 void axisExtractor02::blanquear3(vtkImageData *data )
2501 data->GetExtent(ext);
2505 for(i=ext[0];i<=ext[1];i++){
2506 for(j=ext[2];j<=ext[3];j++){
2507 for(k=ext[4];k<=ext[5];k++){
2508 ptr=(double *) data->GetScalarPointer(i,j,k);
2518 //----------------------------------------------------------------------------
2519 void axisExtractor02::blanquear2(vtkImageData *data )
2525 data->GetExtent(ext);
2527 centro[0]=(ext[1]+ext[0])/2;
2528 centro[1]=(ext[3]+ext[2])/2;
2529 centro[2]=(ext[5]+ext[4])/2;
2533 for(i=ext[0];i<=ext[1];i++){
2534 for(j=ext[2];j<=ext[3];j++){
2535 for(k=ext[4];k<=ext[5];k++){
2536 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2537 if(i==centro[0] && j==centro[1] && k==centro[2]){
2551 //----------------------------------------------------------------------------
2552 void axisExtractor02::cilindro(vtkImageData *data, double vector[3] )
2565 data->GetExtent(ext);
2569 radio=(ext[1]-ext[0])/2;
2570 radiof=espprin[0]*(((double)ext[1]-(double)ext[0])/4);
2571 centro[0]=(ext[1]+ext[0])/2;
2572 centro[1]=(ext[3]+ext[2])/2;
2573 centro[2]=(ext[5]+ext[4])/2;
2575 indextoreal(centro, centrof );
2577 if(vector[0]==0 && vector[1]==0 && vector[2]==0){
2581 puntof2[0]=centrof[0]+vector[0];
2582 puntof2[1]=centrof[1]+vector[1];
2583 puntof2[2]=centrof[2]+vector[2];
2586 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2587 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2588 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2589 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2593 indextoreal(punto, puntof );
2594 tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2);
2595 if( radio<sqrt(tmpsqrt) ){
2599 else if(radiof<distanciaejepunto(puntof, centrof, puntof2)){
2615 //----------------------------------------------------------------------------
2616 void axisExtractor02::modelo(vtkImageData *data, unsigned char cantidad, unsigned long vector[50][4], int candit[10][3], double radioactual, double minis[10])
2619 int i, j, k, radio, radio2;
2630 double dist, prop, rest;
2633 data->GetExtent(ext);
2637 radio=(ext[1]-ext[0])/2;
2638 radio2=(int)radioactual;
2640 centro[0]=(ext[1]+ext[0])/2;
2641 centro[1]=(ext[3]+ext[2])/2;
2642 centro[2]=(ext[5]+ext[4])/2;
2644 indextoreal(centro, centrof );
2648 for(i=ext[0];i<=ext[1];i++){
2649 for(j=ext[2];j<=ext[3];j++){
2650 for(k=ext[4];k<=ext[5];k++){
2651 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2656 if( radioactual<distancia(punto, centro) ){
2667 for(t=0;t<cantidad;t++){
2668 radiof=sqrt(minis[t]);
2670 punto2[0]=candit[t][0];
2671 punto2[1]=candit[t][1];
2672 punto2[2]=candit[t][2];
2674 indextoreal(punto2, puntof2 );
2677 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2678 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2679 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2680 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2684 indextoreal(punto, puntof );
2686 dist=distanciaejepunto(puntof, centrof, puntof2);
2687 prop=proporcioejepunto(puntof, centrof, puntof2);
2689 if(prop>=0 && prop<=1){
2690 rest=(radiof-radioactual)*prop+radioactual;
2692 rest=rest*espprin[0];
2697 tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2);
2698 if( radio>sqrt(tmpsqrt) && rest>dist){
2710 //----------------------------------------------------------------------------
2711 void axisExtractor02::comparacion(vtkImageData *data, vtkImageData *data2, unsigned long minis[4])
2717 data->GetExtent(ext);
2720 unsigned char *ptr2;
2731 radio=((double)ext[1]-(double)ext[0])/2;
2733 centro[0]=(ext[1]+ext[0])/2;
2734 centro[1]=(ext[3]+ext[2])/2;
2735 centro[2]=(ext[5]+ext[4])/2;
2737 for(i=ext[0];i<=ext[1];i++){
2738 for(j=ext[2];j<=ext[3];j++){
2739 for(k=ext[4];k<=ext[5];k++){
2740 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2741 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
2747 if( radio>distancia(punto, centro) ){
2749 if( *ptr==0 && *ptr2==0 ){
2752 else if( *ptr!=0 && *ptr2!=0 ){
2755 else if(*ptr==0 && *ptr2!=0){
2769 //----------------------------------------------------------------------------
2770 void axisExtractor02::copiar(vtkImageData *data, vtkImageData *data2 )
2775 data->GetExtent(ext);
2777 unsigned char *ptr, *ptr2;
2779 for(i=ext[0];i<=ext[1];i++){
2780 for(j=ext[2];j<=ext[3];j++){
2781 for(k=ext[4];k<=ext[5];k++){
2782 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2783 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
2784 if(*ptr!=0 && *ptr2==0){
2797 //----------------------------------------------------------------------------
2798 void axisExtractor02::copiar2(vtkImageData *data, vtkImageData *data2 )
2803 data->GetExtent(ext);
2807 for(i=ext[0];i<=ext[1];i++){
2808 for(j=ext[2];j<=ext[3];j++){
2809 for(k=ext[4];k<=ext[5];k++){
2810 ptr=(double *) data->GetScalarPointer(i,j,k);
2811 ptr2=(double *) data2->GetScalarPointer(i,j,k);
2812 if(*ptr!=0 && *ptr2==0){
2826 //----------------------------------------------------------------------------
2827 double axisExtractor02::angulo(double a[3], double b[3] )
2829 double m1=sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2]));
2830 double m2=sqrt((b[0]*b[0])+(b[1]*b[1])+(b[2]*b[2]));
2831 double res=(acos(((a[0]*b[0])+(a[1]*b[1])+(a[2]*b[2]))/(m1*m2))*180)/3.1415;
2837 if(res>=0 && res<90){
2840 else if(res>=90 && res<180){
2843 else if(res>=180 && res<270){
2856 //----------------------------------------------------------------------------
2857 double axisExtractor02::angulo(double i1, double j1, double k1, double i2, double j2, double k2 )
2859 double m1=sqrt((i1*i1)+(j1*j1)+(k1*k1));
2860 double m2=sqrt((i2*i2)+(j2*j2)+(k2*k2));
2861 double res=(acos(((i1*i2)+(j1*j2)+(k1*k2))/(m1*m2))*180)/3.1415;
2867 if(res>=0 && res<90){
2870 else if(res>=90 && res<180){
2873 else if(res>=180 && res<270){
2887 //----------------------------------------------------------------------------
2888 int axisExtractor02::envolumen(int a[3], vtkImageData *datae )
2893 unsigned short *ptr;
2895 ptr=(unsigned short *) datae->GetScalarPointer(a);
2910 //----------------------------------------------------------------------------
2911 int axisExtractor02::mincandit(int candit[10][3], int cantidad, double puntoanterior[3])
2917 double canditreal[3];
2924 for(i=0;i<cantidad && i<10;i++){
2925 indextoreal(candit[i], canditreal);
2926 distentpu=distancia(canditreal, puntoanterior);
2927 if(distentpu<distmientrp){
2928 distmientrp=distentpu;
2939 //----------------------------------------------------------------------------
2940 int axisExtractor02::maxareacandit(unsigned long vector[50][4], int cantidad)
2950 for(i=0;i<cantidad && i<10;i++){
2951 if(max<vector[i][0]){
2964 //----------------------------------------------------------------------------
2965 unsigned long axisExtractor02::totalarea(unsigned long vector[50][4], unsigned long vectorb[50][4], int cantidad, int cantidadb)
2973 for(i=0;i<cantidad && i<10;i++){
2977 for(i=0;i<cantidadb;i++){
2978 area+=vectorb[i][0];
2988 //----------------------------------------------------------------------------
2989 unsigned long axisExtractor02::conecarea(unsigned long vector[50][4], int cantidad)
2997 for(i=0;i<cantidad && i<10;i++){
3009 //----------------------------------------------------------------------------
3010 int axisExtractor02::bruled(int candit[10][3], int cantidad, vtkImageData *data4)
3016 for(i=0;i<cantidad && i<10;i++){
3017 if(envolumen(candit[i], data4)!=0){
3028 //----------------------------------------------------------------------------
3029 double axisExtractor02::correction(int candit[10][3], int cantidad, vtkImageData *data, int indicecorregido[3], double puntocorregido[3], int indiceanterior[3], double radioanterior, int indicepre[3], double radiopre)
3052 data->GetExtent(ext);
3054 radio=(ext[1]-ext[0])/2;
3055 radiof=((double)ext[1]-(double)ext[0])/2;
3056 centro[0]=(ext[1]+ext[0])/2;
3057 centro[1]=(ext[3]+ext[2])/2;
3058 centro[2]=(ext[5]+ext[4])/2;
3060 indicecorregido[0]=centro[0];
3061 indicecorregido[1]=centro[1];
3062 indicecorregido[2]=centro[2];
3063 indextoreal(indicecorregido, puntocorregido);
3064 indextoreal(indiceanterior, centrof);
3065 indextoreal(indicepre, puntof2);
3067 rap= (int)( (double)radio/2 );
3071 for(i=ext[0];i<=ext[1];i++){
3072 for(j=ext[2];j<=ext[3];j++){
3073 for(k=ext[4];k<=ext[5];k++){
3074 ptr=(double *) data->GetScalarPointer(i,j,k);
3078 indextoreal(punto, puntof );
3080 distcorr3=distancia(indiceanterior, indicepre);
3081 distcorr2=distancia(punto, indicepre);
3082 distcorr=distancia(punto, indiceanterior);
3084 dist=distanciaejepunto(puntof, centrof, puntof2);
3085 prop=proporcioejepunto(puntof, centrof, puntof2);
3087 if(prop>=0 && prop<=1){
3088 rest=(radiopre-radioanterior)*prop+radioanterior;
3089 rest=rest*espprin[0];
3095 if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && rest>dist){
3096 //if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && rest>dist){
3097 //if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && distcorr2<=distcorr3-radioanterior){
3101 indicecorregido[0]=i;
3102 indicecorregido[1]=j;
3103 indicecorregido[2]=k;
3104 indextoreal(indicecorregido, puntocorregido);
3106 else if( maxcent==*ptr){
3107 if( mincent>distcorr){
3110 indicecorregido[0]=i;
3111 indicecorregido[1]=j;
3112 indicecorregido[2]=k;
3113 indextoreal(indicecorregido, puntocorregido);
3121 ptr=(double *) data->GetScalarPointer(centro[0],centro[1],centro[2]);
3129 //----------------------------------------------------------------------------
3130 double axisExtractor02::correction2(int candit[10][3], int cantidad, vtkImageData *data, int indicecorregido[3], double puntocorregido[3], int indiceanterior[3], double radioanterior)
3145 data->GetExtent(ext);
3147 radio=(ext[1]-ext[0])/2;
3148 centro[0]=(ext[1]+ext[0])/2;
3149 centro[1]=(ext[3]+ext[2])/2;
3150 centro[2]=(ext[5]+ext[4])/2;
3152 indicecorregido[0]=centro[0];
3153 indicecorregido[1]=centro[1];
3154 indicecorregido[2]=centro[2];
3155 indextoreal(indicecorregido, puntocorregido);
3157 rap= (int)((double)radio/2);
3161 for(i=centro[0]-rap;i<=centro[0]+rap;i++){
3162 for(j=centro[1]-rap;j<=centro[1]+rap;j++){
3163 for(k=centro[2]-rap;k<=centro[2]+rap;k++){
3164 ptr=(double *) data->GetScalarPointer(i,j,k);
3169 distcorr2=distancia(punto, centro);
3170 distcorr=distancia(punto, indiceanterior);
3172 if( rap>distcorr2 ){
3176 indicecorregido[0]=i;
3177 indicecorregido[1]=j;
3178 indicecorregido[2]=k;
3179 indextoreal(indicecorregido, puntocorregido);
3181 else if( maxcent==*ptr){
3182 if( mincent>distcorr){
3185 indicecorregido[0]=i;
3186 indicecorregido[1]=j;
3187 indicecorregido[2]=k;
3188 indextoreal(indicecorregido, puntocorregido);
3198 return sqrt(maxcent);
3205 //----------------------------------------------------------------------------
3206 void axisExtractor02::avanzar()
3209 double puntoactual[3];
3210 double puntocorregido[3];
3211 double puntoanterior[3];
3212 double puntosiguiente[3];
3213 //double puntointer[3];
3216 int indiceactual[3];
3217 int indiceanterior[3];
3218 int indicecorregido[3];
3219 //int indiceinter[3];
3224 unsigned char cantidad;
3225 unsigned char cantidadb;
3227 unsigned long vector[50][4];
3228 unsigned long vectorb[50][4];
3268 double radiotrabajof;
3273 double radioanterior;
3279 double canditreal[3];
3281 unsigned long totalsurf;
3282 unsigned long conecsurf;
3285 unsigned long caf[4];
3294 if(!m_Stack0.empty()){
3297 /* indexp=m_Stack.top();
3299 puntoactual[0]=m_Stack0.top();
3300 puntoactual[1]=m_Stack1.top();
3301 puntoactual[2]=m_Stack2.top();
3303 puntoanterior[0]=m_Stack3.top();
3304 puntoanterior[1]=m_Stack4.top();
3305 puntoanterior[2]=m_Stack5.top();
3307 puntopre[0]=m_Stack6.top();
3308 puntopre[1]=m_Stack7.top();
3309 puntopre[2]=m_Stack8.top();
3313 radiopred=m_Stackr.top();
3314 radioanterior=m_Stackra.top();
3315 radioprinc=m_Stackrp.top();*/
3317 indexp=m_Stack.front();
3319 puntoactual[0]=m_Stack0.front();
3320 puntoactual[1]=m_Stack1.front();
3321 puntoactual[2]=m_Stack2.front();
3323 puntoanterior[0]=m_Stack3.front();
3324 puntoanterior[1]=m_Stack4.front();
3325 puntoanterior[2]=m_Stack5.front();
3327 puntopre[0]=m_Stack6.front();
3328 puntopre[1]=m_Stack7.front();
3329 puntopre[2]=m_Stack8.front();
3331 radiopred=m_Stackr.front();
3332 radioanterior=m_Stackra.front();
3333 radioprinc=m_Stackrp.front();
3335 radiotrabajo= (int) (proprad*radiopred);
3336 radiotrabajof=proprad*radiopred;
3338 if(radiotrabajof-radiotrabajo>0.5){
3339 radiotrabajo=radiotrabajo+1;
3342 // EED 15 Mai 2007 .NET
3343 // radiotrabajo=max(minant,radiotrabajo);
3344 if (minant>radiotrabajo)
3346 radiotrabajo=minant;
3364 m_Stack0.pop_front();
3365 m_Stack1.pop_front();
3366 m_Stack2.pop_front();
3367 m_Stack3.pop_front();
3368 m_Stack4.pop_front();
3369 m_Stack5.pop_front();
3370 m_Stack6.pop_front();
3371 m_Stack7.pop_front();
3372 m_Stack8.pop_front();
3374 m_Stack.pop_front();
3375 m_Stackr.pop_front();
3376 m_Stackra.pop_front();
3377 m_Stackrp.pop_front();
3379 dirant[0]=puntoanterior[0]-puntoactual[0];
3380 dirant[1]=puntoanterior[1]-puntoactual[1];
3381 dirant[2]=puntoanterior[2]-puntoactual[2];
3384 realtoindex(puntoactual, indiceactual);
3385 realtoindex(puntoanterior, indiceanterior);
3386 realtoindex(puntopre, indicepre);
3388 radioactual=radiopred;
3390 if(radiotrabajo<=maxant){
3392 extint[0]=indiceactual[0]-radiotrabajo;
3393 extint[1]=indiceactual[0]+radiotrabajo;
3394 extint[2]=indiceactual[1]-radiotrabajo;
3395 extint[3]=indiceactual[1]+radiotrabajo;
3396 extint[4]=indiceactual[2]-radiotrabajo;
3397 extint[5]=indiceactual[2]+radiotrabajo;
3399 extrac->SetVOI(extint);
3400 extrac->UpdateWholeExtent();
3402 extrac->GetOutput()->GetExtent(extintreal);
3404 if(extint[0]!=extintreal[0] || extint[1]!=extintreal[1] || extint[2]!=extintreal[2] || extint[3]!=extintreal[3] || extint[4]!=extintreal[4] || extint[5]!=extintreal[5]){
3405 fprintf(stream, "salio\t");
3411 data1=vtkImageData::New();
3412 data1->SetScalarTypeToUnsignedChar();
3413 data1->SetExtent(extrac->GetOutput()->GetExtent());
3414 data1->SetSpacing(extrac->GetOutput()->GetSpacing());
3416 cilindro(data1, dirant );
3418 optim(extrac->GetOutput(), data1 );
3420 connect->SetInput(data1);
3421 connect->RemoveAllSeeds();
3422 connect->AddSeed(indiceactual[0],indiceactual[1],indiceactual[2]);
3423 connect->UpdateWholeExtent();
3427 data2=vtkImageData::New();
3428 data2->SetScalarTypeToUnsignedChar();
3429 data2->SetExtent(connect->GetOutput()->GetExtent());
3430 data2->SetSpacing(connect->GetOutput()->GetSpacing());
3434 cantidad=find_components(connect->GetOutput(), data2, 0, vector);
3437 data3=vtkImageData::New();
3438 data3->SetScalarTypeToUnsignedChar();
3439 data3->SetExtent(connect->GetOutput()->GetExtent());
3440 data3->SetSpacing(connect->GetOutput()->GetSpacing());
3444 cantidadb=find_componentsb(connect->GetOutput(), data3, 0, vectorb);
3446 redondear3(connect->GetOutput());
3450 redondear2(distance->GetOutput());
3452 redondear(connect->GetOutput());
3454 costominimo(distance->GetOutput(), data2);
3458 radioactual=correction2(candit, cantidad, distance->GetOutput(), indicecorregido, puntocorregido, indiceanterior, radioanterior);
3461 radioactual=correction(candit, cantidad, distance->GetOutput(), indicecorregido, puntocorregido, indiceanterior, radioanterior, indicepre, radioprinc);
3465 //inpeq(extrac->GetOutput(), indicecorregido, radioactual);
3467 indant=mincandit(candit, cantidad, puntoanterior);
3469 indmaxarea=maxareacandit(vector, cantidad);
3470 totalsurf=totalarea(vector, vectorb, cantidad, cantidadb);
3471 conecsurf=conecarea(vector, cantidad);
3473 indextoreal(candit[indant], canditreal);
3475 burned=bruled(candit, cantidad, data4);
3478 data6=vtkImageData::New();
3479 data6->SetScalarTypeToUnsignedChar();
3480 data6->SetExtent(extrac->GetOutput()->GetExtent());
3481 data6->SetSpacing(extrac->GetOutput()->GetSpacing());
3483 modelo(data6, cantidad, vector, candit, radioactual, minis);
3484 comparacion(connect->GetOutput(), data6, caf);
3486 rel1=((double)caf[0]+(double)caf[1])/((double)caf[0]+(double)caf[1]+(double)caf[2]+(double)caf[3]);
3487 rel2=((double)caf[1])/((double)caf[1]+(double)caf[3]);
3488 rel3=((double)caf[1])/((double)caf[1]+(double)caf[2]);
3489 rel4=(double)conecsurf/(double)totalsurf;
3490 rel5=(double)vector[indmaxarea][0]/(double)totalsurf;
3491 rel6=(double)cant/((double)cant+(double)cant2);
3494 if(rel1<0.5 || rel2<0.5 || rel3<0.5 || (rel1<0.75 && rel2<0.75) || (rel2<0.75 && rel3<0.75)|| (rel1<0.75 && rel3<0.75)){
3498 if(indicecorregido[0]!=indiceactual[0] || indicecorregido[1]!=indiceactual[1] || indicecorregido[2]!=indiceactual[2]){
3502 if(burned>=cantidad ){
3514 if(flag7==1 && flagg>4){
3526 if(flag7==1 && flag18==1){
3527 for(i=0;i<flagg;i++){
3528 if(visited[i][0]==indiceactual[0] && visited[i][1]==indiceactual[1] && visited[i][2]==indiceactual[2] && visitedrad[i]==radiopred){
3534 if((flag19==1 && flagg2<5) || (flag19==1 && flagg2<10 && ( rel4>0.25 || rel5>0.13 || rel6>0.4 || flag17==1 || flag15==1 || flag5==1) ) || (flagg2<10 && ( ((rel4>0.25 || rel5>0.25 || rel6>0.4) && cantidad<3) || flag17==1 || flag15==1 || flag5==1) )){
3542 if(angulo(dirant[0], dirant[1], dirant[2], canditreal[0]-puntoactual[0], canditreal[1]-puntoactual[1], canditreal[2]-puntoactual[2] )>65){
3546 if(angulo(Vp[0][ejeminp], Vp[1][ejeminp], Vp[2][ejeminp], Vi[0][ejemini], Vi[1][ejemini], Vi[2][ejemini] )>60){
3550 if((flag3==1 || flag5==1 || (flag16==1 && flag11==1 && flag19==0) ) && (flag7==0 && flag8==0)){
3552 //fprintf(stream, "malo\t");
3555 if(flag3==1 || flag5==1 || (flag16==1 && flag11==1 && flag19==0) ){
3559 if((mejordst<radioactual || (mejordst==radioactual && mejorrad>radiopred) || (mejordst==radioactual && mejorrad<=radiopred && mejorcant<cantidad )) && flag17==0 && flag20==0 && flag4==0 && flag12==0 && flag15==0 ){
3560 mejordst=radioactual;
3561 mejor[0]=puntoactual[0];
3562 mejor[1]=puntoactual[1];
3563 mejor[2]=puntoactual[2];
3570 fprintf(stream, "grande\t");
3586 if( flag9==1 || flag4==1 || flag12==1 || ((cantidad<2 || burned==cantidad) && (flag7==0 && flag8==0)) ){
3593 if( cantidad>2 && (flag7==0 && flag8==0) ){
3600 if((flag7==0 && flag8==0 && (mejordst>radioactual || (mejordst==radioactual && mejorrad<radiopred)) && flag19==0) || flag2==1){
3606 fprintf(stream, "mejor\t");
3608 m_Stack0.push_front(mejor[0]);
3609 m_Stack1.push_front(mejor[1]);
3610 m_Stack2.push_front(mejor[2]);
3611 m_Stack3.push_front(puntoanterior[0]);
3612 m_Stack4.push_front(puntoanterior[1]);
3613 m_Stack5.push_front(puntoanterior[2]);
3614 m_Stack6.push_front(puntopre[0]);
3615 m_Stack7.push_front(puntopre[1]);
3616 m_Stack8.push_front(puntopre[2]);
3617 m_Stack.push_front(indexp);
3618 m_Stackr.push_front(mejorrad);
3619 m_Stackra.push_front(radioanterior);
3620 m_Stackrp.push_front(radioprinc);
3627 else if((flag17==1 || flag9==1 || flag4==1 || flag12==1 || flag15==1) && flag18==1){
3632 fprintf(stream, "descorrigio\t");
3634 m_Stack0.push_front(mejor[0]);
3635 m_Stack1.push_front(mejor[1]);
3636 m_Stack2.push_front(mejor[2]);
3637 m_Stack3.push_front(puntoanterior[0]);
3638 m_Stack4.push_front(puntoanterior[1]);
3639 m_Stack5.push_front(puntoanterior[2]);
3640 m_Stack6.push_front(puntopre[0]);
3641 m_Stack7.push_front(puntopre[1]);
3642 m_Stack8.push_front(puntopre[2]);
3643 m_Stack.push_front(indexp);
3644 m_Stackr.push_front(mejorrad);
3645 m_Stackra.push_front(radioanterior);
3646 m_Stackrp.push_front(radioprinc);
3653 else if(flag12==0 && flag4==0 && flag9==0){
3654 if(flag7==0 && flag8==0){
3656 if(flag17==0 && flag15==0){
3663 fprintf(stream, "inserto\t");
3666 /*costominimo2(distance->GetOutput(), connect->GetOutput(), indiceactual, candit[indant], indiceinter);
3667 indextoreal(indiceinter, puntointer);
3668 realtoreal2(puntointer, provvc);
3669 points->InsertPoint(buenos,provvc);
3671 lineas->InsertNextCell(2);
3672 lineas->InsertCellPoint(indexp);
3673 lineas->InsertCellPoint(buenos);
3677 realtoreal2(puntoactual, provvc);
3678 points->InsertPoint(buenos,provvc);
3680 lineas->InsertNextCell(2);
3681 lineas->InsertCellPoint(buenos-1);
3682 lineas->InsertCellPoint(buenos);
3686 realtoreal2(puntoactual, provvc);
3687 points->InsertPoint(buenos,provvc);
3689 lineas->InsertNextCell(2);
3690 lineas->InsertCellPoint(indexp);
3691 lineas->InsertCellPoint(buenos);
3697 realtoreal2(puntoactual, provvc);
3698 points->InsertPoint(buenos,provvc);
3711 fprintf(stream, "para\t");
3714 else if(flag7==0 && flag8==1){
3717 fprintf(stream, "agrando\t");
3719 m_Stack0.push_front(puntoactual[0]);
3720 m_Stack1.push_front(puntoactual[1]);
3721 m_Stack2.push_front(puntoactual[2]);
3722 m_Stack3.push_front(puntoanterior[0]);
3723 m_Stack4.push_front(puntoanterior[1]);
3724 m_Stack5.push_front(puntoanterior[2]);
3725 m_Stack6.push_front(puntopre[0]);
3726 m_Stack7.push_front(puntopre[1]);
3727 m_Stack8.push_front(puntopre[2]);
3729 m_Stack.push_front(indexp);
3730 m_Stackr.push_front(((double)radiotrabajo+1)/2);
3731 m_Stackra.push_front(radioanterior);
3732 m_Stackrp.push_front(radioprinc);
3735 else if(flag7==1 && flag8==1){
3736 visited[flagg][0]=indiceactual[0];
3737 visited[flagg][1]=indiceactual[1];
3738 visited[flagg][2]=indiceactual[2];
3739 visitedrad[flagg]=radiopred;
3743 fprintf(stream, "agrando y corrigiendo\t");
3745 m_Stack0.push_front(puntocorregido[0]);
3746 m_Stack1.push_front(puntocorregido[1]);
3747 m_Stack2.push_front(puntocorregido[2]);
3748 m_Stack3.push_front(puntoanterior[0]);
3749 m_Stack4.push_front(puntoanterior[1]);
3750 m_Stack5.push_front(puntoanterior[2]);
3751 m_Stack6.push_front(puntopre[0]);
3752 m_Stack7.push_front(puntopre[1]);
3753 m_Stack8.push_front(puntopre[2]);
3755 m_Stack.push_front(indexp);
3756 m_Stackr.push_front(((double)radiotrabajo+1)/2);
3757 m_Stackra.push_front(radioanterior);
3758 m_Stackrp.push_front(radioprinc);
3761 visited[flagg][0]=indiceactual[0];
3762 visited[flagg][1]=indiceactual[1];
3763 visited[flagg][2]=indiceactual[2];
3764 visitedrad[flagg]=radiopred;
3771 fprintf(stream, "corrigio\t");
3773 m_Stack0.push_front(puntocorregido[0]);
3774 m_Stack1.push_front(puntocorregido[1]);
3775 m_Stack2.push_front(puntocorregido[2]);
3776 m_Stack3.push_front(puntoanterior[0]);
3777 m_Stack4.push_front(puntoanterior[1]);
3778 m_Stack5.push_front(puntoanterior[2]);
3780 m_Stack6.push_front(puntopre[0]);
3781 m_Stack7.push_front(puntopre[1]);
3782 m_Stack8.push_front(puntopre[2]);
3784 m_Stack.push_front(indexp);
3785 m_Stackr.push_front(radiopred);
3786 m_Stackra.push_front(radioanterior);
3787 m_Stackrp.push_front(radioprinc);
3792 if(flag7==0 && flag8==0){
3794 if(flag17==0 && flag15==0){
3796 for(i=0;i<cantidad && i<10;i++){
3797 for(j=i;j<cantidad && i<10;j++){
3798 if(minis[j]>minis[i]){
3799 vectortemp[0]=candit[i][0];
3800 vectortemp[1]=candit[i][1];
3801 vectortemp[2]=candit[i][2];
3804 candit[i][0]=candit[j][0];
3805 candit[i][1]=candit[j][1];
3806 candit[i][2]=candit[j][2];
3809 candit[j][0]=vectortemp[0];
3810 candit[j][1]=vectortemp[1];
3811 candit[j][2]=vectortemp[2];
3819 for(i=0;i<cantidad && i<10;i++){
3820 if(i!=indant || buenos<2){
3821 indextoreal(candit[i], puntosiguiente);
3823 m_Stack0.push_front(puntosiguiente[0]);
3824 m_Stack1.push_front(puntosiguiente[1]);
3825 m_Stack2.push_front(puntosiguiente[2]);
3826 m_Stack3.push_front(puntoactual[0]);
3827 m_Stack4.push_front(puntoactual[1]);
3828 m_Stack5.push_front(puntoactual[2]);
3830 m_Stack6.push_front(puntosiguiente[0]);
3831 m_Stack7.push_front(puntosiguiente[1]);
3832 m_Stack8.push_front(puntosiguiente[2]);
3834 double prov1a=sqrt(minis[i]);
3836 m_Stack.push_front(buenos-1);
3837 m_Stackr.push_front(prov1a);
3838 m_Stackra.push_front(radioactual);
3839 m_Stackrp.push_front(prov1a);
3845 for(i=0;i<cantidad && i<10;i++){
3846 indextoreal(candit[i], puntosiguiente);
3848 if(envolumen(candit[i], data4)==0){
3849 m_Stack0.push_front(puntosiguiente[0]);
3850 m_Stack1.push_front(puntosiguiente[1]);
3851 m_Stack2.push_front(puntosiguiente[2]);
3852 m_Stack3.push_front(puntoactual[0]);
3853 m_Stack4.push_front(puntoactual[1]);
3854 m_Stack5.push_front(puntoactual[2]);
3856 m_Stack6.push_front(puntosiguiente[0]);
3857 m_Stack7.push_front(puntosiguiente[1]);
3858 m_Stack8.push_front(puntosiguiente[2]);
3860 double prov1b=sqrt(minis[i]);
3862 m_Stack.push_front(buenos-1);
3863 m_Stackr.push_front(prov1b);
3864 m_Stackra.push_front(radioactual);
3865 m_Stackrp.push_front(prov1b);
3870 copiar(connect->GetOutput(), data4 );
3871 //copiar2(distance->GetOutput(), data5 );
3876 // fprintf(stream, "algo");
3889 printf("wi[0]: %f\n", wi[0]);
3890 printf("wi[1]: %f\n", wi[1]);
3891 printf("wi[2]: %f\n", wi[2]);
3893 printf("inerciari: %f\n", inerciari);
3894 printf("inerciariy: %f\n", inerciariy);
3895 printf("inerciariz: %f\n", inerciariz);
3897 printf("ri1: %f\n", inerciari/inerciariy);
3898 printf("ri2: %f\n", inerciariy/inerciariz);
3899 printf("ri3: %f\n", inerciari/inerciariz);
3901 printf("inerciarp: %f\n", inerciarp);
3902 printf("inerciarpy: %f\n", inerciarpy);
3903 printf("inerciarpz: %f\n", inerciarpz);
3905 printf("rp1: %f\n", inerciarp/inerciarpy);
3906 printf("rp2: %f\n", inerciarpy/inerciarpz);
3907 printf("rp3: %f\n", inerciarp/inerciarpz);
3909 printf("angulo: %f\n",angulo(dirant[0], dirant[1], dirant[2], canditreal[0]-puntoactual[0], canditreal[1]-puntoactual[1], canditreal[2]-puntoactual[2] ));
3910 printf("angulo2: %f\n", angulo(Vp[0][ejeminp], Vp[1][ejeminp], Vp[2][ejeminp], Vi[0][ejemini], Vi[1][ejemini], Vi[2][ejemini]));
3912 printf("centi: %f\n", centi);
3913 printf("centi2: %f\n", centi2);
3918 fprintf(stream, "%d\t", indexp);
3919 fprintf(stream, "%d\t", buenos);
3920 fprintf(stream, "%d\t", radiotrabajo);
3921 fprintf(stream, "%f\t", radioactual);
3922 fprintf(stream, "%f\t", radiopred);
3923 fprintf(stream, "%f\t", radioanterior);
3924 fprintf(stream, "%f\t", mejordst);
3925 fprintf(stream, "%f\t", mejorrad);
3926 fprintf(stream, "%d\t", mejorcant);
3930 fprintf(stream, "%d\t", cant);
3931 fprintf(stream, "%d\t", cant2);
3932 fprintf(stream, "%f\t", (double)cant/((double)cant+(double)cant2));
3935 fprintf(stream, "%d\t", flag2);
3936 fprintf(stream, "%d\t", flag3);
3937 fprintf(stream, "%d\t", flag4);
3938 fprintf(stream, "%d\t", flag5);
3939 fprintf(stream, "%d\t", flag6);
3940 fprintf(stream, "%d\t", flag7);
3941 fprintf(stream, "%d\t", flag8);
3942 fprintf(stream, "%d\t", flag9);
3943 fprintf(stream, "%d\t", flag10);
3944 fprintf(stream, "%d\t", flag11);
3945 fprintf(stream, "%d\t", flag12);
3946 fprintf(stream, "%d\t", flag13);
3947 fprintf(stream, "%d\t", flag14);
3948 fprintf(stream, "%d\t", flag15);
3949 fprintf(stream, "%d\t", flag16);
3950 fprintf(stream, "%d\t", flag17);
3951 fprintf(stream, "%d\t", flag18);
3952 fprintf(stream, "%d\t", flag19);
3953 fprintf(stream, "%d\t", flag20);
3954 fprintf(stream, "%d\t", flagg);
3955 fprintf(stream, "%d\t", flagg2);
3956 fprintf(stream, "%d\t", cantidad);
3957 fprintf(stream, "%d\t", cantidadb);
3958 fprintf(stream, "%d\t", burned);
3962 fprintf(stream, "caf[0] %d\t", caf[0]);
3963 fprintf(stream, "caf[1] %d\t", caf[1]);
3964 fprintf(stream, "caf[2] %d\t", caf[2]);
3965 fprintf(stream, "caf[3] %d\t", caf[3]);
3966 fprintf(stream, "REL%f\t", ((double)caf[0]+(double)caf[1])/((double)caf[0]+(double)caf[1]+(double)caf[2]+(double)caf[3]));
3967 fprintf(stream, "REL2%f\t", ((double)caf[1])/((double)caf[3]+(double)caf[1]));
3968 fprintf(stream, "REL3%f\t", ((double)caf[1])/((double)caf[1]+(double)caf[2]));
3970 fprintf(stream, "maxareacandit%d\t", vector[indmaxarea][0]);
3971 fprintf(stream, "totalsurf%d\t", totalsurf);
3972 fprintf(stream, "conecsurf%d\t", conecsurf);
3973 fprintf(stream, "ratio%f\t", (double)conecsurf/(double)totalsurf);
3974 fprintf(stream, "ratio2%f\t", (double)vector[indmaxarea][0]/(double)totalsurf);
3978 fprintf(stream, "[%f, %f, %f]\t", puntoactual[0], puntoactual[1], puntoactual[2]);
3979 fprintf(stream, "[%f, %f, %f]\t", puntoanterior[0], puntoanterior[1], puntoanterior[2]);
3980 fprintf(stream, "[%f, %f, %f]\t", puntopre[0], puntopre[1], puntopre[2]);
3981 fprintf(stream, "[%f, %f, %f]\t", mejor[0], mejor[1], mejor[2]);
3983 fprintf(stream, "\n");
3992 //----------------------------------------------------------------------------
3993 void axisExtractor02::todo()
3997 while( !m_Stack0.empty() && i<5000){
4005 //----------------------------------------------------------------------------
4006 void axisExtractor02::paso()
4012 while( !m_Stack0.empty() && inicio==buenos && i<5000){
4021 //----------------------------------------------------------------------------
4022 void axisExtractor02::rama()
4028 while( !m_Stack0.empty() && frama==0 && i<5000){
4037 //----------------------------------------------------------------------------
4038 void axisExtractor02::segmento()
4045 while( !m_Stack0.empty() && frama==0 && fseg==0 && i<5000){