115 lines
2.6 KiB
Matlab
115 lines
2.6 KiB
Matlab
%%PHASE UNWRAP
|
|
%unwraps given data with three different unwrapping methods
|
|
%compares quantitatively to OMME results
|
|
%%
|
|
clear
|
|
clc
|
|
|
|
%these are from the repository for the Loecher 2014 paper
|
|
addpath('4dflow-lapunwrap/utils');
|
|
addpath('4dflow-lapunwrap/unwrap');
|
|
%%
|
|
struct = load('path_to_data_file');
|
|
wrapped = struct.pca_p;
|
|
mask_struct = load('path_to_mask_file');
|
|
mask = mask_struct.labels;
|
|
si = size(wrapped);
|
|
si(3) = 15;
|
|
phi = zeros(si(1), si(2), 1, si(3));
|
|
noise_std_percent = 0.15;
|
|
noise_std = noise_std_percent *pi;
|
|
rng('default')
|
|
rng(100 + file{1}(2))
|
|
noise = noise_std.*randn(si);
|
|
phi(:, :, 1, 1:si(3)) = (wrapped(:, :, 1:si(3)) + noise).*mask(:, :, 1:si(3));
|
|
phi_wn = zeros(si(1), si(2), 1, si(3));
|
|
phi_wn(:, :, 1, 1:si(3)) = (wrapped(:, :, 1:si(3))).*mask(:, :, 1:si(3));
|
|
phi = phi - 2*pi*(phi > pi);
|
|
phi = phi + 2*pi*(phi < -pi);
|
|
masked = zeros(si(1), si(2), si(3));
|
|
masked(:, :, :) = phi(:, :, 1, :);
|
|
|
|
%%
|
|
%laplacian unwrap
|
|
tic
|
|
n_u4 = unwrap_4D(phi);
|
|
toc
|
|
lap4d_phi = phi + double(n_u4).*2*pi;
|
|
n_u3 = zeros(size(n_u4), 'int8');
|
|
tic
|
|
for i=1:si(3)
|
|
n_u3(:, :, :, i) = unwrap_3D(phi(:, :, :, i));
|
|
end
|
|
lap3d_phi = phi + double(n_u3).*2*pi;
|
|
toc
|
|
|
|
%%
|
|
%temporal unwrap
|
|
t_ref = 15;
|
|
tic
|
|
temp_phi = temporal(phi, t_ref);
|
|
toc
|
|
|
|
%%
|
|
%probability unwrap
|
|
w_s = 1.0;
|
|
w_t = 2.5;
|
|
t_low = 0.32;
|
|
t_high = 0.75;
|
|
tic
|
|
prob_phi = probability(phi, t_low, t_high, w_s, w_t);
|
|
toc
|
|
|
|
%%
|
|
%omme
|
|
struct = load('path_to_higher_venc_file');
|
|
high_data = struct.pca_p;
|
|
high_venc = zeros(size(phi));
|
|
high_venc(:, :, 1, :) = (high_data(:, :, 1:si(3))).*mask(:, :, 1:si(3));
|
|
highv = 150;
|
|
lowv = 75;
|
|
high_venc = high_venc.*(highv/pi);
|
|
low_venc = phi.*(lowv/pi);
|
|
tic
|
|
[omme_v, dyn, conf, std] = omme([high_venc(:), low_venc(:)], [highv, lowv]);
|
|
omme_v = reshape(omme_v, [si(1), si(2), 1, si(3)]);
|
|
toc
|
|
|
|
%%
|
|
%evaluate!
|
|
temp_v = temp_phi.*lowv/pi;
|
|
lap3d_v = lap3d_phi.*lowv/pi;
|
|
lap4d_v = lap4d_phi.*lowv/pi;
|
|
prob_v = prob_phi.*lowv/pi;
|
|
v = phi.*lowv/pi;
|
|
|
|
t = 1:si(3);
|
|
base_error = eval_error(v(:, :, 1, t), omme_v(:, :, 1, t))
|
|
t_error = eval_error(temp_v(:, :, 1, t), omme_v(:, :, 1, t))
|
|
l4_error = eval_error(lap4d_v(:, :, 1, t), omme_v(:, :, 1, t))
|
|
l3_error = eval_error(lap3d_v(:, :, 1, t), omme_v(:, :, 1, t))
|
|
prob_error = eval_error(prob_v(:, :, 1, t), omme_v(:, :, 1, t))
|
|
|
|
%%
|
|
%display!
|
|
t = 4;
|
|
slice = 1;
|
|
h = 145:180;
|
|
w = 35:80;
|
|
im0 = v(h, w, slice, t);
|
|
im1 = temp_v(h, w, slice, t);
|
|
im2 = lap4d_v(h, w, slice, t);
|
|
im4 = prob_v(h, w, slice, t);
|
|
im3 = omme_v(h, w, slice, t);
|
|
im5 = lap3d_v(h, w, slice, t);
|
|
figure();
|
|
|
|
imshow([[im0 im3 im1], [im5, im2 im4]], [], 'border', 'tight', 'InitialMagnification', 300)
|
|
axis off
|
|
colormap("jet")
|
|
a = colorbar
|
|
ylabel(a,'velocity (cm/s)')
|
|
caxis([-160, 220])
|
|
axis off
|
|
|