-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathInteractiveZalignment.m
155 lines (137 loc) · 4.22 KB
/
InteractiveZalignment.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
function [pxRatioZ,cent] = InteractiveZalignment(cellType,tifpath,varargin)
% 5/11/20 by RT
% Interactively compare projected skeleton image with intensity discarded
% depth depth image to find Z alignment parameter
% This does not do re-projection of skeleton so camera parameter must be
% fixed first
%% Default parameters
% camera parameters
cameraPosition = [26000,36000,22000]; % looks at the target along the native depth axis
cameraTarget = [26000,26000,21000]; % Center of the ellipsoid body, based on R4m
cameraUpVector = [0,1,0];
% projection to image conversion parameters
defPxRatioZ = 8/1000;
defCenterZ = 80; % this depends on cameraTarget, obviously
XYPxRatio = [1,1] * 8/620;
XYCenter = [605,191]; % this depends on cameraTarget, obviously
if ismac
keySet = [41,79:82,225];
elseif IsWin
keySet = [0,0,0,0,0,0];
elseif IsLinux
keySet = [0,0,0,0,0,0];
else
error('Unknown OS');
end
% register optional arguments
for ii = 1:2:length(varargin)
eval([varargin{ii} '= varargin{' num2str(ii+1) '};']);
end
%% get MIP
disp('Loading tif image...');
dMap = CDMIPtoDepthMap(tifpath);
dMapSize = size(dMap);
%% get projection
disp('Projecting SWC...');
SWC = showswcnodesByCellType(cellType,0);
[rotXYZ,R] = projectSWC(SWC,cameraPosition,cameraTarget,cameraUpVector,0);
rotXYZ(:,2) = -rotXYZ(:,2);
projdMap = XYZtodMap(rotXYZ, R, XYPxRatio, XYCenter, dMapSize);
% Show target depth Map
figure;
% Update Z scaling/positioning, show XZ slice
slice = round(dMapSize(1)/2);
pxRatioZ = defPxRatioZ;
cent = defCenterZ;
while 1
% Show slice location
s1 = subplot(2,2,1);
cla(s1);
hold(s1, 'on');
imagesc(dMap);
colormap(s1,'gray');
caxis([1,255]);
plot([1,dMapSize(2)],slice*[1,1],'r--');
axis ij tight
xlim([1,dMapSize(2)]); ylim([1,dMapSize(1)]);
colorbar;
hold(s1, 'off');
% Show depth map
s2 = subplot(2,2,3);
cla(s2)
hold(s2, 'on');
projImg = projdMap*pxRatioZ + cent;
imagesc(projImg);
colormap(s2,'gray');
caxis([1,255]);
plot([1,dMapSize(2)],slice*[1,1],'r--');
axis ij tight
xlim([1,dMapSize(2)]); ylim([1,dMapSize(1)]);
colorbar;
hold(s2, 'off');
% show slices
s3 = subplot(1,2,2);
cla(s3);
title(['pixel conversion = ',num2str(pxRatioZ),' center = ',num2str(cent)]);
targSlice = dMap(slice,:);
moveSlice = projImg(slice,:);
hold(s3, 'on');
scatter(1:dMapSize(2),targSlice,'ko');
scatter(1:dMapSize(2),moveSlice,'r+');
xlim([1,dMapSize(2)]); ylim([1,255]);
xlabel('x (px)'); ylabel('slice');
hold(s3, 'off');
drawnow;
% move paramters
while 1
[~,~,kc] = KbCheck;
if any(kc(keySet(1:5)))
break
end
end
if kc(keySet(1)); break; end % ESC
if kc(keySet(6)) % scaling
if kc(keySet(2))
slice = slice + 5;
elseif kc(keySet(3))
slice = slice - 5;
end
else % translation
if kc(keySet(2))
cent = cent + 1;
elseif kc(keySet(3))
cent = cent - 1;
elseif kc(keySet(4))
pxRatioZ = pxRatioZ+0.0001;
elseif kc(keySet(5))
pxRatioZ = pxRatioZ-0.0001;
end
end
end
end
function out = XYZtodMap(XYZ, R, pxR, cent, imSize)
out = nan(imSize);
% convert projected XYZ to image coordinate
imgX = round(XYZ(:,1)*pxR(1)+cent(1));
imgY = round(XYZ(:,2)*pxR(2)+cent(2));
Z = XYZ(:,3); % use native Z
% ignore nodes that are out of image
outOfImage = imgX>imSize(2) | imgX<1 | imgY>imSize(1) | imgY<1;
imgX = imgX(~outOfImage);
imgY = imgY(~outOfImage);
R = R(~outOfImage);
Z = Z(~outOfImage);
% Only search through pixels with something, using combinatorial unique
% pixel number search with this hacky algorithm
% maybe faster than going through all the pixels?
pxIdx = imgX*1000 + imgY;
uniquePxIdx = unique(pxIdx);
for ii = 1:length(uniquePxIdx)
maxR = max(R(pxIdx == uniquePxIdx(ii)));
thisZ = Z(pxIdx == uniquePxIdx(ii));
maxZ = thisZ(R(pxIdx == uniquePxIdx(ii))==maxR);
thisX = floor(uniquePxIdx(ii)/1000);
thisY = mod(uniquePxIdx(ii),1000);
out(thisY,thisX) = maxZ(1);
end
end