Homec4science
Diffusion Motion correction paper (master) Browse

R6289/modified_motion_compensate.m
d373ff55171cmaster

Last Change
Raw File

modified_motion_compensate.m
View Options

function motion_compensate(fnInput1,fnInput2,fnMatch,fnDeepMatching,fnOut1,fnOut2,fnColor,fnVector,N,param,tmpdir,outdir)
%% Description
% Motion is estimated with a brigthness constancy data term defined on fnInput1
% and a feature matching similarity term defined on fnMatch. The sequences fnIput1 and
% fnInput2 are warped according to the estimated motion field.
% For more details see the paper: "Imaging neural activity in the ventral
% nerve cord of behaving adult Drosophila", bioRxiv
%
%% Input
% fnInput1: filename of the sequence used for the brightness constancy term, in TIF format
% fnInput2: filename of the sequence warped with the motion field estimated from fnInput1 and fnMatch, in TIF format
% fnMatch: filename of the sequence used for feature matching, in TIF format
% fnDeepMatching: filename of the deepmatching code
% fnOut1: filename used to save the warped version of fnInput1, in TIF format
% fnOut2: filename used to save the warped version of fnInput2, in TIF format
% fnOut2: filename used to save the color visualization of the estimated motion, in TIF format
% N: number of frames to process
% param: parameters of the algorithm (see 'default_parameters.m')
%
% Copyright (C) 2017 D. Fortun, denis.fortun@epfl.ch
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%% Input data
seq=mijread_stack(fnInput1);
seqW=mijread_stack(fnInput2);
seqMatch=mijread_stack(fnMatch);
if N==-1
N=size(seq,3);
end
seq=double(seq(:,:,1:N));
seqRescale=(seq-min(seq(:)))/(max(seq(:))-min(seq(:)))*255;
seqW=double(seqW(:,:,1:N));
seqMatch=double(seqMatch(:,:,1:N));
seqMatch=(seqMatch-min(seqMatch(:)))/(max(seqMatch(:))-min(seqMatch(:)))*255;
% Motion estimation
w=zeros(size(seqRescale,1),size(seqRescale,2),2,size(seqRescale,3)-1);
colorFlow=zeros(size(seqRescale,1),size(seqRescale,2),3,size(seqRescale,3)-1);
vectorFlow=zeros(512,512,size(seqRescale,3)-1);
i1=seqRescale(:,:,1);
i1Match=seqMatch(:,:,1);
parfor t=1:N-1 % Replace parfor by for if you don't want to parallelize
presteps_to_compute_motion(i1Match,seqMatch(:,:,t+1),t,tmpdir);
end
fprintf('Presteps_to_compute_motion done!\n')
return_number_of_running_processes = 'ps r | wc -l';
for t=1:N-1
[status, n_processes] = system(return_number_of_running_processes);
while str2num(n_processes)-1 > 28
pause(1);
[status, n_processes] = system(return_number_of_running_processes);
end
fprintf(['Strating deep of frame ', num2str(t),'\n']);
fnI1=[tmpdir,'/I1_',num2str(t),'.png'];
fnI2=[tmpdir,'/I2_',num2str(t),'.png'];
fnMatch=[tmpdir,'/match_',num2str(t),'.txt'];
command = [fnDeepMatching,'/deepmatching ',fnI1,' ',fnI2,' -out ',fnMatch,' &'];
s = system(command);
end
parfor t = 1:N-1
i2=seqRescale(:,:,t+1);
[i10,i2]=midway(i1,i2);
w(:,:,:,t) = poststeps_to_compute_motion(i10,i2,i1Match,param,t,tmpdir);
fprintf(['poststeps_to_compute_motion done for frame ',num2str(t),'\n']);
colorFlow(:,:,:,t) = flowToColor(w(:,:,:,t));
vectorFlow(:,:,t) = flowToVectorImage(w(:,:,:,t),[18,18],[512,512],t,tmpdir);
end
parfor t = 1:N-1
dlmwrite([outdir,'/wx_frame',num2str(t),'.dat'],w(:,:,1,t));
dlmwrite([outdir,'/wy_frame',num2str(t),'.dat'],w(:,:,2,t));
end
%% Registration
seqWarped=seq;
seqwWarped=seqW;
parfor t=1:N-1
seqWarped(:,:,t+1)=warpImg(seq(:,:,t+1),w(:,:,1,t),w(:,:,2,t));
seqwWarped(:,:,t+1)=warpImg(seqW(:,:,t+1),w(:,:,1,t),w(:,:,2,t));
end
%% Save
mijwrite_stack(single(seqWarped),fnOut1);
mijwrite_stack(single(seqwWarped),fnOut2);
mijwrite_stack(single(colorFlow),fnColor,1);
mijwrite_stack(single(vectorFlow),fnVector);
c4science · Help