-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpsf_in_fdr_space.m
153 lines (136 loc) · 5.5 KB
/
psf_in_fdr_space.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
function [h,hinv] = psf_in_fdr_space(sizes, steps, starts, zoomdata, psf, psf_starts,psf_steps,psf_sizes,zoompsf);
% compute the psf in FDR space where the latter has coordinates of spatial frequency and angular frequency
% sizes : (Nx, Nphi) dimensions of image
% steps : (sx, sphi) voxel size in x and l dimensions
% starts : (cx, cphi) centre of 0,0 voxel
% psf : ArrayVolume specifying psf
% alloc array
h = zeros(sizes(1), sizes(2));
zero_offset = size(psf,2)/2;
% scale the zoomed psf measurement according to eq 2 in the patent
% guessing values for the various parameters
lamda=473e-9;
nbath=1.33;
e=7.4e-6;
% new custom-built opt
zoomedNA=[0.501 6.328e-3;
0.907 0.011;
1.0 0.012;
1.361 0.017;
2.038 0.025;
2.871 0.036;
5.101 0.062;
5.5 0.067];
% former bioptonics
%zoomedNA = [ 0.80 0.0125; 1.00 0.0155;
% 1.30 0.0180;
% 1.60 0.0225;
% 2.00 0.0270;
% 2.50 0.0320;
% 3.20 0.0390;
% 4.00 0.0465;
% 5.00 0.0505;
% 6.30 0.0620;
% 8.00 0.0625;
% 10.0 0.0625 ];
NApsf=interp1(zoomedNA(:,1),zoomedNA(:,2),zoompsf);
NAdata=interp1(zoomedNA(:,1),zoomedNA(:,2),zoomdata);
DOFpsf=(nbath*(lamda/NApsf^2+e/zoompsf/NApsf));
DOFdata=(nbath*(lamda/NAdata^2+e/zoomdata/NAdata));
lscale=DOFdata/DOFpsf;
% since nobody is providing the position of the focal plane either, guess
% a value as 1/2 the DOF of the experiment, converting from metres to mm
pfp=DOFdata/2 * 1000;
% from the shape of the bowtie, ga sample data seem to have pfp near zero
%pfp=0;
% compute fft of psf along detector element dimension (x)
Rx_psf = fftshift(fft(fftshift((psf(:,:)),2),[],2) ,2);
% calculate R (Fourier space) coord
Rx_values = ([-zero_offset+1: ...
psf_sizes(2)-zero_offset]+0.5)/(psf_steps(2)*psf_sizes(2));
Rz_values=Rx_values;
% calculate the l coord, scale it to the brain dataset, and offset by the pfp
l_values = (psf_starts(3) + [1:psf_sizes(3)]*psf_steps(3))*lscale + pfp;
figure(1),subplot(1,3,2),imagesc(Rx_values', l_values',abs(Rx_psf),[0 2]);
xlabel('R (1/mm)');
ylabel('distance from focal plane (mm)');
title('psf in l,R space');
% generate matrices from the coordinate arrays for bilinear interp
[l_values,Rx_values]=meshgrid(l_values,Rx_values);
do3d=0;
if do3d
Rx_values = ([-zero_offset+1: ...
psf_sizes(2)-zero_offset]+0.5)/(psf_steps(2)*psf_sizes(2)); ...
x_r=1:4:psf_sizes(2);
Rx_values=Rx_values(x_r);
Rz_values=Rx_values;
l_values = (psf_starts(3) + [1:psf_sizes(3)]*psf_steps(3))*lscale + ...
pfp;
l_r=[1:10:floor(psf_sizes(3))];
l_values=l_values(l_r);
[l_grid,Rx_grid,Rz_grid]=meshgrid(l_values,Rx_values,Rz_values);
psf3d_grid=psf3d(l_r,x_r,x_r);
psf3d_grid=permute(psf3d_grid,[2 1 3]);
end
for j=1:sizes(2) % for each angular frequency
for i=1:sizes(1) % for each detector spatial frequency
% these coords are calculated in Fourier space, not sinogram space
% it has to be in Fourier space in order to get units of
% mm for l.
Phi(i,j) = ( - sizes(2)/2+j)/(steps(2)*sizes(2));
Rx(i,j) = (-sizes(1)/2+i+0.5)/(steps(1)*sizes(1));
l(i,j) = -Phi(i,j)/Rx(i,j);
end
end
lmax=DOFdata*1000*1.5;
lmin=-DOFdata*1000/1.5;
h = interp2(l_values, Rx_values,Rx_psf.',l, Rx,'linear');
if do3d
% Rx=[-sizes(1)/2+.5:sizes(1)/2]/(steps(1)*sizes(1));
clear Rx Rz Phi;
Rx=[-sizes(1)/2+.5:10:sizes(1)/2]/(steps(1)*sizes(1));
Rz=[-sizes(1)/2+.5:20:sizes(1)/2]/(steps(1)*sizes(1));
Phi=[(-sizes(2)/2:10:sizes(2)/2)/(steps(2)*sizes(2))];
% Rz=Rx;
% Phi0=4;
% px=-Phi0./([-sizes(1)/2+.5:10:sizes(1)/2]/(steps(1)*sizes(1)));
[Phi,Rx,Rz]=meshgrid(Phi,Rx,Rz);
lx=-Phi./Rx;
% interp2 literally won't work with ndgrid, only meshgrid
h3 = interp3(l_grid,Rx_grid,Rz_grid,single(psf3d_grid),lx,Rx,Rz,'linear');
figure(6),imagesc(sq(abs(h3(:,61,:))));
delete(['fdr.raw']);
fid=fopen('fdr.raw','wb','l');
count=fwrite(fid,abs(h3(:)),'float32');
fclose(fid);
unix(['rawtominc -input fdr.raw -float -clobber -xyz -xstart ' ...
num2str(0) ' -xstep ' num2str(steps(1)) ' -ystart ' ...
num2str(starts(1)) ' -ystep ' num2str(steps(2)) ' -zstart ' ...
num2str(starts(2)) ' -zstep ' num2str(steps(3)) ' fdr.mnc ' ...
num2str(size(h3,3)) ' ' num2str(size(h3,2)) ' ' num2str(size(h3,1)) ...
]);
end
hinv = 1./h;
% a filter to smooth off the sharp edges of the nullspace at the lmax/lmin
% cutoff points, while also preserving a small part of the nullspace
% (ie the "vee") very close to the origin where there looks to be more
% than just noise, since we also know the FDR is only an approximation
% to begin with.
% some hardcoded cutoff values for defining the region of origin preservation
% need to parameterize these later.
Pmin=-1; % 1/rad
Rmin=2; % 1/mm
% generate the cutoff line P(Rx) as a piecewise constant line in Phi,R space,
% and smooth the vertices a little bit
clear P;
for i=1:sizes(1)
P(i)=-Rx(i)*lmax ./(1+exp(-(Rx(i)-Rmin))) + Pmin / (1+exp(-(Rmin-Rx(i)))) + Pmin / (1+exp(Rx(i)+Rmin)) + -Rx(i)*lmin ./ (1+exp((Rx(i)+Rmin)));
end
% map that line to a Phi-based sigmoidal roll off function
P=repmat(P,[sizes(2) 1]);
sg=1./(1+exp(-(Phi-P')));
% generate correct values for the other half of Phi,R space by a
% reflection symmetry argument.
sg(:,end/2+1:end)=flipud(fliplr(sg(:,1:end/2)));
h = h.* sg ;
hinv = hinv .* sg;