path_name='d:/OCT/SOURCES/3.DCM';
fid = fopen(path_name, 'r');
dataa = fread(fid,'uint8');
fclose(fid);
[header_dicom,Ls]=OCT_head_read(dataa);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ls=medfilt2(Ls,[7 7]);
Ls=imresize(Ls,[256 512 ]);
figure; imshow(Ls,[]);
L2=Ls;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
przed=1;
L22=imclose(Ls>max(max(Ls))*0.10,ones([3 3]));
L2lab=bwlabel(L22);
L33=zeros(size(L22));
for ir=1:max(L2lab(:))
L2_=L2lab==ir;
if sum(L2_(:))>100
L33=L33|L2_;
end
end
figure; imshow(L33,[]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L22=bwfill(L33,'holes');
figure; imshow(L22,[]);
L55=bwlabel(xor(L22,L33));
for ir=1:max(L55(:))
L5_=L55==ir;
if sum(L5_(:))<100
L33=L33|L5_;
end
end
L22=L33;
figure; imshow(L22,[]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
linie_12=[];
for i=1:size(L22,2)
Lf=L22(:,i);
Lff=bwlabel(Lf);
if sum(Lff)>0
clear Lnr
for yt=1:max(Lff(:))
Lffd=Lff==yt;
if sum(Lffd(:))>10
Lnr=[(1:length(Lffd))',Lffd]; Lnr(Lnr(:,2)==0,:)=[];
break
end
end
if (exist('Lnr')>0)&(~isempty(Lnr))
linie_12=[linie_12; [i Lnr(1,1) Lnr(end,1)]];
end
end
end
hold on;
plot(linie_12(:,1),linie_12(:,2),'r*'); grid on
plot(linie_12(:,1),linie_12(:,3),'g*'); grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
linie_12(:,2)=medfilt2(linie_12(:,2),[5 1]);
linie_12(:,3)=medfilt2(linie_12(:,3),[5 1]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=linie_12(:,1);
y=linie_12(:,2);
ybw=bwlabel(abs([diff(y') 0])<5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rzad=3;
toler=10;
P=polyfit(x,y,rzad);
Y=polyval(P,x);
yyy=Y-y;
pamm=[0 0 sum( abs(yyy)<toler )/length(yyy)];
for ir=1:(max(ybw)-1)
for irr=(ir+1):max(ybw)
y_=[y(ybw==ir); y(ybw==irr)];
x_=[x(ybw==ir); x(ybw==irr)];
P=polyfit(x_,y_,rzad);
Y=polyval(P,x);
hold on; plot(x,Y,'-g*');
yyy=Y-y;
pamm=[pamm; [ir irr sum( abs(yyy)<toler )/length(yyy) ]];
end
end
pamm_=sortrows(pamm,3); ir=pamm_(end,1); irr=pamm_(end,2);
if ir==0;
y_=y;
x_=x;
P=polyfit(x_,y_,rzad);
Y=polyval(P,x);
yyy=Y-y;
y_=y(abs(yyy)<toler);
x_=x(abs(yyy)<toler);
P=polyfit(x_,y_,rzad);
Y=polyval(P,x);
else
y_=[y(ybw==ir); y(ybw==irr)];
x_=[x(ybw==ir); x(ybw==irr)];
P=polyfit(x_,y_,rzad);
Y=polyval(P,x);
yyy=Y-y;
y_=y(abs(yyy)<toler);
x_=x(abs(yyy)<toler);
P=polyfit(x_,y_,rzad);
Y=polyval(P,x);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%
plot(x,Y,'b*'); grid on
y_1=Y;
x=linie_12(:,1);
y=linie_12(:,3);
ybw=bwlabel(abs([diff(y') 0])<5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%
rzad=3; toler=15;
P=polyfit(x,y,rzad); pamm=[];
Y=polyval(P,x);
yyy=Y-y;
if sum( (Y(:)-y_1(:))<0 )==0
pamm=[0 0 sum( abs(yyy)<toler )/length(yyy)];
end
for ir=1:(max(ybw)-1)
for irr=(ir+1):max(ybw)
y_=[y(ybw==ir); y(ybw==irr)];
x_=[x(ybw==ir); x(ybw==irr)];
P=polyfit(x_,y_,rzad);
Y=polyval(P,x);
yyy=Y-y;
if sum( (Y(:)-y_1(:))<0 )==0
pamm=[pamm; [ir irr sum( abs(yyy)<toler )/length(yyy) ]];
end
end
end
if size(pamm,1)>1
pamm_=sortrows(pamm,3);
ir=pamm_(end,1); irr=pamm_(end,2);
if ir==0;
y_=y;
x_=x;
P=polyfit(x_,y_,rzad);
Y=polyval(P,x);
yyy=Y-y;
y_=y(abs(yyy)<toler);
x_=x(abs(yyy)<toler);
P=polyfit(x_,y_,rzad);
Y=polyval(P,x);
else
y_=[y(ybw==ir); y(ybw==irr)];
x_=[x(ybw==ir); x(ybw==irr)];
P=polyfit(x_,y_,rzad);
Y=polyval(P,x);
yyy=Y-y;
y_=y(abs(yyy)<toler);
x_=x(abs(yyy)<toler);
P=polyfit(x_,y_,rzad);
Y=polyval(P,x);
end
else
x=[]; Y=[];P=[];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%
plot(x,Y,'m*'); grid on
y_2=Y;
figure; imshow(L33,[]); hold on
sf=zeros( [ 1 length(Y) ] ); sf_=zeros( [ 1 length(Y) ] );pole_=zeros([1 length(Y)]);
pole_x=zeros([1 length(Y)]); pole_y=zeros([1 length(Y)]);
p_zn=0; zakres_=40;
Lwys=zeros([zakres_ length(Y)-1]);
Lwys_bin=zeros([zakres_ length(Y)-1]);
L_gridXX=[]; L_gridYY=[];
for nb=1:(length(Y)-1)
PP=polyfit(x(nb:nb+1),Y(nb:nb+1),1);
PP2(2)=x(nb)/PP(1)+Y(nb);
PP2(1)=-1/PP(1);
if Y(nb)>Y(nb+1)
XX=x(nb):1:(x(nb)+zakres_);
else
XX=x(nb):-1:(x(nb)-zakres_);
end
YY=polyval(PP2,XX);
if (max(YY)-min(YY))>(zakres_+1)
YY=Y(nb):1:(Y(nb)+zakres_);
PP3(1)=1/PP2(1);
PP3(2)=-PP2(2)/PP2(1);
XX=polyval(PP3,YY);
plot(XX,YY,'r*'); grid on; hold on
XX(round(YY)>size(L2,1))=[]; YY(round(YY)>size(L2,1))=[];
YY(round(XX)>size(L2,2))=[]; XX(round(XX)>size(L2,2))=[];
for vc=1:length(XX); if (round(YY(vc))>0)&(round(XX(vc))>0) ;Lwys(vc,nb)=L2( round(YY(vc)), round(XX(vc)) ); Lwys_bin(vc,nb)=L22( round(YY(vc)), round(XX(vc)) ); end; end
L_gridXX(1:length(XX),nb)=XX; L_gridYY(1:length(YY),nb)=YY;
else
plot(XX,YY,'c*'); grid on; hold on
XX(round(YY)>size(L2,1))=[]; YY(round(YY)>size(L2,1))=[];
YY(round(XX)>size(L2,2))=[]; XX(round(XX)>size(L2,2))=[];
for vc=1:length(XX); if (round(YY(vc))>0)&(round(XX(vc))>0); Lwys(vc,nb)=L2( round(YY(vc)), round(XX(vc)) ); Lwys_bin(vc,nb)=L22( round(YY(vc)), round(XX(vc)) ); end; end
L_gridXX(1:length(XX),nb)=XX; L_gridYY(1:length(YY),nb)=YY;
end
end
figure; imshow(Lwys,[]);
Lwys_bin=imopen(Lwys_bin,ones(3));
Lss=sum(Lwys);
XX=1:length(Lss);
YY=Lss;
Lm=max(Lss(:));
XXl=XX(1:round(length(XX)/2));
XXp=XX(round(length(XX)/2):end);
YYl=YY(1:round(length(YY)/2));
YYp=YY(round(length(YY)/2):end);
YYlb=YYl>(Lm/4);
YYpb=YYp>(Lm/4);
nr_XXl=1:length(XXl);
nr_XXp=1:length(XXp);
XXl_max=nr_XXl(YYl==max(YYl));
XXp_max=nr_XXp(YYp==max(YYp));
figure
plot(XXl,YYl,'-r*'); hold on
plot(XXp,YYp,'-g*'); grid on
xlabel('XXl-red, XXp - green')
ylabel('YYl,YYp')
figure; imshow(Ls,[]); hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LEWA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YYlb_=bwlabel(YYlb);
pam_l=[];
for ty=1:max(YYlb_)
YYt=YYl(YYlb_==ty);
XXt=XXl(YYlb_==ty);
pam_l=[pam_l; [ty sum(YYt) XXt(end)]];
end
if size(pam_l,1)>0
pam_l=pam_l(YYlb_(XXl_max),:);
plot(pam_l(1,3),Y(pam_l(1,3)),'rs','MarkerSize',10)
end
xy_g_l=[];
xy_d_l=[];
for vv=pam_l(1,3):-1:1
pp=Lwys_bin(:,vv);
ppl=bwlabel(pp);