1 /*=========================================================================
4 =========================================================================*/
5 #include "axisExtractor02.h"
6 #include "vtkPolyData.h"
7 #include "vtkObjectFactory.h"
8 #include "vtkImageThreshold.h"
9 #include "vtkImageCast.h"
10 #include "vtkImageSeedConnectivity.h"
11 #include "vtkImageData.h"
12 #include "vtkMarchingCubes.h"
13 #include "vtkDoubleArray.h"
14 #include "vtkPointData.h"
15 #include "vtkExtractVOI.h"
16 #include "vtkImageConstantPad.h"
17 #include "vtkImageTranslateExtent.h"
19 #include "vtkTransformPolyDataFilter.h"
20 #include "vtkTransform.h"
27 vtkStandardNewMacro(axisExtractor02);
30 //----------------------------------------------------------------------------
31 void axisExtractor02::SetParam(double value)
38 //----------------------------------------------------------------------------
39 double axisExtractor02::GetParam()
47 //----------------------------------------------------------------------------
48 void axisExtractor02::SetParam2(double value)
55 //----------------------------------------------------------------------------
56 double axisExtractor02::GetParam2()
64 //----------------------------------------------------------------------------
65 void axisExtractor02::SetParam3(double value)
72 //----------------------------------------------------------------------------
73 double axisExtractor02::GetParam3()
80 //----------------------------------------------------------------------------
81 void axisExtractor02::SetMaxant(int value)
88 //----------------------------------------------------------------------------
89 int axisExtractor02::GetMaxant()
96 //----------------------------------------------------------------------------
97 void axisExtractor02::SetMinant(int value)
104 //----------------------------------------------------------------------------
105 int axisExtractor02::GetMinant()
111 //----------------------------------------------------------------------------
112 void axisExtractor02::SetPoint(double value[3])
118 this->GetInput()->GetSpacing (spa);
120 value[0]=value[0]+maxant*spa[0];
121 value[1]=value[1]+maxant*spa[1];
122 value[2]=value[2]+maxant*spa[2];
124 this->m_Stack0.push_front(value[0]);
125 this->m_Stack1.push_front(value[1]);
126 this->m_Stack2.push_front(value[2]);
127 this->m_Stack3.push_front(value[0]);
128 this->m_Stack4.push_front(value[1]);
129 this->m_Stack5.push_front(value[2]);
130 this->m_Stack6.push_front(value[0]);
131 this->m_Stack7.push_front(value[1]);
132 this->m_Stack8.push_front(value[2]);
133 this->m_Stack.push_front(0);
134 this->m_Stackr.push_front(0);
135 this->m_Stackra.push_front(0);
136 this->m_Stackrp.push_front(0);
141 //----------------------------------------------------------------------------
142 vtkImageData *axisExtractor02::GetVolumen()
145 vtkImageTranslateExtent *trans;
147 trans = vtkImageTranslateExtent::New();
149 trans->SetInput(this->data4);
151 trans->SetTranslation (maxant, maxant, maxant);
155 return trans->GetOutput();
162 //----------------------------------------------------------------------------
163 vtkPolyData *axisExtractor02::GetOutput ()
169 vtkTransform *transL1;
170 vtkTransformPolyDataFilter *trans;
172 transL1 = vtkTransform::New();
173 trans = vtkTransformPolyDataFilter::New();
177 this->GetInput()->GetSpacing (spa);
178 transL1->Translate(-maxant*spa[0], -maxant*spa[1], -maxant*spa[2]);
179 // transL1->Translate(-maxant, -maxant, -maxant);
182 trans->SetInput((vtkPolyData *)(this->Outputs[0]));
184 trans->SetTransform (transL1);
188 return trans->GetOutput();
197 //----------------------------------------------------------------------------
198 axisExtractor02::axisExtractor02() {
201 this->NumberOfRequiredInputs = 1;
207 /*this->resample= vtkImageResample::New();
208 this->resample->SetDimensionality (3);*/
210 this->data4=vtkImageData::New();
212 this->data1=vtkImageData::New();
213 this->data2=vtkImageData::New();
214 this->data3=vtkImageData::New();
215 this->data6=vtkImageData::New();
217 this->extrac = vtkExtractVOI::New();
218 // this->extrac->SetInput(resample->GetOutput());
219 this->extrac->SetSampleRate(1, 1, 1);
221 this->connect = vtkImageSeedConnectivity::New();
222 this->connect->SetInput(this->data1);
223 this->connect->SetInputConnectValue(255);
224 this->connect->SetOutputConnectedValue(255);
225 this->connect->SetOutputUnconnectedValue(0);
227 this->distance= vtkImageEuclideanDistance::New();
228 this->distance->SetInput(this->connect->GetOutput());
229 this->distance->InitializeOn();
230 this->distance->ConsiderAnisotropyOff();
237 //----------------------------------------------------------------------------
238 void axisExtractor02::SetInput(vtkImageData *input)
241 vtkImageConstantPad *pad;
242 vtkImageTranslateExtent *trans;
244 pad = vtkImageConstantPad::New();
245 trans = vtkImageTranslateExtent::New();
247 pad->SetInput(input);
248 trans->SetInput(pad->GetOutput());
253 pad->SetInput(input);
257 input->GetExtent(ext);
259 ext[0]=ext[0]-maxant;
260 ext[1]=ext[1]+maxant;
261 ext[2]=ext[2]-maxant;
262 ext[3]=ext[3]+maxant;
263 ext[4]=ext[4]-maxant;
264 ext[5]=ext[5]+maxant;
266 pad->SetOutputWholeExtent(ext);
268 trans->SetTranslation (maxant, maxant, maxant);
273 //this->vtkProcessObject::SetNthInput(0, input);
274 this->vtkProcessObject::SetNthInput(0, trans->GetOutput());
276 //this->GetInput()->SetOrigin(0,0,0);
277 //this->resample->SetInput(this->GetInput());
278 this->extrac->SetInput(this->GetInput());
284 //----------------------------------------------------------------------------
285 void axisExtractor02::distanciaejes(vtkPolyData *eje1, vtkPolyData *eje2)
287 vtkIdType totalpuntos1=eje1->GetNumberOfPoints();
288 vtkIdType totalpuntos2=eje2->GetNumberOfPoints();
289 vtkIdType totallineas=eje2->GetNumberOfLines();
295 double point[3], point2[3], point3[3], point4[3];
299 double x, dist, min, y;
302 ff=fopen("comparacion.txt","w");
306 for(i=0;i<totalpuntos1;i++){
307 eje1->GetPoint(i, point);
309 for(j=0;j<totallineas;j++){
310 lineas=eje2->GetCell(j);
312 lineas->GetPoints()->GetPoint(0, point3);
313 lineas->GetPoints()->GetPoint(1, point2);
314 v[0]=point3[0]-point2[0];
315 v[1]=point3[1]-point2[1];
316 v[2]=point3[2]-point2[2];
318 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]));
321 point4[0]=point2[0]+(x*v[0]);
322 point4[1]=point2[1]+(x*v[1]);
323 point4[2]=point2[2]+(x*v[2]);
325 v2[0]=point4[0]-point[0];
326 v2[1]=point4[1]-point[1];
327 v2[2]=point4[2]-point[2];
329 y=(v[0]*v2[0])+(v[1]*v2[1])+(v[2]*v2[2]);
332 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])));
337 /*printf("punto: %f %f %f\n", point[0], point[1], point[2]);
338 printf("punto2: %f %f %f\n", point2[0], point2[1], point2[2]);
339 printf("punto3: %f %f %f\n", point3[0], point3[1], point3[2]);
340 printf("punto4: %f %f %f\n", point4[0], point4[1], point4[2]);
341 printf("v: %f %f %f\n", v[0], v[1], v[2]);
342 printf("v2: %f %f %f\n", v2[0], v2[1], v2[2]);
355 fprintf(ff,"%f\n", min);
371 //----------------------------------------------------------------------------
372 void axisExtractor02::Execute()
388 this->GetInput()->GetExtent(extprin0);
389 this->GetInput()->GetExtent(extprin);
390 this->GetInput()->GetSpacing(espprin);
393 // stream = fopen( "resultado.txt", "w" );
394 // 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");
396 /*minspac=min(espprin[0],min(espprin[1],espprin[2]));
398 this->resample->SetAxisOutputSpacing( 0, minspac);
399 this->resample->SetAxisOutputSpacing( 1, minspac);
400 this->resample->SetAxisOutputSpacing( 2, minspac);
401 resample->Update();*/
404 this->data4->SetScalarTypeToUnsignedChar();
405 this->data4->SetExtent(this->GetInput()->GetExtent());
406 this->data4->SetSpacing(this->GetInput()->GetSpacing());
408 this->blanquear(this->data4);
410 this->points = vtkPoints::New();
412 this->lineas = vtkCellArray::New();
421 ((vtkPolyData *)this->Outputs[0])->SetPoints (this->points);
422 ((vtkPolyData *)this->Outputs[0])->SetLines(this->lineas);
424 dif = difftime (end,start);
426 // ff=fopen("time.txt","w");
427 // fprintf(ff,"%d \n",this->buenos);
428 // fprintf(ff,"%.2lf \n",dif);
429 // fprintf(ff,"%.2lf \n",(double)this->buenos/dif);
439 //----------------------------------------------------------------------------
440 vtkImageData *axisExtractor02::GetInput()
442 if (this->NumberOfInputs < 1){
446 return (vtkImageData *)(this->Inputs[0]);
455 //----------------------------------------------------------------------------
456 void axisExtractor02::PrintSelf(ostream& os, vtkIndent indent)
458 this->Superclass::PrintSelf(os,indent);
464 //----------------------------------------------------------------------------
465 void axisExtractor02::realtoreal(double a[3], double b[3] )
468 b[0]=a[0]+(espprin[0]*extprin[0]);
469 b[1]=a[1]+(espprin[1]*extprin[2]);
478 //----------------------------------------------------------------------------
479 void axisExtractor02::realtoreal2(double a[3], double b[3] )
482 b[0]=a[0]-(espprin[0]*extprin[0]);
483 b[1]=a[1]-(espprin[1]*extprin[2]);
489 //----------------------------------------------------------------------------
490 void axisExtractor02::realtoindex(double a[3], int b[3] )
496 // EED 15 Mai 2007 .NET
497 // minspac=min(espprin[0],min(espprin[1],espprin[2]));
500 if (espprin[1]<minspac) { minspac=espprin[1]; }
501 if (espprin[2]<minspac) { minspac=espprin[2]; }
508 b[0]=(int)(a[0]/minspac);
509 b[1]=(int)(a[1]/minspac);
510 b[2]=(int)(a[2]/minspac);
531 //----------------------------------------------------------------------------
532 void axisExtractor02::indextoreal(int a[3], double b[3] )
536 // EED 15 Mai 2007 .NET
537 // minspac=min(espprin[0],min(espprin[1],espprin[2]));
539 if (espprin[1]<minspac) { minspac=espprin[1]; }
540 if (espprin[2]<minspac) { minspac=espprin[2]; }
549 //----------------------------------------------------------------------------
550 void axisExtractor02::indextoreal(double a[3], double b[3] )
554 // EED 15 Mai 2007 .NET
555 // minspac=min(espprin[0],min(espprin[1],espprin[2]));
558 if (espprin[1]<minspac) { minspac=espprin[1]; }
559 if (espprin[2]<minspac) { minspac=espprin[2]; }
568 //----------------------------------------------------------------------------
569 double axisExtractor02::distanciaejepunto(double point[3], double point2[3], double point3[3])
575 double v[3], v2[3], v3[3];
580 v[0]=point3[0]-point2[0];
581 v[1]=point3[1]-point2[1];
582 v[2]=point3[2]-point2[2];
584 v3[0]=point[0]-point2[0];
585 v3[1]=point[1]-point2[1];
586 v3[2]=point[2]-point2[2];
588 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]));
591 point4[0]=point2[0]+(x*v[0]);
592 point4[1]=point2[1]+(x*v[1]);
593 point4[2]=point2[2]+(x*v[2]);
595 v2[0]=point4[0]-point[0];
596 v2[1]=point4[1]-point[1];
597 v2[2]=point4[2]-point[2];
599 y=(v[0]*v2[0])+(v[1]*v2[1])+(v[2]*v2[2]);
601 dist=sqrt((v2[0]*v2[0])+(v2[1]*v2[1])+(v2[2]*v2[2]));
609 //----------------------------------------------------------------------------
610 double axisExtractor02::proporcioejepunto(double point[3], double point2[3], double point3[3])
618 v[0]=point3[0]-point2[0];
619 v[1]=point3[1]-point2[1];
620 v[2]=point3[2]-point2[2];
622 v3[0]=point[0]-point2[0];
623 v3[1]=point[1]-point2[1];
624 v3[2]=point[2]-point2[2];
626 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]));
642 //----------------------------------------------------------------------------
643 void axisExtractor02::searc(int i, int j, int k, vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
654 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
655 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
659 data->GetExtent(ext);
661 radio=(ext[1]-ext[0])/2;
668 vector[label-1][0]+=1;
669 vector[label-1][1]+=i;
670 vector[label-1][2]+=j;
671 vector[label-1][3]+=k;
673 for(i3=-1;i3<=1;i3++){
674 for(j3=-1;j3<=1;j3++){
675 for(k3=-1;k3<=1;k3++){
676 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] ){
677 ptr=(unsigned char *) data->GetScalarPointer(i+i3,j+j3,k+k3);
678 ptr2=(unsigned char *) data2->GetScalarPointer(i+i3,j+j3,k+k3);
679 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)))){
680 searc(i+i3, j+j3, k+k3, data, data2, label, vector);
692 //----------------------------------------------------------------------------
693 void axisExtractor02::searcb(int i, int j, int k, vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
704 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
705 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
709 data->GetExtent(ext);
711 radio=(ext[1]-ext[0])/2;
718 vector[label-1][0]+=1;
719 vector[label-1][1]+=i;
720 vector[label-1][2]+=j;
721 vector[label-1][3]+=k;
723 for(i3=-1;i3<=1;i3++){
724 for(j3=-1;j3<=1;j3++){
725 for(k3=-1;k3<=1;k3++){
726 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] ){
727 ptr=(unsigned char *) data->GetScalarPointer(i+i3,j+j3,k+k3);
728 ptr2=(unsigned char *) data2->GetScalarPointer(i+i3,j+j3,k+k3);
729 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)))){
730 searcb(i+i3, j+j3, k+k3, data, data2, label, vector);
742 //----------------------------------------------------------------------------
743 unsigned char axisExtractor02::find_components(vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
749 data->GetExtent(ext);
750 double *ff= data->GetOrigin();
756 radio=(ext[1]-ext[0])/2;
758 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
759 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
760 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
761 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
762 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
764 if(*ptr==255 && *ptr2==0 && (radio*radio)>=((i2*i2)+(j2*j2)+(k2*k2)) && ((radio-1)*(radio-1))<=((i2*i2)+(j2*j2)+(k2*k2)) ){
766 vector[label-1][0]=0;
767 vector[label-1][1]=0;
768 vector[label-1][2]=0;
769 vector[label-1][3]=0;
770 searc(i, j, k, data, data2, label, vector);
783 //----------------------------------------------------------------------------
784 unsigned char axisExtractor02::find_componentsb(vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
790 data->GetExtent(ext);
795 radio=(ext[1]-ext[0])/2;
797 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
798 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
799 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
800 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
801 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
803 if(*ptr==0 && *ptr2==0 && (radio*radio)>=((i2*i2)+(j2*j2)+(k2*k2)) && ((radio-1)*(radio-1))<=((i2*i2)+(j2*j2)+(k2*k2)) ){
805 vector[label-1][0]=0;
806 vector[label-1][1]=0;
807 vector[label-1][2]=0;
808 vector[label-1][3]=0;
809 searcb(i, j, k, data, data2, label, vector);
823 //----------------------------------------------------------------------------
824 int axisExtractor02::proporcion(vtkImageData *data )
831 data->GetExtent(ext);
835 for(i=ext[0];i<=ext[1];i++){
836 for(j=ext[2];j<=ext[3];j++){
837 for(k=ext[4];k<=ext[5];k++){
838 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
847 return (max*100)/total;
857 //----------------------------------------------------------------------------
858 bool axisExtractor02::border(vtkImageData *data, int p1[3] )
869 data->GetExtent(ext);
871 centro[0]=(ext[1]+ext[0])/2;
872 centro[1]=(ext[3]+ext[2])/2;
873 centro[2]=(ext[5]+ext[4])/2;
877 ptr=(unsigned char *) data->GetScalarPointer(p1);
891 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] ){
892 ptr=(unsigned char *) data->GetScalarPointer(p2);
907 /*if(p1[0]==centro[0] && p1[1]==centro[1] && p1[2]==centro[2]){
919 //----------------------------------------------------------------------------
920 void axisExtractor02::optim(vtkImageData *data, vtkImageData *data2 )
930 double inerciaxx, inerciayy, inerciazz, inerciaxy, inerciayz, inerciazx;
931 double sumx, sumy, sumz;
932 double sumxx, sumyy, sumzz, sumxy, sumyz, sumzx;
934 double sumix, sumiy, sumiz;
935 double sumixx, sumiyy, sumizz, sumixy, sumiyz, sumizx;
937 double inerciaixx, inerciaiyy, inerciaizz, inerciaixy, inerciaiyz, inerciaizx;
939 double inerciaxx2, inerciayy2, inerciazz2, inerciaxy2, inerciayz2, inerciazx2;
940 double sumx2, sumy2, sumz2;
941 double sumxx2, sumyy2, sumzz2, sumxy2, sumyz2, sumzx2;
943 double sumix2, sumiy2, sumiz2;
944 double sumixx2, sumiyy2, sumizz2, sumixy2, sumiyz2, sumizx2;
946 //double sumixp, sumiyp, sumizp;
947 //double sumix2p, sumiy2p, sumiz2p;
949 double inerciaxxp, inerciayyp, inerciazzp, inerciaxyp, inerciayzp, inerciazxp;
950 double sumxp, sumyp, sumzp;
951 double sumxxp, sumyyp, sumzzp, sumxyp, sumyzp, sumzxp;
953 double inerciaxx2p, inerciayy2p, inerciazz2p, inerciaxy2p, inerciayz2p, inerciazx2p;
954 double sumx2p, sumy2p, sumz2p;
955 double sumxx2p, sumyy2p, sumzz2p, sumxy2p, sumyz2p, sumzx2p;
957 //double inerciaixx2p, inerciaiyy2p, inerciaizz2p, inerciaixy2p, inerciaiyz2p, inerciaizx2p;
965 unsigned long cant2p;
977 double im, jm, km, lm;
1048 data->GetExtent(ext);
1050 unsigned short *ptr;
1051 unsigned char *ptr2;
1053 radio=(ext[1]-ext[0])/2;
1055 centro[0]=(ext[1]+ext[0])/2;
1056 centro[1]=(ext[3]+ext[2])/2;
1057 centro[2]=(ext[5]+ext[4])/2;
1061 for(i=ext[0];i<=ext[1];i++){
1062 for(j=ext[2];j<=ext[3];j++){
1063 for(k=ext[4];k<=ext[5];k++){
1064 tmpsqrt=((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
1065 if( radio>=sqrt( tmpsqrt) ){
1066 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
1080 for(i=ext[0];i<=ext[1];i++){
1081 for(j=ext[2];j<=ext[3];j++){
1082 for(k=ext[4];k<=ext[5];k++){
1083 tmpsqrt=((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
1084 if( radio>=sqrt(tmpsqrt) ){
1085 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
1086 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
1091 im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
1092 jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
1093 km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
1094 lm=2*(((double)*ptr-(double)min)/((double)max-(double)min));
1155 centi=(double)sumi/(double)cant;
1156 inerciaii= ((double)sumii/(double)cant)-(centi*centi);
1158 centi2=(double)sumi2/(double)cant2;
1159 inerciaii2= ((double)sumii2/(double)cant2)-(centi2*centi2);
1161 centx=(double)sumx/(double)cant;
1162 centy=(double)sumy/(double)cant;
1163 centz=(double)sumz/(double)cant;
1165 centx2=(double)sumx2/(double)cant2;
1166 centy2=(double)sumy2/(double)cant2;
1167 centz2=(double)sumz2/(double)cant2;
1169 inerciaxx= ((double)sumxx/(double)cant)-centx*centx;
1170 inerciayy= ((double)sumyy/(double)cant)-centy*centy;
1171 inerciazz= ((double)sumzz/(double)cant)-centz*centz;
1173 inerciaxx2= ((double)sumxx2/(double)cant2)-centx2*centx2;
1174 inerciayy2= ((double)sumyy2/(double)cant2)-centy2*centy2;
1175 inerciazz2= ((double)sumzz2/(double)cant2)-centz2*centz2;
1177 inerciaxy= ((double)sumxy/(double)cant)-centx*centy;
1178 inerciayz= ((double)sumyz/(double)cant)-centy*centz;
1179 inerciazx= ((double)sumzx/(double)cant)-centz*centx;
1181 inerciaxy2= ((double)sumxy2/(double)cant2)-centx2*centy2;
1182 inerciayz2= ((double)sumyz2/(double)cant2)-centy2*centz2;
1183 inerciazx2= ((double)sumzx2/(double)cant2)-centz2*centx2;
1186 centix=((double)sumix+(double)sumix2)/((double)sumi+(double)sumi2);
1187 centiy=((double)sumiy+(double)sumiy2)/((double)sumi+(double)sumi2);
1188 centiz=((double)sumiz+(double)sumiz2)/((double)sumi+(double)sumi2);
1190 inerciaixx= ( ((double)sumixx+(double)sumixx2) / ((double)sumi+(double)sumi2) )-(centix*centix);
1191 inerciaiyy= ( ((double)sumiyy+(double)sumiyy2) / ((double)sumi+(double)sumi2) )-(centiy*centiy);
1192 inerciaizz= ( ((double)sumizz+(double)sumizz2) / ((double)sumi+(double)sumi2) )-(centiz*centiz);
1195 inerciaixy= ( ((double)sumixy+(double)sumixy2) / ((double)sumi+(double)sumi2) )-(centix*centiy);
1196 inerciaiyz= ( ((double)sumiyz+(double)sumiyz2) / ((double)sumi+(double)sumi2) )-(centiy*centiz);
1197 inerciaizx= ( ((double)sumizx+(double)sumizx2) / ((double)sumi+(double)sumi2) )-(centiz*centix);
1203 A[0][1]=A[1][0]=inerciaxy;
1204 A[0][2]=A[2][0]=inerciazx;
1205 A[1][2]=A[2][1]=inerciayz;
1207 A2[0][0]=inerciaxx2;
1208 A2[1][1]=inerciayy2;
1209 A2[2][2]=inerciazz2;
1210 A2[0][1]=A2[1][0]=inerciaxy2;
1211 A2[0][2]=A2[2][0]=inerciazx2;
1212 A2[1][2]=A2[2][1]=inerciayz2;
1214 Ai[0][0]=inerciaixx;
1215 Ai[1][1]=inerciaiyy;
1216 Ai[2][2]=inerciaizz;
1217 Ai[0][1]=Ai[1][0]=inerciaixy;
1218 Ai[0][2]=Ai[2][0]=inerciaizx;
1219 Ai[1][2]=Ai[2][1]=inerciaiyz;
1221 vtkMath::Diagonalize3x3 (A, w, V);
1222 vtkMath::Diagonalize3x3 (A2, w2, V2);
1223 vtkMath::Diagonalize3x3 (Ai, wi, Vi);
1231 inerciary=w[0]+w[2];
1232 inerciarz=w[0]+w[1];
1235 inerciary=w[0]+w[1];
1236 inerciarz=w[0]+w[2];
1243 inerciary=w[0]+w[2];
1244 inerciarz=w[1]+w[2];
1247 inerciary=w[1]+w[2];
1248 inerciarz=w[0]+w[2];
1257 inerciary=w[1]+w[0];
1258 inerciarz=w[1]+w[2];
1261 inerciary=w[1]+w[2];
1262 inerciarz=w[1]+w[0];
1269 inerciary=w[0]+w[2];
1270 inerciarz=w[1]+w[2];
1273 inerciary=w[1]+w[2];
1274 inerciarz=w[0]+w[2];
1283 inercia2r=w2[1]+w2[2];
1287 inercia2r=w2[1]+w2[0];
1293 inercia2r=w2[2]+w2[0];
1297 inercia2r=w2[1]+w2[0];
1305 inerciari=wi[1]+wi[2];
1307 inerciariy=wi[0]+wi[2];
1308 inerciariz=wi[0]+wi[1];
1311 inerciariy=wi[0]+wi[1];
1312 inerciariz=wi[0]+wi[2];
1317 inerciari=wi[1]+wi[0];
1319 inerciariy=wi[0]+wi[2];
1320 inerciariz=wi[1]+wi[2];
1323 inerciariy=wi[1]+wi[2];
1324 inerciariz=wi[0]+wi[2];
1332 inerciari=wi[2]+wi[0];
1334 inerciariy=wi[1]+wi[0];
1335 inerciariz=wi[1]+wi[2];
1338 inerciariy=wi[1]+wi[2];
1339 inerciariz=wi[1]+wi[0];
1344 inerciari=wi[1]+wi[0];
1346 inerciariy=wi[0]+wi[2];
1347 inerciariz=wi[1]+wi[2];
1350 inerciariy=wi[1]+wi[2];
1351 inerciariz=wi[0]+wi[2];
1358 for(i=ext[0];i<=ext[1];i++){
1359 for(j=ext[2];j<=ext[3];j++){
1360 for(k=ext[4];k<=ext[5];k++){
1361 if( radio>=sqrt(((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]))) ){
1362 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
1363 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
1366 im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
1367 jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
1368 km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
1374 point2[0]=centx+V[0][ejemin];
1375 point2[1]=centy+V[1][ejemin];
1376 point2[2]=centz+V[2][ejemin];
1378 point3[0]=centx-V[0][ejemin];
1379 point3[1]=centy-V[1][ejemin];
1380 point3[2]=centz-V[2][ejemin];
1382 //distcent=1+sqrt(centx*centx+centy*centy+centz*centz);
1387 sumitotal+=(distcent-fabs(distanciaejepunto(point1, point2, point3)))*(distcent-fabs(distanciaejepunto(point1, point2, point3)));
1395 inercia2r=sumitotal/cant2;
1401 /*printf("centx: %f\n", centx);
1402 printf("centy: %f\n", centy);
1403 printf("centz: %f\n", centz);
1405 printf("w[0]: %f\n", w[0]);
1406 printf("w[1]: %f\n", w[1]);
1407 printf("w[2]: %f\n", w[2]);
1408 printf("inerciar: %f\n", inerciar);
1410 printf("w[ejemin]: %f\n", w[ejemin]);
1412 printf("w[ejemin]/inerciar: %f\n", w[ejemin]/inerciar);
1414 printf("V[0][0]: %f\n", V[0][0]);
1415 printf("V[1][0]: %f\n", V[1][0]);
1416 printf("V[2][0]: %f\n", V[2][0]);
1418 printf("V[0][1]: %f\n", V[0][1]);
1419 printf("V[1][1]: %f\n", V[1][1]);
1420 printf("V[2][1]: %f\n", V[2][1]);
1422 printf("V[0][2]: %f\n", V[0][2]);
1423 printf("V[1][2]: %f\n", V[1][2]);
1424 printf("V[2][2]: %f\n", V[2][2]);*/
1426 w0=(double)cant/((double)cant+(double)cant2);
1427 w1=(double)cant2/((double)cant+(double)cant2);
1431 ut=w0*centi+w1*centi2;
1432 intervar=w0*inerciaii+w1*inerciaii2;
1440 //costo=w0*inerciaii+w1*inerciaii2+param*inerciar;
1441 //costo=param*(w0*param3*inerciaii+w1*(1-param3)*inerciaii2)+(1-param)*(w0*param2*inerciar+w1*(1-param2)*inercia2r);
1442 costo=param*w0*inerciaii+w1*param2*inerciaii2+w0*param3*inerciar+w1*param4*inercia2r;
1448 while(changes>3 && total<500){
1451 for(i=ext[0];i<=ext[1];i++){
1452 for(j=ext[2];j<=ext[3];j++){
1453 for(k=ext[4];k<=ext[5];k++){
1454 tmpsqrt = ((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
1455 if( radio>=sqrt(tmpsqrt) ){
1456 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
1457 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
1462 im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
1463 jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
1464 km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
1465 lm=2*(((double)*ptr-(double)min)/((double)max-(double)min));
1467 if(border(data2, point)){
1476 sumii2p=sumii2+lm*lm;
1493 sumxx2p=sumxx2+im*im;
1494 sumyy2p=sumyy2+jm*jm;
1495 sumzz2p=sumzz2+km*km;
1496 sumxy2p=sumxy2+im*jm;
1497 sumyz2p=sumyz2+jm*km;
1498 sumzx2p=sumzx2+km*im;
1500 centip=(double)sumip/(double)cantp;
1501 inerciaiip= ((double)sumiip/(double)cantp)-(centip*centip);
1503 centi2p=(double)sumi2p/(double)cant2p;
1504 inerciaii2p= ((double)sumii2p/(double)cant2p)-(centi2p*centi2p);
1506 centxp=(double)sumxp/(double)cantp;
1507 centyp=(double)sumyp/(double)cantp;
1508 centzp=(double)sumzp/(double)cantp;
1510 centx2p=(double)sumx2p/(double)cant2p;
1511 centy2p=(double)sumy2p/(double)cant2p;
1512 centz2p=(double)sumz2p/(double)cant2p;
1514 inerciaxxp= ((double)sumxxp/(double)cantp)-centxp*centxp;
1515 inerciayyp= ((double)sumyyp/(double)cantp)-centyp*centyp;
1516 inerciazzp= ((double)sumzzp/(double)cantp)-centzp*centzp;
1518 inerciaxx2p= ((double)sumxx2p/(double)cant2p)-centx2p*centx2p;
1519 inerciayy2p= ((double)sumyy2p/(double)cant2p)-centy2p*centy2p;
1520 inerciazz2p= ((double)sumzz2p/(double)cant2p)-centz2p*centz2p;
1522 inerciaxyp= ((double)sumxyp/(double)cantp)-centxp*centyp;
1523 inerciayzp= ((double)sumyzp/(double)cantp)-centyp*centzp;
1524 inerciazxp= ((double)sumzxp/(double)cantp)-centzp*centxp;
1526 inerciaxy2p= ((double)sumxy2p/(double)cant2p)-centx2p*centy2p;
1527 inerciayz2p= ((double)sumyz2p/(double)cant2p)-centy2p*centz2p;
1528 inerciazx2p= ((double)sumzx2p/(double)cant2p)-centz2p*centx2p;
1530 Ap[0][0]=inerciaxxp;
1531 Ap[1][1]=inerciayyp;
1532 Ap[2][2]=inerciazzp;
1533 Ap[0][1]=Ap[1][0]=inerciaxyp;
1534 Ap[0][2]=Ap[2][0]=inerciazxp;
1535 Ap[1][2]=Ap[2][1]=inerciayzp;
1538 A2p[0][0]=inerciaxx2p;
1539 A2p[1][1]=inerciayy2p;
1540 A2p[2][2]=inerciazz2p;
1541 A2p[0][1]=A2p[1][0]=inerciaxy2p;
1542 A2p[0][2]=A2p[2][0]=inerciazx2p;
1543 A2p[1][2]=A2p[2][1]=inerciayz2p;
1547 w0p=(double)cantp/((double)cantp+(double)cant2p);
1548 w1p=(double)cant2p/((double)cantp+(double)cant2p);
1553 vtkMath::Diagonalize3x3 (Ap, wp, Vp);
1554 vtkMath::Diagonalize3x3 (A2p, w2p, V2p);
1560 inerciarp=wp[2]+wp[1];
1562 inerciarpy=wp[0]+wp[2];
1563 inerciarpz=wp[0]+wp[1];
1566 inerciarpy=wp[0]+wp[1];
1567 inerciarpz=wp[0]+wp[2];
1572 inerciarp=wp[1]+wp[0];
1574 inerciarpy=wp[0]+wp[2];
1575 inerciarpz=wp[1]+wp[2];
1578 inerciarpy=wp[1]+wp[2];
1579 inerciarpz=wp[0]+wp[2];
1586 inerciarp=wp[2]+wp[0];
1588 inerciarpy=wp[1]+wp[0];
1589 inerciarpz=wp[1]+wp[2];
1592 inerciarpy=wp[1]+wp[2];
1593 inerciarpz=wp[1]+wp[0];
1598 inerciarp=wp[1]+wp[0];
1600 inerciarpy=wp[0]+wp[2];
1601 inerciarpz=wp[1]+wp[2];
1604 inerciarpy=wp[1]+wp[2];
1605 inerciarpz=wp[0]+wp[2];
1614 inercia2rp=w2p[2]+w2p[1];
1618 inercia2rp=w2p[1]+w2p[0];
1624 inercia2rp=w2p[2]+w2p[0];
1628 inercia2rp=w2p[1]+w2p[0];
1636 point2[0]=centxp+Vp[0][ejeminp];
1637 point2[1]=centyp+Vp[1][ejeminp];
1638 point2[2]=centzp+Vp[2][ejeminp];
1640 point3[0]=centxp-Vp[0][ejeminp];
1641 point3[1]=centyp-Vp[1][ejeminp];
1642 point3[2]=centzp-Vp[2][ejeminp];*/
1644 //distcent=1+sqrt(centx*centx+centy*centy+centz*centz);
1647 sumitotalp=sumitotal+(distcent-fabs(distanciaejepunto(point1, point2, point3)))*(distcent-fabs(distanciaejepunto(point1, point2, point3)));
1649 inercia2rp=sumitotalp/cant2p;*/
1651 //costop=w0p*inerciaiip+w1p*inerciaii2p+param*inerciarp;
1653 //costop=param*(w0p*param3*inerciaiip+w1p*(1-param3)*inerciaii2p)+(1-param)*(w0p*param2*inerciarp+w1p*(1-param2)*inercia2rp);
1656 costop=param*w0p*inerciaiip+w1p*param2*inerciaii2p+w0p*param3*inerciarp+w1p*param4*inercia2rp;
1658 if(costop<costo && centip>centi2p){
1692 inerciaii= inerciaiip;
1695 inerciaii2p= inerciaii2p;
1705 inerciaxx= inerciaxxp;
1706 inerciayy= inerciayyp;
1707 inerciazz= inerciazzp;
1709 inerciaxx2= inerciaxx2p;
1710 inerciayy2= inerciayy2p;
1711 inerciazz2= inerciazz2p;
1713 inerciaxy= inerciaxyp;
1714 inerciayz= inerciayzp;
1715 inerciazx= inerciazxp;
1717 inerciaxy2= inerciaxy2p;
1718 inerciayz2= inerciayz2p;
1719 inerciazx2= inerciazx2p;
1721 // sumitotal=sumitotalp;
1734 sumii2p=sumii2-lm*lm;
1751 sumxx2p=sumxx2-im*im;
1752 sumyy2p=sumyy2-jm*jm;
1753 sumzz2p=sumzz2-km*km;
1754 sumxy2p=sumxy2-im*jm;
1755 sumyz2p=sumyz2-jm*km;
1756 sumzx2p=sumzx2-km*im;
1758 centip=(double)sumip/(double)cantp;
1759 inerciaiip= ((double)sumiip/(double)cantp)-(centip*centip);
1761 centi2p=(double)sumi2p/(double)cant2p;
1762 inerciaii2p= ((double)sumii2p/(double)cant2p)-(centi2p*centi2p);
1764 centxp=(double)sumxp/(double)cantp;
1765 centyp=(double)sumyp/(double)cantp;
1766 centzp=(double)sumzp/(double)cantp;
1768 centx2p=(double)sumx2p/(double)cant2p;
1769 centy2p=(double)sumy2p/(double)cant2p;
1770 centz2p=(double)sumz2p/(double)cant2p;
1773 inerciaxxp= ((double)sumxxp/(double)cantp)-centxp*centxp;
1774 inerciayyp= ((double)sumyyp/(double)cantp)-centyp*centyp;
1775 inerciazzp= ((double)sumzzp/(double)cantp)-centzp*centzp;
1777 inerciaxx2p= ((double)sumxx2p/(double)cant2p)-centx2p*centx2p;
1778 inerciayy2p= ((double)sumyy2p/(double)cant2p)-centy2p*centy2p;
1779 inerciazz2p= ((double)sumzz2p/(double)cant2p)-centz2p*centz2p;
1781 inerciaxyp= ((double)sumxyp/(double)cantp)-centxp*centyp;
1782 inerciayzp= ((double)sumyzp/(double)cantp)-centyp*centzp;
1783 inerciazxp= ((double)sumzxp/(double)cantp)-centzp*centxp;
1785 inerciaxy2p= ((double)sumxy2p/(double)cant2p)-centx2p*centy2p;
1786 inerciayz2p= ((double)sumyz2p/(double)cant2p)-centy2p*centz2p;
1787 inerciazx2p= ((double)sumzx2p/(double)cant2p)-centz2p*centx2p;
1789 Ap[0][0]=inerciaxxp;
1790 Ap[1][1]=inerciayyp;
1791 Ap[2][2]=inerciazzp;
1792 Ap[0][1]=Ap[1][0]=inerciaxyp;
1793 Ap[0][2]=Ap[2][0]=inerciazxp;
1794 Ap[1][2]=Ap[2][1]=inerciayzp;
1796 A2p[0][0]=inerciaxx2p;
1797 A2p[1][1]=inerciayy2p;
1798 A2p[2][2]=inerciazz2p;
1799 A2p[0][1]=A2p[1][0]=inerciaxy2p;
1800 A2p[0][2]=A2p[2][0]=inerciazx2p;
1801 A2p[1][2]=A2p[2][1]=inerciayz2p;
1803 w0p=(double)cantp/((double)cantp+(double)cant2p);
1804 w1p=(double)cant2p/((double)cantp+(double)cant2p);
1809 vtkMath::Diagonalize3x3 (Ap, wp, Vp);
1810 vtkMath::Diagonalize3x3 (A2p, w2p, V2p);
1815 inerciarp=wp[2]+wp[1];
1817 inerciarpy=wp[0]+wp[2];
1818 inerciarpz=wp[0]+wp[1];
1821 inerciarpy=wp[0]+wp[1];
1822 inerciarpz=wp[0]+wp[2];
1827 inerciarp=wp[1]+wp[0];
1829 inerciarpy=wp[0]+wp[2];
1830 inerciarpz=wp[1]+wp[2];
1833 inerciarpy=wp[1]+wp[2];
1834 inerciarpz=wp[0]+wp[2];
1841 inerciarp=wp[2]+wp[0];
1843 inerciarpy=wp[1]+wp[0];
1844 inerciarpz=wp[1]+wp[2];
1847 inerciarpy=wp[1]+wp[2];
1848 inerciarpz=wp[1]+wp[0];
1853 inerciarp=wp[1]+wp[0];
1855 inerciarpy=wp[0]+wp[2];
1856 inerciarpz=wp[1]+wp[2];
1859 inerciarpy=wp[1]+wp[2];
1860 inerciarpz=wp[0]+wp[2];
1869 inercia2rp=w2p[2]+w2p[1];
1873 inercia2rp=w2p[1]+w2p[0];
1879 inercia2rp=w2p[2]+w2p[0];
1883 inercia2rp=w2p[1]+w2p[0];
1892 point2[0]=centxp+Vp[0][ejeminp];
1893 point2[1]=centyp+Vp[1][ejeminp];
1894 point2[2]=centzp+Vp[2][ejeminp];
1896 point3[0]=centxp-Vp[0][ejeminp];
1897 point3[1]=centyp-Vp[1][ejeminp];
1898 point3[2]=centzp-Vp[2][ejeminp];*/
1900 //distcent=1+sqrt(centx*centx+centy*centy+centz*centz);
1903 sumitotalp=sumitotal-(distcent-fabs(distanciaejepunto(point1, point2, point3)))*(distcent-fabs(distanciaejepunto(point1, point2, point3)));
1905 inercia2rp=sumitotalp/cant2p;*/
1908 //costop=w0p*inerciaiip+w1p*inerciaii2p+param*inerciarp;
1909 //costop=param*(w0p*param3*inerciaiip+w1p*(1-param3)*inerciaii2p)+(1-param)*(w0p*param2*inerciarp+w1p*(1-param2)*inercia2rp);
1912 costop=param*w0p*inerciaiip+w1p*param2*inerciaii2p+w0p*param3*inerciarp+w1p*param4*inercia2rp;
1914 if(costop<costo && centip>centi2p){
1948 inerciaii= inerciaiip;
1951 inerciaii2p= inerciaii2p;
1961 inerciaxx= inerciaxxp;
1962 inerciayy= inerciayyp;
1963 inerciazz= inerciazzp;
1965 inerciaxx2= inerciaxx2p;
1966 inerciayy2= inerciayy2p;
1967 inerciazz2= inerciazz2p;
1969 inerciaxy= inerciaxyp;
1970 inerciayz= inerciayzp;
1971 inerciazx= inerciazxp;
1973 inerciaxy2= inerciaxy2p;
1974 inerciayz2= inerciayz2p;
1975 inerciazx2= inerciazx2p;
1977 // sumitotal=sumitotalp;
1990 // printf("changes: %d\n", changes);
1995 //printf("total: %d\n", total);
2004 //----------------------------------------------------------------------------
2005 void axisExtractor02::costominimo(vtkImageData *data, vtkImageData *data2 )
2012 data->GetExtent(ext);
2015 unsigned char *ptr2;
2019 radio=(ext[1]-ext[0])/2;
2021 centro[0]=(ext[1]+ext[0])/2;
2022 centro[1]=(ext[3]+ext[2])/2;
2023 centro[2]=(ext[5]+ext[4])/2;
2030 for(i=ext[0];i<=ext[1];i++){
2031 for(j=ext[2];j<=ext[3];j++){
2032 for(k=ext[4];k<=ext[5];k++){
2033 ptr=(double *)data->GetScalarPointer(i,j,k);
2034 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
2039 if(minis[ind-1]<=*ptr && *ptr!=0 && ind<11){
2055 //----------------------------------------------------------------------------
2056 void axisExtractor02::costominimo2(vtkImageData *data, vtkImageData *data3, int p1[3], int p2[3], int p3[3])
2069 unsigned char *ptr3;
2070 unsigned char *ptr4;
2072 unsigned char *ptr6;
2076 vtkImageData *data4;
2078 data4=vtkImageData::New();
2079 data4->SetScalarTypeToUnsignedChar();
2080 data4->SetExtent(data->GetExtent());
2081 data4->SetSpacing(data->GetSpacing());
2083 vtkImageData *data6;
2085 data6=vtkImageData::New();
2086 data6->SetScalarTypeToUnsignedChar();
2087 data6->SetExtent(data->GetExtent());
2088 data6->SetSpacing(data->GetSpacing());
2091 vtkImageData *data2;
2093 data2=vtkImageData::New();
2094 data2->SetScalarTypeToDouble();
2095 data2->SetExtent(data->GetExtent());
2096 data2->SetSpacing(data->GetSpacing());
2098 vtkImageData *data8;
2100 data8=vtkImageData::New();
2101 data8->SetScalarTypeToDouble();
2102 data8->SetExtent(data->GetExtent());
2103 data8->SetSpacing(data->GetSpacing());
2106 double mincst, mincstb;
2112 data->GetExtent(ext);
2114 radio=(ext[1]-ext[0])/2;
2116 for(i=ext[0];i<=ext[1];i++){
2117 for(j=ext[2];j<=ext[3];j++){
2118 for(k=ext[4];k<=ext[5];k++){
2119 ptr2=(double *)data2->GetScalarPointer(i,j,k);
2120 ptr4=(unsigned char *)data4->GetScalarPointer(i,j,k);
2121 ptr6=(unsigned char *)data6->GetScalarPointer(i,j,k);
2122 ptr8=(double *)data8->GetScalarPointer(i,j,k);
2142 ptr2=(double *)data2->GetScalarPointer(p1[0],p1[1],p1[2]);
2145 ptr4=(unsigned char *)data4->GetScalarPointer(p1[0],p1[1],p1[2]);
2148 ptr8=(double *)data8->GetScalarPointer(p2[0],p2[1],p2[2]);
2151 ptr6=(unsigned char *)data6->GetScalarPointer(p2[0],p2[1],p2[2]);
2161 ptr5=(double *)data2->GetScalarPointer(punto[0],punto[1],punto[2]);
2162 ptr7=(double *)data8->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
2167 if(i!=0 || j!=0 || k!=0){
2168 punto2[0]=punto[0]+i;
2169 punto2[1]=punto[1]+j;
2170 punto2[2]=punto[2]+k;
2172 puntob2[0]=puntob[0]+i;
2173 puntob2[1]=puntob[1]+j;
2174 puntob2[2]=puntob[2]+k;
2176 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] ){
2178 ptr=(double *)data->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
2179 ptr2=(double *)data2->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
2180 ptr3=(unsigned char *)data3->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
2181 ptr4=(unsigned char *)data4->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
2184 if(*ptr3!=0 && *ptr4!=2){
2185 costoactual=*ptr5+(1/(*ptr));
2186 if(*ptr2>costoactual){
2193 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] ){
2195 ptr=(double *)data->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
2196 ptr8=(double *)data8->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
2197 ptr3=(unsigned char *)data3->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
2198 ptr6=(unsigned char *)data6->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
2200 if(*ptr3!=0 && *ptr6!=2){
2201 costoactual=*ptr7+(1/(*ptr));
2202 if(*ptr8>costoactual){
2213 ptr4=(unsigned char *)data4->GetScalarPointer(punto[0],punto[1],punto[2]);
2216 ptr6=(unsigned char *)data6->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
2222 for(i=ext[0];i<=ext[1];i++){
2223 for(j=ext[2];j<=ext[3];j++){
2224 for(k=ext[4];k<=ext[5];k++){
2225 ptr2=(double *)data2->GetScalarPointer(i,j,k);
2226 ptr8=(double *)data8->GetScalarPointer(i,j,k);
2227 ptr4=(unsigned char *)data4->GetScalarPointer(i,j,k);
2228 ptr6=(unsigned char *)data6->GetScalarPointer(i,j,k);
2251 punto[0]=puntomin[0];
2252 punto[1]=puntomin[1];
2253 punto[2]=puntomin[2];
2255 puntob[0]=puntobmin[0];
2256 puntob[1]=puntobmin[1];
2257 puntob[2]=puntobmin[2];
2259 ptr4=(unsigned char *)data4->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
2260 ptr6=(unsigned char *)data6->GetScalarPointer(punto[0],punto[1],punto[2]);
2288 //----------------------------------------------------------------------------
2289 void axisExtractor02::invertir(vtkImageData *data )
2295 data->GetExtent(ext);
2299 radio=(ext[1]-ext[0])/2;
2303 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2304 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2305 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2306 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2307 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
2308 if( radio<sqrt(tmpsqrt) ){
2323 //----------------------------------------------------------------------------
2324 void axisExtractor02::redondear(vtkImageData *data )
2330 data->GetExtent(ext);
2334 radio=(ext[1]-ext[0])/2;
2338 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2339 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2340 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2341 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2342 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
2343 if( radio<sqrt(tmpsqrt) ){
2346 /*if( i2==0 && j2==0 && k2==0 ){
2358 //----------------------------------------------------------------------------
2359 void axisExtractor02::redondear2(vtkImageData *data )
2365 data->GetExtent(ext);
2369 radio=(ext[1]-ext[0])/2;
2373 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2374 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2375 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2376 ptr=(double *) data->GetScalarPointer(i,j,k);
2377 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
2378 if( radio<sqrt(tmpsqrt) ){
2391 //----------------------------------------------------------------------------
2392 void axisExtractor02::redondear3(vtkImageData *data )
2398 data->GetExtent(ext);
2402 radio=(ext[1]-ext[0])/2;
2406 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2407 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2408 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2409 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2410 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
2411 if( radio<sqrt(tmpsqrt) ){
2424 //----------------------------------------------------------------------------
2425 double axisExtractor02::distancia(double a[3], double b[3] )
2427 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]));
2428 return sqrt(tmpsqrt);
2435 //----------------------------------------------------------------------------
2436 double axisExtractor02::distancia(int a[3], int b[3] )
2438 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]));
2439 return sqrt(tmpsqrt);
2448 //----------------------------------------------------------------------------
2449 void axisExtractor02::blanquear(vtkImageData *data )
2454 data->GetExtent(ext);
2458 for(i=ext[0];i<=ext[1];i++){
2459 for(j=ext[2];j<=ext[3];j++){
2460 for(k=ext[4];k<=ext[5];k++){
2461 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2470 //----------------------------------------------------------------------------
2471 void axisExtractor02::blanquear3(vtkImageData *data )
2476 data->GetExtent(ext);
2480 for(i=ext[0];i<=ext[1];i++){
2481 for(j=ext[2];j<=ext[3];j++){
2482 for(k=ext[4];k<=ext[5];k++){
2483 ptr=(double *) data->GetScalarPointer(i,j,k);
2493 //----------------------------------------------------------------------------
2494 void axisExtractor02::blanquear2(vtkImageData *data )
2500 data->GetExtent(ext);
2502 centro[0]=(ext[1]+ext[0])/2;
2503 centro[1]=(ext[3]+ext[2])/2;
2504 centro[2]=(ext[5]+ext[4])/2;
2508 for(i=ext[0];i<=ext[1];i++){
2509 for(j=ext[2];j<=ext[3];j++){
2510 for(k=ext[4];k<=ext[5];k++){
2511 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2512 if(i==centro[0] && j==centro[1] && k==centro[2]){
2526 //----------------------------------------------------------------------------
2527 void axisExtractor02::cilindro(vtkImageData *data, double vector[3] )
2540 data->GetExtent(ext);
2544 radio=(ext[1]-ext[0])/2;
2545 radiof=espprin[0]*(((double)ext[1]-(double)ext[0])/4);
2546 centro[0]=(ext[1]+ext[0])/2;
2547 centro[1]=(ext[3]+ext[2])/2;
2548 centro[2]=(ext[5]+ext[4])/2;
2550 indextoreal(centro, centrof );
2552 if(vector[0]==0 && vector[1]==0 && vector[2]==0){
2556 puntof2[0]=centrof[0]+vector[0];
2557 puntof2[1]=centrof[1]+vector[1];
2558 puntof2[2]=centrof[2]+vector[2];
2561 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2562 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2563 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2564 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2568 indextoreal(punto, puntof );
2569 tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2);
2570 if( radio<sqrt(tmpsqrt) ){
2574 else if(radiof<distanciaejepunto(puntof, centrof, puntof2)){
2590 //----------------------------------------------------------------------------
2591 void axisExtractor02::modelo(vtkImageData *data, unsigned char cantidad, unsigned long vector[50][4], int candit[10][3], double radioactual, double minis[10])
2594 int i, j, k, radio, radio2;
2605 double dist, prop, rest;
2608 data->GetExtent(ext);
2612 radio=(ext[1]-ext[0])/2;
2613 radio2=(int)radioactual;
2615 centro[0]=(ext[1]+ext[0])/2;
2616 centro[1]=(ext[3]+ext[2])/2;
2617 centro[2]=(ext[5]+ext[4])/2;
2619 indextoreal(centro, centrof );
2623 for(i=ext[0];i<=ext[1];i++){
2624 for(j=ext[2];j<=ext[3];j++){
2625 for(k=ext[4];k<=ext[5];k++){
2626 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2631 if( radioactual<distancia(punto, centro) ){
2642 for(t=0;t<cantidad;t++){
2643 radiof=sqrt(minis[t]);
2645 punto2[0]=candit[t][0];
2646 punto2[1]=candit[t][1];
2647 punto2[2]=candit[t][2];
2649 indextoreal(punto2, puntof2 );
2652 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2653 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2654 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2655 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2659 indextoreal(punto, puntof );
2661 dist=distanciaejepunto(puntof, centrof, puntof2);
2662 prop=proporcioejepunto(puntof, centrof, puntof2);
2664 if(prop>=0 && prop<=1){
2665 rest=(radiof-radioactual)*prop+radioactual;
2667 rest=rest*espprin[0];
2672 tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2);
2673 if( radio>sqrt(tmpsqrt) && rest>dist){
2685 //----------------------------------------------------------------------------
2686 void axisExtractor02::comparacion(vtkImageData *data, vtkImageData *data2, unsigned long minis[4])
2692 data->GetExtent(ext);
2695 unsigned char *ptr2;
2706 radio=((double)ext[1]-(double)ext[0])/2;
2708 centro[0]=(ext[1]+ext[0])/2;
2709 centro[1]=(ext[3]+ext[2])/2;
2710 centro[2]=(ext[5]+ext[4])/2;
2712 for(i=ext[0];i<=ext[1];i++){
2713 for(j=ext[2];j<=ext[3];j++){
2714 for(k=ext[4];k<=ext[5];k++){
2715 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2716 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
2722 if( radio>distancia(punto, centro) ){
2724 if( *ptr==0 && *ptr2==0 ){
2727 else if( *ptr!=0 && *ptr2!=0 ){
2730 else if(*ptr==0 && *ptr2!=0){
2744 //----------------------------------------------------------------------------
2745 void axisExtractor02::copiar(vtkImageData *data, vtkImageData *data2 )
2750 data->GetExtent(ext);
2752 unsigned char *ptr, *ptr2;
2754 for(i=ext[0];i<=ext[1];i++){
2755 for(j=ext[2];j<=ext[3];j++){
2756 for(k=ext[4];k<=ext[5];k++){
2757 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
2758 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
2759 if(*ptr!=0 && *ptr2==0){
2772 //----------------------------------------------------------------------------
2773 void axisExtractor02::copiar2(vtkImageData *data, vtkImageData *data2 )
2778 data->GetExtent(ext);
2782 for(i=ext[0];i<=ext[1];i++){
2783 for(j=ext[2];j<=ext[3];j++){
2784 for(k=ext[4];k<=ext[5];k++){
2785 ptr=(double *) data->GetScalarPointer(i,j,k);
2786 ptr2=(double *) data2->GetScalarPointer(i,j,k);
2787 if(*ptr!=0 && *ptr2==0){
2801 //----------------------------------------------------------------------------
2802 double axisExtractor02::angulo(double a[3], double b[3] )
2804 double m1=sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2]));
2805 double m2=sqrt((b[0]*b[0])+(b[1]*b[1])+(b[2]*b[2]));
2806 double res=(acos(((a[0]*b[0])+(a[1]*b[1])+(a[2]*b[2]))/(m1*m2))*180)/3.1415;
2812 if(res>=0 && res<90){
2815 else if(res>=90 && res<180){
2818 else if(res>=180 && res<270){
2831 //----------------------------------------------------------------------------
2832 double axisExtractor02::angulo(double i1, double j1, double k1, double i2, double j2, double k2 )
2834 double m1=sqrt((i1*i1)+(j1*j1)+(k1*k1));
2835 double m2=sqrt((i2*i2)+(j2*j2)+(k2*k2));
2836 double res=(acos(((i1*i2)+(j1*j2)+(k1*k2))/(m1*m2))*180)/3.1415;
2842 if(res>=0 && res<90){
2845 else if(res>=90 && res<180){
2848 else if(res>=180 && res<270){
2862 //----------------------------------------------------------------------------
2863 int axisExtractor02::envolumen(int a[3], vtkImageData *datae )
2868 unsigned short *ptr;
2870 ptr=(unsigned short *) datae->GetScalarPointer(a);
2885 //----------------------------------------------------------------------------
2886 int axisExtractor02::mincandit(int candit[10][3], int cantidad, double puntoanterior[3])
2892 double canditreal[3];
2899 for(i=0;i<cantidad && i<10;i++){
2900 indextoreal(candit[i], canditreal);
2901 distentpu=distancia(canditreal, puntoanterior);
2902 if(distentpu<distmientrp){
2903 distmientrp=distentpu;
2914 //----------------------------------------------------------------------------
2915 int axisExtractor02::maxareacandit(unsigned long vector[50][4], int cantidad)
2925 for(i=0;i<cantidad && i<10;i++){
2926 if(max<vector[i][0]){
2939 //----------------------------------------------------------------------------
2940 unsigned long axisExtractor02::totalarea(unsigned long vector[50][4], unsigned long vectorb[50][4], int cantidad, int cantidadb)
2948 for(i=0;i<cantidad && i<10;i++){
2952 for(i=0;i<cantidadb;i++){
2953 area+=vectorb[i][0];
2963 //----------------------------------------------------------------------------
2964 unsigned long axisExtractor02::conecarea(unsigned long vector[50][4], int cantidad)
2972 for(i=0;i<cantidad && i<10;i++){
2984 //----------------------------------------------------------------------------
2985 int axisExtractor02::bruled(int candit[10][3], int cantidad, vtkImageData *data4)
2991 for(i=0;i<cantidad && i<10;i++){
2992 if(envolumen(candit[i], data4)!=0){
3003 //----------------------------------------------------------------------------
3004 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)
3027 data->GetExtent(ext);
3029 radio=(ext[1]-ext[0])/2;
3030 radiof=((double)ext[1]-(double)ext[0])/2;
3031 centro[0]=(ext[1]+ext[0])/2;
3032 centro[1]=(ext[3]+ext[2])/2;
3033 centro[2]=(ext[5]+ext[4])/2;
3035 indicecorregido[0]=centro[0];
3036 indicecorregido[1]=centro[1];
3037 indicecorregido[2]=centro[2];
3038 indextoreal(indicecorregido, puntocorregido);
3039 indextoreal(indiceanterior, centrof);
3040 indextoreal(indicepre, puntof2);
3042 rap= (int)( (double)radio/2 );
3046 for(i=ext[0];i<=ext[1];i++){
3047 for(j=ext[2];j<=ext[3];j++){
3048 for(k=ext[4];k<=ext[5];k++){
3049 ptr=(double *) data->GetScalarPointer(i,j,k);
3053 indextoreal(punto, puntof );
3055 distcorr3=distancia(indiceanterior, indicepre);
3056 distcorr2=distancia(punto, indicepre);
3057 distcorr=distancia(punto, indiceanterior);
3059 dist=distanciaejepunto(puntof, centrof, puntof2);
3060 prop=proporcioejepunto(puntof, centrof, puntof2);
3062 if(prop>=0 && prop<=1){
3063 rest=(radiopre-radioanterior)*prop+radioanterior;
3064 rest=rest*espprin[0];
3070 if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && rest>dist){
3071 //if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && rest>dist){
3072 //if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && distcorr2<=distcorr3-radioanterior){
3076 indicecorregido[0]=i;
3077 indicecorregido[1]=j;
3078 indicecorregido[2]=k;
3079 indextoreal(indicecorregido, puntocorregido);
3081 else if( maxcent==*ptr){
3082 if( mincent>distcorr){
3085 indicecorregido[0]=i;
3086 indicecorregido[1]=j;
3087 indicecorregido[2]=k;
3088 indextoreal(indicecorregido, puntocorregido);
3096 ptr=(double *) data->GetScalarPointer(centro[0],centro[1],centro[2]);
3104 //----------------------------------------------------------------------------
3105 double axisExtractor02::correction2(int candit[10][3], int cantidad, vtkImageData *data, int indicecorregido[3], double puntocorregido[3], int indiceanterior[3], double radioanterior)
3120 data->GetExtent(ext);
3122 radio=(ext[1]-ext[0])/2;
3123 centro[0]=(ext[1]+ext[0])/2;
3124 centro[1]=(ext[3]+ext[2])/2;
3125 centro[2]=(ext[5]+ext[4])/2;
3127 indicecorregido[0]=centro[0];
3128 indicecorregido[1]=centro[1];
3129 indicecorregido[2]=centro[2];
3130 indextoreal(indicecorregido, puntocorregido);
3132 rap= (int)((double)radio/2);
3136 for(i=centro[0]-rap;i<=centro[0]+rap;i++){
3137 for(j=centro[1]-rap;j<=centro[1]+rap;j++){
3138 for(k=centro[2]-rap;k<=centro[2]+rap;k++){
3139 ptr=(double *) data->GetScalarPointer(i,j,k);
3144 distcorr2=distancia(punto, centro);
3145 distcorr=distancia(punto, indiceanterior);
3147 if( rap>distcorr2 ){
3151 indicecorregido[0]=i;
3152 indicecorregido[1]=j;
3153 indicecorregido[2]=k;
3154 indextoreal(indicecorregido, puntocorregido);
3156 else if( maxcent==*ptr){
3157 if( mincent>distcorr){
3160 indicecorregido[0]=i;
3161 indicecorregido[1]=j;
3162 indicecorregido[2]=k;
3163 indextoreal(indicecorregido, puntocorregido);
3173 return sqrt(maxcent);
3180 //----------------------------------------------------------------------------
3181 void axisExtractor02::avanzar()
3184 double puntoactual[3];
3185 double puntocorregido[3];
3186 double puntoanterior[3];
3187 double puntosiguiente[3];
3188 //double puntointer[3];
3191 int indiceactual[3];
3192 int indiceanterior[3];
3193 int indicecorregido[3];
3194 //int indiceinter[3];
3199 unsigned char cantidad;
3200 unsigned char cantidadb;
3202 unsigned long vector[50][4];
3203 unsigned long vectorb[50][4];
3243 double radiotrabajof;
3248 double radioanterior;
3254 double canditreal[3];
3256 unsigned long totalsurf;
3257 unsigned long conecsurf;
3260 unsigned long caf[4];
3269 if(!m_Stack0.empty()){
3272 /* indexp=m_Stack.top();
3274 puntoactual[0]=m_Stack0.top();
3275 puntoactual[1]=m_Stack1.top();
3276 puntoactual[2]=m_Stack2.top();
3278 puntoanterior[0]=m_Stack3.top();
3279 puntoanterior[1]=m_Stack4.top();
3280 puntoanterior[2]=m_Stack5.top();
3282 puntopre[0]=m_Stack6.top();
3283 puntopre[1]=m_Stack7.top();
3284 puntopre[2]=m_Stack8.top();
3288 radiopred=m_Stackr.top();
3289 radioanterior=m_Stackra.top();
3290 radioprinc=m_Stackrp.top();*/
3292 indexp=m_Stack.front();
3294 puntoactual[0]=m_Stack0.front();
3295 puntoactual[1]=m_Stack1.front();
3296 puntoactual[2]=m_Stack2.front();
3298 puntoanterior[0]=m_Stack3.front();
3299 puntoanterior[1]=m_Stack4.front();
3300 puntoanterior[2]=m_Stack5.front();
3302 puntopre[0]=m_Stack6.front();
3303 puntopre[1]=m_Stack7.front();
3304 puntopre[2]=m_Stack8.front();
3306 radiopred=m_Stackr.front();
3307 radioanterior=m_Stackra.front();
3308 radioprinc=m_Stackrp.front();
3310 radiotrabajo= (int) (proprad*radiopred);
3311 radiotrabajof=proprad*radiopred;
3313 if(radiotrabajof-radiotrabajo>0.5){
3314 radiotrabajo=radiotrabajo+1;
3317 // EED 15 Mai 2007 .NET
3318 // radiotrabajo=max(minant,radiotrabajo);
3319 if (minant>radiotrabajo)
3321 radiotrabajo=minant;
3339 m_Stack0.pop_front();
3340 m_Stack1.pop_front();
3341 m_Stack2.pop_front();
3342 m_Stack3.pop_front();
3343 m_Stack4.pop_front();
3344 m_Stack5.pop_front();
3345 m_Stack6.pop_front();
3346 m_Stack7.pop_front();
3347 m_Stack8.pop_front();
3349 m_Stack.pop_front();
3350 m_Stackr.pop_front();
3351 m_Stackra.pop_front();
3352 m_Stackrp.pop_front();
3354 dirant[0]=puntoanterior[0]-puntoactual[0];
3355 dirant[1]=puntoanterior[1]-puntoactual[1];
3356 dirant[2]=puntoanterior[2]-puntoactual[2];
3359 realtoindex(puntoactual, indiceactual);
3360 realtoindex(puntoanterior, indiceanterior);
3361 realtoindex(puntopre, indicepre);
3363 radioactual=radiopred;
3365 if(radiotrabajo<=maxant){
3367 extint[0]=indiceactual[0]-radiotrabajo;
3368 extint[1]=indiceactual[0]+radiotrabajo;
3369 extint[2]=indiceactual[1]-radiotrabajo;
3370 extint[3]=indiceactual[1]+radiotrabajo;
3371 extint[4]=indiceactual[2]-radiotrabajo;
3372 extint[5]=indiceactual[2]+radiotrabajo;
3374 extrac->SetVOI(extint);
3375 extrac->UpdateWholeExtent();
3377 extrac->GetOutput()->GetExtent(extintreal);
3379 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]){
3380 fprintf(stream, "salio\t");
3386 data1=vtkImageData::New();
3387 data1->SetScalarTypeToUnsignedChar();
3388 data1->SetExtent(extrac->GetOutput()->GetExtent());
3389 data1->SetSpacing(extrac->GetOutput()->GetSpacing());
3391 cilindro(data1, dirant );
3393 optim(extrac->GetOutput(), data1 );
3395 connect->SetInput(data1);
3396 connect->RemoveAllSeeds();
3397 connect->AddSeed(indiceactual[0],indiceactual[1],indiceactual[2]);
3398 connect->UpdateWholeExtent();
3402 data2=vtkImageData::New();
3403 data2->SetScalarTypeToUnsignedChar();
3404 data2->SetExtent(connect->GetOutput()->GetExtent());
3405 data2->SetSpacing(connect->GetOutput()->GetSpacing());
3409 cantidad=find_components(connect->GetOutput(), data2, 0, vector);
3412 data3=vtkImageData::New();
3413 data3->SetScalarTypeToUnsignedChar();
3414 data3->SetExtent(connect->GetOutput()->GetExtent());
3415 data3->SetSpacing(connect->GetOutput()->GetSpacing());
3419 cantidadb=find_componentsb(connect->GetOutput(), data3, 0, vectorb);
3421 redondear3(connect->GetOutput());
3425 redondear2(distance->GetOutput());
3427 redondear(connect->GetOutput());
3429 costominimo(distance->GetOutput(), data2);
3433 radioactual=correction2(candit, cantidad, distance->GetOutput(), indicecorregido, puntocorregido, indiceanterior, radioanterior);
3436 radioactual=correction(candit, cantidad, distance->GetOutput(), indicecorregido, puntocorregido, indiceanterior, radioanterior, indicepre, radioprinc);
3440 //inpeq(extrac->GetOutput(), indicecorregido, radioactual);
3442 indant=mincandit(candit, cantidad, puntoanterior);
3444 indmaxarea=maxareacandit(vector, cantidad);
3445 totalsurf=totalarea(vector, vectorb, cantidad, cantidadb);
3446 conecsurf=conecarea(vector, cantidad);
3448 indextoreal(candit[indant], canditreal);
3450 burned=bruled(candit, cantidad, data4);
3453 data6=vtkImageData::New();
3454 data6->SetScalarTypeToUnsignedChar();
3455 data6->SetExtent(extrac->GetOutput()->GetExtent());
3456 data6->SetSpacing(extrac->GetOutput()->GetSpacing());
3458 modelo(data6, cantidad, vector, candit, radioactual, minis);
3459 comparacion(connect->GetOutput(), data6, caf);
3461 rel1=((double)caf[0]+(double)caf[1])/((double)caf[0]+(double)caf[1]+(double)caf[2]+(double)caf[3]);
3462 rel2=((double)caf[1])/((double)caf[1]+(double)caf[3]);
3463 rel3=((double)caf[1])/((double)caf[1]+(double)caf[2]);
3464 rel4=(double)conecsurf/(double)totalsurf;
3465 rel5=(double)vector[indmaxarea][0]/(double)totalsurf;
3466 rel6=(double)cant/((double)cant+(double)cant2);
3469 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)){
3473 if(indicecorregido[0]!=indiceactual[0] || indicecorregido[1]!=indiceactual[1] || indicecorregido[2]!=indiceactual[2]){
3477 if(burned>=cantidad ){
3489 if(flag7==1 && flagg>4){
3501 if(flag7==1 && flag18==1){
3502 for(i=0;i<flagg;i++){
3503 if(visited[i][0]==indiceactual[0] && visited[i][1]==indiceactual[1] && visited[i][2]==indiceactual[2] && visitedrad[i]==radiopred){
3509 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) )){
3517 if(angulo(dirant[0], dirant[1], dirant[2], canditreal[0]-puntoactual[0], canditreal[1]-puntoactual[1], canditreal[2]-puntoactual[2] )>65){
3521 if(angulo(Vp[0][ejeminp], Vp[1][ejeminp], Vp[2][ejeminp], Vi[0][ejemini], Vi[1][ejemini], Vi[2][ejemini] )>60){
3525 if((flag3==1 || flag5==1 || (flag16==1 && flag11==1 && flag19==0) ) && (flag7==0 && flag8==0)){
3527 //fprintf(stream, "malo\t");
3530 if(flag3==1 || flag5==1 || (flag16==1 && flag11==1 && flag19==0) ){
3534 if((mejordst<radioactual || (mejordst==radioactual && mejorrad>radiopred) || (mejordst==radioactual && mejorrad<=radiopred && mejorcant<cantidad )) && flag17==0 && flag20==0 && flag4==0 && flag12==0 && flag15==0 ){
3535 mejordst=radioactual;
3536 mejor[0]=puntoactual[0];
3537 mejor[1]=puntoactual[1];
3538 mejor[2]=puntoactual[2];
3545 fprintf(stream, "grande\t");
3561 if( flag9==1 || flag4==1 || flag12==1 || ((cantidad<2 || burned==cantidad) && (flag7==0 && flag8==0)) ){
3568 if( cantidad>2 && (flag7==0 && flag8==0) ){
3575 if((flag7==0 && flag8==0 && (mejordst>radioactual || (mejordst==radioactual && mejorrad<radiopred)) && flag19==0) || flag2==1){
3581 fprintf(stream, "mejor\t");
3583 m_Stack0.push_front(mejor[0]);
3584 m_Stack1.push_front(mejor[1]);
3585 m_Stack2.push_front(mejor[2]);
3586 m_Stack3.push_front(puntoanterior[0]);
3587 m_Stack4.push_front(puntoanterior[1]);
3588 m_Stack5.push_front(puntoanterior[2]);
3589 m_Stack6.push_front(puntopre[0]);
3590 m_Stack7.push_front(puntopre[1]);
3591 m_Stack8.push_front(puntopre[2]);
3592 m_Stack.push_front(indexp);
3593 m_Stackr.push_front(mejorrad);
3594 m_Stackra.push_front(radioanterior);
3595 m_Stackrp.push_front(radioprinc);
3602 else if((flag17==1 || flag9==1 || flag4==1 || flag12==1 || flag15==1) && flag18==1){
3607 fprintf(stream, "descorrigio\t");
3609 m_Stack0.push_front(mejor[0]);
3610 m_Stack1.push_front(mejor[1]);
3611 m_Stack2.push_front(mejor[2]);
3612 m_Stack3.push_front(puntoanterior[0]);
3613 m_Stack4.push_front(puntoanterior[1]);
3614 m_Stack5.push_front(puntoanterior[2]);
3615 m_Stack6.push_front(puntopre[0]);
3616 m_Stack7.push_front(puntopre[1]);
3617 m_Stack8.push_front(puntopre[2]);
3618 m_Stack.push_front(indexp);
3619 m_Stackr.push_front(mejorrad);
3620 m_Stackra.push_front(radioanterior);
3621 m_Stackrp.push_front(radioprinc);
3628 else if(flag12==0 && flag4==0 && flag9==0){
3629 if(flag7==0 && flag8==0){
3631 if(flag17==0 && flag15==0){
3638 fprintf(stream, "inserto\t");
3641 /*costominimo2(distance->GetOutput(), connect->GetOutput(), indiceactual, candit[indant], indiceinter);
3642 indextoreal(indiceinter, puntointer);
3643 realtoreal2(puntointer, provvc);
3644 points->InsertPoint(buenos,provvc);
3646 lineas->InsertNextCell(2);
3647 lineas->InsertCellPoint(indexp);
3648 lineas->InsertCellPoint(buenos);
3652 realtoreal2(puntoactual, provvc);
3653 points->InsertPoint(buenos,provvc);
3655 lineas->InsertNextCell(2);
3656 lineas->InsertCellPoint(buenos-1);
3657 lineas->InsertCellPoint(buenos);
3661 realtoreal2(puntoactual, provvc);
3662 points->InsertPoint(buenos,provvc);
3664 lineas->InsertNextCell(2);
3665 lineas->InsertCellPoint(indexp);
3666 lineas->InsertCellPoint(buenos);
3672 realtoreal2(puntoactual, provvc);
3673 points->InsertPoint(buenos,provvc);
3686 fprintf(stream, "para\t");
3689 else if(flag7==0 && flag8==1){
3692 fprintf(stream, "agrando\t");
3694 m_Stack0.push_front(puntoactual[0]);
3695 m_Stack1.push_front(puntoactual[1]);
3696 m_Stack2.push_front(puntoactual[2]);
3697 m_Stack3.push_front(puntoanterior[0]);
3698 m_Stack4.push_front(puntoanterior[1]);
3699 m_Stack5.push_front(puntoanterior[2]);
3700 m_Stack6.push_front(puntopre[0]);
3701 m_Stack7.push_front(puntopre[1]);
3702 m_Stack8.push_front(puntopre[2]);
3704 m_Stack.push_front(indexp);
3705 m_Stackr.push_front(((double)radiotrabajo+1)/2);
3706 m_Stackra.push_front(radioanterior);
3707 m_Stackrp.push_front(radioprinc);
3710 else if(flag7==1 && flag8==1){
3711 visited[flagg][0]=indiceactual[0];
3712 visited[flagg][1]=indiceactual[1];
3713 visited[flagg][2]=indiceactual[2];
3714 visitedrad[flagg]=radiopred;
3718 fprintf(stream, "agrando y corrigiendo\t");
3720 m_Stack0.push_front(puntocorregido[0]);
3721 m_Stack1.push_front(puntocorregido[1]);
3722 m_Stack2.push_front(puntocorregido[2]);
3723 m_Stack3.push_front(puntoanterior[0]);
3724 m_Stack4.push_front(puntoanterior[1]);
3725 m_Stack5.push_front(puntoanterior[2]);
3726 m_Stack6.push_front(puntopre[0]);
3727 m_Stack7.push_front(puntopre[1]);
3728 m_Stack8.push_front(puntopre[2]);
3730 m_Stack.push_front(indexp);
3731 m_Stackr.push_front(((double)radiotrabajo+1)/2);
3732 m_Stackra.push_front(radioanterior);
3733 m_Stackrp.push_front(radioprinc);
3736 visited[flagg][0]=indiceactual[0];
3737 visited[flagg][1]=indiceactual[1];
3738 visited[flagg][2]=indiceactual[2];
3739 visitedrad[flagg]=radiopred;
3746 fprintf(stream, "corrigio\t");
3748 m_Stack0.push_front(puntocorregido[0]);
3749 m_Stack1.push_front(puntocorregido[1]);
3750 m_Stack2.push_front(puntocorregido[2]);
3751 m_Stack3.push_front(puntoanterior[0]);
3752 m_Stack4.push_front(puntoanterior[1]);
3753 m_Stack5.push_front(puntoanterior[2]);
3755 m_Stack6.push_front(puntopre[0]);
3756 m_Stack7.push_front(puntopre[1]);
3757 m_Stack8.push_front(puntopre[2]);
3759 m_Stack.push_front(indexp);
3760 m_Stackr.push_front(radiopred);
3761 m_Stackra.push_front(radioanterior);
3762 m_Stackrp.push_front(radioprinc);
3767 if(flag7==0 && flag8==0){
3769 if(flag17==0 && flag15==0){
3771 for(i=0;i<cantidad && i<10;i++){
3772 for(j=i;j<cantidad && i<10;j++){
3773 if(minis[j]>minis[i]){
3774 vectortemp[0]=candit[i][0];
3775 vectortemp[1]=candit[i][1];
3776 vectortemp[2]=candit[i][2];
3779 candit[i][0]=candit[j][0];
3780 candit[i][1]=candit[j][1];
3781 candit[i][2]=candit[j][2];
3784 candit[j][0]=vectortemp[0];
3785 candit[j][1]=vectortemp[1];
3786 candit[j][2]=vectortemp[2];
3794 for(i=0;i<cantidad && i<10;i++){
3795 if(i!=indant || buenos<2){
3796 indextoreal(candit[i], puntosiguiente);
3798 m_Stack0.push_front(puntosiguiente[0]);
3799 m_Stack1.push_front(puntosiguiente[1]);
3800 m_Stack2.push_front(puntosiguiente[2]);
3801 m_Stack3.push_front(puntoactual[0]);
3802 m_Stack4.push_front(puntoactual[1]);
3803 m_Stack5.push_front(puntoactual[2]);
3805 m_Stack6.push_front(puntosiguiente[0]);
3806 m_Stack7.push_front(puntosiguiente[1]);
3807 m_Stack8.push_front(puntosiguiente[2]);
3809 double prov1a=sqrt(minis[i]);
3811 m_Stack.push_front(buenos-1);
3812 m_Stackr.push_front(prov1a);
3813 m_Stackra.push_front(radioactual);
3814 m_Stackrp.push_front(prov1a);
3820 for(i=0;i<cantidad && i<10;i++){
3821 indextoreal(candit[i], puntosiguiente);
3823 if(envolumen(candit[i], data4)==0){
3824 m_Stack0.push_front(puntosiguiente[0]);
3825 m_Stack1.push_front(puntosiguiente[1]);
3826 m_Stack2.push_front(puntosiguiente[2]);
3827 m_Stack3.push_front(puntoactual[0]);
3828 m_Stack4.push_front(puntoactual[1]);
3829 m_Stack5.push_front(puntoactual[2]);
3831 m_Stack6.push_front(puntosiguiente[0]);
3832 m_Stack7.push_front(puntosiguiente[1]);
3833 m_Stack8.push_front(puntosiguiente[2]);
3835 double prov1b=sqrt(minis[i]);
3837 m_Stack.push_front(buenos-1);
3838 m_Stackr.push_front(prov1b);
3839 m_Stackra.push_front(radioactual);
3840 m_Stackrp.push_front(prov1b);
3845 copiar(connect->GetOutput(), data4 );
3846 //copiar2(distance->GetOutput(), data5 );
3851 // fprintf(stream, "algo");
3864 printf("wi[0]: %f\n", wi[0]);
3865 printf("wi[1]: %f\n", wi[1]);
3866 printf("wi[2]: %f\n", wi[2]);
3868 printf("inerciari: %f\n", inerciari);
3869 printf("inerciariy: %f\n", inerciariy);
3870 printf("inerciariz: %f\n", inerciariz);
3872 printf("ri1: %f\n", inerciari/inerciariy);
3873 printf("ri2: %f\n", inerciariy/inerciariz);
3874 printf("ri3: %f\n", inerciari/inerciariz);
3876 printf("inerciarp: %f\n", inerciarp);
3877 printf("inerciarpy: %f\n", inerciarpy);
3878 printf("inerciarpz: %f\n", inerciarpz);
3880 printf("rp1: %f\n", inerciarp/inerciarpy);
3881 printf("rp2: %f\n", inerciarpy/inerciarpz);
3882 printf("rp3: %f\n", inerciarp/inerciarpz);
3884 printf("angulo: %f\n",angulo(dirant[0], dirant[1], dirant[2], canditreal[0]-puntoactual[0], canditreal[1]-puntoactual[1], canditreal[2]-puntoactual[2] ));
3885 printf("angulo2: %f\n", angulo(Vp[0][ejeminp], Vp[1][ejeminp], Vp[2][ejeminp], Vi[0][ejemini], Vi[1][ejemini], Vi[2][ejemini]));
3887 printf("centi: %f\n", centi);
3888 printf("centi2: %f\n", centi2);
3893 fprintf(stream, "%d\t", indexp);
3894 fprintf(stream, "%d\t", buenos);
3895 fprintf(stream, "%d\t", radiotrabajo);
3896 fprintf(stream, "%f\t", radioactual);
3897 fprintf(stream, "%f\t", radiopred);
3898 fprintf(stream, "%f\t", radioanterior);
3899 fprintf(stream, "%f\t", mejordst);
3900 fprintf(stream, "%f\t", mejorrad);
3901 fprintf(stream, "%d\t", mejorcant);
3905 fprintf(stream, "%d\t", cant);
3906 fprintf(stream, "%d\t", cant2);
3907 fprintf(stream, "%f\t", (double)cant/((double)cant+(double)cant2));
3910 fprintf(stream, "%d\t", flag2);
3911 fprintf(stream, "%d\t", flag3);
3912 fprintf(stream, "%d\t", flag4);
3913 fprintf(stream, "%d\t", flag5);
3914 fprintf(stream, "%d\t", flag6);
3915 fprintf(stream, "%d\t", flag7);
3916 fprintf(stream, "%d\t", flag8);
3917 fprintf(stream, "%d\t", flag9);
3918 fprintf(stream, "%d\t", flag10);
3919 fprintf(stream, "%d\t", flag11);
3920 fprintf(stream, "%d\t", flag12);
3921 fprintf(stream, "%d\t", flag13);
3922 fprintf(stream, "%d\t", flag14);
3923 fprintf(stream, "%d\t", flag15);
3924 fprintf(stream, "%d\t", flag16);
3925 fprintf(stream, "%d\t", flag17);
3926 fprintf(stream, "%d\t", flag18);
3927 fprintf(stream, "%d\t", flag19);
3928 fprintf(stream, "%d\t", flag20);
3929 fprintf(stream, "%d\t", flagg);
3930 fprintf(stream, "%d\t", flagg2);
3931 fprintf(stream, "%d\t", cantidad);
3932 fprintf(stream, "%d\t", cantidadb);
3933 fprintf(stream, "%d\t", burned);
3937 fprintf(stream, "caf[0] %d\t", caf[0]);
3938 fprintf(stream, "caf[1] %d\t", caf[1]);
3939 fprintf(stream, "caf[2] %d\t", caf[2]);
3940 fprintf(stream, "caf[3] %d\t", caf[3]);
3941 fprintf(stream, "REL%f\t", ((double)caf[0]+(double)caf[1])/((double)caf[0]+(double)caf[1]+(double)caf[2]+(double)caf[3]));
3942 fprintf(stream, "REL2%f\t", ((double)caf[1])/((double)caf[3]+(double)caf[1]));
3943 fprintf(stream, "REL3%f\t", ((double)caf[1])/((double)caf[1]+(double)caf[2]));
3945 fprintf(stream, "maxareacandit%d\t", vector[indmaxarea][0]);
3946 fprintf(stream, "totalsurf%d\t", totalsurf);
3947 fprintf(stream, "conecsurf%d\t", conecsurf);
3948 fprintf(stream, "ratio%f\t", (double)conecsurf/(double)totalsurf);
3949 fprintf(stream, "ratio2%f\t", (double)vector[indmaxarea][0]/(double)totalsurf);
3953 fprintf(stream, "[%f, %f, %f]\t", puntoactual[0], puntoactual[1], puntoactual[2]);
3954 fprintf(stream, "[%f, %f, %f]\t", puntoanterior[0], puntoanterior[1], puntoanterior[2]);
3955 fprintf(stream, "[%f, %f, %f]\t", puntopre[0], puntopre[1], puntopre[2]);
3956 fprintf(stream, "[%f, %f, %f]\t", mejor[0], mejor[1], mejor[2]);
3958 fprintf(stream, "\n");
3967 //----------------------------------------------------------------------------
3968 void axisExtractor02::todo()
3972 while( !m_Stack0.empty() && i<5000){
3980 //----------------------------------------------------------------------------
3981 void axisExtractor02::paso()
3987 while( !m_Stack0.empty() && inicio==buenos && i<5000){
3996 //----------------------------------------------------------------------------
3997 void axisExtractor02::rama()
4003 while( !m_Stack0.empty() && frama==0 && i<5000){
4012 //----------------------------------------------------------------------------
4013 void axisExtractor02::segmento()
4020 while( !m_Stack0.empty() && frama==0 && fseg==0 && i<5000){