Skip to content
Snippets Groups Projects
Commit 47c36a90 authored by Nicola Vigano's avatar Nicola Vigano
Browse files

PSF: added direct-space convolution, OTF estimation from ESRF measurements...

PSF: added direct-space convolution, OTF estimation from ESRF measurements (imported), and some small fixes

Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent a860d1a4
No related branches found
No related tags found
No related merge requests found
function mt = mncf(opt, varargin)
% function mt = mncf(opt, varargin)
%
% Input:
% opt : camera-setup
% 1: july 97, 10um set-up
% 2: april 98, 1.93um pixel set-up: optique 10x + 25um YAG
% 3: august 98, 0.932um pixel set-up: optique 10x + eyepiece + 25um YAG
% 4: september 97, 0.7 um pixel set-up: 20* YAG 5um
% 5: november 98, 0.5 um pixel set-up: optique 20x + eyepiece + LAG 3um
% 6: november 99, 0.394 um pixe set-up: optique 20x + eyepiece + LAG 6mu
% 7: november 2001, 1 micron setup with 0.69um pixelsize (2048 ccd): 10x optics + eyepiece + 25um YAG
% 8: for testing
% 9: for fits
% Optional Input:
% m : number of row samples
% n : number of column samples
% oversampling : samplerate is fixed by oversamp, sigmas expressed in pixels,
% changing n changes the width of the window and sampling rate
% in frequency space
% dimmensions: 1 | {2} create 1-dim mtf plot or 2-dim mtf map
conf = struct( ...
'size', 11, ...
'oversampling', 1, ...
'dimensions', 2, ...
'p', []);
conf = parse_pv_pairs(conf, varargin);
switch (opt)
case 1
% july 97, 10 mu set-up
mi = 0.01;
ma = 1;
p = [0.5884, 0.9093, 0.3570, 2.8537, 0.0352, 49.3503, 748.8534];
case 2
% april 98, 1.93 um pixel set-up: optique 10x + 25um YAG
mi = 0.0117;
ma = 0.9862;
p = [0.4327, 0.9147, 0.2296, 7.9874, 0.1683, 54.6962, 197.9522];
case 3
% august 98, 0.932 mu pixel set-up: optique 10x + eyepiece 2x + 25um YAG
%mi = 0.0149;
%ma = 0.996;
%p = [0.4723 1.3329 0.2154 5.8018 0.1569 56.1836 268.3804];
%best 100 columns:
%mi = 0.0151;
%ma = 0.994;
%p = [0.4530 1.2271 0.2279 5.2286 0.1550 52.5938 255.9305];
%measurement in february 99 - optimised uptilt of GaAs
mi = 0.013;
ma = 0.984;
p = [0.5985, 0.9903, 0.1622, 8.4452, 0.1423, 41.2661, 375.8298];
case 4
% september 97 0.7 um pixel set-up: optique 20x + eyepiece 1.6x + YAG 5um
mi = 0.0214;
ma = 0.9841;
p = [0.3396, 1.2313, 0.2130, 10.4607, 0.1801, 51.0301, 240.2380];
case 5
% november 98 0.5 um pixel set-up: optique 20x + eyepiece 2x + LAG 3um
mi = 0.022;
ma = 1.001;
p = [0.4182, 1.6436, 0.2293, 12.4082, 0.1922, 288.6455, 67.6439];
case 6
% november 99 0.394 um pixel set-up: optique 20x + eyepiece 2x + LAG 6um
mi = 0.005;
ma = 0.996;
p = [0.6023, 1.6800, 0.2545, 5.0455, 0.1051, 29.9829, 235.2297];
case 7
% november 2001, 1 micron setup with 0.69um pixelsize (2048 ccd): 10x optic
% determined with jofit and a defocus series
mi = 0;
ma = 1;
p = [0.3978, 1.0493, 0.4608, 3.4978, 0.0674, 17.8457, 679.3300];
case 8
% for testing
mi = 0;
ma = 1;
p = [0.3246, 1.1299, 0.4810, 3.4915, 0.1187, 15.4062, 729.8940];
case 9
% for fits
mi = 0;
ma = 1;
p = conf.p;
case 10
mi = 0;
ma = 1;
p = [0.5305, 1.5617, 0.4765, 7.4535, -0.0676, 15.7088, 683.5288];
end
% a's are the widths in erfc = sqrt(2)*sigma
a1 = p(2) * conf.oversampling;
a2 = p(4) * conf.oversampling;
a3 = p(6) * conf.oversampling;
a4 = p(7) * conf.oversampling;
% b's are weigths
b1 = p(1);
b2 = p(3);
b3 = p(5);
b4 = ma - mi - p(1) - p(3) - p(5);
n = conf.size;
m = conf.size;
switch (conf.dimensions)
case 1
f = (0:1/n:0.5)';
mt = b1 * exp(-(pi * a1 * f) .^ 2) + b2 * exp(-(pi * a2 * f) .^ 2) + b3 * exp(-(pi * a3 * f) .^2 ) + b4 * exp(-(pi * a4 * f) .^ 2);
mt = [mt; flipud(mt(2:n/2))];
mt(1) = 1;
case 2
[x, y] = meshgrid(0:1/n:0.5, 0:1/m:0.5);
f2 = x .^ 2 + y .^ 2;
mt = b1 * exp(-(pi * a1) ^ 2 * f2) + b2 * exp(-(pi * a2) ^ 2 * f2) + b3 * exp(-(pi * a3) ^ 2 * f2) + b4 * exp(-(pi * a4) ^ 2 * f2);
mt = [ ...
mt, fliplr(mt(:, 2:floor(n/2))); ...
flipud(mt(2:floor(m/2), :)), rot90(mt(2:floor(m/2), 2:floor(n/2)),2)];
mt(1, 1) = 1;
otherwise
error([mfilename ':wrong_argument'], 'Dimemsions must be 1 or 2');
end
end
......@@ -10,10 +10,16 @@ classdef GtPSF < handle
is_projection_invariant = true;
image_size = [];
use_otf = true;
data_type = 'single';
end
methods (Access = public)
function self = GtPSF(varargin)
self = parse_pv_pairs(self, varargin);
end
function set_psf_direct(self, psf_d, image_size)
self.reset();
......@@ -33,7 +39,7 @@ classdef GtPSF < handle
psfs_norm = 1 ./ sum(sum(abs(psf_d), 1), 3);
psf_d = bsxfun(@times, psf_d, psfs_norm);
self.psf_direct = psf_d;
self.psf_direct = cast(psf_d, self.data_type);
if (exist('image_size', 'var') && ~isempty(image_size))
self.image_size = image_size;
end
......@@ -48,7 +54,7 @@ classdef GtPSF < handle
self.is_projection_invariant = size(psf_d, 2) == 1;
% if we know the images sizes, we can already compute the OTF
if (~isempty(image_size))
if (~isempty(self.image_size))
self.init_otfs();
end
end
......@@ -56,13 +62,21 @@ classdef GtPSF < handle
function imgs = apply_psf_direct(self, imgs)
self.check_incoming_images(imgs);
imgs = self.apply_otf(imgs, true);
if (self.use_otf)
imgs = self.apply_otf(imgs, true);
else
imgs = self.apply_psf(imgs, true);
end
end
function imgs = apply_psf_transpose(self, imgs)
self.check_incoming_images(imgs);
imgs = self.apply_otf(imgs, false);
if (self.use_otf)
imgs = self.apply_otf(imgs, false);
else
imgs = self.apply_psf(imgs, false);
end
end
end
......@@ -132,7 +146,7 @@ classdef GtPSF < handle
self.image_size = imgs_size;
end
if (isempty(self.otf_direct))
if (isempty(self.otf_direct) && self.use_otf)
self.init_otfs();
end
end
......@@ -164,5 +178,21 @@ classdef GtPSF < handle
imgs = real(imgs);
end
function imgs = apply_psf(self, imgs, is_direct)
if (is_direct || self.is_symmetric)
psf = self.psf_direct;
else
psf = self.psf_transpose;
end
if (self.is_projection_invariant)
imgs = convn(imgs, psf, 'same');
else
for ii_i = 1:size(imgs, 2)
imgs(:, ii_i, :) = convn(imgs(:, ii_i, :), psf(:, ii_i, :), 'same');
end
end
end
end
end
\ No newline at end of file
function [otf, psf] = gtOTFGetMeasured(opt, varargin)
% function mt = gtOTFGetMeasured(opt, varargin)
%
% Input:
% opt : camera-setup
% 1: july 97, 10um set-up
% 2: april 98, 1.93um pixel set-up: optique 10x + 25um YAG
% 3: august 98, 0.932um pixel set-up: optique 10x + eyepiece + 25um YAG
% 4: september 97, 0.7 um pixel set-up: 20* YAG 5um
% 5: november 98, 0.5 um pixel set-up: optique 20x + eyepiece + LAG 3um
% 6: november 99, 0.394 um pixe set-up: optique 20x + eyepiece + LAG 6mu
% 7: november 2001, 1 micron setup with 0.69um pixelsize (2048 ccd): 10x optics + eyepiece + 25um YAG
% 8: for testing
% 9: for fits
% Optional Input:
% m : number of row samples
% n : number of column samples
% oversampling : samplerate is fixed by oversamp, sigmas expressed in pixels,
% changing n changes the width of the window and sampling rate
% in frequency space
%
% Derived from a modified version of the function: mncf from
% Wolfgang/Peter's matlab functions
conf = struct( ...
'size', 11, ...
'oversampling', 1, ...
'p', []);
conf = parse_pv_pairs(conf, varargin);
% if (mod(conf.size, 2) ~= 1)
% error([mfilename ':wrong_argument'], ...
% 'The OTF should have an odd size');
% end
switch (opt)
case 1
% july 97, 10 mu set-up
mi = 0.01;
ma = 1;
p = [0.5884, 0.9093, 0.3570, 2.8537, 0.0352, 49.3503, 748.8534];
case 2
% april 98, 1.93 um pixel set-up: optique 10x + 25um YAG
mi = 0.0117;
ma = 0.9862;
p = [0.4327, 0.9147, 0.2296, 7.9874, 0.1683, 54.6962, 197.9522];
case 3
% august 98, 0.932 mu pixel set-up: optique 10x + eyepiece 2x + 25um YAG
%mi = 0.0149;
%ma = 0.996;
%p = [0.4723 1.3329 0.2154 5.8018 0.1569 56.1836 268.3804];
%best 100 columns:
%mi = 0.0151;
%ma = 0.994;
%p = [0.4530 1.2271 0.2279 5.2286 0.1550 52.5938 255.9305];
%measurement in february 99 - optimised uptilt of GaAs
mi = 0.013;
ma = 0.984;
p = [0.5985, 0.9903, 0.1622, 8.4452, 0.1423, 41.2661, 375.8298];
case 4
% september 97 0.7 um pixel set-up: optique 20x + eyepiece 1.6x + YAG 5um
mi = 0.0214;
ma = 0.9841;
p = [0.3396, 1.2313, 0.2130, 10.4607, 0.1801, 51.0301, 240.2380];
case 5
% november 98 0.5 um pixel set-up: optique 20x + eyepiece 2x + LAG 3um
mi = 0.022;
ma = 1.001;
p = [0.4182, 1.6436, 0.2293, 12.4082, 0.1922, 288.6455, 67.6439];
case 6
% november 99 0.394 um pixel set-up: optique 20x + eyepiece 2x + LAG 6um
mi = 0.005;
ma = 0.996;
p = [0.6023, 1.6800, 0.2545, 5.0455, 0.1051, 29.9829, 235.2297];
case 7
% november 2001, 1 micron setup with 0.69um pixelsize (2048 ccd): 10x optic
% determined with jofit and a defocus series
mi = 0;
ma = 1;
p = [0.3978, 1.0493, 0.4608, 3.4978, 0.0674, 17.8457, 679.3300];
case 8
% for testing
mi = 0;
ma = 1;
p = [0.3246, 1.1299, 0.4810, 3.4915, 0.1187, 15.4062, 729.8940];
case 9
% for fits
mi = 0;
ma = 1;
p = conf.p;
case 10
mi = 0;
ma = 1;
p = [0.5305, 1.5617, 0.4765, 7.4535, -0.0676, 15.7088, 683.5288];
otherwise
error([mfilename ':wrong_argument'], ...
'Not recognized option: %d', opt)
end
% a's are the widths in erfc = sqrt(2)*sigma
a = (pi .* p([2 4 6 7]) .* conf.oversampling) .^ 2;
a = reshape(a, 1, 1, []);
% b's are weigths
b = [p([1 3 5]), (ma - mi - sum(p([1 3 5])))];
b = reshape(b, 1, 1, []);
grid_points = linspace(-0.5, 0.5, conf.size);
[x, y] = meshgrid(grid_points, grid_points);
otf = x .^ 2 + y .^ 2;
otf = bsxfun(@times, otf, a);
otf = exp(-otf);
otf = sum(bsxfun(@times, otf, b), 3);
center_pos = floor(conf.size / 2) + 1;
otf(center_pos, center_pos) = 1;
otf = ifftshift(otf);
psf = real(ifft2(otf));
psf = fftshift(psf);
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment