Skip to content
Snippets Groups Projects

Resolve "Implement appropriate filter for axial velocities"

Merged Resolve "Implement appropriate filter for axial velocities"
1-implement-appropriate-filter-for-axial-velocities into main
Merged Sebastian Hofmann requested to merge 1-implement-appropriate-filter-for-axial-velocities into main
Compare
and
5 files
+ 384
− 47
Compare changes
  • Side-by-side
  • Inline
Files
5
code/f_prep_mix_time_mov.m 0 → 100644
+ 226
− 0
%% Prepare 3DPTV images for DAVIS
% Date: 20210307
% Author: Alexandra von Kameke (alexandra.vonkameke@tuhh.de)
% Editor: Sebastian Hofmann (sebastian.hofmann@tuhh.de)
close all
clear all
%% load images with readb16
%dirPath = 'D:\20201204\HSA\350rpm\'
dirPath = strcat(pwd,"\" );
mkdir("results")
files = dir(fullfile(dirPath,'*.b16'));
files = strcat(dirPath,{files.name}');
fname = char(files(3));
img1 = readB16(fname);
img1 = floor(img1./16);
[M, N] = size(img1);
img1(1030:1050,:) = 0;
maxi = max(img1(:));
mini= min(img1(:));
%% Create mask in reactor
close all
clims = [mini maxi]
figure;
imagesc((img1),[20 130])
roi = drawpolygon
mask = createMask(roi);
%% Define steps and minimum image conditions
nst = 50; %stepimgs
nn = length(files);
anf = 3;
intensity = zeros(1,nn);
minimgsimp = img1;
minimg = img1;
img_tot = zeros(round(nn/nst),size(img1,1),size(img1,2));
%% Create minimum image according steps nst
j = 0
for i = anf:nst:nn
j =j+1;
fname = char(files{i});
img = readB16(fname);
img = floor(img./16);
img(1030:1050,:) = 0;
%imgtot(j,:,:) = img;
B = minimgsimp>img;
minimgsimp(B) = img(B);
end
%minimg = squeeze(min(imgtot(:,:,:),[],1));
%% Show created images
close all
% figure
% imagesc((minimgsimp))
% figure;
% imagesc(img1,clims)
% figure;
% imagesc((img1-minimgsimp).*mask)
%% save results
save('results', 'minimgsimp','mask')
%load('results', 'minimgsimp','mask')
%% Design Butterworth high-pass filter (Sharpening)
% Assign the order value
n = 2; % one can change this value accordingly
% Assign Cut-off Frequency
D0 = 70 ; % one can change this value accordingly
% Designing filter
u = 0:(M-1);
v = 0:(N-1);
idx = find(u > M/2);
u(idx) = u(idx) - M;
idy = find(v > N/2);
v(idy) = v(idy) - N;
% MATLAB library function meshgrid(v, u) returns
% 2D grid which contains the coordinates of vectors
% v and u. Matrix V with each row is a copy of v
% and matrix U with each column is a copy of u
[V, U] = meshgrid(v, u);
% Calculating Euclidean Distance
D = sqrt(U.^2 + V.^2);
% determining the filtering mask
hh = 1./(1 + (D./D0).^(2*n));
hh = ones(size(hh,1),size(hh,2))-hh;
% figure
% surf((hh),'LineStyle','None')
%% Design Butterworth low-pass filter (Denoising)
% Assign the order value
ord = 2; % one can change this value accordingly
% Assign Cut-off Frequency
D0 = 320; % one can change this value accordingly, low value: much smoothening (50), high value: no smoothing (500)
% Designing filter
u = 0:(M-1);
v = 0:(N-1);
idx = find(u > M/2);
u(idx) = u(idx) - M;
idy = find(v > N/2);
v(idy) = v(idy) - N;
% MATLAB library function meshgrid(v, u) returns
% 2D grid which contains the coordinates of vectors
% v and u. Matrix V with each row is a copy of v
% and matrix U with each column is a copy of u
[V, U] = meshgrid(v, u);
% Calculating Euclidean Distance
D = sqrt(U.^2 + V.^2);
% determining the filtering mask
h = 1./(1 + (D./D0).^(2*ord));
% figure
% surf((h),'LineStyle','None')
%% subtract minimum image and filter (low pass high pass)
nn = length(files);
anf = 3
close all
multi = 650; %If multi too high, then the entire brighntess "room" won't be used - image too dark later,
% multi should be so high, that the center of the brightest particles is at
% approx. value of 4096 (see line 220) -> Perfect value is, if
% there are still brightness differences seen in one particle
% towards the edge
j = 0
for i = anf:nn
j =j+1
fname = char(files{i});
img = readB16(fname);
img = floor(img./16);
img(1030:1050,:) = 0;
img = double(img - minimgsimp);
% figure
% surf(img(590:670,100:200),'LineStyle','None')
% image(img(590:670,100:200)) % 1ST IMAGE ORIGINAL
% Getting Fourier Transform of the input_image
% using MATLAB library function fft2 (2D fast fourier transform)
FT_img = fft2(double(img));
% figure
% imagesc((log(abs((FT_img)))))
% Convolution between the Fourier Transformed
% image and the mask
G = hh.*FT_img;
% Getting the resultant image by Inverse Fourier Transform
% of the convoluted image using MATLAB library function
% ifft2 (2D inverse fast fourier transform)
img = real(ifft2(double(G)));
% figure
% surf(img(100:500,100:600),'LineStyle','None')
% image(img(590:670,100:200)) % 2ND IMAGE HIGH PASS
% Getting Fourier Transform of the input_image
% using MATLAB library function fft2 (2D fast fourier transform)
FT_img = fft2(double(img));
% Convolution between the Fourier Transformed
% image and the mask
G = h.*FT_img;
% Getting the resultant image by Inverse Fourier Transform
% of the convoluted image using MATLAB library function
% ifft2 (2D inverse fast fourier transform)
img = real(ifft2(double(G)));
% figure
%surf(img(100:500,100:600),'LineStyle','None')
% image(img(590:670,100:200)) % 3RD IMAGE LOW PASS
%
%imge = img./max(max(img));
imge = img./multi;
% figure
% imagesc(imge) % IMAGE TOTAL divided by "multi"
imge(imge>1)=1;
imge(imge<0)=0;
imge = imge.*4096;
% img_tot(j,:,:) = img;
imge = (double(imge.*mask));
% figure;
% %image(img(400:500,400:500)) %
% imagesc(imge,[0 4096]) % 4TH IMAGE TOTAL FINAL
imwrite(uint16(fliplr(imrotate(imge,90))),sprintf('results//Prep_%04d.tiff',i));
fclose('all');
end
% %%
% figure
% imagesc(imrotate(img,90))