1 /*=========================================================================
4 =========================================================================*/
5 #include "carotidaBifurcacion.h"
7 #include "vtkPolyData.h"
8 #include "vtkObjectFactory.h"
9 #include "vtkImageThreshold.h"
10 #include "vtkImageCast.h"
11 #include "vtkImageSeedConnectivity.h"
12 #include "vtkImageData.h"
13 #include "vtkMarchingCubes.h"
14 #include "vtkDoubleArray.h"
15 #include "vtkPointData.h"
16 #include "vtkExtractVOI.h"
21 vtkStandardNewMacro(carotidaBifurcacion);
24 carotidaBifurcacion::carotidaBifurcacion() {
26 this->NumberOfRequiredInputs = 2;
31 this->minpropmasa=0.001;
33 this->extrac = vtkExtractVOI::New();
34 this->extrac->SetSampleRate(1, 1, 1);
36 this->thresh = vtkImageThreshold::New();
37 this->thresh->SetInput(this->extrac->GetOutput());
38 this->thresh->SetInValue(255);
39 this->thresh->SetOutValue(0);
40 this->thresh->ReleaseDataFlagOff();
42 this->cast = vtkImageCast::New();
43 this->cast->SetInput(this->thresh->GetOutput());
44 this->cast->SetOutputScalarTypeToUnsignedChar();
46 this->connect = vtkImageSeedConnectivity::New();
47 this->connect->SetInput(this->cast->GetOutput());
48 this->connect->SetInputConnectValue(255);
49 this->connect->SetOutputConnectedValue(255);
50 this->connect->SetOutputUnconnectedValue(0);
52 this->dataprov=vtkImageData::New();
53 this->dataprov->SetScalarTypeToUnsignedChar();
59 //----------------------------------------------------------------------------
60 // Specify the input data or filter.
61 void carotidaBifurcacion::SetInput(vtkPolyData *input)
63 this->vtkProcessObject::SetNthInput(0, input);
67 //----------------------------------------------------------------------------
68 // Specify the input data or filter.
69 void carotidaBifurcacion::SetInput2(vtkImageData *input)
71 this->vtkProcessObject::SetNthInput(1, input);
72 this->extrac->SetInput(this->GetInput2());
77 //----------------------------------------------------------------------------
78 // Specify the input data or filter.
80 void carotidaBifurcacion::Execute()
82 vtkIdType totalpuntos=this->GetInput()->GetNumberOfPoints();
83 this->points = vtkPoints::New();
84 this->points2 = vtkPoints::New();
85 this->points3 = vtkPoints::New();
86 this->points4 = vtkPoints::New();
89 /*this->lineas = vtkCellArray::New();
90 this->salidas= vtkDoubleArray::New();
91 this->salidas->SetNumberOfComponents(3);*/
93 for(this->iter=0;this->iter<totalpuntos;this->iter++){
94 this->avanzar(this->iter);
97 this->limpiar(points, points2, points3, points4 );
99 this->GetOutput()->SetPoints (this->points);
100 /*this->GetOutput()->SetLines(this->lineas);
101 this->GetOutput()->GetPointData()->SetVectors(this->salidas);*/
110 //----------------------------------------------------------------------------
111 // Specify the input data or filter.
112 vtkPolyData *carotidaBifurcacion::GetInput(){
113 if (this->NumberOfInputs < 1){
117 return (vtkPolyData *)(this->Inputs[0]);
126 //----------------------------------------------------------------------------
127 // Specify the input data or filter.
128 vtkImageData *carotidaBifurcacion::GetInput2()
130 if (this->NumberOfInputs < 2){
134 return (vtkImageData *)(this->Inputs[1]);
143 //----------------------------------------------------------------------------
144 void carotidaBifurcacion::PrintSelf(ostream& os, vtkIndent indent)
146 this->Superclass::PrintSelf(os,indent);
156 //----------------------------------------------------------------------------
157 void carotidaBifurcacion::SetMaxPropRadio(double value)
159 this->maxpropradio=value;
164 //----------------------------------------------------------------------------
165 double carotidaBifurcacion::GetMaxPropRadio()
167 return this->maxpropradio;
172 //----------------------------------------------------------------------------
174 void carotidaBifurcacion::SetHumbral(double value)
181 //----------------------------------------------------------------------------
183 double carotidaBifurcacion::GetHumbral()
185 return this->humbral;
189 //----------------------------------------------------------------------------
191 void carotidaBifurcacion::SetMaxPropMasa(double value)
193 this->maxpropmasa=value;
198 //----------------------------------------------------------------------------
200 double carotidaBifurcacion::GetMaxPropMasa()
202 return this->maxpropmasa;
206 //----------------------------------------------------------------------------
208 void carotidaBifurcacion::SetMinPropMasa(double value)
210 this->minpropmasa=value;
215 //----------------------------------------------------------------------------
217 double carotidaBifurcacion::GetMinPropMasa()
219 return this->minpropmasa;
228 void carotidaBifurcacion::searc(int i, int j, int k)
239 ptr=(unsigned char *) this->connect->GetOutput()->GetScalarPointer(i,j,k);
240 ptr2=(unsigned char *) this->dataprov->GetScalarPointer(i,j,k);
244 this->connect->GetOutput()->GetExtent(ext);
246 radio=(ext[1]-ext[0])/2;
254 this->vector[this->label-1][0]+=1;
255 this->vector[this->label-1][1]+=i;
256 this->vector[this->label-1][2]+=j;
257 this->vector[this->label-1][3]+=k;
260 for(i3=-1;i3<=1;i3++){
261 for(j3=-1;j3<=1;j3++){
262 for(k3=-1;k3<=1;k3++){
263 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] ){
264 ptr=(unsigned char *) this->connect->GetOutput()->GetScalarPointer(i+i3,j+j3,k+k3);
265 ptr2=(unsigned char *) this->dataprov->GetScalarPointer(i+i3,j+j3,k+k3);
266 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))){
267 this->searc(i+i3, j+j3, k+k3);
277 void carotidaBifurcacion::find_components( )
287 this->connect->GetOutput()->GetExtent(ext);
290 radio=(ext[1]-ext[0])/2;
296 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
297 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
298 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
299 ptr=(unsigned char *) this->connect->GetOutput()->GetScalarPointer(i,j,k);
300 ptr2=(unsigned char *) this->dataprov->GetScalarPointer(i,j,k);
302 if(*ptr==255 && *ptr2==0 && (radio*radio)>(i2*i2)+(j2*j2)+(k2*k2) && ((radio-2)*(radio-2))<(i2*i2)+(j2*j2)+(k2*k2) ){
303 this->label=this->label+1;
304 this->vector[this->label-1][0]=0;
305 this->vector[this->label-1][1]=0;
306 this->vector[this->label-1][2]=0;
307 this->vector[this->label-1][3]=0;
308 this->searc(i, j, k);
321 unsigned short carotidaBifurcacion::maximo( )
325 unsigned short max=0;
327 this->extrac->GetOutput()->GetExtent(ext);
332 for(i=ext[0];i<=ext[1];i++){
333 for(j=ext[2];j<=ext[3];j++){
334 for(k=ext[4];k<=ext[5];k++){
335 ptr=(unsigned short *)this->extrac->GetOutput()->GetScalarPointer(i,j,k);
351 void carotidaBifurcacion::blanquear()
356 this->dataprov->GetExtent(ext);
362 for(i=ext[0];i<=ext[1];i++){
363 for(j=ext[2];j<=ext[3];j++){
364 for(k=ext[4];k<=ext[5];k++){
365 ptr=(unsigned char *) this->dataprov->GetScalarPointer(i,j,k);
376 double carotidaBifurcacion::angulo(double i1, double j1, double k1, double i2, double j2, double k2 )
378 double m1=sqrt((i1*i1)+(j1*j1)+(k1*k1));
379 double m2=sqrt((i2*i2)+(j2*j2)+(k2*k2));
380 double res=(acos(((i1*i2)+(j1*j2)+(k1*k2))/(m1*m2))*180)/3.1415;
386 if(res>=0 && res<180){
398 void carotidaBifurcacion::corte(double punto1[3], double punto2[3], double punto3[3], double centro[3], double radio )
408 m1=punto2[0]-punto1[0];
409 m2=punto2[1]-punto1[1];
410 m3=punto2[2]-punto1[2];
412 b1=punto1[0]-centro[0];
413 b2=punto1[1]-centro[1];
414 b3=punto1[2]-centro[2];
416 c0=m1*m1+m2*m2+m3*m3;
417 c1=2*m1*b1+2*m2*b2+2*m3*b3;
418 c2=b1*b1+b2*b2+b3*b3-radio*radio;
420 roots=vtkMath::SolveQuadratic (c0, c1, c2);
422 if(roots[1]>=0 && roots[1]<=1){
429 punto3[0]=punto1[0]+root*(punto2[0]-punto1[0]);
430 punto3[1]=punto1[1]+root*(punto2[1]-punto1[1]);
431 punto3[2]=punto1[2]+root*(punto2[2]-punto1[2]);
439 void carotidaBifurcacion::direcciones(vtkPolyData *profile, int iter, double radio, double puntocortea[3], double puntocortes[3] )
441 double puntoactual[3];
449 vtkIdType totalpuntos=profile->GetNumberOfPoints();
451 if(iter<totalpuntos){
453 profile->GetPoint(iter, puntoactual);
456 profile->GetPoint(iter+1, puntosig);
457 profile->GetPoint(iter, puntosig2);
460 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
464 profile->GetPoint(iter+1+i, puntosig);
465 profile->GetPoint(iter+i, puntosig2);
466 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
471 this->corte(puntosig, puntosig2, puntocortes, puntoactual, radio );
473 puntocortea[0]=2*puntoactual[0]-puntocortes[0];
474 puntocortea[1]=2*puntoactual[1]-puntocortes[1];
475 puntocortea[2]=2*puntoactual[2]-puntocortes[2];
478 else if(iter==totalpuntos-1){
479 profile->GetPoint(iter-1, puntoant);
480 profile->GetPoint(iter, puntoant2);
483 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
487 profile->GetPoint(iter-1-j, puntoant);
488 profile->GetPoint(iter-j, puntoant2);
489 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
494 this->corte(puntoant, puntoant2, puntocortea, puntoactual, radio );
496 puntocortes[0]=2*puntoactual[0]-puntocortea[0];
497 puntocortes[1]=2*puntoactual[1]-puntocortea[1];
498 puntocortes[2]=2*puntoactual[2]-puntocortea[2];
501 profile->GetPoint(iter+1, puntosig);
502 profile->GetPoint(iter, puntosig2);
505 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
507 while(dist< radio && iter+1+i<totalpuntos){
509 profile->GetPoint(iter+1+i, puntosig);
510 profile->GetPoint(iter+i, puntosig2);
511 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
516 profile->GetPoint(iter-1, puntoant);
517 profile->GetPoint(iter, puntoant2);
520 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
522 while(dist< radio && iter-1-j>=0){
524 profile->GetPoint(iter-1-j, puntoant);
525 profile->GetPoint(iter-j, puntoant2);
526 dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
531 if(iter+1+i<totalpuntos){
532 this->corte(puntosig, puntosig2, puntocortes, puntoactual, radio );
536 this->corte(puntoant, puntoant2, puntocortea, puntoactual, radio );
539 if(iter+1+i>=totalpuntos){
540 puntocortes[0]=2*puntoactual[0]-puntocortea[0];
541 puntocortes[1]=2*puntoactual[1]-puntocortea[1];
542 puntocortes[2]=2*puntoactual[2]-puntocortea[2];
546 puntocortea[0]=2*puntoactual[0]-puntocortes[0];
547 puntocortea[1]=2*puntoactual[1]-puntocortes[1];
548 puntocortea[2]=2*puntoactual[2]-puntocortes[2];
555 int carotidaBifurcacion::igual(double puntoactual[3], double puntoactualdis[3], double puntoantguar[3], double puntoantguardis[3] )
570 disact=sqrt(vtkMath::Distance2BetweenPoints (puntoactual,puntoactualdis));
571 disant=sqrt(vtkMath::Distance2BetweenPoints (puntoantguar,puntoantguardis));
572 discru=sqrt(vtkMath::Distance2BetweenPoints (puntoactual,puntoantguar));
575 if(discru<disact+disant){
576 mat[0][0]=puntoantguardis[0]-puntoantguar[0];
577 mat[0][1]=-(puntoactualdis[0]-puntoactual[0]);
580 mat[1][0]=puntoantguardis[1]-puntoantguar[1];
581 mat[1][1]=-(puntoactualdis[1]-puntoactual[1]);
588 b[0]=puntoactual[0]-puntoantguar[0];
589 b[1]=puntoactual[1]-puntoantguar[1];
593 vtkMath::Invert3x3(mat,matinv);
594 vtkMath::Multiply3x3 (matinv, b, t);
598 if(t[0]>0 && t[1]>0 ){
610 int carotidaBifurcacion::proporcion(vtkImageData *data )
617 data->GetExtent(ext);
621 for(i=ext[0];i<=ext[1];i++){
622 for(j=ext[2];j<=ext[3];j++){
623 for(k=ext[4];k<=ext[5];k++){
624 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
634 return (max*100)/total;
641 void carotidaBifurcacion::redondear(vtkImageData *data )
647 data->GetExtent(ext);
651 radio=(ext[1]-ext[0])/2;
654 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
655 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
656 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
657 ptr=(unsigned char *) data->GetScalarPointer(i,j,k);
658 tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2) ;
659 if( radio<sqrt(tmpsqrt) ){
671 void carotidaBifurcacion::limpiar(vtkPoints *profile, vtkPoints *profile2, vtkPoints *profile3, vtkPoints *profile4 )
674 double puntoactual[3];
675 double puntoactualdis[3];
676 double puntoantguar[3];
677 double puntoantguardis[3];
684 vtkIdType totalpuntos=profile2->GetNumberOfPoints();
688 for(i=0;i<totalpuntos;i++){
690 profile2->GetPoint(i, puntoactual);
691 profile3->GetPoint(i, puntoactualdis);
693 for(j=i+1;j<totalpuntos;j++){
694 profile2->GetPoint(j, puntoantguar);
695 profile3->GetPoint(j, puntoantguardis);
696 if(puntoactual[0]!=puntoantguar[0] && puntoactual[1]!=puntoantguar[1] && puntoactual[2]!=puntoantguar[2] && 1==this->igual(puntoactual, puntoactualdis, puntoantguar, puntoantguardis)){
729 total[l][m]=pairs[i][0];
733 total[l][m]=pairs[i][1];
736 if(pairs[i][0]!=-1 || pairs[i][1]!=-1){
746 if(total[l-1][n]==pairs[j][0]){
749 if(total[l-1][p]==pairs[j][1]){
754 total[l-1][m]=pairs[j][1];
762 else if(total[l-1][n]==pairs[j][1]){
765 if(total[l-1][p]==pairs[j][0]){
770 total[l-1][m]=pairs[j][0];
793 profile4->GetPoint(total[i][m], puntoactual);
794 promedio[0]+=puntoactual[0];
795 promedio[1]+=puntoactual[1];
796 promedio[2]+=puntoactual[2];
804 promedio[0]=promedio[0]/cant;
805 promedio[1]=promedio[1]/cant;
806 promedio[2]=promedio[2]/cant;
808 profile->InsertPoint(i,promedio);
817 void carotidaBifurcacion::avanzar(vtkIdType iter)
823 double puntoactual[3];
824 double puntoactualdis[3];
825 int puntoactualarr[3];
832 double vector2[50][3];
835 unsigned long vectortemp[4];
836 double vector2temp[3];
851 this->GetInput2()->GetExtent(extprin);
852 this->GetInput2()->GetSpacing(espprin);
854 vtkIdType totalpuntos=this->GetInput()->GetNumberOfPoints();
856 if(iter<totalpuntos){
858 vtkDoubleArray *radio = (vtkDoubleArray *)this->GetInput()->GetPointData()->GetScalars("radio");
859 //vtkDoubleArray *dir = (vtkDoubleArray *)this->GetInput()->GetPointData()->GetVectors("norm");
860 this->GetInput()->GetPoint(iter, puntoactual);
861 radioactual=radio->GetValue(iter);
864 for(this->label=0, prop=100;radiotrabajo<radioactual*this->maxpropradio && (this->label<2 || prop>this->maxpropmasa);radiotrabajo++){
866 puntoactualarr[0]= (int) ( (puntoactual[0]/espprin[0])+extprin[0] );
867 puntoactualarr[1]= (int) ( (puntoactual[1]/espprin[1])+extprin[2] );
868 puntoactualarr[2]= (int) ( (puntoactual[2]/espprin[2]) );
870 extint[0]=puntoactualarr[0]-radiotrabajo;
871 extint[1]=puntoactualarr[0]+radiotrabajo;
872 extint[2]=puntoactualarr[1]-radiotrabajo;
873 extint[3]=puntoactualarr[1]+radiotrabajo;
874 extint[4]=puntoactualarr[2]-radiotrabajo;
875 extint[5]=puntoactualarr[2]+radiotrabajo;
877 extrac->SetVOI(extint);
878 extrac->UpdateWholeExtent();
880 extrac->GetOutput()->GetExtent(extintreal);
882 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] ){
887 this->thresh->ThresholdByUpper(this->maximo()*this->humbral);
888 this->thresh->UpdateWholeExtent();
889 this->thresh->Update();
891 prop=this->proporcion(this->thresh->GetOutput());
892 this->redondear(this->thresh->GetOutput());
895 this->cast->UpdateWholeExtent();
896 this->cast->Update();
899 this->connect->RemoveAllSeeds();
900 this->connect->AddSeed(puntoactualarr[0],puntoactualarr[1],puntoactualarr[2]);
901 this->connect->UpdateWholeExtent();
902 this->connect->Update();
904 this->dataprov->Delete();
906 this->dataprov=vtkImageData::New();
907 this->dataprov->SetScalarTypeToUnsignedChar();
909 this->dataprov->SetExtent(this->connect->GetOutput()->GetExtent());
913 this->find_components();
923 this->direcciones(this->GetInput(), this->iter, radiotrabajo-1, dirant, dirsig);
927 dirant[0]=dirant[0]-puntoactual[0];
928 dirant[1]=dirant[1]-puntoactual[1];
929 dirant[2]=dirant[2]-puntoactual[2];
930 dirsig[0]=dirsig[0]-puntoactual[0];
931 dirsig[1]=dirsig[1]-puntoactual[1];
932 dirsig[2]=dirsig[2]-puntoactual[2];
941 for(i=0;i<this->label;i++){
942 vector2[i][0]=((double)this->vector[i][1]/(double)this->vector[i][0])-(double)puntoactualarr[0];
943 vector2[i][1]=((double)this->vector[i][2]/(double)this->vector[i][0])-(double)puntoactualarr[1];
944 vector2[i][2]=((double)this->vector[i][3]/(double)this->vector[i][0])-(double)puntoactualarr[2];
945 double tmp_EED_A = this->angulo(dirsig[0], dirsig[1], dirsig[2], vector2[i][0], vector2[i][1], vector2[i][2]);
946 double tmp_EED_B = this->angulo(dirant[0], dirant[1], dirant[2], vector2[i][0], vector2[i][1], vector2[i][2]);
947 // angulos[i] =min( tmp_EED_A , tmp_EED_B );
948 angulos[i] = tmp_EED_A ;
949 if (tmp_EED_B<angulos[i] )
951 angulos[i]=tmp_EED_B;
953 if(maxmasa<this->vector[i][0]){
954 maxmasa=this->vector[i][0];
959 for(i=0;i<this->label;i++){
960 for(j=i;j<this->label;j++){
961 if(angulos[j]>angulos[i]){
962 vectortemp[0]=this->vector[i][0];
963 vectortemp[1]=this->vector[i][1];
964 vectortemp[2]=this->vector[i][2];
965 vectortemp[3]=this->vector[i][3];
966 vector2temp[0]=vector2[i][0];
967 vector2temp[1]=vector2[i][1];
968 vector2temp[2]=vector2[i][2];
969 angulostemp=angulos[i];
971 this->vector[i][0]=this->vector[j][0];
972 this->vector[i][1]=this->vector[j][1];
973 this->vector[i][2]=this->vector[j][2];
974 this->vector[i][3]=this->vector[j][3];
975 vector2[i][0]=vector2[j][0];
976 vector2[i][1]=vector2[j][1];
977 vector2[i][2]=vector2[j][2];
978 angulos[i]=angulos[j];
980 this->vector[j][0]=vectortemp[0];
981 this->vector[j][1]=vectortemp[1];
982 this->vector[j][2]=vectortemp[2];
983 this->vector[j][3]=vectortemp[3];
984 vector2[j][0]=vector2temp[0];
985 vector2[j][1]=vector2temp[1];
986 vector2[j][2]=vector2temp[2];
987 angulos[j]=angulostemp;
994 for(i=this->label-2;i<this->label;i++){
995 this->vector[i][0]=0;
1000 for(i=0;i<this->label;i++){
1001 for(j=i;j<this->label;j++){
1002 if(this->vector[j][0]>this->vector[i][0]){
1003 vectortemp[0]=this->vector[i][0];
1004 vectortemp[1]=this->vector[i][1];
1005 vectortemp[2]=this->vector[i][2];
1006 vectortemp[3]=this->vector[i][3];
1007 vector2temp[0]=vector2[i][0];
1008 vector2temp[1]=vector2[i][1];
1009 vector2temp[2]=vector2[i][2];
1010 angulostemp=angulos[i];
1012 this->vector[i][0]=this->vector[j][0];
1013 this->vector[i][1]=this->vector[j][1];
1014 this->vector[i][2]=this->vector[j][2];
1015 this->vector[i][3]=this->vector[j][3];
1016 vector2[i][0]=vector2[j][0];
1017 vector2[i][1]=vector2[j][1];
1018 vector2[i][2]=vector2[j][2];
1019 angulos[i]=angulos[j];
1021 this->vector[j][0]=vectortemp[0];
1022 this->vector[j][1]=vectortemp[1];
1023 this->vector[j][2]=vectortemp[2];
1024 this->vector[j][3]=vectortemp[3];
1025 vector2[j][0]=vector2temp[0];
1026 vector2[j][1]=vector2temp[1];
1027 vector2[j][2]=vector2temp[2];
1028 angulos[j]=angulostemp;
1035 for(i=0;i<this->label;i++){
1036 if(maxmasa*this->minpropmasa<this->vector[i][0]){
1037 puntoactualdis[0]=(puntoactualarr[0]+vector2[i][0]-extprin[0])*espprin[0];
1038 puntoactualdis[1]=(puntoactualarr[1]+vector2[i][1]-extprin[2])*espprin[1];
1039 puntoactualdis[2]=(puntoactualarr[2]+vector2[i][2])*espprin[2];
1041 points2->InsertPoint(buenos2,puntoactual);
1042 points3->InsertPoint(buenos2,puntoactualdis);
1043 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]));