diff --git a/.gitignore b/.gitignore index 427e368..02f00ec 100644 --- a/.gitignore +++ b/.gitignore @@ -1,124 +1,121 @@ #swap files **.swp # executable for deepmatching code/external/deepmatching_1.2.2_c++_mac/deepmatching # temporary folder created by program tmp/ # macOS folder attributes **/.DS_Store -# matlab files -**/*.m - **.o data/ results/ # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] *$py.class # C extensions *.so # Distribution / packaging .Python build/ develop-eggs/ dist/ downloads/ eggs/ .eggs/ lib/ lib64/ parts/ sdist/ var/ wheels/ *.egg-info/ .installed.cfg *.egg MANIFEST # PyInstaller # Usually these files are written by a python script from a template # before PyInstaller builds the exe, so as to inject date/other infos into it. *.manifest *.spec # Installer logs pip-log.txt pip-delete-this-directory.txt # Unit test / coverage reports htmlcov/ .tox/ .coverage .coverage.* .cache nosetests.xml coverage.xml *.cover .hypothesis/ .pytest_cache/ # Translations *.mo *.pot # Django stuff: *.log local_settings.py db.sqlite3 # Flask stuff: instance/ .webassets-cache # Scrapy stuff: .scrapy # Sphinx documentation docs/_build/ # PyBuilder target/ # Jupyter Notebook .ipynb_checkpoints # pyenv .python-version # celery beat schedule file celerybeat-schedule # SageMath parsed files *.sage.py # Environments .env .venv env/ venv/ ENV/ env.bak/ venv.bak/ # Spyder project settings .spyderproject .spyproject # Rope project settings .ropeproject # mkdocs documentation /site # mypy .mypy_cache/ diff --git a/code/flowToVectorImage.m b/code/flowToVectorImage.m new file mode 100644 index 0000000..7b775be --- /dev/null +++ b/code/flowToVectorImage.m @@ -0,0 +1,33 @@ +function A = flowToVectorImage(w, n_sampling_points, resolution, frame_number, tmpdir) +% flowToVectorImage(w, sampling_rate, resolution) +% w is flow of form n x n x 2 +% n_sampling_points is the number of sampling point in x and y direction +% resolution is an array of length 2 which specifies the resolution of the resulting vector image +% frame_number is needed for parallel execution to separate tmp files + sampling_rate_1 = round(size(w,1)/n_sampling_points(1)); + sampling_rate_2 = round(size(w,1)/n_sampling_points(2)); + u = w(1:sampling_rate_1:end,1:sampling_rate_1:end,1); + v = w(1:sampling_rate_2:end,1:sampling_rate_2:end,2); + [x,y] = meshgrid(1:sampling_rate_2:size(w,2),1:sampling_rate_1:size(w,1)); + figure; + quiver(x,y,u,v,'color','k','AutoScale','off','LineWidth', 8); + set(gca,'XTick',[]); % Remove the ticks in the x axis! + set(gca,'YTick',[]); % Remove the ticks in the y axis + set(gca,'XLim',[1,size(w,1)]); + set(gca,'YLim',[1,size(w,2)]); + set(gca,'Visible','off') + set(gcf,'PaperUnits','inches','PaperPosition',[0 0 resolution(1) resolution(2)]); + ax = gca; + outerpos = ax.OuterPosition; + ti = ax.TightInset; + left = outerpos(1) + ti(1); + bottom = outerpos(2) + ti(2); + ax_width = outerpos(3) - ti(1) - ti(3); + ax_height = outerpos(4) - ti(2) - ti(4); + ax.Position = [left bottom ax_width ax_height]; + + print([tmpdir,'tmpvector',num2str(frame_number),'.tiff'],'-dtiff','-r1'); + im = imread([tmpdir,'tmpvector',num2str(frame_number),'.tiff']); + %A = imerode(imadjust(rgb2gray(im)), strel('square',1)); + + delete([tmpdir,'tmpvector',num2str(frame_number),'.tiff']); diff --git a/code/parse_json.m b/code/parse_json.m new file mode 100644 index 0000000..77cd7bd --- /dev/null +++ b/code/parse_json.m @@ -0,0 +1,205 @@ +function data = parse_json(string) +% DATA = PARSE_JSON(string) +% This function parses a JSON string and returns a cell array with the +% parsed data. JSON objects are converted to structures and JSON arrays are +% converted to cell arrays. + +% F. Glineur, 2009 +% (inspired by the JSON parser by Joel Feenstra on MATLAB File Exchange +% (http://www.mathworks.com/matlabcentral/fileexchange/20565) but with +% faster handling of strings) + +pos = 1; +len = length(string); +% String delimiters and escape characters are identified beforehand to improve speed +esc = regexp(string, '["\\]'); index_esc = 1; len_esc = length(esc); + +if pos <= len + switch(next_char) + case '{' + data = parse_object; + case '[' + data = parse_array; + otherwise + error_pos('Outer level structure must be an object or an array'); + end +end + + function object = parse_object + parse_char('{'); + object = []; + if next_char ~= '}' + while 1 + str = parse_string; + if isempty(str) + error_pos('Name of value at position %d cannot be empty'); + end + parse_char(':'); + val = parse_value; + object.(valid_field(str)) = val; + if next_char == '}' + break; + end + parse_char(','); + end + end + parse_char('}'); + end + + function object = parse_array + parse_char('['); + object = cell(0, 1); + if next_char ~= ']' + while 1 + val = parse_value; + object{end+1} = val; + if next_char == ']' + break; + end + parse_char(','); + end + end + parse_char(']'); + end + + function parse_char(c) + skip_whitespace; + if pos > len || string(pos) ~= c + error_pos(sprintf('Expected %c at position %%d', c)); + else + pos = pos + 1; + skip_whitespace; + end + end + + function c = next_char + skip_whitespace; + if pos > len + c = []; + else + c = string(pos); + end + end + + function skip_whitespace + while pos <= len && isspace(string(pos)) + pos = pos + 1; + end + end + + function str = parse_string + if string(pos) ~= '"' + error_pos('String starting with " expected at position %d'); + else + pos = pos + 1; + end + str = ''; + while pos <= len + while index_esc <= len_esc && esc(index_esc) < pos + index_esc = index_esc + 1; + end + if index_esc > len_esc + str = [str string(pos:end)]; + pos = len + 1; + break; + else + str = [str string(pos:esc(index_esc)-1)]; + pos = esc(index_esc); + end + switch string(pos) + case '"' + pos = pos + 1; + return; + case '\' + if pos+1 > len + error_pos('End of file reached right after escape character'); + end + pos = pos + 1; + switch string(pos) + case {'"' '\' '/'} + str(end+1) = string(pos); + pos = pos + 1; + case {'b' 'f' 'n' 'r' 't'} + str(end+1) = sprintf(['\' string(pos)]); + pos = pos + 1; + case 'u' + if pos+4 > len + error_pos('End of file reached in escaped unicode character'); + end + str(end+1:end+6) = string(pos-1:pos+4); + pos = pos + 5; + end + otherwise % should never happen + str(end+1) = string(pos); + pos = pos + 1; + end + end + error_pos('End of file while expecting end of string'); + end + + function num = parse_number + [num, one, err, delta] = sscanf(string(pos:min(len,pos+20)), '%f', 1); % TODO : compare with json(pos:end) + if ~isempty(err) + error_pos('Error reading number at position %d'); + end + pos = pos + delta-1; + end + + function val = parse_value + switch(string(pos)) + case '"' + val = parse_string; + return; + case '[' + val = parse_array; + return; + case '{' + val = parse_object; + return; + case {'-','0','1','2','3','4','5','6','7','8','9'} + val = parse_number; + return; + case 't' + if pos+3 <= len && strcmpi(string(pos:pos+3), 'true') + val = true; + pos = pos + 4; + return; + end + case 'f' + if pos+4 <= len && strcmpi(string(pos:pos+4), 'false') + val = false; + pos = pos + 5; + return; + end + case 'n' + if pos+3 <= len && strcmpi(string(pos:pos+3), 'null') + val = []; + pos = pos + 4; + return; + end + end + error_pos('Value expected at position %d'); + end + + function error_pos(msg) + poss = max(min([pos-15 pos-1 pos pos+20],len),1); + if poss(3) == poss(2) + poss(3:4) = poss(2)+[0 -1]; % display nothing after + end + msg = [sprintf(msg, pos) ' : ... ' string(poss(1):poss(2)) '' string(poss(3):poss(4)) ' ... ']; + ME = MException('JSONparser:invalidFormat', msg); + throw(ME); + end + + function str = valid_field(str) + % From MATLAB doc: field names must begin with a letter, which may be + % followed by any combination of letters, digits, and underscores. + % Invalid characters will be converted to underscores, and the prefix + % "alpha_" will be added if first character is not a letter. + if ~isletter(str(1)) + str = ['alpha_' str]; + end + str(~isletter(str) & ~('0' <= str & str <= '9')) = '_'; + end + +end \ No newline at end of file diff --git a/code/poststeps_to_compute_motion.m b/code/poststeps_to_compute_motion.m new file mode 100755 index 0000000..62203d5 --- /dev/null +++ b/code/poststeps_to_compute_motion.m @@ -0,0 +1,84 @@ +function w = poststeps_to_compute_motion(I1,I2,I1Match,param,t,tmpdir) +%% Description +% Compute motion with a brigthness constancy data term defined on (I1,I2) +% and a feature matching similarity term defined on (I1Match,I2Match). The +% feature matches are computed by the DeepMatching method [Revaud et al. 2016], with the +% source code provided by the authors +% For more details see the paper: "Imaging neural activity in the ventral +% nerve cord of behaving adult Drosophila", bioRxiv +% +% [Revaud et al. 2016] J. Revaud, P. Weinzaepfel, Z. Harchaoui and C. Schmid (2016) "DeepMatching: +% Hierarchical Deformable Dense Matching. Int J Comput Vis 120:300–323 +%% Input +% I1,I2: Images used for the brightness constancy term +% I1Match,I2Match: Images used by the feature matching similarity term +% fnDeepMatching: filename of the deepmatching code +% param: parameters of the algorithm (see 'default_parameters.m') +% t: number of the frame in the image sequence +% +% 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. + +sz0 = size(I1); +I1=padarray(I1,[15,15],'replicate'); +I2=padarray(I2,[15,15],'replicate'); + +fnI1= [tmpdir,'/I1_',num2str(t),'.png']; +fnI2= [tmpdir,'/I2_',num2str(t),'.png']; +fnMatch=[tmpdir,'/match_',num2str(t),'.txt']; + +formatSpec = '%u %u %u %u %f %u'; +sizeCorresp = [6 Inf]; + +f=fopen(fnMatch); +while f==-1 + f=fopen(fnMatch); +end +corresp = fscanf(f,formatSpec,sizeCorresp); +while size(corresp,1)==0 + corresp = fscanf(f,formatSpec,sizeCorresp); +end +command = ['rm ',fnI1]; s = system(command); +command = ['rm ',fnI2]; s = system(command); +command = ['rm ',fnMatch]; s = system(command); + +thresh=param.threshMatch; +Iseg=segment(I1Match,'variance',thresh); + +matches = zeros(5,1); +k=0; +for i=1:size(corresp,2) + if corresp(1,i)>0 && corresp(1,i)<=sz0(2) && corresp(2,i)>0 && corresp(2,i)<=sz0(1) + if Iseg(corresp(2,i),corresp(1,i)) == 1 + k = k+1; + matches(1,k) = corresp(1,i); + matches(2,k) = corresp(2,i); + matches(4,k) = corresp(4,i)-corresp(2,i); + matches(3,k) = corresp(3,i)-corresp(1,i); + end + end +end + +scores = corresp(5,:); +idx = find(corresp(5,:)); +scores = scores(scores~=0); +scores = (scores-min(scores))/(max(scores)-min(scores)); +for j=1:length(idx) + matches(5,idx(j)) = scores(j); +end + +%% Optical flow +print = 0; % if print =1:verbose mode +disp = 0; % if disp=1; displays the results at each iteration +w = of_l1_l2_fm_admm(I1,I2,sz0,matches,param,print,disp); + +w = crop_fit_size_center(w,[sz0(1),sz0(2),2]); diff --git a/code/presteps_to_compute_motion.m b/code/presteps_to_compute_motion.m new file mode 100755 index 0000000..e392160 --- /dev/null +++ b/code/presteps_to_compute_motion.m @@ -0,0 +1,34 @@ +function presteps_to_compute_motion(I1Match,I2Match,t,tmpdir) +%% Description +% Compute motion with a brigthness constancy data term defined on (I1,I2) +% and a feature matching similarity term defined on (I1Match,I2Match). The +% feature matches are computed by the DeepMatching method [Revaud et al. 2016], with the +% source code provided by the authors +% For more details see the paper: "Imaging neural activity in the ventral +% nerve cord of behaving adult Drosophila", bioRxiv +% +% [Revaud et al. 2016] J. Revaud, P. Weinzaepfel, Z. Harchaoui and C. Schmid (2016) "DeepMatching: +% Hierarchical Deformable Dense Matching. Int J Comput Vis 120:300–323 +%% Input +% I1,I2: Images used for the brightness constancy term +% I1Match,I2Match: Images used by the feature matching similarity term +% fnDeepMatching: filename of the deepmatching code +% param: parameters of the algorithm (see 'default_parameters.m') +% t: number of the frame in the image sequence +% +% 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. + +fnI1=[tmpdir,'/I1_',num2str(t),'.png']; +fnI2=[tmpdir,'/I2_',num2str(t),'.png']; +imwrite(uint8(I1Match),fnI1); +imwrite(uint8(I2Match),fnI2); diff --git a/modified_motion_compensate.m b/modified_motion_compensate.m new file mode 100755 index 0000000..7ea22c4 --- /dev/null +++ b/modified_motion_compensate.m @@ -0,0 +1,97 @@ +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);