basic-isp/MTF50.m at main · pranav-avn/basic-isp · GitHub

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
% =========================================================================
% Project: Robust Slanted Edge MTF
% Author: Pranav Aravindhan V
% =========================================================================

clc; clear; close all;

img_raw = imread('test_sfr2.jpg');

I = rgb2gray(img_raw);
I = im2double(I);

[h, w] = size(I);
crop_size = 300;
r_start = floor(h/2 - crop_size/2);
c_start = floor(w/2 - crop_size/2);
I = I(r_start:r_start+crop_size, c_start:c_start+crop_size);
[rows, cols] = size(I);

I_smooth = imgaussfilt(I, 2);
edge_centers = zeros(rows, 1);
for r = 1:rows
    row_signal = I_smooth(r, :);
    diff_signal = diff(row_signal);
    [~, max_idx] = max(abs(diff_signal));

    window = 3;
    idx_start = max(1, max_idx - window);
    idx_end = min(length(diff_signal), max_idx + window);

    numerator = 0; denominator = 0;
    for k = idx_start:idx_end
        weight = abs(diff_signal(k));
        numerator = numerator + k * weight;
        denominator = denominator + weight;
    end

    if denominator == 0, edge_centers(r) = max_idx;
    else, edge_centers(r) = numerator / denominator; end
end

row_indices = (1:rows)';
p = polyfit(row_indices, edge_centers, 1);
fitted_edge = polyval(p, row_indices);

oversample_factor = 4;
esf_bins = zeros(cols * oversample_factor, 1);
esf_counts = zeros(cols * oversample_factor, 1);
center_target = cols / 2;

for r = 1:rows
    shift_amount = center_target - fitted_edge(r);
    original_indices = 1:cols;
    shifted_indices = original_indices + shift_amount;
    bin_indices = round(shifted_indices * oversample_factor);

    valid_mask = bin_indices > 0 & bin_indices <= length(esf_bins);

    esf_bins(bin_indices(valid_mask)) = esf_bins(bin_indices(valid_mask)) + I(r, valid_mask)';
    esf_counts(bin_indices(valid_mask)) = esf_counts(bin_indices(valid_mask)) + 1;
end

esf_counts(esf_counts == 0) = 1;
esf_profile = esf_bins ./ esf_counts;

valid_idx = find(esf_profile > 0);
esf_profile = esf_profile(valid_idx(1):valid_idx(end));

if mean(esf_profile(1:10)) > mean(esf_profile(end-10:end))
    esf_profile = 1 - esf_profile;
end

esf_profile = smoothdata(esf_profile, 'gaussian', 10);
lsf_profile = diff(esf_profile);
window = hamming(length(lsf_profile));
lsf_windowed = lsf_profile .* window;

fft_result = fft(lsf_windowed, 4096);
mtf_raw = abs(fft_result);
mtf_norm = mtf_raw / mtf_raw(1);

nyquist = 0.5 * oversample_factor;
freq_axis = linspace(0, nyquist, length(mtf_norm));

% MTF50 Calculation
target = 0.5;
% Find index where it first drops below 0.5
idx_drop = find(mtf_norm < 0.5, 1, 'first');
if isempty(idx_drop)
    mtf50 = 0;
else
    x1 = freq_axis(idx_drop-1); x2 = freq_axis(idx_drop);
    y1 = mtf_norm(idx_drop-1);  y2 = mtf_norm(idx_drop);
    mtf50 = x1 + (target - y1) * (x2 - x1) / (y2 - y1);
end

fprintf('Calculated MTF50: %.4f cycles/pixel\n', mtf50);

figure('Name', 'Robust MTF Analysis');
subplot(1,2,1); imshow(I); title('Cropped ROI');
hold on; plot(fitted_edge, row_indices, 'r');

subplot(1,2,2); plot(freq_axis, mtf_norm, 'k', 'LineWidth', 2);
grid on; xlim([0 1]);
yline(0.5, '--b'); xline(mtf50, '--b');
title(sprintf('MTF50 = %.3f cy/px', mtf50));