-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathtest_3DR_topadjfilt.m~
132 lines (122 loc) · 4.72 KB
/
test_3DR_topadjfilt.m~
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
clc; clear; close all; warning('off');
dataset='150811_DeepPIV_BathoStyg_3DR_xlarge';
display(dataset);
vel=10; %mm/s %UNKNOWN
inc=1;
%retrieving data set-specific parameters
[dir,start,finish,fps,fstop,shutter,calib,red,aspectratio]=videoinfo(dataset,vel);
nFrames=finish-start;
display('Loading video')
viddir=[dir,'clips/'];
vid=VideoReader([viddir,dataset,'.mov']);
%% arranging planar cross sections into "volume" matrix
if red==1
display('Arranging images into a stack using the red color channel only');
else
display('Arranging images into a stack using grayscale of images');
end
im_file=[viddir,dataset,'_IMAGE.mat'];
if exist(im_file,'file')==2
load(im_file); %if exists, load original image stack
else
im=read(vid,1)*0;
IMAGE(:,:,(finish-start)/inc)=im;
for i=start:inc:finish
im=read(vid,i);
if red==1
IMAGE(:,:,(i-start)/inc+1)=im(:,:,1);
else
imgray=rgb2gray(im);
IMAGE(:,:,(i-start)/inc+1)=imgray;
end
end
save(im_file,'IMAGE','-v7.3'); %save post-processing parameters
end
close all; clear vid
%% retrieving/determining data set-specific post-processing parameters
display('Post-processing images')
par_file=[viddir,dataset,'_',num2str(vel),'vel_params.mat'];
if exist(par_file,'file')==2
load(par_file); %if exists, load post-proc parameters
else
hfig=crop_new(IMAGE); %define cropped boundaries on image
waitfor(hfig);
hfig=manipulate_new(IMAGE,output); %define threshold and filter settings
waitfor(hfig);
save(par_file,'output'); %save post-processing parameters
end
close all
%% post-process all cross section images
display('Saving post-processed image stack')
im_filenew=[viddir,dataset,'_IMAGEcrop.mat'];
if exist(im_filenew,'file')==2
load(im_filenew);
else
IMAGEcrop=postprocIM(IMAGE,output); %crop, threshold, and filter images
save(im_filenew,'IMAGEcrop','-v7.3'); %save post-processing parameters
end
close all; clear IMAGE
%%
display('Saving post-processed image stack')
im_filenew=[viddir,dataset,'_IMAGEsm.mat'];
if exist(im_filenew,'file')==2
load(im_filenew);
else
im=IMAGEcrop(:,:,1)*0;imtop=im;imdil=im;imadj=im;imbin=im;imero=im;imnew=im;
IMAGEnew=IMAGEcrop*0;
for i=1:1:size(IMAGEcrop,3)
im=IMAGEcrop(:,:,i);
imtop = imtophat(im,strel('disk',15));
imadj=imadjust(imtop,[0.08,0.39],[]);%using low and high limits provides absolute values for image adjustments
imfilt=medfilt2(imadj,[5,5]);
figure(1)
imshow([im,imtop;imadj,imfilt])
width=size(im,2);height=size(im,1);
text(width*0.5,20,'im','Color','w','FontSize',16);text(width*1.5,20,'imtop','Color','w','FontSize',16);
text(width*0.5,height,'imadj','Color','w','FontSize',16);text(width*1.5,height,'imfilt','Color','w','FontSize',16);
imdil=imdilate(imadj,strel('disk',10));
imbin=im2bw(imdil,0.01);
% figure(2)
% imshow([im;imdil;imbin*255])
% imdil=imdilate(imbin,strel('rectangle',[2,4]));
imfilled=imfill(imbin,'holes');
imopen=bwareaopen(imfilled,500,4);
figure(2)
imshow([imbin*255,imdil*255;imfilled*255,imopen*255])
text(width*0.5,20,'imbin','Color','w','FontSize',16);text(width*1.5,20,'imdil','Color','w','FontSize',16);
text(width*0.5,height,'imfilled','Color','w','FontSize',16);text(width*1.5,height,'imopen','Color','w','FontSize',16);
imnew=imfilt.*uint8(imopen);
figure(3)
imshow([im,imnew])
text(width*0.5,20,'im','Color','w','FontSize',16);text(width*1.5,20,'imnew','Color','w','FontSize',16);
IMAGEnew(:,:,i)=imnew;
end
IMAGEsm=smooth3(IMAGEnew); %smooth images
save(im_filenew,'IMAGEsm','-v7.3'); %save post-processing parameters
end
close all; clear('IMAGEcrop','IMAGEnew')
%% plotting 3D reconstruction
display('Plotting image stack using isosurfaces')
%define physical domain
%%NEED TO CORRECT THIS SOON!
% length=(finish-start)/fps*vel*calib;
% [X,Y,Z]=meshgrid((1:size(IMAGEsm,2))/calib,(1:size(IMAGEsm,1))/calib,1:length/size(IMAGEsm,3):length);
[X,Y,Z]=meshgrid(1:size(IMAGEsm,2),1:size(IMAGEsm,1),1:size(IMAGEsm,3));
xx=single(X);yy=single(Y);zz=single(Z);
clear X Y Z
%plot isosurfaces
figure(3)
val=10;
display(['Plotting isosurface with val = ',num2str(val)]);
hiso=patch(isosurface(xx,yy,zz,IMAGEsm,val,'verbose'),'FaceColor',[0.5,0.5,0.65],'EdgeColor','none');
isonormals(IMAGEsm,hiso)
alpha(0.2) %sets level of transparency
view(3)
camlight
lighting gouraud
grid on
daspect(aspectratio);
% axis equal
outfile=[viddir,dataset,'_3DR_topadjfilt_',num2str(val),'.mp4'];
display(['Saving 3DR movie with val = ',num2str(val)]);
makevideo(outfile);