-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathorientation.m
162 lines (115 loc) · 4.15 KB
/
orientation.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
function [imgs_o,ov,cm_sum_mean_vert,vari_vert,ocheckhorz,vari_horz]=orientation(imgs,tiporf)
%Calcula la orientacion de la imagen img de tamaño XxYxZ.
%devuelve un vector de 3 componentes con valores comprendidos entre -3 y 3
%y la imagen orientada correctamente, de forma que:
A=size(imgs);
posdim=[1 2 3];
BW=zeros(A);
level=graythresh(imgs(:,:));
if level==0;
level=0.1;
end
for i=1:A(3)
imint=im2bw(imgs(:,:,i),level);
se = strel('disk',1);
imintop=imopen(imint,se);
BW(:,:,i)=imfill(squeeze(imintop),'holes');
end
%%%%%%%%%%%
% Encuentra las lineas:
% dim1 = la linea oreja-oreja
% dim2 = la linea frente-cogote
% dim3 = la linea coronilla-cuello
%%%%%%%%%%%%
% Path para error cuando las imágenes tienen mucho espacio en negro.
zmin = find(sum(sum(BW,1),2)>0,1,'first');
zmax = find(sum(sum(BW,1),2)>0,1,'last');
ymin = find(sum(sum(BW,1),3)>0,1,'first');
ymax = find(sum(sum(BW,1),3)>0,1,'last');
xmin = find(sum(sum(BW,2),3)>0,1,'first');
xmax = find(sum(sum(BW,2),3)>0,1,'last');
BW = BW(xmin:xmax, ymin:ymax, zmin:zmax);
imgs_mod = imgs(xmin:xmax, ymin:ymax, zmin:zmax);
cc=bwconncomp(BW);
stats=regionprops(cc);fprintf('.')
[~,rg]=max([stats.Area]);
reBoundingBox(1)=stats(rg).BoundingBox(2);
reBoundingBox(2)=stats(rg).BoundingBox(1);
reBoundingBox(3)=stats(rg).BoundingBox(3);
reBoundingBox(4)=stats(rg).BoundingBox(5);
reBoundingBox(5)=stats(rg).BoundingBox(4);
reBoundingBox(6)=stats(rg).BoundingBox(6);
%% Encuentra el eje mayor (coincidirá con la linea frente-cogote)
%%ax,ax1,vga1,vga2,ccoerf1,ccoerf2,pes,
[~,x]=max(reBoundingBox);
dim2=(x-3==posdim);
%% Encuentra el eje con la linea oreja-oreja)
%%
dimes=find(~dim2);
Ct=squeeze(sum(imgs_mod,find(dim2)));
windowt=Ct(ceil(reBoundingBox(dimes(1))):floor(reBoundingBox(dimes(1))+reBoundingBox(dimes(1)+3)),ceil(reBoundingBox(dimes(2))):floor(reBoundingBox(dimes(2))+(reBoundingBox(dimes(2)+3))));
ccr1=normxcorr2(flipdim(windowt,1),Ct);
ccr2=normxcorr2(flipdim(windowt,2),Ct);
ccoerf1=max(ccr1(:));
ccoerf2=max(ccr2(:));
if ccoerf1>ccoerf2
canddim=1;
else
canddim=2;
end
dim1=[false false false];
dim1(dimes(canddim))=true;
%% Encuentra el eje restante (coincidirá con la linea coronilla-cuello)
%%
dim3=~or(dim1,dim2);
% Perumta la imagen de modo que la direccion Z sea la 1 la oreja-oreja la 3
ori_vect(3)=posdim(dim1);
ori_vect(2)=posdim(dim2);
ori_vect(1)=posdim(dim3);
fprintf('.')
ov=ori_vect;
if ~all(ori_vect==posdim)
imgs_op=permute(imgs,ori_vect);
BWN=permute(BW,ori_vect);
else
imgs_op=imgs;
BWN=BW;
end
% ORIENTACION AXIAL: Determina si la orientacion del eje frente-cogote es correcta
BWS2=squeeze(sum(BWN,1));
s2=sum(BWS2,2);vari_horz=s2(s2>0);
N=5;
x=-numel(vari_horz):numel(vari_horz); y=gaussmf(x,[numel(vari_horz)/N 0]);
corcruzf=xcorr(flipud(vari_horz),y);
corcruz =xcorr(vari_horz,y);
mpic=max(corcruz);
[~,maxdif]=max(corcruz(corcruz>mpic/100)-corcruzf(corcruzf>mpic/100));%gradient(xcorr(flipud(vari_horz),vari_horz)));%max(vari_horz-flipud(vari_horz));
[~,mindif]=min(corcruz(corcruz>mpic/100)-corcruzf(corcruzf>mpic/100));%gradient(xcorr(flipud(vari_horz),vari_horz)));%min(vari_horz-flipud(vari_horz));
BWS2p=squeeze(sum(imgs_op,1));
s2o=max(BWS2p,[],2);
ax= find(s2==max(s2),1)-find(s2o==max(s2o),1); %diferencia basada en el maximo central comparado con el máximo en los estriados
if tiporf == 303 %SPECT DaTSCAN
ocheckhorz=(ax<=0);
else
ocheckhorz=((maxdif-mindif)<=0);
end
% ORIENTACION SAGITAL: Determina la orientacion del eje coronilla-cuello
for i=ceil(reBoundingBox(find(dim3))):floor(reBoundingBox(find(dim3)))+reBoundingBox(find(dim3)+3)-1
gradient_vert(i,:,:)=BWN(i+1,:,:)-BWN(i,:,:);
vari_vert(i)=var(gradient_vert(i,abs(gradient_vert(i,:))>0));
end
cm_sum_mean_vert= (mean(cumsum(fliplr(vari_vert)))>=mean(cumsum(vari_vert)));%(mean(find(area==max(area)))>(numel(area)/2));
ocheckvert=(~cm_sum_mean_vert);
if ocheckvert % si es positivo quiere decir que la linea es coronilla-cuello
imgs_op=flipdim(imgs_op,1);
imgs_op=flipdim(imgs_op,3);
ov(2)=-ov(2);
end
if ocheckhorz
imgs_op=flipdim(imgs_op,2);
imgs_op=flipdim(imgs_op,3);
ov(3)=-ov(3);
end
imgs_o=imgs_op;
fprintf('. \n')
end