close all clear all % Cargo la Imagen Inicial Imagen_Inicial=imread('01_original.bmp'); figure(1) imshow(Imagen_Inicial) % Obtenemos la componente verde Verde=double(Imagen_Inicial(:,:,2))/255; figure(2) imshow(Verde) %Obtenemos el filtro de media aritmetica grande y se lo aplicamos a la %componente verde Filtro_Media=fspecial('average',[60,60]); Imagen_Filtrada=imfilter(Verde,Filtro_Media); figure(3) imshow(Imagen_Filtrada) Imagen_Normalizada=Imagen_Filtrada-Verde; figure(4) imshow(Imagen_Normalizada) [Filas,Columnas]=size(Imagen_Normalizada); %Eliminamos el borde de la imagen a traves de la mascara Mascara=imread('01_Mask.gif'); for i=1:Filas for j=1:Columnas if Mascara(i,j)==0 Imagen_Normalizada(i,j)=0; end end end figure(5) imshow(Imagen_Normalizada) %Thin Vessel enhancement H1=1/6*[-1 -1 -1; 2 2 2; -1 -1 -1]; H2=1/6*[2 -1 -1; -1 2 -1; -1 -1 2]; H3=1/6*[-1 2 -1; -1 2 -1; -1 2 -1]; H4=1/6*[-1 -1 2; -1 2 -1; 2 -1 -1]; Orientacion_0=imfilter(Imagen_Normalizada,H1); Orientacion_45=imfilter(Imagen_Normalizada,H2); Orientacion_90=imfilter(Imagen_Normalizada,H3); Orientacion_135=imfilter(Imagen_Normalizada,H4); for i=1:Filas for j=1:Columnas m=1; X(m)=Orientacion_0(i,j); m=m+1; X(m)=Orientacion_45(i,j); m=m+1; X(m)=Orientacion_90(i,j); m=m+1; X(m)=Orientacion_135(i,j); Realzado(i,j)=max(X); end end Imagen_Realzada=Imagen_Normalizada+Realzado; So8=strel('disk',8); Sc=strel('disk',1); TopHat8=Imagen_Realzada-min(imopen(imclose(Imagen_Realzada,Sc),So8),Imagen_Realzada); figure(6) imshow(Imagen_Realzada) % Obtenemos gradientes en direccion x e y [GradienteX,GradienteY]=gradient(TopHat8); % Obtenemos la magnitud del gradiente Mag_Gradiente=sqrt(GradienteX.^2+GradienteY.^2); % Obtenemos la orientacion del gradiente forzando que atan(0/0)=0 y asi no % de error for f=1:Filas for c=1:Columnas if GradienteX(f,c)==0 && GradienteY(f,c)==0 Orient_Gradiente(f,c)=0; else Orient_Gradiente(f,c)=atan(GradienteY(f,c)/GradienteX(f,c)); end end end % N*N tamaño de la vecindad alrededor de (i,j) N=31; % Obtenemos el valor M que se utilizara posteriormente para hacer los % sumatorios y delimitarlos M=(N-1)/2; % Pasamos a Calcular el angulo de anisotropia, utilizando variables % desconocidas hasta ahora que se explicaran a medida que vayan apareciendo Angulo_Anisotropia=[]; % Analizamos primero las M primeras filas Contador_F=0; for i=1:M Contador_C=0; for j=1:M resultado_num=0; resultado_den=0; for m=(M-Contador_F):(N-1) for n=(M-Contador_C):(N-1) resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_den==0 && resultado_num==0 Angulo_Anisotropia(i,j)=0; else Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den); end Contador_C=Contador_C+1; end Contador_F=Contador_F+1; end Contador_F=0; for i=1:M for j=(M+1):(Columnas-M) resultado_num=0; resultado_den=0; for m=(M-Contador_F):(N-1) for n=0:(N-1) resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_den==0 && resultado_num==0 Angulo_Anisotropia(i,j)=0; else Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den); end end Contador_F=Contador_F+1; end Contador_F=0; for i=1:M Contador_C=0; for j=(Columnas-M+1):Columnas resultado_num=0; resultado_den=0; for m=(M-Contador_F):(N-1) for n=0:(N-2-Contador_C) resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_den==0 && resultado_num==0 Angulo_Anisotropia(i,j)=0; else Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den); end Contador_C=Contador_C+1; end Contador_F=Contador_F+1; end %Analizamos desde la fila M+1 hasta la Fila-M for i=(M+1):(Filas-M) Contador_C=0; for j=1:M resultado_num=0; resultado_den=0; for m=0:(N-1) % Recorremos todas las filas for n=(M-Contador_C):(N-1) % Recorremos todas las columnas dentro de cada fila resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_den==0 && resultado_num==0 Angulo_Anisotropia(i,j)=0; else Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den); end Contador_C=Contador_C+1; end for j=(M+1):(Columnas-M) resultado_num=0; resultado_den=0; for m=0:(N-1) % Recorremos todas las filas for n=0:(N-1) % Recorremos todas las columnas dentro de cada fila resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_den==0 && resultado_num==0 Angulo_Anisotropia(i,j)=0; else Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den); end end Contador_C=0; for j=(Columnas-M+1):Columnas resultado_num=0; resultado_den=0; for m=0:(N-1) for n=0:(N-2-Contador_C) resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_den==0 && resultado_num==0 Angulo_Anisotropia(i,j)=0; else Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den); end Contador_C=Contador_C+1; end end % Analizamos las ultimas M Filas Contador_F=0; for i=(Filas-M+1):Filas Contador_C=0; for j=1:M resultado_num=0; resultado_den=0; for m=0:(N-2-Contador_F) % Recorremos todas las filas for n=(M-Contador_C):(N-1) % Recorremos todas las columnas dentro de cada fila resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_den==0 && resultado_num==0 Angulo_Anisotropia(i,j)=0; else Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den); end Contador_C=Contador_C+1; end Contador_F=Contador_F+1; end Contador_F=0; for i=(Filas-M+1):Filas for j=(M+1):(Columnas-M) resultado_num=0; resultado_den=0; for m=0:(N-2-Contador_F) for n=0:(N-1) resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_den==0 && resultado_num==0 Angulo_Anisotropia(i,j)=0; else Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den); end end Contador_F=Contador_F+1; end Contador_F=0; for i=(Filas-M+1):Filas Contador_C=0; for j=(Columnas-M+1):Columnas resultado_num=0; resultado_den=0; for m=0:(N-2-Contador_F) for n=0:(N-2-Contador_C) resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n))); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_den==0 && resultado_num==0 Angulo_Anisotropia(i,j)=0; else Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den); end Contador_C=Contador_C+1; end Contador_F=Contador_F+1; end figure(7) imshow(Angulo_Anisotropia) Angulo_Orientacion=Angulo_Anisotropia+(pi/2); % Deben estar todos comprendidos entre 0 y pi Angulo_Anisotropia=abs(Angulo_Anisotropia); figure(8) imshow(Angulo_Anisotropia) figure(9) imshow(Angulo_Orientacion) % Calculamos la coherencia Coherencia=[]; Contador_F=0; for i=1:M Contador_C=0; for j=1:M resultado_num=0; resultado_den=0; for m=(M-Contador_F):(N-1) % Recorremos todas las filas for n=(M-Contador_C):(N-1) % Recorremos todas las columnas dentro de cada fila resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j)); resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_num==0&resultado_den==0 Coherencia(i,j)=0; else Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den); end Contador_C=Contador_C+1; end Contador_F=Contador_F+1; end Contador_F=0; for i=1:M for j=(M+1):(Columnas-M) resultado_num=0; resultado_den=0; for m=(M-Contador_F):(N-1) for n=0:(N-1) resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j)); resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_num==0&resultado_den==0 Coherencia(i,j)=0; else Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den); end end Contador_F=Contador_F+1; end Contador_F=0; for i=1:M Contador_C=0; for j=(Columnas-M+1):Columnas resultado_num=0; resultado_den=0; for m=(M-Contador_F):(N-1) for n=0:(N-2-Contador_C) resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j)); resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_num==0&resultado_den==0 Coherencia(i,j)=0; else Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den); end Contador_C=Contador_C+1; end Contador_F=Contador_F+1; end %Analizamos desde la fila M+1 hasta la Fila-M for i=(M+1):(Filas-M) Contador_C=0; for j=1:M resultado_num=0; resultado_den=0; for m=0:(N-1) % Recorremos todas las filas for n=(M-Contador_C):(N-1) % Recorremos todas las columnas dentro de cada fila resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j)); resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_num==0&resultado_den==0 Coherencia(i,j)=0; else Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den); end Contador_C=Contador_C+1; end for j=(M+1):(Columnas-M) resultado_num=0; resultado_den=0; for m=0:(N-1) % Recorremos todas las filas for n=0:(N-1) % Recorremos todas las columnas dentro de cada fila resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j)); resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_num==0&resultado_den==0 Coherencia(i,j)=0; else Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den); end end Contador_C=0; for j=(Columnas-M+1):Columnas resultado_num=0; resultado_den=0; for m=0:(N-1) for n=0:(N-2-Contador_C) resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j)); resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_num==0&resultado_den==0 Coherencia(i,j)=0; else Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den); end Contador_C=Contador_C+1; end end % Analizamos las ultimas M Filas Contador_F=0; for i=(Filas-M+1):Filas Contador_C=0; for j=1:M resultado_num=0; resultado_den=0; for m=0:(N-2-Contador_F) % Recorremos todas las filas for n=(M-Contador_C):(N-1) % Recorremos todas las columnas dentro de cada fila resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j)); resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_num==0&resultado_den==0 Coherencia(i,j)=0; else Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den); end Contador_C=Contador_C+1; end Contador_F=Contador_F+1; end Contador_F=0; for i=(Filas-M+1):Filas for j=(M+1):(Columnas-M) resultado_num=0; resultado_den=0; for m=0:(N-2-Contador_F) for n=0:(N-1) resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j)); resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_num==0&resultado_den==0 Coherencia(i,j)=0; else Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den); end end Contador_F=Contador_F+1; end Contador_F=0; for i=(Filas-M+1):Filas Contador_C=0; for j=(Columnas-M+1):Columnas resultado_num=0; resultado_den=0; for m=0:(N-2-Contador_F) for n=0:(N-2-Contador_C) resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j)); resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n); resultado_num=resultado_num+resultado_provisional_num; resultado_den=resultado_den+resultado_provisional_den; end end if resultado_num==0&resultado_den==0 Coherencia(i,j)=0; else Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den); end Contador_C=Contador_C+1; end Contador_F=Contador_F+1; end figure(10) imshow(Coherencia) Maximo_Valor=max(max(Coherencia)); Coherencia_Maxima=Coherencia/Maximo_Valor; figure(11) imshow(Coherencia_Maxima)