1 /*# ---------------------------------------------------------------------
3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
5 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
6 # Previous Authors : Laurent Guigues, Jean-Pierre Roux
7 # CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
9 # This software is governed by the CeCILL-B license under French law and
10 # abiding by the rules of distribution of free software. You can use,
11 # modify and/ or redistribute the software under the terms of the CeCILL-B
12 # license as circulated by CEA, CNRS and INRIA at the following URL
13 # http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
14 # or in the file LICENSE.txt.
16 # As a counterpart to the access to the source code and rights to copy,
17 # modify and redistribute granted by the license, users are provided only
18 # with a limited warranty and the software's author, the holder of the
19 # economic rights, and the successive licensors have only limited
22 # The fact that you are presently reading this means that you have had
23 # knowledge of the CeCILL-B license and that you accept its terms.
24 # ------------------------------------------------------------------------ */
26 /*=========================================================================
29 =========================================================================*/
30 #include "carotidaBifurcacion.h"
32 #include "vtkPolyData.h"
33 #include "vtkObjectFactory.h"
34 #include "vtkImageThreshold.h"
35 #include "vtkImageCast.h"
36 #include "vtkImageSeedConnectivity.h"
37 #include "vtkImageData.h"
38 #include "vtkMarchingCubes.h"
39 #include "vtkDoubleArray.h"
40 #include "vtkPointData.h"
41 #include "vtkExtractVOI.h"
46 vtkStandardNewMacro(carotidaBifurcacion);
49 carotidaBifurcacion::carotidaBifurcacion() {
51 this->NumberOfRequiredInputs = 2;
56 this->minpropmasa=0.001;
58 this->extrac = vtkExtractVOI::New();
59 this->extrac->SetSampleRate(1, 1, 1);
61 this->thresh = vtkImageThreshold::New();
62 this->thresh->SetInput(this->extrac->GetOutput());
63 this->thresh->SetInValue(255);
64 this->thresh->SetOutValue(0);
65 this->thresh->ReleaseDataFlagOff();
67 this->cast = vtkImageCast::New();
68 this->cast->SetInput(this->thresh->GetOutput());
69 this->cast->SetOutputScalarTypeToUnsignedChar();
71 this->connect = vtkImageSeedConnectivity::New();
72 this->connect->SetInput(this->cast->GetOutput());
73 this->connect->SetInputConnectValue(255);
74 this->connect->SetOutputConnectedValue(255);
75 this->connect->SetOutputUnconnectedValue(0);
77 this->dataprov=vtkImageData::New();
78 this->dataprov->SetScalarTypeToUnsignedChar();
84 //----------------------------------------------------------------------------
85 // Specify the input data or filter.
86 void carotidaBifurcacion::SetInput(vtkPolyData *input)
88 this->vtkProcessObject::SetNthInput(0, input);
92 //----------------------------------------------------------------------------
93 // Specify the input data or filter.
94 void carotidaBifurcacion::SetInput2(vtkImageData *input)
96 this->vtkProcessObject::SetNthInput(1, input);
97 this->extrac->SetInput(this->GetInput2());
102 //----------------------------------------------------------------------------
103 // Specify the input data or filter.
105 void carotidaBifurcacion::Execute()
107 vtkIdType totalpuntos=this->GetInput()->GetNumberOfPoints();
108 this->points = vtkPoints::New();
109 this->points2 = vtkPoints::New();
110 this->points3 = vtkPoints::New();
111 this->points4 = vtkPoints::New();
114 /*this->lineas = vtkCellArray::New();
115 this->salidas= vtkDoubleArray::New();
116 this->salidas->SetNumberOfComponents(3);*/
118 for(this->iter=0;this->iter<totalpuntos;this->iter++){
119 this->avanzar(this->iter);
122 this->limpiar(points, points2, points3, points4 );
124 this->GetOutput()->SetPoints (this->points);
125 /*this->GetOutput()->SetLines(this->lineas);
126 this->GetOutput()->GetPointData()->SetVectors(this->salidas);*/
135 //----------------------------------------------------------------------------
136 // Specify the input data or filter.
137 vtkPolyData *carotidaBifurcacion::GetInput(){
138 if (this->NumberOfInputs < 1){
142 return (vtkPolyData *)(this->Inputs[0]);
151 //----------------------------------------------------------------------------
152 // Specify the input data or filter.
153 vtkImageData *carotidaBifurcacion::GetInput2()
155 if (this->NumberOfInputs < 2){
159 return (vtkImageData *)(this->Inputs[1]);
168 //----------------------------------------------------------------------------
169 void carotidaBifurcacion::PrintSelf(ostream& os, vtkIndent indent)
171 this->Superclass::PrintSelf(os,indent);
181 //----------------------------------------------------------------------------
182 void carotidaBifurcacion::SetMaxPropRadio(double value)
184 this->maxpropradio=value;
189 //----------------------------------------------------------------------------
190 double carotidaBifurcacion::GetMaxPropRadio()
192 return this->maxpropradio;
197 //----------------------------------------------------------------------------
199 void carotidaBifurcacion::SetHumbral(double value)
206 //----------------------------------------------------------------------------
208 double carotidaBifurcacion::GetHumbral()
210 return this->humbral;
214 //----------------------------------------------------------------------------
216 void carotidaBifurcacion::SetMaxPropMasa(double value)
218 this->maxpropmasa=value;
223 //----------------------------------------------------------------------------
225 double carotidaBifurcacion::GetMaxPropMasa()
227 return this->maxpropmasa;
231 //----------------------------------------------------------------------------
233 void carotidaBifurcacion::SetMinPropMasa(double value)
235 this->minpropmasa=value;
240 //----------------------------------------------------------------------------
242 double carotidaBifurcacion::GetMinPropMasa()
244 return this->minpropmasa;
253 void carotidaBifurcacion::searc(int i, int j, int k)
264 ptr=(unsigned char *) this->connect->GetOutput()->GetScalarPointer(i,j,k);
265 ptr2=(unsigned char *) this->dataprov->GetScalarPointer(i,j,k);
269 this->connect->GetOutput()->GetExtent(ext);
271 radio=(ext[1]-ext[0])/2;
279 this->vector[this->label-1][0]+=1;
280 this->vector[this->label-1][1]+=i;
281 this->vector[this->label-1][2]+=j;
282 this->vector[this->label-1][3]+=k;
285 for(i3=-1;i3<=1;i3++){
286 for(j3=-1;j3<=1;j3++){
287 for(k3=-1;k3<=1;k3++){
288 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] ){
289 ptr=(unsigned char *) this->connect->GetOutput()->GetScalarPointer(i+i3,j+j3,k+k3);
290 ptr2=(unsigned char *) this->dataprov->GetScalarPointer(i+i3,j+j3,k+k3);
291 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))){
292 this->searc(i+i3, j+j3, k+k3);
302 void carotidaBifurcacion::find_components( )
312 this->connect->GetOutput()->GetExtent(ext);
315 radio=(ext[1]-ext[0])/2;
321 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
322 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
323 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
324 ptr=(unsigned char *) this->connect->GetOutput()->GetScalarPointer(i,j,k);
325 ptr2=(unsigned char *) this->dataprov->GetScalarPointer(i,j,k);
327 if(*ptr==255 && *ptr2==0 && (radio*radio)>(i2*i2)+(j2*j2)+(k2*k2) && ((radio-2)*(radio-2))<(i2*i2)+(j2*j2)+(k2*k2) ){
328 this->label=this->label+1;
329 this->vector[this->label-1][0]=0;
330 this->vector[this->label-1][1]=0;
331 this->vector[this->label-1][2]=0;
332 this->vector[this->label-1][3]=0;
333 this->searc(i, j, k);
346 unsigned short carotidaBifurcacion::maximo( )
350 unsigned short max=0;
352 this->extrac->GetOutput()->GetExtent(ext);
357 for(i=ext[0];i<=ext[1];i++){
358 for(j=ext[2];j<=ext[3];j++){
359 for(k=ext[4];k<=ext[5];k++){
360 ptr=(unsigned short *)this->extrac->GetOutput()->GetScalarPointer(i,j,k);
376 void carotidaBifurcacion::blanquear()
381 this->dataprov->GetExtent(ext);
387 for(i=ext[0];i<=ext[1];i++){
388 for(j=ext[2];j<=ext[3];j++){
389 for(k=ext[4];k<=ext[5];k++){
390 ptr=(unsigned char *) this->dataprov->GetScalarPointer(i,j,k);
401 double carotidaBifurcacion::angulo(double i1, double j1, double k1, double i2, double j2, double k2 )
403 double m1=sqrt((i1*i1)+(j1*j1)+(k1*k1));
404 double m2=sqrt((i2*i2)+(j2*j2)+(k2*k2));
405 double res=(acos(((i1*i2)+(j1*j2)+(k1*k2))/(m1*m2))*180)/3.1415;
411 if(res>=0 && res<180){
423 void carotidaBifurcacion::corte(double punto1[3], double punto2[3], double punto3[3], double centro[3], double radio )
433 m1=punto2[0]-punto1[0];
434 m2=punto2[1]-punto1[1];
435 m3=punto2[2]-punto1[2];
437 b1=punto1[0]-centro[0];
438 b2=punto1[1]-centro[1];
439 b3=punto1[2]-centro[2];
441 c0=m1*m1+m2*m2+m3*m3;
442 c1=2*m1*b1+2*m2*b2+2*m3*b3;
443 c2=b1*b1+b2*b2+b3*b3-radio*radio;
445 roots=vtkMath::SolveQuadratic (c0, c1, c2);
447 if(roots[1]>=0 && roots[1]<=1){
454 punto3[0]=punto1[0]+root*(punto2[0]-punto1[0]);
455 punto3[1]=punto1[1]+root*(punto2[1]-punto1[1]);
456 punto3[2]=punto1[2]+root*(punto2[2]-punto1[2]);
464 void carotidaBifurcacion::direcciones(vtkPolyData *profile, int iter, double radio, double puntocortea[3], double puntocortes[3] )
466 double puntoactual[3];
474 vtkIdType totalpuntos=profile->GetNumberOfPoints();
476 if(iter<totalpuntos){
478 profile->GetPoint(iter, puntoactual);
481 profile->GetPoint(iter+1, puntosig);
482 profile->GetPoint(iter, puntosig2);
485 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
489 profile->GetPoint(iter+1+i, puntosig);
490 profile->GetPoint(iter+i, puntosig2);
491 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
496 this->corte(puntosig, puntosig2, puntocortes, puntoactual, radio );
498 puntocortea[0]=2*puntoactual[0]-puntocortes[0];
499 puntocortea[1]=2*puntoactual[1]-puntocortes[1];
500 puntocortea[2]=2*puntoactual[2]-puntocortes[2];
503 else if(iter==totalpuntos-1){
504 profile->GetPoint(iter-1, puntoant);
505 profile->GetPoint(iter, puntoant2);
508 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
512 profile->GetPoint(iter-1-j, puntoant);
513 profile->GetPoint(iter-j, puntoant2);
514 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
519 this->corte(puntoant, puntoant2, puntocortea, puntoactual, radio );
521 puntocortes[0]=2*puntoactual[0]-puntocortea[0];
522 puntocortes[1]=2*puntoactual[1]-puntocortea[1];
523 puntocortes[2]=2*puntoactual[2]-puntocortea[2];
526 profile->GetPoint(iter+1, puntosig);
527 profile->GetPoint(iter, puntosig2);
530 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
532 while(dist< radio && iter+1+i<totalpuntos){
534 profile->GetPoint(iter+1+i, puntosig);
535 profile->GetPoint(iter+i, puntosig2);
536 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
541 profile->GetPoint(iter-1, puntoant);
542 profile->GetPoint(iter, puntoant2);
545 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
547 while(dist< radio && iter-1-j>=0){
549 profile->GetPoint(iter-1-j, puntoant);
550 profile->GetPoint(iter-j, puntoant2);
551 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
556 if(iter+1+i<totalpuntos){
557 this->corte(puntosig, puntosig2, puntocortes, puntoactual, radio );
561 this->corte(puntoant, puntoant2, puntocortea, puntoactual, radio );
564 if(iter+1+i>=totalpuntos){
565 puntocortes[0]=2*puntoactual[0]-puntocortea[0];
566 puntocortes[1]=2*puntoactual[1]-puntocortea[1];
567 puntocortes[2]=2*puntoactual[2]-puntocortea[2];
571 puntocortea[0]=2*puntoactual[0]-puntocortes[0];
572 puntocortea[1]=2*puntoactual[1]-puntocortes[1];
573 puntocortea[2]=2*puntoactual[2]-puntocortes[2];
580 int carotidaBifurcacion::igual(double puntoactual[3], double puntoactualdis[3], double puntoantguar[3], double puntoantguardis[3] )
595 disact=sqrt(vtkMath::Distance2BetweenPoints (puntoactual,puntoactualdis));
596 disant=sqrt(vtkMath::Distance2BetweenPoints (puntoantguar,puntoantguardis));
597 discru=sqrt(vtkMath::Distance2BetweenPoints (puntoactual,puntoantguar));
600 if(discru<disact+disant){
601 mat[0][0]=puntoantguardis[0]-puntoantguar[0];
602 mat[0][1]=-(puntoactualdis[0]-puntoactual[0]);
605 mat[1][0]=puntoantguardis[1]-puntoantguar[1];
606 mat[1][1]=-(puntoactualdis[1]-puntoactual[1]);
613 b[0]=puntoactual[0]-puntoantguar[0];
614 b[1]=puntoactual[1]-puntoantguar[1];
618 vtkMath::Invert3x3(mat,matinv);
619 vtkMath::Multiply3x3 (matinv, b, t);
623 if(t[0]>0 && t[1]>0 ){
635 int carotidaBifurcacion::proporcion(vtkImageData *data )
642 data->GetExtent(ext);
646 for(i=ext[0];i<=ext[1];i++){
647 for(j=ext[2];j<=ext[3];j++){
648 for(k=ext[4];k<=ext[5];k++){
649 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
659 return (max*100)/total;
666 void carotidaBifurcacion::redondear(vtkImageData *data )
672 data->GetExtent(ext);
676 radio=(ext[1]-ext[0])/2;
679 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
680 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
681 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
682 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
683 tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2) ;
684 if( radio<sqrt(tmpsqrt) ){
696 void carotidaBifurcacion::limpiar(vtkPoints *profile, vtkPoints *profile2, vtkPoints *profile3, vtkPoints *profile4 )
699 double puntoactual[3];
700 double puntoactualdis[3];
701 double puntoantguar[3];
702 double puntoantguardis[3];
709 vtkIdType totalpuntos=profile2->GetNumberOfPoints();
713 for(i=0;i<totalpuntos;i++){
715 profile2->GetPoint(i, puntoactual);
716 profile3->GetPoint(i, puntoactualdis);
718 for(j=i+1;j<totalpuntos;j++){
719 profile2->GetPoint(j, puntoantguar);
720 profile3->GetPoint(j, puntoantguardis);
721 if(puntoactual[0]!=puntoantguar[0] && puntoactual[1]!=puntoantguar[1] && puntoactual[2]!=puntoantguar[2] && 1==this->igual(puntoactual, puntoactualdis, puntoantguar, puntoantguardis)){
754 total[l][m]=pairs[i][0];
758 total[l][m]=pairs[i][1];
761 if(pairs[i][0]!=-1 || pairs[i][1]!=-1){
771 if(total[l-1][n]==pairs[j][0]){
774 if(total[l-1][p]==pairs[j][1]){
779 total[l-1][m]=pairs[j][1];
787 else if(total[l-1][n]==pairs[j][1]){
790 if(total[l-1][p]==pairs[j][0]){
795 total[l-1][m]=pairs[j][0];
818 profile4->GetPoint(total[i][m], puntoactual);
819 promedio[0]+=puntoactual[0];
820 promedio[1]+=puntoactual[1];
821 promedio[2]+=puntoactual[2];
829 promedio[0]=promedio[0]/cant;
830 promedio[1]=promedio[1]/cant;
831 promedio[2]=promedio[2]/cant;
833 profile->InsertPoint(i,promedio);
842 void carotidaBifurcacion::avanzar(vtkIdType iter)
848 double puntoactual[3];
849 double puntoactualdis[3];
850 int puntoactualarr[3];
857 double vector2[50][3];
860 unsigned long vectortemp[4];
861 double vector2temp[3];
876 this->GetInput2()->GetExtent(extprin);
877 this->GetInput2()->GetSpacing(espprin);
879 vtkIdType totalpuntos=this->GetInput()->GetNumberOfPoints();
881 if(iter<totalpuntos){
883 vtkDoubleArray *radio = (vtkDoubleArray *)this->GetInput()->GetPointData()->GetScalars("radio");
884 //vtkDoubleArray *dir = (vtkDoubleArray *)this->GetInput()->GetPointData()->GetVectors("norm");
885 this->GetInput()->GetPoint(iter, puntoactual);
886 radioactual=radio->GetValue(iter);
889 for(this->label=0, prop=100;radiotrabajo<radioactual*this->maxpropradio && (this->label<2 || prop>this->maxpropmasa);radiotrabajo++){
891 puntoactualarr[0]= (int) ( (puntoactual[0]/espprin[0])+extprin[0] );
892 puntoactualarr[1]= (int) ( (puntoactual[1]/espprin[1])+extprin[2] );
893 puntoactualarr[2]= (int) ( (puntoactual[2]/espprin[2]) );
895 extint[0]=puntoactualarr[0]-radiotrabajo;
896 extint[1]=puntoactualarr[0]+radiotrabajo;
897 extint[2]=puntoactualarr[1]-radiotrabajo;
898 extint[3]=puntoactualarr[1]+radiotrabajo;
899 extint[4]=puntoactualarr[2]-radiotrabajo;
900 extint[5]=puntoactualarr[2]+radiotrabajo;
902 extrac->SetVOI(extint);
903 extrac->UpdateWholeExtent();
905 extrac->GetOutput()->GetExtent(extintreal);
907 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] ){
912 this->thresh->ThresholdByUpper(this->maximo()*this->humbral);
913 this->thresh->UpdateWholeExtent();
914 this->thresh->Update();
916 prop=this->proporcion(this->thresh->GetOutput());
917 this->redondear(this->thresh->GetOutput());
920 this->cast->UpdateWholeExtent();
921 this->cast->Update();
924 this->connect->RemoveAllSeeds();
925 this->connect->AddSeed(puntoactualarr[0],puntoactualarr[1],puntoactualarr[2]);
926 this->connect->UpdateWholeExtent();
927 this->connect->Update();
929 this->dataprov->Delete();
931 this->dataprov=vtkImageData::New();
932 this->dataprov->SetScalarTypeToUnsignedChar();
934 this->dataprov->SetExtent(this->connect->GetOutput()->GetExtent());
938 this->find_components();
948 this->direcciones(this->GetInput(), this->iter, radiotrabajo-1, dirant, dirsig);
952 dirant[0]=dirant[0]-puntoactual[0];
953 dirant[1]=dirant[1]-puntoactual[1];
954 dirant[2]=dirant[2]-puntoactual[2];
955 dirsig[0]=dirsig[0]-puntoactual[0];
956 dirsig[1]=dirsig[1]-puntoactual[1];
957 dirsig[2]=dirsig[2]-puntoactual[2];
966 for(i=0;i<this->label;i++){
967 vector2[i][0]=((double)this->vector[i][1]/(double)this->vector[i][0])-(double)puntoactualarr[0];
968 vector2[i][1]=((double)this->vector[i][2]/(double)this->vector[i][0])-(double)puntoactualarr[1];
969 vector2[i][2]=((double)this->vector[i][3]/(double)this->vector[i][0])-(double)puntoactualarr[2];
970 double tmp_EED_A = this->angulo(dirsig[0], dirsig[1], dirsig[2], vector2[i][0], vector2[i][1], vector2[i][2]);
971 double tmp_EED_B = this->angulo(dirant[0], dirant[1], dirant[2], vector2[i][0], vector2[i][1], vector2[i][2]);
972 // angulos[i] =min( tmp_EED_A , tmp_EED_B );
973 angulos[i] = tmp_EED_A ;
974 if (tmp_EED_B<angulos[i] )
976 angulos[i]=tmp_EED_B;
978 if(maxmasa<this->vector[i][0]){
979 maxmasa=this->vector[i][0];
984 for(i=0;i<this->label;i++){
985 for(j=i;j<this->label;j++){
986 if(angulos[j]>angulos[i]){
987 vectortemp[0]=this->vector[i][0];
988 vectortemp[1]=this->vector[i][1];
989 vectortemp[2]=this->vector[i][2];
990 vectortemp[3]=this->vector[i][3];
991 vector2temp[0]=vector2[i][0];
992 vector2temp[1]=vector2[i][1];
993 vector2temp[2]=vector2[i][2];
994 angulostemp=angulos[i];
996 this->vector[i][0]=this->vector[j][0];
997 this->vector[i][1]=this->vector[j][1];
998 this->vector[i][2]=this->vector[j][2];
999 this->vector[i][3]=this->vector[j][3];
1000 vector2[i][0]=vector2[j][0];
1001 vector2[i][1]=vector2[j][1];
1002 vector2[i][2]=vector2[j][2];
1003 angulos[i]=angulos[j];
1005 this->vector[j][0]=vectortemp[0];
1006 this->vector[j][1]=vectortemp[1];
1007 this->vector[j][2]=vectortemp[2];
1008 this->vector[j][3]=vectortemp[3];
1009 vector2[j][0]=vector2temp[0];
1010 vector2[j][1]=vector2temp[1];
1011 vector2[j][2]=vector2temp[2];
1012 angulos[j]=angulostemp;
1019 for(i=this->label-2;i<this->label;i++){
1020 this->vector[i][0]=0;
1025 for(i=0;i<this->label;i++){
1026 for(j=i;j<this->label;j++){
1027 if(this->vector[j][0]>this->vector[i][0]){
1028 vectortemp[0]=this->vector[i][0];
1029 vectortemp[1]=this->vector[i][1];
1030 vectortemp[2]=this->vector[i][2];
1031 vectortemp[3]=this->vector[i][3];
1032 vector2temp[0]=vector2[i][0];
1033 vector2temp[1]=vector2[i][1];
1034 vector2temp[2]=vector2[i][2];
1035 angulostemp=angulos[i];
1037 this->vector[i][0]=this->vector[j][0];
1038 this->vector[i][1]=this->vector[j][1];
1039 this->vector[i][2]=this->vector[j][2];
1040 this->vector[i][3]=this->vector[j][3];
1041 vector2[i][0]=vector2[j][0];
1042 vector2[i][1]=vector2[j][1];
1043 vector2[i][2]=vector2[j][2];
1044 angulos[i]=angulos[j];
1046 this->vector[j][0]=vectortemp[0];
1047 this->vector[j][1]=vectortemp[1];
1048 this->vector[j][2]=vectortemp[2];
1049 this->vector[j][3]=vectortemp[3];
1050 vector2[j][0]=vector2temp[0];
1051 vector2[j][1]=vector2temp[1];
1052 vector2[j][2]=vector2temp[2];
1053 angulos[j]=angulostemp;
1060 for(i=0;i<this->label;i++){
1061 if(maxmasa*this->minpropmasa<this->vector[i][0]){
1062 puntoactualdis[0]=(puntoactualarr[0]+vector2[i][0]-extprin[0])*espprin[0];
1063 puntoactualdis[1]=(puntoactualarr[1]+vector2[i][1]-extprin[2])*espprin[1];
1064 puntoactualdis[2]=(puntoactualarr[2]+vector2[i][2])*espprin[2];
1066 points2->InsertPoint(buenos2,puntoactual);
1067 points3->InsertPoint(buenos2,puntoactualdis);
1068 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]));