00001
00002
00003
00004
00005 #include "axisExtractor02.h"
00006 #include "vtkPolyData.h"
00007 #include "vtkObjectFactory.h"
00008 #include "vtkImageThreshold.h"
00009 #include "vtkImageCast.h"
00010 #include "vtkImageSeedConnectivity.h"
00011 #include "vtkImageData.h"
00012 #include "vtkMarchingCubes.h"
00013 #include "vtkDoubleArray.h"
00014 #include "vtkPointData.h"
00015 #include "vtkExtractVOI.h"
00016 #include "vtkImageConstantPad.h"
00017 #include "vtkImageTranslateExtent.h"
00018 #include "vtkMath.h"
00019 #include "vtkTransformPolyDataFilter.h"
00020 #include "vtkTransform.h"
00021 #include <time.h>
00022
00023
00024
00025
00026
00027 vtkStandardNewMacro(axisExtractor02);
00028
00029
00030
00031 void axisExtractor02::SetParam(double value)
00032 {
00033 this->param=value;
00034 }
00035
00036
00037
00038
00039 double axisExtractor02::GetParam()
00040 {
00041 return this->param;
00042 }
00043
00044
00045
00046
00047
00048 void axisExtractor02::SetParam2(double value)
00049 {
00050 this->param2=value;
00051 }
00052
00053
00054
00055
00056 double axisExtractor02::GetParam2()
00057 {
00058 return this->param2;
00059 }
00060
00061
00062
00063
00064
00065 void axisExtractor02::SetParam3(double value)
00066 {
00067 this->param3=value;
00068 }
00069
00070
00071
00072
00073 double axisExtractor02::GetParam3()
00074 {
00075 return this->param3;
00076 }
00077
00078
00079
00080
00081 void axisExtractor02::SetMaxant(int value)
00082 {
00083 this->maxant=value;
00084 }
00085
00086
00087
00088
00089 int axisExtractor02::GetMaxant()
00090 {
00091 return this->maxant;
00092 }
00093
00094
00095
00096
00097 void axisExtractor02::SetMinant(int value)
00098 {
00099 this->minant=value;
00100 }
00101
00102
00103
00104
00105 int axisExtractor02::GetMinant()
00106 {
00107 return this->minant;
00108 }
00109
00110
00111
00112 void axisExtractor02::SetPoint(double value[3])
00113 {
00114
00115
00116 double spa[3];
00117
00118 this->GetInput()->GetSpacing (spa);
00119
00120 value[0]=value[0]+maxant*spa[0];
00121 value[1]=value[1]+maxant*spa[1];
00122 value[2]=value[2]+maxant*spa[2];
00123
00124 this->m_Stack0.push_front(value[0]);
00125 this->m_Stack1.push_front(value[1]);
00126 this->m_Stack2.push_front(value[2]);
00127 this->m_Stack3.push_front(value[0]);
00128 this->m_Stack4.push_front(value[1]);
00129 this->m_Stack5.push_front(value[2]);
00130 this->m_Stack6.push_front(value[0]);
00131 this->m_Stack7.push_front(value[1]);
00132 this->m_Stack8.push_front(value[2]);
00133 this->m_Stack.push_front(0);
00134 this->m_Stackr.push_front(0);
00135 this->m_Stackra.push_front(0);
00136 this->m_Stackrp.push_front(0);
00137 }
00138
00139
00140
00141
00142 vtkImageData *axisExtractor02::GetVolumen()
00143 {
00144
00145 vtkImageTranslateExtent *trans;
00146
00147 trans = vtkImageTranslateExtent::New();
00148
00149 trans->SetInput(this->data4);
00150
00151 trans->SetTranslation (maxant, maxant, maxant);
00152
00153 trans->Update();
00154
00155 return trans->GetOutput();
00156 }
00157
00158
00159
00160
00161
00162
00163 vtkPolyData *axisExtractor02::GetOutput ()
00164 {
00165
00166
00167
00168
00169 vtkTransform *transL1;
00170 vtkTransformPolyDataFilter *trans;
00171
00172 transL1 = vtkTransform::New();
00173 trans = vtkTransformPolyDataFilter::New();
00174
00175
00176 double spa[3];
00177 this->GetInput()->GetSpacing (spa);
00178 transL1->Translate(-maxant*spa[0], -maxant*spa[1], -maxant*spa[2]);
00179
00180
00181
00182 trans->SetInput((vtkPolyData *)(this->Outputs[0]));
00183
00184 trans->SetTransform (transL1);
00185
00186 trans->Update();
00187
00188 return trans->GetOutput();
00189 }
00190
00191
00192
00193
00194
00195
00196
00197
00198 axisExtractor02::axisExtractor02() {
00199
00200
00201 this->NumberOfRequiredInputs = 1;
00202 this->param=1;
00203 this->param2=1;
00204 this->param3=0.5;
00205 this->param4=0;
00206
00207
00208
00209
00210 this->data4=vtkImageData::New();
00211
00212 this->data1=vtkImageData::New();
00213 this->data2=vtkImageData::New();
00214 this->data3=vtkImageData::New();
00215 this->data6=vtkImageData::New();
00216
00217 this->extrac = vtkExtractVOI::New();
00218
00219 this->extrac->SetSampleRate(1, 1, 1);
00220
00221 this->connect = vtkImageSeedConnectivity::New();
00222 this->connect->SetInput(this->data1);
00223 this->connect->SetInputConnectValue(255);
00224 this->connect->SetOutputConnectedValue(255);
00225 this->connect->SetOutputUnconnectedValue(0);
00226
00227 this->distance= vtkImageEuclideanDistance::New();
00228 this->distance->SetInput(this->connect->GetOutput());
00229 this->distance->InitializeOn();
00230 this->distance->ConsiderAnisotropyOff();
00231
00232
00233 }
00234
00235
00236
00237
00238 void axisExtractor02::SetInput(vtkImageData *input)
00239 {
00240
00241 vtkImageConstantPad *pad;
00242 vtkImageTranslateExtent *trans;
00243
00244 pad = vtkImageConstantPad::New();
00245 trans = vtkImageTranslateExtent::New();
00246
00247 pad->SetInput(input);
00248 trans->SetInput(pad->GetOutput());
00249
00250
00251 pad->SetConstant(0);
00252
00253 pad->SetInput(input);
00254
00255 int ext[6];
00256
00257 input->GetExtent(ext);
00258
00259 ext[0]=ext[0]-maxant;
00260 ext[1]=ext[1]+maxant;
00261 ext[2]=ext[2]-maxant;
00262 ext[3]=ext[3]+maxant;
00263 ext[4]=ext[4]-maxant;
00264 ext[5]=ext[5]+maxant;
00265
00266 pad->SetOutputWholeExtent(ext);
00267
00268 trans->SetTranslation (maxant, maxant, maxant);
00269
00270 trans->Update();
00271
00272
00273
00274 this->vtkProcessObject::SetNthInput(0, trans->GetOutput());
00275
00276
00277
00278 this->extrac->SetInput(this->GetInput());
00279
00280 }
00281
00282
00283
00284
00285 void axisExtractor02::distanciaejes(vtkPolyData *eje1, vtkPolyData *eje2)
00286 {
00287 vtkIdType totalpuntos1=eje1->GetNumberOfPoints();
00288 vtkIdType totalpuntos2=eje2->GetNumberOfPoints();
00289 vtkIdType totallineas=eje2->GetNumberOfLines();
00290
00291 vtkCell * lineas;
00292
00293 int i, j;
00294
00295 double point[3], point2[3], point3[3], point4[3];
00296
00297 double v[3], v2[3];
00298
00299 double x, dist, min, y;
00300 FILE *ff;
00301
00302 ff=fopen("comparacion.txt","w");
00303
00304
00305
00306 for(i=0;i<totalpuntos1;i++){
00307 eje1->GetPoint(i, point);
00308 min =5000;
00309 for(j=0;j<totallineas;j++){
00310 lineas=eje2->GetCell(j);
00311
00312 lineas->GetPoints()->GetPoint(0, point3);
00313 lineas->GetPoints()->GetPoint(1, point2);
00314 v[0]=point3[0]-point2[0];
00315 v[1]=point3[1]-point2[1];
00316 v[2]=point3[2]-point2[2];
00317
00318 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]));
00319
00320
00321 point4[0]=point2[0]+(x*v[0]);
00322 point4[1]=point2[1]+(x*v[1]);
00323 point4[2]=point2[2]+(x*v[2]);
00324
00325 v2[0]=point4[0]-point[0];
00326 v2[1]=point4[1]-point[1];
00327 v2[2]=point4[2]-point[2];
00328
00329 y=(v[0]*v2[0])+(v[1]*v2[1])+(v[2]*v2[2]);
00330
00331 if(x>=0 && x<=1){
00332 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])));
00333 if(dist<min){
00334 min=dist;
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349 }
00350 }
00351
00352
00353 }
00354 if(min<5){
00355 fprintf(ff,"%f\n", min);
00356 }
00357
00358
00359
00360 }
00361
00362 fclose(ff);
00363
00364
00365
00366 }
00367
00368
00369
00370
00371
00372 void axisExtractor02::Execute()
00373 {
00374 flagg=0;
00375 flagg2=0;
00376 frama=0;
00377 fseg=0;
00378 mejordst=0;
00379 mejorrad=0;
00380 mejorcant=0;
00381
00382
00383
00384 time_t start,end;
00385 double dif;
00386
00387
00388 this->GetInput()->GetExtent(extprin0);
00389 this->GetInput()->GetExtent(extprin);
00390 this->GetInput()->GetSpacing(espprin);
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404 this->data4->SetScalarTypeToUnsignedChar();
00405 this->data4->SetExtent(this->GetInput()->GetExtent());
00406 this->data4->SetSpacing(this->GetInput()->GetSpacing());
00407
00408 this->blanquear(this->data4);
00409
00410 this->points = vtkPoints::New();
00411 this->buenos=0;
00412 this->lineas = vtkCellArray::New();
00413
00414
00415 time (&start);
00416 this->todo();
00417
00418
00419
00420
00421 ((vtkPolyData *)this->Outputs[0])->SetPoints (this->points);
00422 ((vtkPolyData *)this->Outputs[0])->SetLines(this->lineas);
00423 time (&end);
00424 dif = difftime (end,start);
00425
00426
00427
00428
00429
00430
00431
00432 }
00433
00434
00435
00436
00437
00438
00439
00440 vtkImageData *axisExtractor02::GetInput()
00441 {
00442 if (this->NumberOfInputs < 1){
00443 return NULL;
00444 }
00445
00446 return (vtkImageData *)(this->Inputs[0]);
00447 }
00448
00449
00450
00451
00452
00453
00454
00455
00456 void axisExtractor02::PrintSelf(ostream& os, vtkIndent indent)
00457 {
00458 this->Superclass::PrintSelf(os,indent);
00459 }
00460
00461
00462
00463
00464
00465 void axisExtractor02::realtoreal(double a[3], double b[3] )
00466 {
00467
00468 b[0]=a[0]+(espprin[0]*extprin[0]);
00469 b[1]=a[1]+(espprin[1]*extprin[2]);
00470 b[2]=a[2];
00471
00472
00473
00474 }
00475
00476
00477
00478
00479 void axisExtractor02::realtoreal2(double a[3], double b[3] )
00480 {
00481
00482 b[0]=a[0]-(espprin[0]*extprin[0]);
00483 b[1]=a[1]-(espprin[1]*extprin[2]);
00484 b[2]=a[2];
00485
00486 }
00487
00488
00489
00490 void axisExtractor02::realtoindex(double a[3], int b[3] )
00491 {
00492 double c[3];
00493
00494 double minspac;
00495
00496
00497
00498
00499 minspac=espprin[0];
00500 if (espprin[1]<minspac) { minspac=espprin[1]; }
00501 if (espprin[2]<minspac) { minspac=espprin[2]; }
00502
00503
00504 c[0]=(a[0]/minspac);
00505 c[1]=(a[1]/minspac);
00506 c[2]=(a[2]/minspac);
00507
00508 b[0]=(int)(a[0]/minspac);
00509 b[1]=(int)(a[1]/minspac);
00510 b[2]=(int)(a[2]/minspac);
00511
00512 if(c[0]-b[0]>0.5){
00513 b[0]+=1;
00514 }
00515 if(c[1]-b[1]>0.5){
00516 b[1]+=1;
00517 }
00518 if(c[2]-b[2]>0.5){
00519 b[2]+=1;
00520 }
00521
00522
00523 }
00524
00525
00526
00527
00528
00529
00530
00531
00532 void axisExtractor02::indextoreal(int a[3], double b[3] )
00533 {
00534 double minspac;
00535
00536
00537
00538 minspac=espprin[0];
00539 if (espprin[1]<minspac) { minspac=espprin[1]; }
00540 if (espprin[2]<minspac) { minspac=espprin[2]; }
00541
00542 b[0]=(a[0])*minspac;
00543 b[1]=(a[1])*minspac;
00544 b[2]=a[2]*minspac;
00545 }
00546
00547
00548
00549
00550 void axisExtractor02::indextoreal(double a[3], double b[3] )
00551 {
00552 double minspac;
00553
00554
00555
00556
00557 minspac=espprin[0];
00558 if (espprin[1]<minspac) { minspac=espprin[1]; }
00559 if (espprin[2]<minspac) { minspac=espprin[2]; }
00560
00561 b[0]=(a[0])*minspac;
00562 b[1]=(a[1])*minspac;
00563 b[2]=a[2]*minspac;
00564 }
00565
00566
00567
00568
00569 double axisExtractor02::distanciaejepunto(double point[3], double point2[3], double point3[3])
00570 {
00571
00572
00573 double point4[3];
00574
00575 double v[3], v2[3], v3[3];
00576
00577 double x, dist=0, y;
00578
00579
00580 v[0]=point3[0]-point2[0];
00581 v[1]=point3[1]-point2[1];
00582 v[2]=point3[2]-point2[2];
00583
00584 v3[0]=point[0]-point2[0];
00585 v3[1]=point[1]-point2[1];
00586 v3[2]=point[2]-point2[2];
00587
00588 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]));
00589
00590
00591 point4[0]=point2[0]+(x*v[0]);
00592 point4[1]=point2[1]+(x*v[1]);
00593 point4[2]=point2[2]+(x*v[2]);
00594
00595 v2[0]=point4[0]-point[0];
00596 v2[1]=point4[1]-point[1];
00597 v2[2]=point4[2]-point[2];
00598
00599 y=(v[0]*v2[0])+(v[1]*v2[1])+(v[2]*v2[2]);
00600
00601 dist=sqrt((v2[0]*v2[0])+(v2[1]*v2[1])+(v2[2]*v2[2]));
00602
00603 return dist;
00604
00605 }
00606
00607
00608
00609
00610 double axisExtractor02::proporcioejepunto(double point[3], double point2[3], double point3[3])
00611 {
00612
00613
00614 double v[3], v3[3];
00615
00616 double x;
00617
00618 v[0]=point3[0]-point2[0];
00619 v[1]=point3[1]-point2[1];
00620 v[2]=point3[2]-point2[2];
00621
00622 v3[0]=point[0]-point2[0];
00623 v3[1]=point[1]-point2[1];
00624 v3[2]=point[2]-point2[2];
00625
00626 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]));
00627
00628 return x;
00629
00630
00631
00632
00633 }
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643 void axisExtractor02::searc(int i, int j, int k, vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
00644 {
00645
00646 unsigned char *ptr;
00647 unsigned char *ptr2;
00648
00649 int radio;
00650
00651 int i2, j2, k2;
00652 int i3, j3, k3;
00653
00654 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
00655 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
00656
00657 int ext[6];
00658
00659 data->GetExtent(ext);
00660
00661 radio=(ext[1]-ext[0])/2;
00662
00663 i2=i-ext[0]-radio;
00664 j2=j-ext[2]-radio;
00665 k2=k-ext[4]-radio;
00666
00667 *ptr2=label;
00668 vector[label-1][0]+=1;
00669 vector[label-1][1]+=i;
00670 vector[label-1][2]+=j;
00671 vector[label-1][3]+=k;
00672
00673 for(i3=-1;i3<=1;i3++){
00674 for(j3=-1;j3<=1;j3++){
00675 for(k3=-1;k3<=1;k3++){
00676 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] ){
00677 ptr=(unsigned char *) data->GetScalarPointer(i+i3,j+j3,k+k3);
00678 ptr2=(unsigned char *) data2->GetScalarPointer(i+i3,j+j3,k+k3);
00679 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)))){
00680 searc(i+i3, j+j3, k+k3, data, data2, label, vector);
00681 }
00682 }
00683 }
00684 }
00685 }
00686
00687 }
00688
00689
00690
00691
00692
00693 void axisExtractor02::searcb(int i, int j, int k, vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
00694 {
00695
00696 unsigned char *ptr;
00697 unsigned char *ptr2;
00698
00699 int radio;
00700
00701 int i2, j2, k2;
00702 int i3, j3, k3;
00703
00704 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
00705 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
00706
00707 int ext[6];
00708
00709 data->GetExtent(ext);
00710
00711 radio=(ext[1]-ext[0])/2;
00712
00713 i2=i-ext[0]-radio;
00714 j2=j-ext[2]-radio;
00715 k2=k-ext[4]-radio;
00716
00717 *ptr2=label;
00718 vector[label-1][0]+=1;
00719 vector[label-1][1]+=i;
00720 vector[label-1][2]+=j;
00721 vector[label-1][3]+=k;
00722
00723 for(i3=-1;i3<=1;i3++){
00724 for(j3=-1;j3<=1;j3++){
00725 for(k3=-1;k3<=1;k3++){
00726 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] ){
00727 ptr=(unsigned char *) data->GetScalarPointer(i+i3,j+j3,k+k3);
00728 ptr2=(unsigned char *) data2->GetScalarPointer(i+i3,j+j3,k+k3);
00729 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)))){
00730 searcb(i+i3, j+j3, k+k3, data, data2, label, vector);
00731 }
00732 }
00733 }
00734 }
00735 }
00736
00737 }
00738
00739
00740
00741
00742
00743 unsigned char axisExtractor02::find_components(vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
00744 {
00745 int ext[6];
00746 int i, j, k, radio;
00747 int i2, j2, k2;
00748
00749 data->GetExtent(ext);
00750 double *ff= data->GetOrigin();
00751
00752
00753 unsigned char *ptr;
00754 unsigned char *ptr2;
00755
00756 radio=(ext[1]-ext[0])/2;
00757
00758 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
00759 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
00760 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
00761 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
00762 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
00763
00764 if(*ptr==255 && *ptr2==0 && (radio*radio)>=((i2*i2)+(j2*j2)+(k2*k2)) && ((radio-1)*(radio-1))<=((i2*i2)+(j2*j2)+(k2*k2)) ){
00765 label=label+1;
00766 vector[label-1][0]=0;
00767 vector[label-1][1]=0;
00768 vector[label-1][2]=0;
00769 vector[label-1][3]=0;
00770 searc(i, j, k, data, data2, label, vector);
00771
00772 }
00773 }
00774 }
00775 }
00776
00777 return label;
00778 }
00779
00780
00781
00782
00783
00784 unsigned char axisExtractor02::find_componentsb(vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
00785 {
00786 int ext[6];
00787 int i, j, k, radio;
00788 int i2, j2, k2;
00789
00790 data->GetExtent(ext);
00791
00792 unsigned char *ptr;
00793 unsigned char *ptr2;
00794
00795 radio=(ext[1]-ext[0])/2;
00796
00797 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
00798 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
00799 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
00800 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
00801 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
00802
00803 if(*ptr==0 && *ptr2==0 && (radio*radio)>=((i2*i2)+(j2*j2)+(k2*k2)) && ((radio-1)*(radio-1))<=((i2*i2)+(j2*j2)+(k2*k2)) ){
00804 label=label+1;
00805 vector[label-1][0]=0;
00806 vector[label-1][1]=0;
00807 vector[label-1][2]=0;
00808 vector[label-1][3]=0;
00809 searcb(i, j, k, data, data2, label, vector);
00810
00811 }
00812 }
00813 }
00814 }
00815
00816 return label;
00817 }
00818
00819
00820
00821
00822
00823
00824 int axisExtractor02::proporcion(vtkImageData *data )
00825 {
00826 int ext[6];
00827 int i, j, k;
00828 int max=0;
00829 int total=0;
00830
00831 data->GetExtent(ext);
00832
00833 unsigned short *ptr;
00834
00835 for(i=ext[0];i<=ext[1];i++){
00836 for(j=ext[2];j<=ext[3];j++){
00837 for(k=ext[4];k<=ext[5];k++){
00838 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
00839 if(*ptr!=0){
00840 max++;
00841 }
00842 total++;
00843 }
00844 }
00845 }
00846
00847 return (max*100)/total;
00848 }
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858 bool axisExtractor02::border(vtkImageData *data, int p1[3] )
00859 {
00860 int ext[6];
00861 int i, j, k;
00862 short val;
00863 int total1, total2;
00864 int centro[3];
00865
00866
00867 int p2[3];
00868
00869 data->GetExtent(ext);
00870
00871 centro[0]=(ext[1]+ext[0])/2;
00872 centro[1]=(ext[3]+ext[2])/2;
00873 centro[2]=(ext[5]+ext[4])/2;
00874
00875 unsigned char *ptr;
00876
00877 ptr=(unsigned char *) data->GetScalarPointer(p1);
00878
00879 val=*ptr;
00880
00881 bool res=false;
00882 total1=0;
00883 total2=0;
00884
00885 for(i=-1;i<=1;i++){
00886 for(j=-1;j<=1;j++){
00887 for(k=-1;k<=1;k++){
00888 p2[0]=p1[0]+i;
00889 p2[1]=p1[1]+j;
00890 p2[2]=p1[2]+k;
00891 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] ){
00892 ptr=(unsigned char *) data->GetScalarPointer(p2);
00893 if(*ptr!=val){
00894 total1++;
00895 }
00896 else{
00897 total2++;
00898 }
00899 }
00900 }
00901 }
00902 }
00903
00904 if(total1>0){
00905 res=true;
00906 }
00907
00908
00909
00910
00911 return res;
00912 }
00913
00914
00915
00916
00917
00918
00919
00920 void axisExtractor02::optim(vtkImageData *data, vtkImageData *data2 )
00921 {
00922 int ext[6];
00923 int i, j, k, radio;
00924 int centro[3];
00925
00926 double w0, w1;
00927
00928 double w0p, w1p;
00929
00930 double inerciaxx, inerciayy, inerciazz, inerciaxy, inerciayz, inerciazx;
00931 double sumx, sumy, sumz;
00932 double sumxx, sumyy, sumzz, sumxy, sumyz, sumzx;
00933
00934 double sumix, sumiy, sumiz;
00935 double sumixx, sumiyy, sumizz, sumixy, sumiyz, sumizx;
00936
00937 double inerciaixx, inerciaiyy, inerciaizz, inerciaixy, inerciaiyz, inerciaizx;
00938
00939 double inerciaxx2, inerciayy2, inerciazz2, inerciaxy2, inerciayz2, inerciazx2;
00940 double sumx2, sumy2, sumz2;
00941 double sumxx2, sumyy2, sumzz2, sumxy2, sumyz2, sumzx2;
00942
00943 double sumix2, sumiy2, sumiz2;
00944 double sumixx2, sumiyy2, sumizz2, sumixy2, sumiyz2, sumizx2;
00945
00946
00947
00948
00949 double inerciaxxp, inerciayyp, inerciazzp, inerciaxyp, inerciayzp, inerciazxp;
00950 double sumxp, sumyp, sumzp;
00951 double sumxxp, sumyyp, sumzzp, sumxyp, sumyzp, sumzxp;
00952
00953 double inerciaxx2p, inerciayy2p, inerciazz2p, inerciaxy2p, inerciayz2p, inerciazx2p;
00954 double sumx2p, sumy2p, sumz2p;
00955 double sumxx2p, sumyy2p, sumzz2p, sumxy2p, sumyz2p, sumzx2p;
00956
00957
00958
00959 unsigned long cantp;
00960 double sumip;
00961 double sumiip;
00962 double inerciaiip;
00963 double centip;
00964
00965 unsigned long cant2p;
00966 double sumi2p;
00967 double sumii2p;
00968 double inerciaii2p;
00969 double centi2p;
00970
00971 int point[3];
00972
00973
00974 double sumitotal=0;
00975 double distcent=0;
00976
00977 double im, jm, km, lm;
00978
00979
00980 max2=0;
00981 cant=0;
00982 sumi=0;
00983 sumii=0;
00984 des2=0;
00985 min2=3200;
00986 centi=0;
00987
00988 max3=0;
00989 cant2=0;
00990 sumi2=0;
00991 sumii2=0;
00992 des3=0;
00993 min3=3200;
00994 centi2=0;
00995
00996 sumx=0;
00997 sumy=0;
00998 sumz=0;
00999
01000 sumxx=0;
01001 sumyy=0;
01002 sumzz=0;
01003
01004 sumxy=0;
01005 sumyz=0;
01006 sumzx=0;
01007
01008 sumix=0;
01009 sumiy=0;
01010 sumiz=0;
01011
01012 sumixx=0;
01013 sumiyy=0;
01014 sumizz=0;
01015
01016 sumixy=0;
01017 sumiyz=0;
01018 sumizx=0;
01019
01020 sumx2=0;
01021 sumy2=0;
01022 sumz2=0;
01023
01024 sumxx2=0;
01025 sumyy2=0;
01026 sumzz2=0;
01027
01028 sumxy2=0;
01029 sumyz2=0;
01030 sumzx2=0;
01031
01032 sumix2=0;
01033 sumiy2=0;
01034 sumiz2=0;
01035
01036 sumixx2=0;
01037 sumiyy2=0;
01038 sumizz2=0;
01039
01040 sumixy2=0;
01041 sumiyz2=0;
01042 sumizx2=0;
01043
01044 max=0;
01045 min=3200;
01046
01047
01048 data->GetExtent(ext);
01049
01050 unsigned short *ptr;
01051 unsigned char *ptr2;
01052
01053 radio=(ext[1]-ext[0])/2;
01054
01055 centro[0]=(ext[1]+ext[0])/2;
01056 centro[1]=(ext[3]+ext[2])/2;
01057 centro[2]=(ext[5]+ext[4])/2;
01058
01059 double tmpsqrt;
01060
01061 for(i=ext[0];i<=ext[1];i++){
01062 for(j=ext[2];j<=ext[3];j++){
01063 for(k=ext[4];k<=ext[5];k++){
01064 tmpsqrt=((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
01065 if( radio>=sqrt( tmpsqrt) ){
01066 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
01067 if(*ptr>max){
01068 max=*ptr;
01069
01070 }
01071 if(*ptr<min){
01072 min=*ptr;
01073
01074 }
01075 }
01076 }
01077 }
01078 }
01079
01080 for(i=ext[0];i<=ext[1];i++){
01081 for(j=ext[2];j<=ext[3];j++){
01082 for(k=ext[4];k<=ext[5];k++){
01083 tmpsqrt=((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
01084 if( radio>=sqrt(tmpsqrt) ){
01085 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
01086 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
01087 point[0]=i;
01088 point[1]=j;
01089 point[2]=k;
01090
01091 im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
01092 jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
01093 km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
01094 lm=2*(((double)*ptr-(double)min)/((double)max-(double)min));
01095
01096 if(*ptr2!=0){
01097 sumi+=lm;
01098 sumii+=lm*lm;
01099 cant+=1;
01100
01101 sumx+=im;
01102 sumy+=jm;
01103 sumz+=km;
01104
01105 sumxx+=im*im;
01106 sumyy+=jm*jm;
01107 sumzz+=km*km;
01108 sumxy+=im*jm;
01109 sumyz+=jm*km;
01110 sumzx+=km*im;
01111
01112 sumix+=lm*im;
01113 sumiy+=lm*jm;
01114 sumiz+=lm*km;
01115
01116 sumixx+=lm*im*im;
01117 sumiyy+=lm*jm*jm;
01118 sumizz+=lm*km*km;
01119 sumixy+=lm*im*jm;
01120 sumiyz+=lm*jm*km;
01121 sumizx+=lm*km*im;
01122 }
01123 else{
01124 sumi2+=lm;
01125 sumii2+=lm*lm;
01126 cant2+=1;
01127
01128 sumx2+=im;
01129 sumy2+=jm;
01130 sumz2+=km;
01131
01132 sumxx2+=im*im;
01133 sumyy2+=jm*jm;
01134 sumzz2+=km*km;
01135 sumxy2+=im*jm;
01136 sumyz2+=jm*km;
01137 sumzx2+=km*im;
01138
01139 sumix2+=lm*im;
01140 sumiy2+=lm*jm;
01141 sumiz2+=lm*km;
01142
01143 sumixx2+=lm*im*im;
01144 sumiyy2+=lm*jm*jm;
01145 sumizz2+=lm*km*km;
01146 sumixy2+=lm*im*jm;
01147 sumiyz2+=lm*jm*km;
01148 sumizx2+=lm*km*im;
01149 }
01150 }
01151 }
01152 }
01153 }
01154
01155 centi=(double)sumi/(double)cant;
01156 inerciaii= ((double)sumii/(double)cant)-(centi*centi);
01157
01158 centi2=(double)sumi2/(double)cant2;
01159 inerciaii2= ((double)sumii2/(double)cant2)-(centi2*centi2);
01160
01161 centx=(double)sumx/(double)cant;
01162 centy=(double)sumy/(double)cant;
01163 centz=(double)sumz/(double)cant;
01164
01165 centx2=(double)sumx2/(double)cant2;
01166 centy2=(double)sumy2/(double)cant2;
01167 centz2=(double)sumz2/(double)cant2;
01168
01169 inerciaxx= ((double)sumxx/(double)cant)-centx*centx;
01170 inerciayy= ((double)sumyy/(double)cant)-centy*centy;
01171 inerciazz= ((double)sumzz/(double)cant)-centz*centz;
01172
01173 inerciaxx2= ((double)sumxx2/(double)cant2)-centx2*centx2;
01174 inerciayy2= ((double)sumyy2/(double)cant2)-centy2*centy2;
01175 inerciazz2= ((double)sumzz2/(double)cant2)-centz2*centz2;
01176
01177 inerciaxy= ((double)sumxy/(double)cant)-centx*centy;
01178 inerciayz= ((double)sumyz/(double)cant)-centy*centz;
01179 inerciazx= ((double)sumzx/(double)cant)-centz*centx;
01180
01181 inerciaxy2= ((double)sumxy2/(double)cant2)-centx2*centy2;
01182 inerciayz2= ((double)sumyz2/(double)cant2)-centy2*centz2;
01183 inerciazx2= ((double)sumzx2/(double)cant2)-centz2*centx2;
01184
01185
01186 centix=((double)sumix+(double)sumix2)/((double)sumi+(double)sumi2);
01187 centiy=((double)sumiy+(double)sumiy2)/((double)sumi+(double)sumi2);
01188 centiz=((double)sumiz+(double)sumiz2)/((double)sumi+(double)sumi2);
01189
01190 inerciaixx= ( ((double)sumixx+(double)sumixx2) / ((double)sumi+(double)sumi2) )-(centix*centix);
01191 inerciaiyy= ( ((double)sumiyy+(double)sumiyy2) / ((double)sumi+(double)sumi2) )-(centiy*centiy);
01192 inerciaizz= ( ((double)sumizz+(double)sumizz2) / ((double)sumi+(double)sumi2) )-(centiz*centiz);
01193
01194
01195 inerciaixy= ( ((double)sumixy+(double)sumixy2) / ((double)sumi+(double)sumi2) )-(centix*centiy);
01196 inerciaiyz= ( ((double)sumiyz+(double)sumiyz2) / ((double)sumi+(double)sumi2) )-(centiy*centiz);
01197 inerciaizx= ( ((double)sumizx+(double)sumizx2) / ((double)sumi+(double)sumi2) )-(centiz*centix);
01198
01199
01200 A[0][0]=inerciaxx;
01201 A[1][1]=inerciayy;
01202 A[2][2]=inerciazz;
01203 A[0][1]=A[1][0]=inerciaxy;
01204 A[0][2]=A[2][0]=inerciazx;
01205 A[1][2]=A[2][1]=inerciayz;
01206
01207 A2[0][0]=inerciaxx2;
01208 A2[1][1]=inerciayy2;
01209 A2[2][2]=inerciazz2;
01210 A2[0][1]=A2[1][0]=inerciaxy2;
01211 A2[0][2]=A2[2][0]=inerciazx2;
01212 A2[1][2]=A2[2][1]=inerciayz2;
01213
01214 Ai[0][0]=inerciaixx;
01215 Ai[1][1]=inerciaiyy;
01216 Ai[2][2]=inerciaizz;
01217 Ai[0][1]=Ai[1][0]=inerciaixy;
01218 Ai[0][2]=Ai[2][0]=inerciaizx;
01219 Ai[1][2]=Ai[2][1]=inerciaiyz;
01220
01221 vtkMath::Diagonalize3x3 (A, w, V);
01222 vtkMath::Diagonalize3x3 (A2, w2, V2);
01223 vtkMath::Diagonalize3x3 (Ai, wi, Vi);
01224
01225
01226 if(w[0]>w[1]){
01227 if(w[0]>w[2]){
01228 ejemin=0;
01229 inerciar=w[1]+w[2];
01230 if(w[1]>w[2]){
01231 inerciary=w[0]+w[2];
01232 inerciarz=w[0]+w[1];
01233 }
01234 else{
01235 inerciary=w[0]+w[1];
01236 inerciarz=w[0]+w[2];
01237 }
01238 }
01239 else{
01240 ejemin=2;
01241 inerciar=w[0]+w[1];
01242 if(w[1]>w[0]){
01243 inerciary=w[0]+w[2];
01244 inerciarz=w[1]+w[2];
01245 }
01246 else{
01247 inerciary=w[1]+w[2];
01248 inerciarz=w[0]+w[2];
01249 }
01250 }
01251 }
01252 else{
01253 if(w[1]>w[2]){
01254 ejemin=1;
01255 inerciar=w[2]+w[0];
01256 if(w[2]>w[0]){
01257 inerciary=w[1]+w[0];
01258 inerciarz=w[1]+w[2];
01259 }
01260 else{
01261 inerciary=w[1]+w[2];
01262 inerciarz=w[1]+w[0];
01263 }
01264 }
01265 else{
01266 ejemin=2;
01267 inerciar=w[1]+w[0];
01268 if(w[1]>w[0]){
01269 inerciary=w[0]+w[2];
01270 inerciarz=w[1]+w[2];
01271 }
01272 else{
01273 inerciary=w[1]+w[2];
01274 inerciarz=w[0]+w[2];
01275 }
01276 }
01277 }
01278
01279
01280 if(w2[0]>w2[1]){
01281 if(w2[0]>w2[2]){
01282 ejemin2=0;
01283 inercia2r=w2[1]+w2[2];
01284 }
01285 else{
01286 ejemin2=2;
01287 inercia2r=w2[1]+w2[0];
01288 }
01289 }
01290 else{
01291 if(w2[1]>w2[2]){
01292 ejemin2=1;
01293 inercia2r=w2[2]+w2[0];
01294 }
01295 else{
01296 ejemin2=2;
01297 inercia2r=w2[1]+w2[0];
01298 }
01299 }
01300
01301
01302 if(wi[0]>wi[1]){
01303 if(wi[0]>wi[2]){
01304 ejemini=0;
01305 inerciari=wi[1]+wi[2];
01306 if(wi[1]>wi[2]){
01307 inerciariy=wi[0]+wi[2];
01308 inerciariz=wi[0]+wi[1];
01309 }
01310 else{
01311 inerciariy=wi[0]+wi[1];
01312 inerciariz=wi[0]+wi[2];
01313 }
01314 }
01315 else{
01316 ejemini=2;
01317 inerciari=wi[1]+wi[0];
01318 if(wi[1]>wi[0]){
01319 inerciariy=wi[0]+wi[2];
01320 inerciariz=wi[1]+wi[2];
01321 }
01322 else{
01323 inerciariy=wi[1]+wi[2];
01324 inerciariz=wi[0]+wi[2];
01325 }
01326
01327 }
01328 }
01329 else{
01330 if(wi[1]>wi[2]){
01331 ejemini=1;
01332 inerciari=wi[2]+wi[0];
01333 if(wi[2]>wi[0]){
01334 inerciariy=wi[1]+wi[0];
01335 inerciariz=wi[1]+wi[2];
01336 }
01337 else{
01338 inerciariy=wi[1]+wi[2];
01339 inerciariz=wi[1]+wi[0];
01340 }
01341 }
01342 else{
01343 ejemini=2;
01344 inerciari=wi[1]+wi[0];
01345 if(wi[1]>wi[0]){
01346 inerciariy=wi[0]+wi[2];
01347 inerciariz=wi[1]+wi[2];
01348 }
01349 else{
01350 inerciariy=wi[1]+wi[2];
01351 inerciariz=wi[0]+wi[2];
01352 }
01353 }
01354 }
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426 w0=(double)cant/((double)cant+(double)cant2);
01427 w1=(double)cant2/((double)cant+(double)cant2);
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442 costo=param*w0*inerciaii+w1*param2*inerciaii2+w0*param3*inerciar+w1*param4*inercia2r;
01443
01444 int changes=4;
01445
01446 int total=0;
01447
01448 while(changes>3 && total<500){
01449 changes=0;
01450 total++;
01451 for(i=ext[0];i<=ext[1];i++){
01452 for(j=ext[2];j<=ext[3];j++){
01453 for(k=ext[4];k<=ext[5];k++){
01454 tmpsqrt = ((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
01455 if( radio>=sqrt(tmpsqrt) ){
01456 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
01457 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
01458 point[0]=i;
01459 point[1]=j;
01460 point[2]=k;
01461
01462 im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
01463 jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
01464 km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
01465 lm=2*(((double)*ptr-(double)min)/((double)max-(double)min));
01466
01467 if(border(data2, point)){
01468 if(*ptr2!=0){
01469
01470 sumip=sumi-lm;
01471 sumiip=sumii-lm*lm;
01472
01473 cantp=cant-1;
01474
01475 sumi2p=sumi2+lm;
01476 sumii2p=sumii2+lm*lm;
01477
01478 cant2p=cant2+1;
01479
01480 sumxp=sumx-im;
01481 sumyp=sumy-jm;
01482 sumzp=sumz-km;
01483 sumxxp=sumxx-im*im;
01484 sumyyp=sumyy-jm*jm;
01485 sumzzp=sumzz-km*km;
01486 sumxyp=sumxy-im*jm;
01487 sumyzp=sumyz-jm*km;
01488 sumzxp=sumzx-km*im;
01489
01490 sumx2p=sumx2+im;
01491 sumy2p=sumy2+jm;
01492 sumz2p=sumz2+km;
01493 sumxx2p=sumxx2+im*im;
01494 sumyy2p=sumyy2+jm*jm;
01495 sumzz2p=sumzz2+km*km;
01496 sumxy2p=sumxy2+im*jm;
01497 sumyz2p=sumyz2+jm*km;
01498 sumzx2p=sumzx2+km*im;
01499
01500 centip=(double)sumip/(double)cantp;
01501 inerciaiip= ((double)sumiip/(double)cantp)-(centip*centip);
01502
01503 centi2p=(double)sumi2p/(double)cant2p;
01504 inerciaii2p= ((double)sumii2p/(double)cant2p)-(centi2p*centi2p);
01505
01506 centxp=(double)sumxp/(double)cantp;
01507 centyp=(double)sumyp/(double)cantp;
01508 centzp=(double)sumzp/(double)cantp;
01509
01510 centx2p=(double)sumx2p/(double)cant2p;
01511 centy2p=(double)sumy2p/(double)cant2p;
01512 centz2p=(double)sumz2p/(double)cant2p;
01513
01514 inerciaxxp= ((double)sumxxp/(double)cantp)-centxp*centxp;
01515 inerciayyp= ((double)sumyyp/(double)cantp)-centyp*centyp;
01516 inerciazzp= ((double)sumzzp/(double)cantp)-centzp*centzp;
01517
01518 inerciaxx2p= ((double)sumxx2p/(double)cant2p)-centx2p*centx2p;
01519 inerciayy2p= ((double)sumyy2p/(double)cant2p)-centy2p*centy2p;
01520 inerciazz2p= ((double)sumzz2p/(double)cant2p)-centz2p*centz2p;
01521
01522 inerciaxyp= ((double)sumxyp/(double)cantp)-centxp*centyp;
01523 inerciayzp= ((double)sumyzp/(double)cantp)-centyp*centzp;
01524 inerciazxp= ((double)sumzxp/(double)cantp)-centzp*centxp;
01525
01526 inerciaxy2p= ((double)sumxy2p/(double)cant2p)-centx2p*centy2p;
01527 inerciayz2p= ((double)sumyz2p/(double)cant2p)-centy2p*centz2p;
01528 inerciazx2p= ((double)sumzx2p/(double)cant2p)-centz2p*centx2p;
01529
01530 Ap[0][0]=inerciaxxp;
01531 Ap[1][1]=inerciayyp;
01532 Ap[2][2]=inerciazzp;
01533 Ap[0][1]=Ap[1][0]=inerciaxyp;
01534 Ap[0][2]=Ap[2][0]=inerciazxp;
01535 Ap[1][2]=Ap[2][1]=inerciayzp;
01536
01537
01538 A2p[0][0]=inerciaxx2p;
01539 A2p[1][1]=inerciayy2p;
01540 A2p[2][2]=inerciazz2p;
01541 A2p[0][1]=A2p[1][0]=inerciaxy2p;
01542 A2p[0][2]=A2p[2][0]=inerciazx2p;
01543 A2p[1][2]=A2p[2][1]=inerciayz2p;
01544
01545
01546
01547 w0p=(double)cantp/((double)cantp+(double)cant2p);
01548 w1p=(double)cant2p/((double)cantp+(double)cant2p);
01549
01550
01551
01552
01553 vtkMath::Diagonalize3x3 (Ap, wp, Vp);
01554 vtkMath::Diagonalize3x3 (A2p, w2p, V2p);
01555
01556
01557 if(wp[0]>wp[1]){
01558 if(wp[0]>wp[2]){
01559 ejeminp=0;
01560 inerciarp=wp[2]+wp[1];
01561 if(wp[1]>wp[2]){
01562 inerciarpy=wp[0]+wp[2];
01563 inerciarpz=wp[0]+wp[1];
01564 }
01565 else{
01566 inerciarpy=wp[0]+wp[1];
01567 inerciarpz=wp[0]+wp[2];
01568 }
01569 }
01570 else{
01571 ejeminp=2;
01572 inerciarp=wp[1]+wp[0];
01573 if(wp[1]>wp[0]){
01574 inerciarpy=wp[0]+wp[2];
01575 inerciarpz=wp[1]+wp[2];
01576 }
01577 else{
01578 inerciarpy=wp[1]+wp[2];
01579 inerciarpz=wp[0]+wp[2];
01580 }
01581 }
01582 }
01583 else{
01584 if(wp[1]>wp[2]){
01585 ejeminp=1;
01586 inerciarp=wp[2]+wp[0];
01587 if(wp[2]>wp[0]){
01588 inerciarpy=wp[1]+wp[0];
01589 inerciarpz=wp[1]+wp[2];
01590 }
01591 else{
01592 inerciarpy=wp[1]+wp[2];
01593 inerciarpz=wp[1]+wp[0];
01594 }
01595 }
01596 else{
01597 ejeminp=2;
01598 inerciarp=wp[1]+wp[0];
01599 if(wp[1]>wp[0]){
01600 inerciarpy=wp[0]+wp[2];
01601 inerciarpz=wp[1]+wp[2];
01602 }
01603 else{
01604 inerciarpy=wp[1]+wp[2];
01605 inerciarpz=wp[0]+wp[2];
01606 }
01607 }
01608 }
01609
01610
01611 if(w2p[0]>w2p[1]){
01612 if(w2p[0]>w2p[2]){
01613 ejemin2p=0;
01614 inercia2rp=w2p[2]+w2p[1];
01615 }
01616 else{
01617 ejemin2p=2;
01618 inercia2rp=w2p[1]+w2p[0];
01619 }
01620 }
01621 else{
01622 if(w2p[1]>w2p[2]){
01623 ejemin2p=1;
01624 inercia2rp=w2p[2]+w2p[0];
01625 }
01626 else{
01627 ejemin2p=2;
01628 inercia2rp=w2p[1]+w2p[0];
01629 }
01630 }
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656 costop=param*w0p*inerciaiip+w1p*param2*inerciaii2p+w0p*param3*inerciarp+w1p*param4*inercia2rp;
01657
01658 if(costop<costo && centip>centi2p){
01659
01660 *ptr2=0;
01661 changes++;
01662
01663 sumi=sumip;
01664 sumii=sumiip;
01665 cant=cantp;
01666
01667 sumi2=sumi2p;
01668 sumii2=sumii2p;
01669 cant2=cant2p;
01670
01671 sumx=sumxp;
01672 sumy=sumyp;
01673 sumz=sumzp;
01674 sumxx=sumxxp;
01675 sumyy=sumyyp;
01676 sumzz=sumzzp;
01677 sumxy=sumxyp;
01678 sumyz=sumyzp;
01679 sumzx=sumzxp;
01680
01681 sumx2=sumx2p;
01682 sumy2=sumy2p;
01683 sumz2=sumz2p;
01684 sumxx2=sumxx2p;
01685 sumyy2=sumyy2p;
01686 sumzz2=sumzz2p;
01687 sumxy2=sumxy2p;
01688 sumyz2=sumyz2p;
01689 sumzx2=sumzx2p;
01690
01691 centi=centip;
01692 inerciaii= inerciaiip;
01693
01694 centi2=centi2p;
01695 inerciaii2p= inerciaii2p;
01696
01697 centx=centxp;
01698 centy=centyp;
01699 centz=centzp;
01700
01701 centx2=centx2p;
01702 centy2=centy2p;
01703 centz2=centz2p;
01704
01705 inerciaxx= inerciaxxp;
01706 inerciayy= inerciayyp;
01707 inerciazz= inerciazzp;
01708
01709 inerciaxx2= inerciaxx2p;
01710 inerciayy2= inerciayy2p;
01711 inerciazz2= inerciazz2p;
01712
01713 inerciaxy= inerciaxyp;
01714 inerciayz= inerciayzp;
01715 inerciazx= inerciazxp;
01716
01717 inerciaxy2= inerciaxy2p;
01718 inerciayz2= inerciayz2p;
01719 inerciazx2= inerciazx2p;
01720
01721
01722
01723 costo=costop;
01724
01725 }
01726 }
01727 else{
01728 sumip=sumi+lm;
01729 sumiip=sumii+lm*lm;
01730
01731 cantp=cant+1;
01732
01733 sumi2p=sumi2-lm;
01734 sumii2p=sumii2-lm*lm;
01735
01736 cant2p=cant2-1;
01737
01738 sumxp=sumx+im;
01739 sumyp=sumy+jm;
01740 sumzp=sumz+km;
01741 sumxxp=sumxx+im*im;
01742 sumyyp=sumyy+jm*jm;
01743 sumzzp=sumzz+km*km;
01744 sumxyp=sumxy+im*jm;
01745 sumyzp=sumyz+jm*km;
01746 sumzxp=sumzx+km*im;
01747
01748 sumx2p=sumx2-im;
01749 sumy2p=sumy2-jm;
01750 sumz2p=sumz2-km;
01751 sumxx2p=sumxx2-im*im;
01752 sumyy2p=sumyy2-jm*jm;
01753 sumzz2p=sumzz2-km*km;
01754 sumxy2p=sumxy2-im*jm;
01755 sumyz2p=sumyz2-jm*km;
01756 sumzx2p=sumzx2-km*im;
01757
01758 centip=(double)sumip/(double)cantp;
01759 inerciaiip= ((double)sumiip/(double)cantp)-(centip*centip);
01760
01761 centi2p=(double)sumi2p/(double)cant2p;
01762 inerciaii2p= ((double)sumii2p/(double)cant2p)-(centi2p*centi2p);
01763
01764 centxp=(double)sumxp/(double)cantp;
01765 centyp=(double)sumyp/(double)cantp;
01766 centzp=(double)sumzp/(double)cantp;
01767
01768 centx2p=(double)sumx2p/(double)cant2p;
01769 centy2p=(double)sumy2p/(double)cant2p;
01770 centz2p=(double)sumz2p/(double)cant2p;
01771
01772
01773 inerciaxxp= ((double)sumxxp/(double)cantp)-centxp*centxp;
01774 inerciayyp= ((double)sumyyp/(double)cantp)-centyp*centyp;
01775 inerciazzp= ((double)sumzzp/(double)cantp)-centzp*centzp;
01776
01777 inerciaxx2p= ((double)sumxx2p/(double)cant2p)-centx2p*centx2p;
01778 inerciayy2p= ((double)sumyy2p/(double)cant2p)-centy2p*centy2p;
01779 inerciazz2p= ((double)sumzz2p/(double)cant2p)-centz2p*centz2p;
01780
01781 inerciaxyp= ((double)sumxyp/(double)cantp)-centxp*centyp;
01782 inerciayzp= ((double)sumyzp/(double)cantp)-centyp*centzp;
01783 inerciazxp= ((double)sumzxp/(double)cantp)-centzp*centxp;
01784
01785 inerciaxy2p= ((double)sumxy2p/(double)cant2p)-centx2p*centy2p;
01786 inerciayz2p= ((double)sumyz2p/(double)cant2p)-centy2p*centz2p;
01787 inerciazx2p= ((double)sumzx2p/(double)cant2p)-centz2p*centx2p;
01788
01789 Ap[0][0]=inerciaxxp;
01790 Ap[1][1]=inerciayyp;
01791 Ap[2][2]=inerciazzp;
01792 Ap[0][1]=Ap[1][0]=inerciaxyp;
01793 Ap[0][2]=Ap[2][0]=inerciazxp;
01794 Ap[1][2]=Ap[2][1]=inerciayzp;
01795
01796 A2p[0][0]=inerciaxx2p;
01797 A2p[1][1]=inerciayy2p;
01798 A2p[2][2]=inerciazz2p;
01799 A2p[0][1]=A2p[1][0]=inerciaxy2p;
01800 A2p[0][2]=A2p[2][0]=inerciazx2p;
01801 A2p[1][2]=A2p[2][1]=inerciayz2p;
01802
01803 w0p=(double)cantp/((double)cantp+(double)cant2p);
01804 w1p=(double)cant2p/((double)cantp+(double)cant2p);
01805
01806
01807
01808
01809 vtkMath::Diagonalize3x3 (Ap, wp, Vp);
01810 vtkMath::Diagonalize3x3 (A2p, w2p, V2p);
01811
01812 if(wp[0]>wp[1]){
01813 if(wp[0]>wp[2]){
01814 ejeminp=0;
01815 inerciarp=wp[2]+wp[1];
01816 if(wp[1]>wp[2]){
01817 inerciarpy=wp[0]+wp[2];
01818 inerciarpz=wp[0]+wp[1];
01819 }
01820 else{
01821 inerciarpy=wp[0]+wp[1];
01822 inerciarpz=wp[0]+wp[2];
01823 }
01824 }
01825 else{
01826 ejeminp=2;
01827 inerciarp=wp[1]+wp[0];
01828 if(wp[1]>wp[0]){
01829 inerciarpy=wp[0]+wp[2];
01830 inerciarpz=wp[1]+wp[2];
01831 }
01832 else{
01833 inerciarpy=wp[1]+wp[2];
01834 inerciarpz=wp[0]+wp[2];
01835 }
01836 }
01837 }
01838 else{
01839 if(wp[1]>wp[2]){
01840 ejeminp=1;
01841 inerciarp=wp[2]+wp[0];
01842 if(wp[2]>wp[0]){
01843 inerciarpy=wp[1]+wp[0];
01844 inerciarpz=wp[1]+wp[2];
01845 }
01846 else{
01847 inerciarpy=wp[1]+wp[2];
01848 inerciarpz=wp[1]+wp[0];
01849 }
01850 }
01851 else{
01852 ejeminp=2;
01853 inerciarp=wp[1]+wp[0];
01854 if(wp[1]>wp[0]){
01855 inerciarpy=wp[0]+wp[2];
01856 inerciarpz=wp[1]+wp[2];
01857 }
01858 else{
01859 inerciarpy=wp[1]+wp[2];
01860 inerciarpz=wp[0]+wp[2];
01861 }
01862 }
01863 }
01864
01865
01866 if(w2p[0]>w2p[1]){
01867 if(w2p[0]>w2p[2]){
01868 ejemin2p=0;
01869 inercia2rp=w2p[2]+w2p[1];
01870 }
01871 else{
01872 ejemin2p=2;
01873 inercia2rp=w2p[1]+w2p[0];
01874 }
01875 }
01876 else{
01877 if(w2p[1]>w2p[2]){
01878 ejemin2p=1;
01879 inercia2rp=w2p[2]+w2p[0];
01880 }
01881 else{
01882 ejemin2p=2;
01883 inercia2rp=w2p[1]+w2p[0];
01884 }
01885 }
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912 costop=param*w0p*inerciaiip+w1p*param2*inerciaii2p+w0p*param3*inerciarp+w1p*param4*inercia2rp;
01913
01914 if(costop<costo && centip>centi2p){
01915
01916 *ptr2=255;
01917 changes++;
01918
01919 sumi=sumip;
01920 sumii=sumiip;
01921 cant=cantp;
01922
01923 sumi2=sumi2p;
01924 sumii2=sumii2p;
01925 cant2=cant2p;
01926
01927 sumx=sumxp;
01928 sumy=sumyp;
01929 sumz=sumzp;
01930 sumxx=sumxxp;
01931 sumyy=sumyyp;
01932 sumzz=sumzzp;
01933 sumxy=sumxyp;
01934 sumyz=sumyzp;
01935 sumzx=sumzxp;
01936
01937 sumx2=sumx2p;
01938 sumy2=sumy2p;
01939 sumz2=sumz2p;
01940 sumxx2=sumxx2p;
01941 sumyy2=sumyy2p;
01942 sumzz2=sumzz2p;
01943 sumxy2=sumxy2p;
01944 sumyz2=sumyz2p;
01945 sumzx2=sumzx2p;
01946
01947 centi=centip;
01948 inerciaii= inerciaiip;
01949
01950 centi2=centi2p;
01951 inerciaii2p= inerciaii2p;
01952
01953 centx=centxp;
01954 centy=centyp;
01955 centz=centzp;
01956
01957 centx2=centx2p;
01958 centy2=centy2p;
01959 centz2=centz2p;
01960
01961 inerciaxx= inerciaxxp;
01962 inerciayy= inerciayyp;
01963 inerciazz= inerciazzp;
01964
01965 inerciaxx2= inerciaxx2p;
01966 inerciayy2= inerciayy2p;
01967 inerciazz2= inerciazz2p;
01968
01969 inerciaxy= inerciaxyp;
01970 inerciayz= inerciayzp;
01971 inerciazx= inerciazxp;
01972
01973 inerciaxy2= inerciaxy2p;
01974 inerciayz2= inerciayz2p;
01975 inerciazx2= inerciazx2p;
01976
01977
01978
01979 costo=costop;
01980
01981 }
01982
01983 }
01984 }
01985 }
01986 }
01987 }
01988 }
01989
01990
01991
01992
01993 }
01994
01995
01996
01997
01998 }
01999
02000
02001
02002
02003
02004
02005 void axisExtractor02::costominimo(vtkImageData *data, vtkImageData *data2 )
02006 {
02007 int ext[6];
02008 int i, j, k, radio;
02009 int centro[3];
02010
02011
02012 data->GetExtent(ext);
02013
02014 double *ptr;
02015 unsigned char *ptr2;
02016
02017 int ind;
02018
02019 radio=(ext[1]-ext[0])/2;
02020
02021 centro[0]=(ext[1]+ext[0])/2;
02022 centro[1]=(ext[3]+ext[2])/2;
02023 centro[2]=(ext[5]+ext[4])/2;
02024
02025
02026 for(i=0;i<=10;i++){
02027 minis[i]=0;
02028 }
02029
02030 for(i=ext[0];i<=ext[1];i++){
02031 for(j=ext[2];j<=ext[3];j++){
02032 for(k=ext[4];k<=ext[5];k++){
02033 ptr=(double *)data->GetScalarPointer(i,j,k);
02034 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
02035
02036 ind=(int)*ptr2;
02037
02038 if(*ptr2!=0){
02039 if(minis[ind-1]<=*ptr && *ptr!=0 && ind<11){
02040 minis[ind-1]=*ptr;
02041 candit[ind-1][0]=i;
02042 candit[ind-1][1]=j;
02043 candit[ind-1][2]=k;
02044 }
02045 }
02046 }
02047 }
02048 }
02049 }
02050
02051
02052
02053
02054
02055
02056 void axisExtractor02::costominimo2(vtkImageData *data, vtkImageData *data3, int p1[3], int p2[3], int p3[3])
02057 {
02058 int ext[6];
02059 int i, j, k, radio;
02060 int punto[3];
02061 int puntob[3];
02062 int punto2[3];
02063 int puntob2[3];
02064 int puntomin[3];
02065 int puntobmin[3];
02066
02067 double *ptr;
02068 double *ptr2;
02069 unsigned char *ptr3;
02070 unsigned char *ptr4;
02071 double *ptr5;
02072 unsigned char *ptr6;
02073 double *ptr7;
02074 double *ptr8;
02075
02076 vtkImageData *data4;
02077
02078 data4=vtkImageData::New();
02079 data4->SetScalarTypeToUnsignedChar();
02080 data4->SetExtent(data->GetExtent());
02081 data4->SetSpacing(data->GetSpacing());
02082
02083 vtkImageData *data6;
02084
02085 data6=vtkImageData::New();
02086 data6->SetScalarTypeToUnsignedChar();
02087 data6->SetExtent(data->GetExtent());
02088 data6->SetSpacing(data->GetSpacing());
02089
02090
02091 vtkImageData *data2;
02092
02093 data2=vtkImageData::New();
02094 data2->SetScalarTypeToDouble();
02095 data2->SetExtent(data->GetExtent());
02096 data2->SetSpacing(data->GetSpacing());
02097
02098 vtkImageData *data8;
02099
02100 data8=vtkImageData::New();
02101 data8->SetScalarTypeToDouble();
02102 data8->SetExtent(data->GetExtent());
02103 data8->SetSpacing(data->GetSpacing());
02104
02105
02106 double mincst, mincstb;
02107
02108 bool flag=true;
02109
02110 double costoactual;
02111
02112 data->GetExtent(ext);
02113
02114 radio=(ext[1]-ext[0])/2;
02115
02116 for(i=ext[0];i<=ext[1];i++){
02117 for(j=ext[2];j<=ext[3];j++){
02118 for(k=ext[4];k<=ext[5];k++){
02119 ptr2=(double *)data2->GetScalarPointer(i,j,k);
02120 ptr4=(unsigned char *)data4->GetScalarPointer(i,j,k);
02121 ptr6=(unsigned char *)data6->GetScalarPointer(i,j,k);
02122 ptr8=(double *)data8->GetScalarPointer(i,j,k);
02123
02124 *ptr2 = 1000000000;
02125 *ptr4 = 0;
02126 *ptr6 = 0;
02127 *ptr8 = 1000000000;
02128 }
02129 }
02130 }
02131
02132
02133
02134 punto[0]=p1[0];
02135 punto[1]=p1[1];
02136 punto[2]=p1[2];
02137 puntob[0]=p2[0];
02138 puntob[1]=p2[1];
02139 puntob[2]=p2[2];
02140
02141
02142 ptr2=(double *)data2->GetScalarPointer(p1[0],p1[1],p1[2]);
02143 *ptr2=0;
02144
02145 ptr4=(unsigned char *)data4->GetScalarPointer(p1[0],p1[1],p1[2]);
02146 *ptr4=1;
02147
02148 ptr8=(double *)data8->GetScalarPointer(p2[0],p2[1],p2[2]);
02149 *ptr8=0;
02150
02151 ptr6=(unsigned char *)data6->GetScalarPointer(p2[0],p2[1],p2[2]);
02152 *ptr6=1;
02153
02154
02155
02156
02157
02158
02159 while(flag){
02160
02161 ptr5=(double *)data2->GetScalarPointer(punto[0],punto[1],punto[2]);
02162 ptr7=(double *)data8->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
02163
02164 for(i=-1;i<=1;i++){
02165 for(j=-1;j<=1;j++){
02166 for(k=-1;k<=1;k++){
02167 if(i!=0 || j!=0 || k!=0){
02168 punto2[0]=punto[0]+i;
02169 punto2[1]=punto[1]+j;
02170 punto2[2]=punto[2]+k;
02171
02172 puntob2[0]=puntob[0]+i;
02173 puntob2[1]=puntob[1]+j;
02174 puntob2[2]=puntob[2]+k;
02175
02176 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] ){
02177
02178 ptr=(double *)data->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
02179 ptr2=(double *)data2->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
02180 ptr3=(unsigned char *)data3->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
02181 ptr4=(unsigned char *)data4->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
02182
02183
02184 if(*ptr3!=0 && *ptr4!=2){
02185 costoactual=*ptr5+(1/(*ptr));
02186 if(*ptr2>costoactual){
02187 *ptr2=costoactual;
02188 *ptr4=1;
02189 }
02190 }
02191 }
02192
02193 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] ){
02194
02195 ptr=(double *)data->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
02196 ptr8=(double *)data8->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
02197 ptr3=(unsigned char *)data3->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
02198 ptr6=(unsigned char *)data6->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
02199
02200 if(*ptr3!=0 && *ptr6!=2){
02201 costoactual=*ptr7+(1/(*ptr));
02202 if(*ptr8>costoactual){
02203 *ptr8=costoactual;
02204 *ptr6=1;
02205 }
02206 }
02207 }
02208 }
02209 }
02210 }
02211 }
02212
02213 ptr4=(unsigned char *)data4->GetScalarPointer(punto[0],punto[1],punto[2]);
02214 *ptr4=2;
02215
02216 ptr6=(unsigned char *)data6->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
02217 *ptr6=2;
02218
02219 mincst=1000000000;
02220 mincstb=1000000000;
02221
02222 for(i=ext[0];i<=ext[1];i++){
02223 for(j=ext[2];j<=ext[3];j++){
02224 for(k=ext[4];k<=ext[5];k++){
02225 ptr2=(double *)data2->GetScalarPointer(i,j,k);
02226 ptr8=(double *)data8->GetScalarPointer(i,j,k);
02227 ptr4=(unsigned char *)data4->GetScalarPointer(i,j,k);
02228 ptr6=(unsigned char *)data6->GetScalarPointer(i,j,k);
02229
02230 if(*ptr4==1){
02231 if(mincst>=*ptr2){
02232 mincst=*ptr2;
02233 puntomin[0]=i;
02234 puntomin[1]=j;
02235 puntomin[2]=k;
02236 }
02237 }
02238
02239 if(*ptr6==1){
02240 if(mincstb>=*ptr8){
02241 mincstb=*ptr8;
02242 puntobmin[0]=i;
02243 puntobmin[1]=j;
02244 puntobmin[2]=k;
02245 }
02246 }
02247 }
02248 }
02249 }
02250
02251 punto[0]=puntomin[0];
02252 punto[1]=puntomin[1];
02253 punto[2]=puntomin[2];
02254
02255 puntob[0]=puntobmin[0];
02256 puntob[1]=puntobmin[1];
02257 puntob[2]=puntobmin[2];
02258
02259 ptr4=(unsigned char *)data4->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
02260 ptr6=(unsigned char *)data6->GetScalarPointer(punto[0],punto[1],punto[2]);
02261
02262
02263 if(*ptr4==2){
02264 p3[0]=puntob[0];
02265 p3[1]=puntob[1];
02266 p3[2]=puntob[2];
02267 flag=false;
02268 }
02269
02270 if(*ptr6==2){
02271 p3[0]=punto[0];
02272 p3[1]=punto[1];
02273 p3[2]=punto[2];
02274 flag=false;
02275 }
02276 }
02277
02278 data2->Delete();
02279 data4->Delete();
02280 data6->Delete();
02281 data8->Delete();
02282 }
02283
02284
02285
02286
02287
02288
02289 void axisExtractor02::invertir(vtkImageData *data )
02290 {
02291 int ext[6];
02292 int i, j, k, radio;
02293 int i2, j2, k2;
02294
02295 data->GetExtent(ext);
02296
02297 unsigned char *ptr;
02298
02299 radio=(ext[1]-ext[0])/2;
02300
02301 double tmpsqrt;
02302
02303 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
02304 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
02305 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
02306 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
02307 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
02308 if( radio<sqrt(tmpsqrt) ){
02309 if(*ptr==0){
02310 *ptr=255;
02311 }
02312 else{
02313 *ptr=0;
02314 }
02315 }
02316 }
02317 }
02318 }
02319 }
02320
02321
02322
02323
02324 void axisExtractor02::redondear(vtkImageData *data )
02325 {
02326 int ext[6];
02327 int i, j, k, radio;
02328 int i2, j2, k2;
02329
02330 data->GetExtent(ext);
02331
02332 unsigned char *ptr;
02333
02334 radio=(ext[1]-ext[0])/2;
02335
02336 double tmpsqrt;
02337
02338 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
02339 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
02340 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
02341 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
02342 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
02343 if( radio<sqrt(tmpsqrt) ){
02344 *ptr=0;
02345 }
02346
02347
02348
02349 }
02350 }
02351 }
02352 }
02353
02354
02355
02356
02357
02358
02359 void axisExtractor02::redondear2(vtkImageData *data )
02360 {
02361 int ext[6];
02362 int i, j, k, radio;
02363 int i2, j2, k2;
02364
02365 data->GetExtent(ext);
02366
02367 double *ptr;
02368
02369 radio=(ext[1]-ext[0])/2;
02370
02371 double tmpsqrt;
02372
02373 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
02374 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
02375 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
02376 ptr=(double *) data->GetScalarPointer(i,j,k);
02377 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
02378 if( radio<sqrt(tmpsqrt) ){
02379 *ptr=0;
02380
02381 }
02382
02383 }
02384 }
02385 }
02386 }
02387
02388
02389
02390
02391
02392 void axisExtractor02::redondear3(vtkImageData *data )
02393 {
02394 int ext[6];
02395 int i, j, k, radio;
02396 int i2, j2, k2;
02397
02398 data->GetExtent(ext);
02399
02400 unsigned char *ptr;
02401
02402 radio=(ext[1]-ext[0])/2;
02403
02404 double tmpsqrt;
02405
02406 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
02407 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
02408 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
02409 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
02410 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
02411 if( radio<sqrt(tmpsqrt) ){
02412 *ptr=255;
02413
02414 }
02415
02416 }
02417 }
02418 }
02419 }
02420
02421
02422
02423
02424
02425 double axisExtractor02::distancia(double a[3], double b[3] )
02426 {
02427 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]));
02428 return sqrt(tmpsqrt);
02429 }
02430
02431
02432
02433
02434
02435
02436 double axisExtractor02::distancia(int a[3], int b[3] )
02437 {
02438 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]));
02439 return sqrt(tmpsqrt);
02440 }
02441
02442
02443
02444
02445
02446
02447
02448
02449 void axisExtractor02::blanquear(vtkImageData *data )
02450 {
02451 int ext[6];
02452 int i, j, k;
02453
02454 data->GetExtent(ext);
02455
02456 unsigned char *ptr;
02457
02458 for(i=ext[0];i<=ext[1];i++){
02459 for(j=ext[2];j<=ext[3];j++){
02460 for(k=ext[4];k<=ext[5];k++){
02461 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
02462 *ptr=0;
02463 }
02464 }
02465 }
02466 }
02467
02468
02469
02470
02471 void axisExtractor02::blanquear3(vtkImageData *data )
02472 {
02473 int ext[6];
02474 int i, j, k;
02475
02476 data->GetExtent(ext);
02477
02478 double *ptr;
02479
02480 for(i=ext[0];i<=ext[1];i++){
02481 for(j=ext[2];j<=ext[3];j++){
02482 for(k=ext[4];k<=ext[5];k++){
02483 ptr=(double *) data->GetScalarPointer(i,j,k);
02484 *ptr=0;
02485 }
02486 }
02487 }
02488 }
02489
02490
02491
02492
02493
02494 void axisExtractor02::blanquear2(vtkImageData *data )
02495 {
02496 int ext[6];
02497 int centro[3];
02498 int i, j, k;
02499
02500 data->GetExtent(ext);
02501
02502 centro[0]=(ext[1]+ext[0])/2;
02503 centro[1]=(ext[3]+ext[2])/2;
02504 centro[2]=(ext[5]+ext[4])/2;
02505
02506 unsigned char *ptr;
02507
02508 for(i=ext[0];i<=ext[1];i++){
02509 for(j=ext[2];j<=ext[3];j++){
02510 for(k=ext[4];k<=ext[5];k++){
02511 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
02512 if(i==centro[0] && j==centro[1] && k==centro[2]){
02513 *ptr=255;
02514 }
02515 else{
02516 *ptr=0;
02517 }
02518 }
02519 }
02520 }
02521 }
02522
02523
02524
02525
02526
02527 void axisExtractor02::cilindro(vtkImageData *data, double vector[3] )
02528 {
02529 int ext[6];
02530 int i, j, k, radio;
02531 int i2, j2, k2;
02532 int centro[3];
02533 int punto[3];
02534 double centrof[3];
02535 double puntof[3];
02536 double radiof;
02537 double puntof2[3];
02538
02539
02540 data->GetExtent(ext);
02541
02542 unsigned char *ptr;
02543
02544 radio=(ext[1]-ext[0])/2;
02545 radiof=espprin[0]*(((double)ext[1]-(double)ext[0])/4);
02546 centro[0]=(ext[1]+ext[0])/2;
02547 centro[1]=(ext[3]+ext[2])/2;
02548 centro[2]=(ext[5]+ext[4])/2;
02549
02550 indextoreal(centro, centrof );
02551
02552 if(vector[0]==0 && vector[1]==0 && vector[2]==0){
02553 blanquear2(data);
02554 }
02555 else{
02556 puntof2[0]=centrof[0]+vector[0];
02557 puntof2[1]=centrof[1]+vector[1];
02558 puntof2[2]=centrof[2]+vector[2];
02559
02560 double tmpsqrt;
02561 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
02562 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
02563 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
02564 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
02565 punto[0]=i;
02566 punto[1]=j;
02567 punto[2]=k;
02568 indextoreal(punto, puntof );
02569 tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2);
02570 if( radio<sqrt(tmpsqrt) ){
02571 *ptr=0;
02572
02573 }
02574 else if(radiof<distanciaejepunto(puntof, centrof, puntof2)){
02575 *ptr=0;
02576 }
02577 else{
02578 *ptr=255;
02579 }
02580 }
02581 }
02582 }
02583 }
02584 }
02585
02586
02587
02588
02589
02590
02591 void axisExtractor02::modelo(vtkImageData *data, unsigned char cantidad, unsigned long vector[50][4], int candit[10][3], double radioactual, double minis[10])
02592 {
02593 int ext[6];
02594 int i, j, k, radio, radio2;
02595 int i2, j2, k2;
02596 int t;
02597 int centro[3];
02598 int punto[3];
02599 int punto2[3];
02600 double centrof[3];
02601 double puntof[3];
02602 double radiof;
02603 double puntof2[3];
02604
02605 double dist, prop, rest;
02606
02607
02608 data->GetExtent(ext);
02609
02610 unsigned char *ptr;
02611
02612 radio=(ext[1]-ext[0])/2;
02613 radio2=(int)radioactual;
02614
02615 centro[0]=(ext[1]+ext[0])/2;
02616 centro[1]=(ext[3]+ext[2])/2;
02617 centro[2]=(ext[5]+ext[4])/2;
02618
02619 indextoreal(centro, centrof );
02620
02621 double tmpsqrt;
02622
02623 for(i=ext[0];i<=ext[1];i++){
02624 for(j=ext[2];j<=ext[3];j++){
02625 for(k=ext[4];k<=ext[5];k++){
02626 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
02627 punto[0]=i;
02628 punto[1]=j;
02629 punto[2]=k;
02630
02631 if( radioactual<distancia(punto, centro) ){
02632 *ptr=0;
02633 }
02634 else{
02635 *ptr=255;
02636
02637 }
02638 }
02639 }
02640 }
02641
02642 for(t=0;t<cantidad;t++){
02643 radiof=sqrt(minis[t]);
02644
02645 punto2[0]=candit[t][0];
02646 punto2[1]=candit[t][1];
02647 punto2[2]=candit[t][2];
02648
02649 indextoreal(punto2, puntof2 );
02650
02651
02652 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
02653 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
02654 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
02655 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
02656 punto[0]=i;
02657 punto[1]=j;
02658 punto[2]=k;
02659 indextoreal(punto, puntof );
02660
02661 dist=distanciaejepunto(puntof, centrof, puntof2);
02662 prop=proporcioejepunto(puntof, centrof, puntof2);
02663
02664 if(prop>=0 && prop<=1){
02665 rest=(radiof-radioactual)*prop+radioactual;
02666 rest=rest+0.25;
02667 rest=rest*espprin[0];
02668 }
02669 else{
02670 rest=dist-1;
02671 }
02672 tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2);
02673 if( radio>sqrt(tmpsqrt) && rest>dist){
02674 *ptr=255;
02675 }
02676 }
02677 }
02678 }
02679 }
02680 }
02681
02682
02683
02684
02685
02686 void axisExtractor02::comparacion(vtkImageData *data, vtkImageData *data2, unsigned long minis[4])
02687 {
02688 int ext[6];
02689 int i, j, k;
02690 double radio;
02691
02692 data->GetExtent(ext);
02693
02694 unsigned char *ptr;
02695 unsigned char *ptr2;
02696
02697
02698 minis[0]=0;
02699 minis[1]=0;
02700 minis[2]=0;
02701 minis[3]=0;
02702
02703 int centro[3];
02704 int punto[3];
02705
02706 radio=((double)ext[1]-(double)ext[0])/2;
02707
02708 centro[0]=(ext[1]+ext[0])/2;
02709 centro[1]=(ext[3]+ext[2])/2;
02710 centro[2]=(ext[5]+ext[4])/2;
02711
02712 for(i=ext[0];i<=ext[1];i++){
02713 for(j=ext[2];j<=ext[3];j++){
02714 for(k=ext[4];k<=ext[5];k++){
02715 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
02716 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
02717
02718 punto[0]=i;
02719 punto[1]=j;
02720 punto[2]=k;
02721
02722 if( radio>distancia(punto, centro) ){
02723
02724 if( *ptr==0 && *ptr2==0 ){
02725 minis[0]++;
02726 }
02727 else if( *ptr!=0 && *ptr2!=0 ){
02728 minis[1]++;
02729 }
02730 else if(*ptr==0 && *ptr2!=0){
02731 minis[2]++;
02732 }
02733 else{
02734 minis[3]++;
02735 }
02736 }
02737 }
02738 }
02739 }
02740 }
02741
02742
02743
02744
02745 void axisExtractor02::copiar(vtkImageData *data, vtkImageData *data2 )
02746 {
02747 int ext[6];
02748 int i, j, k;
02749
02750 data->GetExtent(ext);
02751
02752 unsigned char *ptr, *ptr2;
02753
02754 for(i=ext[0];i<=ext[1];i++){
02755 for(j=ext[2];j<=ext[3];j++){
02756 for(k=ext[4];k<=ext[5];k++){
02757 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
02758 ptr2=(unsigned char *) data2->GetScalarPointer(i,j,k);
02759 if(*ptr!=0 && *ptr2==0){
02760 *ptr2=*ptr;
02761 }
02762 }
02763 }
02764 }
02765
02766 data2->Modified();
02767 }
02768
02769
02770
02771
02772
02773 void axisExtractor02::copiar2(vtkImageData *data, vtkImageData *data2 )
02774 {
02775 int ext[6];
02776 int i, j, k;
02777
02778 data->GetExtent(ext);
02779
02780 double *ptr, *ptr2;
02781
02782 for(i=ext[0];i<=ext[1];i++){
02783 for(j=ext[2];j<=ext[3];j++){
02784 for(k=ext[4];k<=ext[5];k++){
02785 ptr=(double *) data->GetScalarPointer(i,j,k);
02786 ptr2=(double *) data2->GetScalarPointer(i,j,k);
02787 if(*ptr!=0 && *ptr2==0){
02788 *ptr2=*ptr;
02789 }
02790 }
02791 }
02792 }
02793
02794 data2->Modified();
02795 }
02796
02797
02798
02799
02800
02801
02802 double axisExtractor02::angulo(double a[3], double b[3] )
02803 {
02804 double m1=sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2]));
02805 double m2=sqrt((b[0]*b[0])+(b[1]*b[1])+(b[2]*b[2]));
02806 double res=(acos(((a[0]*b[0])+(a[1]*b[1])+(a[2]*b[2]))/(m1*m2))*180)/3.1415;
02807
02808 if(res<0){
02809 res=-res;
02810 }
02811
02812 if(res>=0 && res<90){
02813 res=res;
02814 }
02815 else if(res>=90 && res<180){
02816 res=180-res;
02817 }
02818 else if(res>=180 && res<270){
02819 res=res-180;
02820 }
02821 else{
02822 res=360-res;
02823 }
02824
02825 return res;
02826 }
02827
02828
02829
02830
02831
02832 double axisExtractor02::angulo(double i1, double j1, double k1, double i2, double j2, double k2 )
02833 {
02834 double m1=sqrt((i1*i1)+(j1*j1)+(k1*k1));
02835 double m2=sqrt((i2*i2)+(j2*j2)+(k2*k2));
02836 double res=(acos(((i1*i2)+(j1*j2)+(k1*k2))/(m1*m2))*180)/3.1415;
02837
02838 if(res<0){
02839 res=-res;
02840 }
02841
02842 if(res>=0 && res<90){
02843 res=res;
02844 }
02845 else if(res>=90 && res<180){
02846 res=180-res;
02847 }
02848 else if(res>=180 && res<270){
02849 res=res-180;
02850 }
02851 else{
02852 res=360-res;
02853 }
02854
02855 return res;
02856 }
02857
02858
02859
02860
02861
02862
02863 int axisExtractor02::envolumen(int a[3], vtkImageData *datae )
02864 {
02865
02866 int res;
02867
02868 unsigned short *ptr;
02869
02870 ptr=(unsigned short *) datae->GetScalarPointer(a);
02871
02872 if(*ptr==0){
02873 res=0;
02874 }
02875 else{
02876 res=1;
02877 }
02878
02879 return res;
02880 }
02881
02882
02883
02884
02885
02886 int axisExtractor02::mincandit(int candit[10][3], int cantidad, double puntoanterior[3])
02887 {
02888 double distentpu;
02889 double distmientrp;
02890 int inminenpu=0;
02891
02892 double canditreal[3];
02893 int i;
02894
02895
02896 distmientrp=64000;
02897 inminenpu=0;
02898
02899 for(i=0;i<cantidad && i<10;i++){
02900 indextoreal(candit[i], canditreal);
02901 distentpu=distancia(canditreal, puntoanterior);
02902 if(distentpu<distmientrp){
02903 distmientrp=distentpu;
02904 inminenpu=i;
02905 }
02906 }
02907
02908 return inminenpu;
02909
02910 }
02911
02912
02913
02914
02915 int axisExtractor02::maxareacandit(unsigned long vector[50][4], int cantidad)
02916 {
02917
02918 unsigned long max;
02919 int i;
02920 int indmax=0;
02921
02922
02923 max=0;
02924
02925 for(i=0;i<cantidad && i<10;i++){
02926 if(max<vector[i][0]){
02927 max=vector[i][0];
02928 indmax=i;
02929 }
02930 }
02931
02932 return indmax;
02933
02934 }
02935
02936
02937
02938
02939
02940 unsigned long axisExtractor02::totalarea(unsigned long vector[50][4], unsigned long vectorb[50][4], int cantidad, int cantidadb)
02941 {
02942
02943 unsigned long area;
02944 int i;
02945
02946 area=0;
02947
02948 for(i=0;i<cantidad && i<10;i++){
02949 area+=vector[i][0];
02950 }
02951
02952 for(i=0;i<cantidadb;i++){
02953 area+=vectorb[i][0];
02954 }
02955
02956 return area;
02957 }
02958
02959
02960
02961
02962
02963
02964 unsigned long axisExtractor02::conecarea(unsigned long vector[50][4], int cantidad)
02965 {
02966
02967 unsigned long area;
02968 int i;
02969
02970 area=0;
02971
02972 for(i=0;i<cantidad && i<10;i++){
02973 area+=vector[i][0];
02974 }
02975
02976 return area;
02977 }
02978
02979
02980
02981
02982
02983
02984
02985 int axisExtractor02::bruled(int candit[10][3], int cantidad, vtkImageData *data4)
02986 {
02987 int i;
02988 int burned;
02989
02990 burned=0;
02991 for(i=0;i<cantidad && i<10;i++){
02992 if(envolumen(candit[i], data4)!=0){
02993 burned++;
02994 }
02995 }
02996
02997 return burned;
02998 }
02999
03000
03001
03002
03003
03004 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)
03005 {
03006
03007
03008 int rap;
03009 int punto[3];
03010 double puntof[3];
03011 double puntof2[3];
03012 double centrof[3];
03013 double radiof;
03014 double maxcent;
03015 double mincent;
03016 double distcorr;
03017 double distcorr2;
03018 double distcorr3;
03019 double dist;
03020 double prop;
03021 double rest;
03022 double *ptr;
03023 int ext[6];
03024 int i, j, k, radio;
03025 int centro[3];
03026
03027 data->GetExtent(ext);
03028
03029 radio=(ext[1]-ext[0])/2;
03030 radiof=((double)ext[1]-(double)ext[0])/2;
03031 centro[0]=(ext[1]+ext[0])/2;
03032 centro[1]=(ext[3]+ext[2])/2;
03033 centro[2]=(ext[5]+ext[4])/2;
03034
03035 indicecorregido[0]=centro[0];
03036 indicecorregido[1]=centro[1];
03037 indicecorregido[2]=centro[2];
03038 indextoreal(indicecorregido, puntocorregido);
03039 indextoreal(indiceanterior, centrof);
03040 indextoreal(indicepre, puntof2);
03041
03042 rap= (int)( (double)radio/2 );
03043 maxcent=0;
03044 mincent=64000;
03045
03046 for(i=ext[0];i<=ext[1];i++){
03047 for(j=ext[2];j<=ext[3];j++){
03048 for(k=ext[4];k<=ext[5];k++){
03049 ptr=(double *) data->GetScalarPointer(i,j,k);
03050 punto[0]=i;
03051 punto[1]=j;
03052 punto[2]=k;
03053 indextoreal(punto, puntof );
03054
03055 distcorr3=distancia(indiceanterior, indicepre);
03056 distcorr2=distancia(punto, indicepre);
03057 distcorr=distancia(punto, indiceanterior);
03058
03059 dist=distanciaejepunto(puntof, centrof, puntof2);
03060 prop=proporcioejepunto(puntof, centrof, puntof2);
03061
03062 if(prop>=0 && prop<=1){
03063 rest=(radiopre-radioanterior)*prop+radioanterior;
03064 rest=rest*espprin[0];
03065 }
03066 else{
03067 rest=dist-1;
03068 }
03069
03070 if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && rest>dist){
03071
03072
03073 if( maxcent<*ptr){
03074 mincent=distcorr;
03075 maxcent=*ptr;
03076 indicecorregido[0]=i;
03077 indicecorregido[1]=j;
03078 indicecorregido[2]=k;
03079 indextoreal(indicecorregido, puntocorregido);
03080 }
03081 else if( maxcent==*ptr){
03082 if( mincent>distcorr){
03083 mincent=distcorr;
03084 maxcent=*ptr;
03085 indicecorregido[0]=i;
03086 indicecorregido[1]=j;
03087 indicecorregido[2]=k;
03088 indextoreal(indicecorregido, puntocorregido);
03089 }
03090 }
03091 }
03092 }
03093 }
03094 }
03095
03096 ptr=(double *) data->GetScalarPointer(centro[0],centro[1],centro[2]);
03097
03098 return sqrt(*ptr);
03099 }
03100
03101
03102
03103
03104
03105 double axisExtractor02::correction2(int candit[10][3], int cantidad, vtkImageData *data, int indicecorregido[3], double puntocorregido[3], int indiceanterior[3], double radioanterior)
03106 {
03107
03108
03109 int rap;
03110 int punto[3];
03111 double maxcent;
03112 double mincent;
03113 double distcorr;
03114 double distcorr2;
03115 double *ptr;
03116 int ext[6];
03117 int i, j, k, radio;
03118 int centro[3];
03119
03120 data->GetExtent(ext);
03121
03122 radio=(ext[1]-ext[0])/2;
03123 centro[0]=(ext[1]+ext[0])/2;
03124 centro[1]=(ext[3]+ext[2])/2;
03125 centro[2]=(ext[5]+ext[4])/2;
03126
03127 indicecorregido[0]=centro[0];
03128 indicecorregido[1]=centro[1];
03129 indicecorregido[2]=centro[2];
03130 indextoreal(indicecorregido, puntocorregido);
03131
03132 rap= (int)((double)radio/2);
03133 maxcent=0;
03134 mincent=64000;
03135
03136 for(i=centro[0]-rap;i<=centro[0]+rap;i++){
03137 for(j=centro[1]-rap;j<=centro[1]+rap;j++){
03138 for(k=centro[2]-rap;k<=centro[2]+rap;k++){
03139 ptr=(double *) data->GetScalarPointer(i,j,k);
03140 punto[0]=i;
03141 punto[1]=j;
03142 punto[2]=k;
03143
03144 distcorr2=distancia(punto, centro);
03145 distcorr=distancia(punto, indiceanterior);
03146
03147 if( rap>distcorr2 ){
03148 if( maxcent<*ptr){
03149 mincent=distcorr;
03150 maxcent=*ptr;
03151 indicecorregido[0]=i;
03152 indicecorregido[1]=j;
03153 indicecorregido[2]=k;
03154 indextoreal(indicecorregido, puntocorregido);
03155 }
03156 else if( maxcent==*ptr){
03157 if( mincent>distcorr){
03158 mincent=distcorr;
03159 maxcent=*ptr;
03160 indicecorregido[0]=i;
03161 indicecorregido[1]=j;
03162 indicecorregido[2]=k;
03163 indextoreal(indicecorregido, puntocorregido);
03164 }
03165
03166 }
03167 }
03168
03169 }
03170 }
03171 }
03172
03173 return sqrt(maxcent);
03174 }
03175
03176
03177
03178
03179
03180
03181 void axisExtractor02::avanzar()
03182 {
03183
03184 double puntoactual[3];
03185 double puntocorregido[3];
03186 double puntoanterior[3];
03187 double puntosiguiente[3];
03188
03189 double puntopre[3];
03190
03191 int indiceactual[3];
03192 int indiceanterior[3];
03193 int indicecorregido[3];
03194
03195 int indicepre[3];
03196 double radioactual;
03197 double dirant[3];
03198
03199 unsigned char cantidad;
03200 unsigned char cantidadb;
03201
03202 unsigned long vector[50][4];
03203 unsigned long vectorb[50][4];
03204
03205 int vectortemp[3];
03206 double tempc;
03207
03208 int extint[6];
03209
03210 int extintreal[6];
03211
03212 int indexp;
03213
03214 double provvc[3];
03215
03216 int radiotrabajo;
03217
03218 int flag =0;
03219 int flag2 =0;
03220 int flag3 =0;
03221 int flag4 =0;
03222 int flag5 =0;
03223 int flag6 =0;
03224 int flag7 =0;
03225 int flag8 =0;
03226 int flag9 =0;
03227 int flag10 =0;
03228 int flag11 =0;
03229 int flag12 =0;
03230 int flag13 =0;
03231 int flag14 =0;
03232 int flag15 =0;
03233 int flag16 =0;
03234 int flag17 =0;
03235 int flag18 =0;
03236 int flag19 =0;
03237 int flag20 =0;
03238 int flag21 =0;
03239 int flag22 =0;
03240
03241 int i, j;
03242
03243 double radiotrabajof;
03244
03245 double proprad=2;
03246
03247 double radiopred;
03248 double radioanterior;
03249 double radioprinc;
03250
03251 int burned;
03252
03253 int indant;
03254 double canditreal[3];
03255 int indmaxarea;
03256 unsigned long totalsurf;
03257 unsigned long conecsurf;
03258
03259 indmaxarea=0;
03260 unsigned long caf[4];
03261
03262 double rel1;
03263 double rel2;
03264 double rel3;
03265 double rel4;
03266 double rel5;
03267 double rel6;
03268
03269 if(!m_Stack0.empty()){
03270
03271
03272
03273
03274
03275
03276
03277
03278
03279
03280
03281
03282
03283
03284
03285
03286
03287
03288
03289
03290
03291
03292 indexp=m_Stack.front();
03293
03294 puntoactual[0]=m_Stack0.front();
03295 puntoactual[1]=m_Stack1.front();
03296 puntoactual[2]=m_Stack2.front();
03297
03298 puntoanterior[0]=m_Stack3.front();
03299 puntoanterior[1]=m_Stack4.front();
03300 puntoanterior[2]=m_Stack5.front();
03301
03302 puntopre[0]=m_Stack6.front();
03303 puntopre[1]=m_Stack7.front();
03304 puntopre[2]=m_Stack8.front();
03305
03306 radiopred=m_Stackr.front();
03307 radioanterior=m_Stackra.front();
03308 radioprinc=m_Stackrp.front();
03309
03310 radiotrabajo= (int) (proprad*radiopred);
03311 radiotrabajof=proprad*radiopred;
03312
03313 if(radiotrabajof-radiotrabajo>0.5){
03314 radiotrabajo=radiotrabajo+1;
03315 }
03316
03317
03318
03319 if (minant>radiotrabajo)
03320 {
03321 radiotrabajo=minant;
03322 }
03323
03324
03325
03326
03327
03328
03329
03330
03331
03332
03333
03334
03335
03336
03337
03338
03339 m_Stack0.pop_front();
03340 m_Stack1.pop_front();
03341 m_Stack2.pop_front();
03342 m_Stack3.pop_front();
03343 m_Stack4.pop_front();
03344 m_Stack5.pop_front();
03345 m_Stack6.pop_front();
03346 m_Stack7.pop_front();
03347 m_Stack8.pop_front();
03348
03349 m_Stack.pop_front();
03350 m_Stackr.pop_front();
03351 m_Stackra.pop_front();
03352 m_Stackrp.pop_front();
03353
03354 dirant[0]=puntoanterior[0]-puntoactual[0];
03355 dirant[1]=puntoanterior[1]-puntoactual[1];
03356 dirant[2]=puntoanterior[2]-puntoactual[2];
03357
03358
03359 realtoindex(puntoactual, indiceactual);
03360 realtoindex(puntoanterior, indiceanterior);
03361 realtoindex(puntopre, indicepre);
03362
03363 radioactual=radiopred;
03364
03365 if(radiotrabajo<=maxant){
03366
03367 extint[0]=indiceactual[0]-radiotrabajo;
03368 extint[1]=indiceactual[0]+radiotrabajo;
03369 extint[2]=indiceactual[1]-radiotrabajo;
03370 extint[3]=indiceactual[1]+radiotrabajo;
03371 extint[4]=indiceactual[2]-radiotrabajo;
03372 extint[5]=indiceactual[2]+radiotrabajo;
03373
03374 extrac->SetVOI(extint);
03375 extrac->UpdateWholeExtent();
03376 extrac->Update();
03377 extrac->GetOutput()->GetExtent(extintreal);
03378
03379 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]){
03380 fprintf(stream, "salio\t");
03381 flag4 =1;
03382 }
03383 else{
03384
03385 data1->Delete();
03386 data1=vtkImageData::New();
03387 data1->SetScalarTypeToUnsignedChar();
03388 data1->SetExtent(extrac->GetOutput()->GetExtent());
03389 data1->SetSpacing(extrac->GetOutput()->GetSpacing());
03390
03391 cilindro(data1, dirant );
03392
03393 optim(extrac->GetOutput(), data1 );
03394
03395 connect->SetInput(data1);
03396 connect->RemoveAllSeeds();
03397 connect->AddSeed(indiceactual[0],indiceactual[1],indiceactual[2]);
03398 connect->UpdateWholeExtent();
03399 connect->Update();
03400
03401 data2->Delete();
03402 data2=vtkImageData::New();
03403 data2->SetScalarTypeToUnsignedChar();
03404 data2->SetExtent(connect->GetOutput()->GetExtent());
03405 data2->SetSpacing(connect->GetOutput()->GetSpacing());
03406
03407 blanquear(data2);
03408
03409 cantidad=find_components(connect->GetOutput(), data2, 0, vector);
03410
03411 data3->Delete();
03412 data3=vtkImageData::New();
03413 data3->SetScalarTypeToUnsignedChar();
03414 data3->SetExtent(connect->GetOutput()->GetExtent());
03415 data3->SetSpacing(connect->GetOutput()->GetSpacing());
03416
03417 blanquear(data3);
03418
03419 cantidadb=find_componentsb(connect->GetOutput(), data3, 0, vectorb);
03420
03421 redondear3(connect->GetOutput());
03422
03423 distance->Update();
03424
03425 redondear2(distance->GetOutput());
03426
03427 redondear(connect->GetOutput());
03428
03429 costominimo(distance->GetOutput(), data2);
03430
03431
03432 if(buenos==0){
03433 radioactual=correction2(candit, cantidad, distance->GetOutput(), indicecorregido, puntocorregido, indiceanterior, radioanterior);
03434 }
03435 else{
03436 radioactual=correction(candit, cantidad, distance->GetOutput(), indicecorregido, puntocorregido, indiceanterior, radioanterior, indicepre, radioprinc);
03437 }
03438
03439
03440
03441
03442 indant=mincandit(candit, cantidad, puntoanterior);
03443
03444 indmaxarea=maxareacandit(vector, cantidad);
03445 totalsurf=totalarea(vector, vectorb, cantidad, cantidadb);
03446 conecsurf=conecarea(vector, cantidad);
03447
03448 indextoreal(candit[indant], canditreal);
03449
03450 burned=bruled(candit, cantidad, data4);
03451
03452 data6->Delete();
03453 data6=vtkImageData::New();
03454 data6->SetScalarTypeToUnsignedChar();
03455 data6->SetExtent(extrac->GetOutput()->GetExtent());
03456 data6->SetSpacing(extrac->GetOutput()->GetSpacing());
03457
03458 modelo(data6, cantidad, vector, candit, radioactual, minis);
03459 comparacion(connect->GetOutput(), data6, caf);
03460
03461 rel1=((double)caf[0]+(double)caf[1])/((double)caf[0]+(double)caf[1]+(double)caf[2]+(double)caf[3]);
03462 rel2=((double)caf[1])/((double)caf[1]+(double)caf[3]);
03463 rel3=((double)caf[1])/((double)caf[1]+(double)caf[2]);
03464 rel4=(double)conecsurf/(double)totalsurf;
03465 rel5=(double)vector[indmaxarea][0]/(double)totalsurf;
03466 rel6=(double)cant/((double)cant+(double)cant2);
03467
03468
03469 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)){
03470 flag5=1;
03471 }
03472
03473 if(indicecorregido[0]!=indiceactual[0] || indicecorregido[1]!=indiceactual[1] || indicecorregido[2]!=indiceactual[2]){
03474 flag7=1;
03475 }
03476
03477 if(burned>=cantidad ){
03478 flag15=1;
03479 }
03480
03481 if(burned==0 ){
03482 flag16=1;
03483 }
03484
03485 if(cantidad<2 ){
03486 flag17=1;
03487 }
03488
03489 if(flag7==1 && flagg>4){
03490 flag7=0;
03491 }
03492
03493 if(buenos==0){
03494 flag19=1;
03495 }
03496
03497 if(mejordst!=0){
03498 flag18=1;
03499 }
03500
03501 if(flag7==1 && flag18==1){
03502 for(i=0;i<flagg;i++){
03503 if(visited[i][0]==indiceactual[0] && visited[i][1]==indiceactual[1] && visited[i][2]==indiceactual[2] && visitedrad[i]==radiopred){
03504 flag2=1;
03505 }
03506 }
03507 }
03508
03509 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) )){
03510 flag8=1;
03511 }
03512
03513 if(centi<centi2){
03514 flag3=1;
03515 }
03516
03517 if(angulo(dirant[0], dirant[1], dirant[2], canditreal[0]-puntoactual[0], canditreal[1]-puntoactual[1], canditreal[2]-puntoactual[2] )>65){
03518 flag11=1;
03519 }
03520
03521 if(angulo(Vp[0][ejeminp], Vp[1][ejeminp], Vp[2][ejeminp], Vi[0][ejemini], Vi[1][ejemini], Vi[2][ejemini] )>60){
03522 flag6=1;
03523 }
03524
03525 if((flag3==1 || flag5==1 || (flag16==1 && flag11==1 && flag19==0) ) && (flag7==0 && flag8==0)){
03526 flag9=1;
03527
03528 }
03529
03530 if(flag3==1 || flag5==1 || (flag16==1 && flag11==1 && flag19==0) ){
03531 flag20=1;
03532 }
03533
03534 if((mejordst<radioactual || (mejordst==radioactual && mejorrad>radiopred) || (mejordst==radioactual && mejorrad<=radiopred && mejorcant<cantidad )) && flag17==0 && flag20==0 && flag4==0 && flag12==0 && flag15==0 ){
03535 mejordst=radioactual;
03536 mejor[0]=puntoactual[0];
03537 mejor[1]=puntoactual[1];
03538 mejor[2]=puntoactual[2];
03539 mejorrad=radiopred;
03540 mejorcant=cantidad;
03541 }
03542 }
03543 }
03544 else{
03545 fprintf(stream, "grande\t");
03546 flag12 =1;
03547 }
03548
03549 if(mejordst!=0){
03550 flag18=1;
03551 }
03552
03553 if(flagg>4){
03554 flag13=1;
03555 }
03556
03557 if(flagg2>4){
03558 flag14=1;
03559 }
03560
03561 if( flag9==1 || flag4==1 || flag12==1 || ((cantidad<2 || burned==cantidad) && (flag7==0 && flag8==0)) ){
03562 frama=1;
03563 }
03564 else{
03565 frama=0;
03566 }
03567
03568 if( cantidad>2 && (flag7==0 && flag8==0) ){
03569 fseg=1;
03570 }
03571 else{
03572 fseg=0;
03573 }
03574
03575 if((flag7==0 && flag8==0 && (mejordst>radioactual || (mejordst==radioactual && mejorrad<radiopred)) && flag19==0) || flag2==1){
03576
03577 flagg=10;
03578 flagg2=10;
03579 mejordst=0;
03580
03581 fprintf(stream, "mejor\t");
03582
03583 m_Stack0.push_front(mejor[0]);
03584 m_Stack1.push_front(mejor[1]);
03585 m_Stack2.push_front(mejor[2]);
03586 m_Stack3.push_front(puntoanterior[0]);
03587 m_Stack4.push_front(puntoanterior[1]);
03588 m_Stack5.push_front(puntoanterior[2]);
03589 m_Stack6.push_front(puntopre[0]);
03590 m_Stack7.push_front(puntopre[1]);
03591 m_Stack8.push_front(puntopre[2]);
03592 m_Stack.push_front(indexp);
03593 m_Stackr.push_front(mejorrad);
03594 m_Stackra.push_front(radioanterior);
03595 m_Stackrp.push_front(radioprinc);
03596
03597 mejorrad=0;
03598 mejorcant=0;
03599
03600 }
03601
03602 else if((flag17==1 || flag9==1 || flag4==1 || flag12==1 || flag15==1) && flag18==1){
03603 flagg=10;
03604 flagg2=10;
03605 mejordst=0;
03606
03607 fprintf(stream, "descorrigio\t");
03608
03609 m_Stack0.push_front(mejor[0]);
03610 m_Stack1.push_front(mejor[1]);
03611 m_Stack2.push_front(mejor[2]);
03612 m_Stack3.push_front(puntoanterior[0]);
03613 m_Stack4.push_front(puntoanterior[1]);
03614 m_Stack5.push_front(puntoanterior[2]);
03615 m_Stack6.push_front(puntopre[0]);
03616 m_Stack7.push_front(puntopre[1]);
03617 m_Stack8.push_front(puntopre[2]);
03618 m_Stack.push_front(indexp);
03619 m_Stackr.push_front(mejorrad);
03620 m_Stackra.push_front(radioanterior);
03621 m_Stackrp.push_front(radioprinc);
03622
03623 mejorrad=0;
03624 mejorcant=0;
03625
03626 }
03627
03628 else if(flag12==0 && flag4==0 && flag9==0){
03629 if(flag7==0 && flag8==0){
03630
03631 if(flag17==0 && flag15==0){
03632 flagg=0;
03633 flagg2=0;
03634 mejordst=0;
03635 mejorrad=0;
03636 mejorcant=0;
03637
03638 fprintf(stream, "inserto\t");
03639
03640 if(flag19==0){
03641
03642
03643
03644
03645
03646
03647
03648
03649
03650
03651
03652
03653
03654
03655
03656
03657
03658
03659
03660
03661 realtoreal2(puntoactual, provvc);
03662 points->InsertPoint(buenos,provvc);
03663
03664 lineas->InsertNextCell(2);
03665 lineas->InsertCellPoint(indexp);
03666 lineas->InsertCellPoint(buenos);
03667
03668 buenos++;
03669
03670 }
03671 else{
03672 realtoreal2(puntoactual, provvc);
03673 points->InsertPoint(buenos,provvc);
03674
03675 buenos++;
03676 }
03677
03678
03679 }
03680 else{
03681 flagg=0;
03682 flagg2=0;
03683 mejordst=0;
03684 mejorrad=0;
03685 mejorcant=0;
03686 fprintf(stream, "para\t");
03687 }
03688 }
03689 else if(flag7==0 && flag8==1){
03690 flagg2++;
03691
03692 fprintf(stream, "agrando\t");
03693
03694 m_Stack0.push_front(puntoactual[0]);
03695 m_Stack1.push_front(puntoactual[1]);
03696 m_Stack2.push_front(puntoactual[2]);
03697 m_Stack3.push_front(puntoanterior[0]);
03698 m_Stack4.push_front(puntoanterior[1]);
03699 m_Stack5.push_front(puntoanterior[2]);
03700 m_Stack6.push_front(puntopre[0]);
03701 m_Stack7.push_front(puntopre[1]);
03702 m_Stack8.push_front(puntopre[2]);
03703
03704 m_Stack.push_front(indexp);
03705 m_Stackr.push_front(((double)radiotrabajo+1)/2);
03706 m_Stackra.push_front(radioanterior);
03707 m_Stackrp.push_front(radioprinc);
03708
03709 }
03710 else if(flag7==1 && flag8==1){
03711 visited[flagg][0]=indiceactual[0];
03712 visited[flagg][1]=indiceactual[1];
03713 visited[flagg][2]=indiceactual[2];
03714 visitedrad[flagg]=radiopred;
03715 flagg++;
03716 flagg2++;
03717
03718 fprintf(stream, "agrando y corrigiendo\t");
03719
03720 m_Stack0.push_front(puntocorregido[0]);
03721 m_Stack1.push_front(puntocorregido[1]);
03722 m_Stack2.push_front(puntocorregido[2]);
03723 m_Stack3.push_front(puntoanterior[0]);
03724 m_Stack4.push_front(puntoanterior[1]);
03725 m_Stack5.push_front(puntoanterior[2]);
03726 m_Stack6.push_front(puntopre[0]);
03727 m_Stack7.push_front(puntopre[1]);
03728 m_Stack8.push_front(puntopre[2]);
03729
03730 m_Stack.push_front(indexp);
03731 m_Stackr.push_front(((double)radiotrabajo+1)/2);
03732 m_Stackra.push_front(radioanterior);
03733 m_Stackrp.push_front(radioprinc);
03734 }
03735 else{
03736 visited[flagg][0]=indiceactual[0];
03737 visited[flagg][1]=indiceactual[1];
03738 visited[flagg][2]=indiceactual[2];
03739 visitedrad[flagg]=radiopred;
03740 flagg++;
03741
03742 if(flag19==1){
03743 flagg=10;
03744 }
03745
03746 fprintf(stream, "corrigio\t");
03747
03748 m_Stack0.push_front(puntocorregido[0]);
03749 m_Stack1.push_front(puntocorregido[1]);
03750 m_Stack2.push_front(puntocorregido[2]);
03751 m_Stack3.push_front(puntoanterior[0]);
03752 m_Stack4.push_front(puntoanterior[1]);
03753 m_Stack5.push_front(puntoanterior[2]);
03754
03755 m_Stack6.push_front(puntopre[0]);
03756 m_Stack7.push_front(puntopre[1]);
03757 m_Stack8.push_front(puntopre[2]);
03758
03759 m_Stack.push_front(indexp);
03760 m_Stackr.push_front(radiopred);
03761 m_Stackra.push_front(radioanterior);
03762 m_Stackrp.push_front(radioprinc);
03763
03764
03765 }
03766
03767 if(flag7==0 && flag8==0){
03768
03769 if(flag17==0 && flag15==0){
03770
03771 for(i=0;i<cantidad && i<10;i++){
03772 for(j=i;j<cantidad && i<10;j++){
03773 if(minis[j]>minis[i]){
03774 vectortemp[0]=candit[i][0];
03775 vectortemp[1]=candit[i][1];
03776 vectortemp[2]=candit[i][2];
03777 tempc=minis[i];
03778
03779 candit[i][0]=candit[j][0];
03780 candit[i][1]=candit[j][1];
03781 candit[i][2]=candit[j][2];
03782 minis[i]=minis[j];
03783
03784 candit[j][0]=vectortemp[0];
03785 candit[j][1]=vectortemp[1];
03786 candit[j][2]=vectortemp[2];
03787 minis[j]=tempc;
03788 }
03789 }
03790 }
03791
03792
03793 if(flag16==1){
03794 for(i=0;i<cantidad && i<10;i++){
03795 if(i!=indant || buenos<2){
03796 indextoreal(candit[i], puntosiguiente);
03797
03798 m_Stack0.push_front(puntosiguiente[0]);
03799 m_Stack1.push_front(puntosiguiente[1]);
03800 m_Stack2.push_front(puntosiguiente[2]);
03801 m_Stack3.push_front(puntoactual[0]);
03802 m_Stack4.push_front(puntoactual[1]);
03803 m_Stack5.push_front(puntoactual[2]);
03804
03805 m_Stack6.push_front(puntosiguiente[0]);
03806 m_Stack7.push_front(puntosiguiente[1]);
03807 m_Stack8.push_front(puntosiguiente[2]);
03808
03809 double prov1a=sqrt(minis[i]);
03810
03811 m_Stack.push_front(buenos-1);
03812 m_Stackr.push_front(prov1a);
03813 m_Stackra.push_front(radioactual);
03814 m_Stackrp.push_front(prov1a);
03815
03816 }
03817 }
03818 }
03819 else{
03820 for(i=0;i<cantidad && i<10;i++){
03821 indextoreal(candit[i], puntosiguiente);
03822
03823 if(envolumen(candit[i], data4)==0){
03824 m_Stack0.push_front(puntosiguiente[0]);
03825 m_Stack1.push_front(puntosiguiente[1]);
03826 m_Stack2.push_front(puntosiguiente[2]);
03827 m_Stack3.push_front(puntoactual[0]);
03828 m_Stack4.push_front(puntoactual[1]);
03829 m_Stack5.push_front(puntoactual[2]);
03830
03831 m_Stack6.push_front(puntosiguiente[0]);
03832 m_Stack7.push_front(puntosiguiente[1]);
03833 m_Stack8.push_front(puntosiguiente[2]);
03834
03835 double prov1b=sqrt(minis[i]);
03836
03837 m_Stack.push_front(buenos-1);
03838 m_Stackr.push_front(prov1b);
03839 m_Stackra.push_front(radioactual);
03840 m_Stackrp.push_front(prov1b);
03841 }
03842 }
03843 }
03844
03845 copiar(connect->GetOutput(), data4 );
03846
03847 }
03848 }
03849 }
03850 else{
03851
03852 mejordst=0;
03853 mejorrad=0;
03854 mejorcant=0;
03855 flagg=0;
03856 flagg2=0;
03857
03858 }
03859
03860
03861
03862
03863
03864
03865
03866
03867
03868
03869
03870
03871
03872
03873
03874
03875
03876
03877
03878
03879
03880
03881
03882
03883
03884
03885
03886
03887
03888
03889
03890
03891
03892
03893 fprintf(stream, "%d\t", indexp);
03894 fprintf(stream, "%d\t", buenos);
03895 fprintf(stream, "%d\t", radiotrabajo);
03896 fprintf(stream, "%f\t", radioactual);
03897 fprintf(stream, "%f\t", radiopred);
03898 fprintf(stream, "%f\t", radioanterior);
03899 fprintf(stream, "%f\t", mejordst);
03900 fprintf(stream, "%f\t", mejorrad);
03901 fprintf(stream, "%d\t", mejorcant);
03902
03903
03904
03905
03906
03907
03908
03909
03910 fprintf(stream, "%d\t", flag2);
03911 fprintf(stream, "%d\t", flag3);
03912 fprintf(stream, "%d\t", flag4);
03913 fprintf(stream, "%d\t", flag5);
03914 fprintf(stream, "%d\t", flag6);
03915 fprintf(stream, "%d\t", flag7);
03916 fprintf(stream, "%d\t", flag8);
03917 fprintf(stream, "%d\t", flag9);
03918 fprintf(stream, "%d\t", flag10);
03919 fprintf(stream, "%d\t", flag11);
03920 fprintf(stream, "%d\t", flag12);
03921 fprintf(stream, "%d\t", flag13);
03922 fprintf(stream, "%d\t", flag14);
03923 fprintf(stream, "%d\t", flag15);
03924 fprintf(stream, "%d\t", flag16);
03925 fprintf(stream, "%d\t", flag17);
03926 fprintf(stream, "%d\t", flag18);
03927 fprintf(stream, "%d\t", flag19);
03928 fprintf(stream, "%d\t", flag20);
03929 fprintf(stream, "%d\t", flagg);
03930 fprintf(stream, "%d\t", flagg2);
03931 fprintf(stream, "%d\t", cantidad);
03932 fprintf(stream, "%d\t", cantidadb);
03933 fprintf(stream, "%d\t", burned);
03934
03935
03936
03937
03938
03939
03940
03941
03942
03943
03944
03945
03946
03947
03948
03949
03950
03951
03952
03953 fprintf(stream, "[%f, %f, %f]\t", puntoactual[0], puntoactual[1], puntoactual[2]);
03954 fprintf(stream, "[%f, %f, %f]\t", puntoanterior[0], puntoanterior[1], puntoanterior[2]);
03955 fprintf(stream, "[%f, %f, %f]\t", puntopre[0], puntopre[1], puntopre[2]);
03956 fprintf(stream, "[%f, %f, %f]\t", mejor[0], mejor[1], mejor[2]);
03957
03958 fprintf(stream, "\n");
03959
03960 }
03961
03962 }
03963
03964
03965
03966
03967
03968 void axisExtractor02::todo()
03969 {
03970 int i=0;
03971
03972 while( !m_Stack0.empty() && i<5000){
03973 avanzar();
03974 i++;
03975 }
03976 }
03977
03978
03979
03980
03981 void axisExtractor02::paso()
03982 {
03983
03984 int i=0;
03985 int inicio=buenos;
03986
03987 while( !m_Stack0.empty() && inicio==buenos && i<5000){
03988 avanzar();
03989 i++;
03990 }
03991 }
03992
03993
03994
03995
03996
03997 void axisExtractor02::rama()
03998 {
03999
04000 int i=0;
04001 frama=0;
04002
04003 while( !m_Stack0.empty() && frama==0 && i<5000){
04004 avanzar();
04005 i++;
04006 }
04007 }
04008
04009
04010
04011
04012
04013 void axisExtractor02::segmento()
04014 {
04015
04016 int i=0;
04017 fseg=0;
04018 frama=0;
04019
04020 while( !m_Stack0.empty() && frama==0 && fseg==0 && i<5000){
04021 avanzar();
04022 i++;
04023 }
04024 }
04025
04026
