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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
|
function [xGridMinutes,yGridMinutes,psf] = OtfToPsf(xSfGridCyclesDeg,ySfGridCyclesDeg,otf,varargin)
%OTFTOPSF Convert a 2D optical transfer fucntion to a 2D point spread function.
% [xSfGridCyclesDeg,ySfGridCyclesDeg,otf] = PsfToOtf([xGridMinutes,yGridMinutes],psf)
%
% Converts a optical transfer function specified over two-dimensional
% spatial frequency in cycles per degree to a point spread function
% specified over positions in minutes. For human vision, these are each
% natural units.
%
% The input spatial frequencies should be specified in matlab's grid
% matrix format and sf for x and y should be specified over the same
% spatialfrequency extent and with the same number of evenly spaced
% samples. Spatial frequency (0,0) should be at location floor(n/2)+1 in
% each dimension. The PSF is returned with position (0,0) at location
% floor(n/2)+1 in each dimension. This is the form we mostly want to
% look at and use.
%
% The OTF is assummed to have the DC term in the center poistion,
% floor(n/2)+1 of the passed matrix. If you have your hands on an OTF
% with the DC term in the upper left position (1,1), apply fftshift to
% it before passing to this routine. The DC in upper left is the Matlab
% native format for applying the ifft, and is also the format stored by
% isetbio in its optics structure.
%
% Positions are returned using the same conventions.
%
% No normalization is performed. The psf should be real, and we
% complain (throw an error) if it is not, to reasonable numerial
% precision. If it seems OK, we make sure it is real.
%
% We also make sure that the returned psf is all postive and sums to
% 1. In some cases, we found that there were small negative values
% and after setting these to zero renormalization was needed.
%
% We wrote this rather than simply relying on Matlab otf2psf/psf2otf
% because we want a routine that deals with the conversion of spatial
% frequency to spatial support.
%
% If you pass the both sf args as empty, both position grids are
% returned as empty and just the conversion on the OTF is performed.
%
% PsychOpticsTest shows that this works very well when we go back and
% forth for diffraction limited OTF/PSF. But not exactly exactly
% perfectly. A signal processing maven might be able to track down
% whether this is just a numerical thing or whether some is some small
% error, for example in how position is converted to sf or back again in
% the PsfToOtf.
%
% Optional key/value pairs:
% 'warningInsteadOfErrorForNegativeValuedPSF' - Set to 1 (default
% 0) to get a warning
% not an error if the psf
% values are too negative
% before they are forced
% to zero. Set to 2 for
% no warning msg or
% error.
% 'negativeFractionalTolerance' - The error/warning is
% thrown if the magnitude
% of the most negative
% value is more than this
% fraction of the maximum
% (positve) value.
% (Default 1e-3).
%
% See also: PsfToOtf, PsychOpticsTest.
% History:
% 01/26/18 dhb
% 03/31/18 dhb Document key/value pair added by someone else.
% dhb Add key/value pair for negative value tolerance.
% This is now 1e-3 rather than 1e-10
% 10/17/18 dhb Now 10-2.
%% Parse input
p = inputParser;
p.addParameter('warningInsteadOfErrorForNegativeValuedPSF', 0, @isnumeric);
p.addParameter('negativeFractionalTolerance', 1e-2, @isnumeric);
p.parse(varargin{:});
%% Handle sf args and converstion
if (~isempty(xSfGridCyclesDeg) & ~isempty(ySfGridCyclesDeg))
% They can both be passed as non-empty, in which case we do a set of sanity
% checks and then do the conversion.
[m,n] = size(xSfGridCyclesDeg);
centerPosition = floor(n/2) + 1;
if (m ~= n)
error('psf must be passed on a square array');
end
[m1,n1] = size(ySfGridCyclesDeg);
if (m1 ~= m || n1 ~= n)
error('x and y positions are not consistent');
end
[m2,n2] = size(otf);
if (m2 ~= m || n2 ~= n)
error('x and y positions are not consistent');
end
if (~all(xSfGridCyclesDeg(:,centerPosition) == 0))
error('Zero spatial frequency is not in right place in the passed xGrid');
end
if (~all(ySfGridCyclesDeg(centerPosition,:) == 0))
error('Zero spatial frequency is not in right place in the passed yGrid');
end
if (xSfGridCyclesDeg(1,centerPosition) ~= ySfGridCyclesDeg(centerPosition,1))
error('Spatial frequency extent of x and y grids does not match');
end
diffX = diff(xSfGridCyclesDeg(:,centerPosition));
if (any(diffX ~= diffX(1)))
error('X positions not evenly spaced');
end
diffY = diff(ySfGridCyclesDeg(centerPosition,:));
if (any(diffY ~= diffY(1)))
error('Y positions not evenly spaced');
end
if (diffX(1) ~= diffY(1))
error('Spatial frequency sampling in x and y not matched');e
end
%% Generate position grids
%
% Samples are evenly spaced and the same for both x and y (checked above).
% Handle even versus odd dimension properly for fft conventions.
[xGridMinutes,yGridMinutes] = SfGridCyclesDegToPositionGridMinutes(xSfGridCyclesDeg,ySfGridCyclesDeg);
elseif (isempty(xSfGridCyclesDeg) & isempty(ySfGridCyclesDeg))
% This case is OK, we set the output grids to empty
xGridMinutes = [];
yGridMinutes = [];
else
% This case is not allowable
error('Either both sf grids must be empty, or neither');
end
%% Compute otf and put spatial center in the middle of the grid
psf = fftshift(ifft2(ifftshift(otf)));
%% See if there is stray imaginary stuff, get rid of it if so.
%
% Throw an error if the returned psf isn't in essence real valued.
% Then set residual imaginary parts to 0.
if (any(abs(imag(psf(:))) > 1e-10))
error('Computed psf is not sufficiently real');
end
if (any(imag(psf(:))) ~= 0)
psf = psf - imag(psf)*1i;
end
% Check for large negative psf values, and then set any small
% negative values to zero. For some cases (e.g. Marimont-Wandell
% OTF, we do get negative values because the way that was constructed
% is an approximation to measurements that does not absolutely guarantee an
% all positive OTF.
if (max(psf(:)) <= 0)
error('Computed PSF has no positive values. This is not good.');
end
if (min(psf(:)) < 0 && abs(min(psf(:))) > p.Results.negativeFractionalTolerance*max(psf(:)))
if (p.Results.warningInsteadOfErrorForNegativeValuedPSF == 1)
fprintf(2,'Mysteriously large negative psf values, min value is %g, relative to max of %g, fraction %g\n',min(psf(:)),max(psf(:)),abs(min(psf(:)))/max(psf(:)));
elseif (p.Results.warningInsteadOfErrorForNegativeValuedPSF == 0)
fprintf(2,'Mysteriously large negative psf values: min value is %g, relative to max of %g, fraction %g\n',min(psf(:)),max(psf(:)),abs(min(psf(:)))/max(psf(:)));
error('Mysteriously large negative psf values');
end
end
psf(psf < 0) = 0;
% Make sure return is real
psf = abs(psf);
% Make sure return sums to 1. It might not be because of the
% above fussing with imaginary and negative values.
psf = psf/sum(psf(:));
end
|
