carotidaBifurcacion.cxx

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003 
00004 =========================================================================*/
00005 #include "carotidaBifurcacion.h"
00006 
00007 #include "vtkPolyData.h"
00008 #include "vtkObjectFactory.h"
00009 #include "vtkImageThreshold.h"
00010 #include "vtkImageCast.h"
00011 #include "vtkImageSeedConnectivity.h"
00012 #include "vtkImageData.h"
00013 #include "vtkMarchingCubes.h"
00014 #include "vtkDoubleArray.h"
00015 #include "vtkPointData.h"
00016 #include "vtkExtractVOI.h"
00017 #include "vtkMath.h"
00018 
00019 
00020 
00021 vtkStandardNewMacro(carotidaBifurcacion);
00022 
00023 
00024 carotidaBifurcacion::carotidaBifurcacion() {
00025         
00026         this->NumberOfRequiredInputs = 2;
00027 
00028         this->humbral=0.45 ;
00029         this->maxpropradio=5;
00030         this->maxpropmasa=10;
00031         this->minpropmasa=0.001;
00032         
00033         this->extrac = vtkExtractVOI::New();
00034         this->extrac->SetSampleRate(1, 1, 1);
00035         
00036         this->thresh = vtkImageThreshold::New();
00037         this->thresh->SetInput(this->extrac->GetOutput());
00038         this->thresh->SetInValue(255);
00039         this->thresh->SetOutValue(0);
00040         this->thresh->ReleaseDataFlagOff();
00041 
00042         this->cast = vtkImageCast::New();
00043         this->cast->SetInput(this->thresh->GetOutput());
00044         this->cast->SetOutputScalarTypeToUnsignedChar();
00045 
00046         this->connect = vtkImageSeedConnectivity::New();
00047         this->connect->SetInput(this->cast->GetOutput());
00048         this->connect->SetInputConnectValue(255);
00049         this->connect->SetOutputConnectedValue(255);
00050         this->connect->SetOutputUnconnectedValue(0);
00051 
00052         this->dataprov=vtkImageData::New();
00053         this->dataprov->SetScalarTypeToUnsignedChar();
00054 
00055         
00056                         
00057 
00058 }
00059 //----------------------------------------------------------------------------
00060 // Specify the input data or filter.
00061 void carotidaBifurcacion::SetInput(vtkPolyData *input)
00062 {
00063         this->vtkProcessObject::SetNthInput(0, input);
00064 }
00065 
00066 
00067 //----------------------------------------------------------------------------
00068 // Specify the input data or filter.
00069 void carotidaBifurcacion::SetInput2(vtkImageData *input)
00070 {
00071         this->vtkProcessObject::SetNthInput(1, input);
00072         this->extrac->SetInput(this->GetInput2());
00073   
00074 }
00075 
00076 
00077 //----------------------------------------------------------------------------
00078 // Specify the input data or filter.
00079 
00080 void carotidaBifurcacion::Execute()
00081 {
00082         vtkIdType totalpuntos=this->GetInput()->GetNumberOfPoints();
00083         this->points = vtkPoints::New();
00084         this->points2 = vtkPoints::New();
00085         this->points3 = vtkPoints::New();
00086         this->points4 = vtkPoints::New();
00087         this->buenos2=0;
00088         this->iter=0;
00089         /*this->lineas = vtkCellArray::New();
00090         this->salidas= vtkDoubleArray::New();
00091         this->salidas->SetNumberOfComponents(3);*/
00092 
00093         for(this->iter=0;this->iter<totalpuntos;this->iter++){
00094                 this->avanzar(this->iter);
00095         }
00096 
00097         this->limpiar(points, points2, points3, points4 );
00098    
00099         this->GetOutput()->SetPoints (this->points);
00100         /*this->GetOutput()->SetLines(this->lineas);
00101         this->GetOutput()->GetPointData()->SetVectors(this->salidas);*/
00102 }
00103 
00104 
00105 
00106 
00107 
00108 
00109 
00110 //----------------------------------------------------------------------------
00111 // Specify the input data or filter.
00112 vtkPolyData *carotidaBifurcacion::GetInput(){
00113         if (this->NumberOfInputs < 1){
00114                 return NULL;
00115         }
00116   
00117         return (vtkPolyData *)(this->Inputs[0]);
00118 }
00119 
00120 
00121 
00122 
00123 
00124 
00125 
00126 //----------------------------------------------------------------------------
00127 // Specify the input data or filter.
00128 vtkImageData *carotidaBifurcacion::GetInput2()
00129 {
00130         if (this->NumberOfInputs < 2){
00131                 return NULL;
00132     }
00133   
00134         return (vtkImageData *)(this->Inputs[1]);
00135 }
00136 
00137 
00138 
00139 
00140 
00141 
00142 
00143 //----------------------------------------------------------------------------
00144 void carotidaBifurcacion::PrintSelf(ostream& os, vtkIndent indent)
00145 {
00146         this->Superclass::PrintSelf(os,indent);
00147 }
00148 
00149 
00150 
00151 
00152 
00153 
00154 
00155 
00156 //----------------------------------------------------------------------------
00157 void carotidaBifurcacion::SetMaxPropRadio(double value)
00158 {
00159         this->maxpropradio=value;
00160 }
00161 
00162 
00163 
00164 //----------------------------------------------------------------------------
00165 double carotidaBifurcacion::GetMaxPropRadio()
00166 {
00167         return this->maxpropradio;
00168 }
00169 
00170 
00171   
00172 //----------------------------------------------------------------------------
00173 
00174 void carotidaBifurcacion::SetHumbral(double value)
00175 {
00176         this->humbral=value;
00177 }
00178 
00179   
00180   
00181   //----------------------------------------------------------------------------
00182 
00183 double carotidaBifurcacion::GetHumbral()
00184 {
00185         return this->humbral;
00186 }
00187 
00188 
00189 //----------------------------------------------------------------------------
00190 
00191 void carotidaBifurcacion::SetMaxPropMasa(double value)
00192 {
00193         this->maxpropmasa=value;
00194 }
00195 
00196   
00197   
00198   //----------------------------------------------------------------------------
00199 
00200 double carotidaBifurcacion::GetMaxPropMasa()
00201 {
00202         return this->maxpropmasa;
00203 }
00204 
00205 
00206 //----------------------------------------------------------------------------
00207 
00208 void carotidaBifurcacion::SetMinPropMasa(double value)
00209 {
00210         this->minpropmasa=value;
00211 }
00212 
00213   
00214   
00215   //----------------------------------------------------------------------------
00216 
00217 double carotidaBifurcacion::GetMinPropMasa()
00218 {
00219         return this->minpropmasa;
00220 }
00221 
00222 
00223 
00224 
00225 
00226 
00227 
00228 void carotidaBifurcacion::searc(int i, int j, int k)
00229 {
00230 
00231         unsigned char *ptr;
00232         unsigned char *ptr2;
00233         
00234         int radio;
00235 
00236         int i2, j2, k2;
00237         int i3, j3, k3;
00238                 
00239         ptr=(unsigned char *)   this->connect->GetOutput()->GetScalarPointer(i,j,k);
00240         ptr2=(unsigned char *)  this->dataprov->GetScalarPointer(i,j,k);
00241 
00242         int ext[6];
00243         
00244         this->connect->GetOutput()->GetExtent(ext);
00245 
00246         radio=(ext[1]-ext[0])/2;
00247 
00248         i2=i-ext[0]-radio;
00249         j2=j-ext[2]-radio;
00250         k2=k-ext[4]-radio;
00251 
00252 
00253         *ptr2=this->label;
00254         this->vector[this->label-1][0]+=1;
00255         this->vector[this->label-1][1]+=i;
00256         this->vector[this->label-1][2]+=j;
00257         this->vector[this->label-1][3]+=k;
00258 
00259 
00260         for(i3=-1;i3<=1;i3++){
00261                         for(j3=-1;j3<=1;j3++){
00262                                 for(k3=-1;k3<=1;k3++){
00263                                         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] ){
00264                                                 ptr=(unsigned char *)   this->connect->GetOutput()->GetScalarPointer(i+i3,j+j3,k+k3);
00265                                                 ptr2=(unsigned char *)  this->dataprov->GetScalarPointer(i+i3,j+j3,k+k3);
00266                                                 if(*ptr==255 && *ptr2==0 && (radio*radio)>((i2+i3)*(i2+i3))+((j2+j3)*(j2+j3))+((k2+k3)*(k2+k3)) && ((radio-2)*(radio-2))<((i2+i3)*(i2+i3))+((j2+j3)*(j2+j3))+((k2+k3)*(k2+k3))){
00267                                                         this->searc(i+i3, j+j3, k+k3);
00268                                                 }
00269                                         }
00270                                 }
00271                         }
00272                 }
00273 }
00274 
00275 
00276 
00277 void carotidaBifurcacion::find_components( )
00278 {
00279         int ext[6];
00280         int i, j, k, radio;
00281         int i2, j2, k2;
00282 
00283         
00284         unsigned char *ptr;
00285         unsigned char *ptr2;
00286 
00287         this->connect->GetOutput()->GetExtent(ext);
00288 
00289 
00290         radio=(ext[1]-ext[0])/2;
00291 
00292         this->label=0;
00293                 
00294 
00295 
00296         for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
00297                         for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
00298                                 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
00299                                         ptr=(unsigned char *)   this->connect->GetOutput()->GetScalarPointer(i,j,k);
00300                                         ptr2=(unsigned char *)  this->dataprov->GetScalarPointer(i,j,k);
00301 
00302                                         if(*ptr==255 && *ptr2==0 && (radio*radio)>(i2*i2)+(j2*j2)+(k2*k2) && ((radio-2)*(radio-2))<(i2*i2)+(j2*j2)+(k2*k2)  ){
00303                                                 this->label=this->label+1;
00304                                                 this->vector[this->label-1][0]=0;
00305                                                 this->vector[this->label-1][1]=0;
00306                                                 this->vector[this->label-1][2]=0;
00307                                                 this->vector[this->label-1][3]=0;
00308                                                 this->searc(i, j, k);
00309                                                 
00310                                         }
00311                                 }
00312                         }
00313                 }
00314         
00315 }
00316 
00317 
00318 
00319 
00320 
00321 unsigned short carotidaBifurcacion::maximo( )
00322 {
00323         int ext[6];
00324         int i, j, k;
00325         unsigned short max=0;
00326 
00327         this->extrac->GetOutput()->GetExtent(ext);
00328 
00329         unsigned short *ptr;
00330         
00331 
00332         for(i=ext[0];i<=ext[1];i++){
00333                         for(j=ext[2];j<=ext[3];j++){
00334                                 for(k=ext[4];k<=ext[5];k++){
00335                                         ptr=(unsigned short *)this->extrac->GetOutput()->GetScalarPointer(i,j,k);
00336                                         if(*ptr>max){
00337                                                 max=*ptr;
00338                                                 
00339                                         }
00340                                 }
00341                         }
00342                 }
00343 
00344                 return max;
00345 }
00346 
00347 
00348 
00349 
00350 
00351 void carotidaBifurcacion::blanquear()
00352 {
00353         int ext[6];
00354         int i, j, k;
00355         
00356         this->dataprov->GetExtent(ext);
00357 
00358         unsigned char *ptr;
00359 
00360 
00361 
00362         for(i=ext[0];i<=ext[1];i++){
00363                         for(j=ext[2];j<=ext[3];j++){
00364                                 for(k=ext[4];k<=ext[5];k++){
00365                                         ptr=(unsigned char *)   this->dataprov->GetScalarPointer(i,j,k);
00366                                         *ptr=0;
00367                                 }
00368                         }
00369                 }
00370 }
00371 
00372 
00373 
00374 
00375 
00376 double carotidaBifurcacion::angulo(double i1, double j1, double k1, double i2, double j2, double k2 )
00377 {
00378         double m1=sqrt((i1*i1)+(j1*j1)+(k1*k1));
00379         double m2=sqrt((i2*i2)+(j2*j2)+(k2*k2));
00380         double res=(acos(((i1*i2)+(j1*j2)+(k1*k2))/(m1*m2))*180)/3.1415;
00381 
00382         if(res<0){
00383                 res=-res;
00384         }
00385 
00386         if(res>=0 && res<180){
00387                 res=res;
00388         }
00389         else{
00390                 res=360-res;
00391         }
00392         
00393         return res;
00394 }
00395 
00396 
00397 
00398 void carotidaBifurcacion::corte(double punto1[3], double punto2[3], double punto3[3], double centro[3],  double radio )
00399 {
00400         
00401         double m1, m2, m3;
00402         double b1, b2, b3;
00403         double c0, c1, c2;
00404         double* roots;
00405         double root;
00406         
00407 
00408         m1=punto2[0]-punto1[0];
00409         m2=punto2[1]-punto1[1];
00410         m3=punto2[2]-punto1[2];
00411 
00412         b1=punto1[0]-centro[0];
00413         b2=punto1[1]-centro[1];
00414         b3=punto1[2]-centro[2];
00415 
00416         c0=m1*m1+m2*m2+m3*m3;
00417         c1=2*m1*b1+2*m2*b2+2*m3*b3;
00418         c2=b1*b1+b2*b2+b3*b3-radio*radio;
00419 
00420         roots=vtkMath::SolveQuadratic  (c0, c1, c2);
00421 
00422         if(roots[1]>=0 && roots[1]<=1){
00423                 root=roots[1];
00424         }
00425         else{
00426                 root=roots[2];
00427         }
00428 
00429         punto3[0]=punto1[0]+root*(punto2[0]-punto1[0]);
00430         punto3[1]=punto1[1]+root*(punto2[1]-punto1[1]);
00431         punto3[2]=punto1[2]+root*(punto2[2]-punto1[2]);
00432                 
00433 }
00434 
00435 
00436 
00437 
00438 
00439 void carotidaBifurcacion::direcciones(vtkPolyData *profile, int iter, double radio, double puntocortea[3], double puntocortes[3] )
00440 {
00441         double puntoactual[3];
00442         double puntoant[3];
00443         double puntosig[3];
00444         double puntoant2[3];
00445         double puntosig2[3];
00446         double dist;
00447         int i,j;
00448         
00449         vtkIdType totalpuntos=profile->GetNumberOfPoints();
00450 
00451         if(iter<totalpuntos){
00452 
00453                 profile->GetPoint(iter, puntoactual);
00454         
00455                 if(iter==0){
00456                         profile->GetPoint(iter+1, puntosig);
00457                         profile->GetPoint(iter, puntosig2);
00458                         i=1;
00459 
00460                         dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
00461 
00462                         while(dist< radio){
00463                         
00464                                 profile->GetPoint(iter+1+i, puntosig);
00465                                 profile->GetPoint(iter+i, puntosig2);
00466                                 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
00467                                 i++;
00468 
00469                         }
00470 
00471                         this->corte(puntosig, puntosig2, puntocortes, puntoactual,  radio );
00472 
00473                         puntocortea[0]=2*puntoactual[0]-puntocortes[0];
00474                         puntocortea[1]=2*puntoactual[1]-puntocortes[1];
00475                         puntocortea[2]=2*puntoactual[2]-puntocortes[2];
00476                         
00477                 }
00478                 else if(iter==totalpuntos-1){
00479                         profile->GetPoint(iter-1, puntoant);
00480                         profile->GetPoint(iter, puntoant2);
00481                         j=1;
00482 
00483                         dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
00484 
00485                         while(dist< radio){
00486                         
00487                                 profile->GetPoint(iter-1-j, puntoant);
00488                                 profile->GetPoint(iter-j, puntoant2);
00489                                 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
00490                                 j++;
00491 
00492                         }
00493 
00494                         this->corte(puntoant, puntoant2, puntocortea, puntoactual,  radio );
00495 
00496                         puntocortes[0]=2*puntoactual[0]-puntocortea[0];
00497                         puntocortes[1]=2*puntoactual[1]-puntocortea[1];
00498                         puntocortes[2]=2*puntoactual[2]-puntocortea[2];
00499                 }
00500                 else{
00501                         profile->GetPoint(iter+1, puntosig);
00502                         profile->GetPoint(iter, puntosig2);
00503                         i=1;
00504 
00505                         dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
00506 
00507                         while(dist< radio && iter+1+i<totalpuntos){
00508                         
00509                                 profile->GetPoint(iter+1+i, puntosig);
00510                                 profile->GetPoint(iter+i, puntosig2);
00511                                 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
00512                                 i++;
00513 
00514                         }
00515 
00516                         profile->GetPoint(iter-1, puntoant);
00517                         profile->GetPoint(iter, puntoant2);
00518                         j=1;
00519 
00520                         dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
00521 
00522                         while(dist< radio && iter-1-j>=0){
00523                         
00524                                 profile->GetPoint(iter-1-j, puntoant);
00525                                 profile->GetPoint(iter-j, puntoant2);
00526                                 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
00527                                 j++;
00528                         }
00529 
00530 
00531                         if(iter+1+i<totalpuntos){
00532                                 this->corte(puntosig, puntosig2, puntocortes, puntoactual,  radio );                            
00533                         }
00534 
00535                         if(iter-1-j>=0){
00536                                 this->corte(puntoant, puntoant2, puntocortea, puntoactual,  radio );    
00537                         }
00538 
00539                         if(iter+1+i>=totalpuntos){
00540                                 puntocortes[0]=2*puntoactual[0]-puntocortea[0];
00541                                 puntocortes[1]=2*puntoactual[1]-puntocortea[1];
00542                                 puntocortes[2]=2*puntoactual[2]-puntocortea[2];
00543                         }
00544                         
00545                         if(iter-1-j<0){
00546                                 puntocortea[0]=2*puntoactual[0]-puntocortes[0];
00547                                 puntocortea[1]=2*puntoactual[1]-puntocortes[1];
00548                                 puntocortea[2]=2*puntoactual[2]-puntocortes[2];
00549                         }
00550                 }
00551         }
00552         
00553 }
00554 
00555 int carotidaBifurcacion::igual(double puntoactual[3], double puntoactualdis[3], double puntoantguar[3], double puntoantguardis[3] )
00556 {
00557         int res=0;
00558         double mat[3][3];
00559         double matinv[3][3];
00560                                 
00561         double t[3];
00562         double b[3];
00563 
00564         double disact;
00565         double disant;
00566         double discru;
00567 
00568                                                                 
00569 
00570         disact=sqrt(vtkMath::Distance2BetweenPoints (puntoactual,puntoactualdis));
00571         disant=sqrt(vtkMath::Distance2BetweenPoints (puntoantguar,puntoantguardis));
00572         discru=sqrt(vtkMath::Distance2BetweenPoints (puntoactual,puntoantguar));
00573 
00574                                                         
00575         if(discru<disact+disant){
00576                 mat[0][0]=puntoantguardis[0]-puntoantguar[0];
00577                 mat[0][1]=-(puntoactualdis[0]-puntoactual[0]);
00578                 mat[0][2]=0;
00579                                                                                 
00580                 mat[1][0]=puntoantguardis[1]-puntoantguar[1];
00581                 mat[1][1]=-(puntoactualdis[1]-puntoactual[1]);
00582                 mat[1][2]=0;
00583 
00584                 mat[2][0]=0;
00585                 mat[2][1]=0;
00586                 mat[2][2]=1;
00587 
00588                 b[0]=puntoactual[0]-puntoantguar[0];
00589                 b[1]=puntoactual[1]-puntoantguar[1];
00590                 b[2]=0;
00591 
00592                                                                         
00593                 vtkMath::Invert3x3(mat,matinv);
00594                 vtkMath::Multiply3x3 (matinv, b, t);
00595 
00596 
00597                                                                 
00598                 if(t[0]>0 && t[1]>0 ){
00599                         res=1;
00600                 }
00601 
00602                                                                         
00603         }
00604 
00605         return res;
00606 }
00607 
00608 
00609 
00610 int carotidaBifurcacion::proporcion(vtkImageData *data )
00611 {
00612         int ext[6];
00613         int i, j, k;
00614         int max=0;
00615         int total=0;
00616 
00617         data->GetExtent(ext);
00618 
00619         unsigned short *ptr;
00620         
00621         for(i=ext[0];i<=ext[1];i++){
00622                         for(j=ext[2];j<=ext[3];j++){
00623                                 for(k=ext[4];k<=ext[5];k++){
00624                                         ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
00625                                         if(*ptr!=0){
00626                                                 max++;
00627                                                 
00628                                         }
00629                                         total++;
00630                                 }
00631                         }
00632                 }
00633 
00634                 return (max*100)/total;
00635 }
00636 
00637 
00638 
00639 
00640 
00641 void carotidaBifurcacion::redondear(vtkImageData *data )
00642 {
00643         int ext[6];
00644         int i, j, k, radio;
00645         int i2, j2, k2;
00646 
00647         data->GetExtent(ext);
00648 
00649         unsigned char *ptr;
00650 
00651         radio=(ext[1]-ext[0])/2;
00652 
00653         double tmpsqrt;
00654         for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
00655                         for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
00656                                 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
00657                                         ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
00658                                     tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2) ;                                 
00659                                         if( radio<sqrt(tmpsqrt) ){
00660                                                 *ptr=0;
00661                                                 
00662                                         }
00663                                 }
00664                         }
00665                 }
00666 }
00667 
00668 
00669 
00670 
00671 void carotidaBifurcacion::limpiar(vtkPoints *profile, vtkPoints *profile2, vtkPoints *profile3, vtkPoints *profile4 )
00672 {
00673         int i,j,k,l,m, n, p;
00674         double puntoactual[3];
00675         double puntoactualdis[3];
00676         double puntoantguar[3];
00677         double puntoantguardis[3];
00678         double promedio[3];
00679         int cant;
00680         int total[100][10];
00681         int pairs[100][2];
00682         int flag;
00683         
00684         vtkIdType totalpuntos=profile2->GetNumberOfPoints();
00685 
00686         k=0;
00687 
00688         for(i=0;i<totalpuntos;i++){
00689                 flag=0;
00690                 profile2->GetPoint(i, puntoactual);
00691                 profile3->GetPoint(i, puntoactualdis);
00692                         
00693                 for(j=i+1;j<totalpuntos;j++){
00694                         profile2->GetPoint(j, puntoantguar);
00695                         profile3->GetPoint(j, puntoantguardis);
00696                         if(puntoactual[0]!=puntoantguar[0] && puntoactual[1]!=puntoantguar[1] && puntoactual[2]!=puntoantguar[2] && 1==this->igual(puntoactual, puntoactualdis, puntoantguar, puntoantguardis)){
00697                                 pairs[k][0]=i;
00698                                 pairs[k][1]=j;
00699                                 k++;
00700                                 flag=1;
00701                         }
00702 
00703                 }
00704                 for(n=0;n<k;n++){
00705                         if(pairs[n][1]==i){
00706                                 flag=1;
00707                         }
00708 
00709                 }
00710                 if(flag==0){
00711                         pairs[k][0]=i;
00712                         pairs[k][1]=-1;
00713                         k++;
00714                 }
00715         }
00716 
00717         for(i=0;i<10;i++){
00718                 for(j=0;j<100;j++){
00719                         total[j][i]=-1;
00720                 }
00721         }
00722 
00723 
00724         l=0;
00725 
00726         for(i=0;i<k;i++){
00727                 m=0;
00728                 if(pairs[i][0]!=-1){
00729                         total[l][m]=pairs[i][0];
00730                         m++;
00731                 }
00732                 if(pairs[i][1]!=-1){
00733                         total[l][m]=pairs[i][1];
00734                         m++;
00735                 }
00736                 if(pairs[i][0]!=-1 || pairs[i][1]!=-1){
00737                         pairs[i][0]=-1;
00738                         pairs[i][1]=-1;
00739                         l++;
00740                         
00741                 
00742                         for(j=i+1;j<k;j++){
00743 
00744                                 for(n=0;n<m;n++){
00745                                 
00746                                         if(total[l-1][n]==pairs[j][0]){
00747                                                 flag=0;
00748                                                 for(p=0;p<m;p++){
00749                                                         if(total[l-1][p]==pairs[j][1]){
00750                                                                 flag=1;
00751                                                         }
00752                                                 }
00753                                                 if(flag==0){
00754                                                         total[l-1][m]=pairs[j][1];
00755                                                         m++;
00756                                                 }
00757                                                 pairs[j][0]=-1;
00758                                                 pairs[j][1]=-1;
00759                                                         
00760 
00761                                         }
00762                                         else if(total[l-1][n]==pairs[j][1]){
00763                                                 flag=0;
00764                                                 for(p=0;p<m;p++){
00765                                                         if(total[l-1][p]==pairs[j][0]){
00766                                                                 flag=1;
00767                                                         }
00768                                                 }
00769                                                 if(flag==0){
00770                                                         total[l-1][m]=pairs[j][0];
00771                                                         m++;
00772                                                 }
00773                                                 pairs[j][0]=-1;
00774                                                 pairs[j][1]=-1;
00775                                         }
00776                                         
00777                                 }
00778                         }
00779                 }
00780         
00781         }
00782 
00783         for(i=0;i<l;i++){
00784                 promedio[0]=0;
00785                 promedio[1]=0;
00786                 promedio[2]=0;
00787                 cant=0;
00788                 flag=0;
00789                 m=0;
00790 
00791                 while(flag==0){
00792                         if(total[i][m]!=-1){
00793                                 profile4->GetPoint(total[i][m], puntoactual);
00794                                 promedio[0]+=puntoactual[0];
00795                                 promedio[1]+=puntoactual[1];
00796                                 promedio[2]+=puntoactual[2];
00797                                 cant++;
00798                         }
00799                         else{
00800                                 flag=1;
00801                         }
00802                         m++;
00803                 }
00804                 promedio[0]=promedio[0]/cant;
00805                 promedio[1]=promedio[1]/cant;
00806                 promedio[2]=promedio[2]/cant;
00807                 
00808                 profile->InsertPoint(i,promedio);
00809 
00810 
00811         }
00812 
00813 
00814 }
00815 
00816 
00817 void carotidaBifurcacion::avanzar(vtkIdType iter)
00818 {
00819 
00820 
00821         double espprin[3];
00822         int extprin[6];
00823         double puntoactual[3];
00824         double puntoactualdis[3];
00825         int puntoactualarr[3];
00826         double radioactual;
00827         double dirant[3];
00828         double dirsig[3];
00829         
00830         
00831         
00832         double vector2[50][3];
00833         double angulos[50];
00834 
00835         unsigned long vectortemp[4];
00836         double vector2temp[3];
00837         double angulostemp;
00838 
00839         double maxmasa;
00840                 
00841         int extint[6];
00842         int extintreal[6];
00843         int prop;
00844 
00845         int i, j;
00846         int radiotrabajo;
00847 
00848         
00849                                 
00850         
00851         this->GetInput2()->GetExtent(extprin);
00852         this->GetInput2()->GetSpacing(espprin);
00853 
00854         vtkIdType totalpuntos=this->GetInput()->GetNumberOfPoints();
00855 
00856         if(iter<totalpuntos){
00857 
00858                 vtkDoubleArray *radio = (vtkDoubleArray *)this->GetInput()->GetPointData()->GetScalars("radio");
00859                 //vtkDoubleArray *dir = (vtkDoubleArray *)this->GetInput()->GetPointData()->GetVectors("norm");
00860                 this->GetInput()->GetPoint(iter, puntoactual);
00861                 radioactual=radio->GetValue(iter);
00862                 radiotrabajo=1;
00863 
00864                 for(this->label=0, prop=100;radiotrabajo<radioactual*this->maxpropradio && (this->label<2 || prop>this->maxpropmasa);radiotrabajo++){
00865 
00866                         puntoactualarr[0]= (int) ( (puntoactual[0]/espprin[0])+extprin[0] );
00867                         puntoactualarr[1]= (int) ( (puntoactual[1]/espprin[1])+extprin[2] );
00868                         puntoactualarr[2]= (int) ( (puntoactual[2]/espprin[2])            );
00869 
00870                         extint[0]=puntoactualarr[0]-radiotrabajo;
00871                         extint[1]=puntoactualarr[0]+radiotrabajo;
00872                         extint[2]=puntoactualarr[1]-radiotrabajo;
00873                         extint[3]=puntoactualarr[1]+radiotrabajo;
00874                         extint[4]=puntoactualarr[2]-radiotrabajo;
00875                         extint[5]=puntoactualarr[2]+radiotrabajo;
00876 
00877                         extrac->SetVOI(extint);
00878                         extrac->UpdateWholeExtent();
00879                         extrac->Update();
00880                         extrac->GetOutput()->GetExtent(extintreal);
00881                         
00882                         if(puntoactualarr[0]>=extintreal[0] && puntoactualarr[0]<=extintreal[1] && puntoactualarr[1]>=extintreal[2] && puntoactualarr[1]<=extintreal[3] && puntoactualarr[2]>=extintreal[4] && puntoactualarr[2]<=extintreal[5]  ){
00883 
00884                                 
00885                                 
00886                                 
00887                         this->thresh->ThresholdByUpper(this->maximo()*this->humbral);
00888                         this->thresh->UpdateWholeExtent();
00889                         this->thresh->Update();
00890 
00891                         prop=this->proporcion(this->thresh->GetOutput());
00892                         this->redondear(this->thresh->GetOutput());
00893 
00894 
00895                         this->cast->UpdateWholeExtent();
00896                         this->cast->Update();
00897 
00898 
00899                         this->connect->RemoveAllSeeds();
00900                         this->connect->AddSeed(puntoactualarr[0],puntoactualarr[1],puntoactualarr[2]);
00901                         this->connect->UpdateWholeExtent();
00902                         this->connect->Update();
00903 
00904                         this->dataprov->Delete();
00905 
00906                         this->dataprov=vtkImageData::New();
00907                         this->dataprov->SetScalarTypeToUnsignedChar();
00908 
00909                         this->dataprov->SetExtent(this->connect->GetOutput()->GetExtent());
00910                 
00911                         this->blanquear();
00912 
00913                         this->find_components();
00914                         
00915                                 
00916 
00917                                 
00918                                 
00919                         }
00920                 }
00921 
00922 
00923                 this->direcciones(this->GetInput(), this->iter, radiotrabajo-1, dirant, dirsig);
00924 
00925                         
00926 
00927                 dirant[0]=dirant[0]-puntoactual[0];
00928                 dirant[1]=dirant[1]-puntoactual[1];
00929                 dirant[2]=dirant[2]-puntoactual[2];
00930                 dirsig[0]=dirsig[0]-puntoactual[0];
00931                 dirsig[1]=dirsig[1]-puntoactual[1];
00932                 dirsig[2]=dirsig[2]-puntoactual[2];
00933 
00934 
00935 
00936 
00937                         
00938                         
00939                         maxmasa=0;
00940                         if(this->label>2){
00941                                 for(i=0;i<this->label;i++){
00942                                         vector2[i][0]=((double)this->vector[i][1]/(double)this->vector[i][0])-(double)puntoactualarr[0];
00943                                         vector2[i][1]=((double)this->vector[i][2]/(double)this->vector[i][0])-(double)puntoactualarr[1];
00944                                         vector2[i][2]=((double)this->vector[i][3]/(double)this->vector[i][0])-(double)puntoactualarr[2];
00945                                         double tmp_EED_A = this->angulo(dirsig[0], dirsig[1], dirsig[2], vector2[i][0], vector2[i][1], vector2[i][2]);
00946                                         double tmp_EED_B = this->angulo(dirant[0], dirant[1], dirant[2], vector2[i][0], vector2[i][1], vector2[i][2]);
00947 //                                      angulos[i] =min( tmp_EED_A , tmp_EED_B );
00948                                         angulos[i] = tmp_EED_A ;
00949                                         if (tmp_EED_B<angulos[i] ) 
00950                                         {
00951                                                 angulos[i]=tmp_EED_B;
00952                                         }
00953                                         if(maxmasa<this->vector[i][0]){
00954                                                 maxmasa=this->vector[i][0];
00955                                         }
00956                                         
00957                                 }
00958 
00959                                 for(i=0;i<this->label;i++){
00960                                         for(j=i;j<this->label;j++){
00961                                                 if(angulos[j]>angulos[i]){
00962                                                         vectortemp[0]=this->vector[i][0];
00963                                                         vectortemp[1]=this->vector[i][1];
00964                                                         vectortemp[2]=this->vector[i][2];
00965                                                         vectortemp[3]=this->vector[i][3];
00966                                                         vector2temp[0]=vector2[i][0];
00967                                                         vector2temp[1]=vector2[i][1];
00968                                                         vector2temp[2]=vector2[i][2];
00969                                                         angulostemp=angulos[i];
00970 
00971                                                         this->vector[i][0]=this->vector[j][0];
00972                                                         this->vector[i][1]=this->vector[j][1];
00973                                                         this->vector[i][2]=this->vector[j][2];
00974                                                         this->vector[i][3]=this->vector[j][3];
00975                                                         vector2[i][0]=vector2[j][0];
00976                                                         vector2[i][1]=vector2[j][1];
00977                                                         vector2[i][2]=vector2[j][2];
00978                                                         angulos[i]=angulos[j];
00979 
00980                                                         this->vector[j][0]=vectortemp[0];
00981                                                         this->vector[j][1]=vectortemp[1];
00982                                                         this->vector[j][2]=vectortemp[2];
00983                                                         this->vector[j][3]=vectortemp[3];
00984                                                         vector2[j][0]=vector2temp[0];
00985                                                         vector2[j][1]=vector2temp[1];
00986                                                         vector2[j][2]=vector2temp[2];
00987                                                         angulos[j]=angulostemp;
00988                                                 }
00989                                         }
00990                                         
00991                                 }
00992 
00993 
00994                                 for(i=this->label-2;i<this->label;i++){
00995                                         this->vector[i][0]=0;
00996                                 }
00997 
00998 
00999 
01000                                 for(i=0;i<this->label;i++){
01001                                         for(j=i;j<this->label;j++){
01002                                                 if(this->vector[j][0]>this->vector[i][0]){
01003                                                         vectortemp[0]=this->vector[i][0];
01004                                                         vectortemp[1]=this->vector[i][1];
01005                                                         vectortemp[2]=this->vector[i][2];
01006                                                         vectortemp[3]=this->vector[i][3];
01007                                                         vector2temp[0]=vector2[i][0];
01008                                                         vector2temp[1]=vector2[i][1];
01009                                                         vector2temp[2]=vector2[i][2];
01010                                                         angulostemp=angulos[i];
01011 
01012                                                         this->vector[i][0]=this->vector[j][0];
01013                                                         this->vector[i][1]=this->vector[j][1];
01014                                                         this->vector[i][2]=this->vector[j][2];
01015                                                         this->vector[i][3]=this->vector[j][3];
01016                                                         vector2[i][0]=vector2[j][0];
01017                                                         vector2[i][1]=vector2[j][1];
01018                                                         vector2[i][2]=vector2[j][2];
01019                                                         angulos[i]=angulos[j];
01020 
01021                                                         this->vector[j][0]=vectortemp[0];
01022                                                         this->vector[j][1]=vectortemp[1];
01023                                                         this->vector[j][2]=vectortemp[2];
01024                                                         this->vector[j][3]=vectortemp[3];
01025                                                         vector2[j][0]=vector2temp[0];
01026                                                         vector2[j][1]=vector2temp[1];
01027                                                         vector2[j][2]=vector2temp[2];
01028                                                         angulos[j]=angulostemp;
01029                                                 }
01030                                         }
01031                                         
01032                                 }
01033 
01034                                 
01035                                 for(i=0;i<this->label;i++){
01036                                         if(maxmasa*this->minpropmasa<this->vector[i][0]){
01037                                                 puntoactualdis[0]=(puntoactualarr[0]+vector2[i][0]-extprin[0])*espprin[0];
01038                                                 puntoactualdis[1]=(puntoactualarr[1]+vector2[i][1]-extprin[2])*espprin[1];
01039                                                 puntoactualdis[2]=(puntoactualarr[2]+vector2[i][2])*espprin[2];
01040 
01041                                                 points2->InsertPoint(buenos2,puntoactual);
01042                                                 points3->InsertPoint(buenos2,puntoactualdis);
01043                                                 points4->InsertPoint(buenos2,puntoactual[0]+(radioactual/radiotrabajo)*(puntoactualdis[0]-puntoactual[0]), puntoactual[1]+(radioactual/radiotrabajo)*(puntoactualdis[1]-puntoactual[1]), puntoactual[2]+(radioactual/radiotrabajo)*(puntoactualdis[2]-puntoactual[2]));
01044                                                 buenos2++;
01045                                         }
01046                                         
01047                                 }
01048                                         
01049                         }
01050         
01051         }
01052 
01053 
01054 }
Generated on 18 Mar 2010 for creaMaracasVisu_lib by  doxygen 1.6.1