Index: trunk/src/series/particle_detect.m =================================================================== --- trunk/src/series/particle_detect.m (revision 1088) +++ trunk/src/series/particle_detect.m (revision 1088) @@ -0,0 +1,629 @@ +%'merge_proj': concatene several fields from series, can project them on a regular grid in phys coordinates +%------------------------------------------------------------------------ +% function ParamOut=merge_proj(Param) +%------------------------------------------------------------------------ +%%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%% +% +%OUTPUT +% ParamOut: sets options in the GUI series.fig needed for the function +% +%INPUT: +% In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series. +% In batch mode, Param is the name of the corresponding xml file containing the same information +% when Param.Action.RUN=0 (as activated when the current Action is selected +% in series), the function ouput paramOut set the activation of the needed GUI elements +% +% Param contains the elements:(use the menu bar command 'export/GUI config' in series to +% see the current structure Param) +% .InputTable: cell of input file names, (several lines for multiple input) +% each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension} +% .OutputSubDir: name of the subdirectory for data outputs +% .OutputDirExt: directory extension for data outputs +% .Action: .ActionName: name of the current activated function +% .ActionPath: path of the current activated function +% .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled Matlab fct +% .RUN =0 for GUI input, =1 for function activation +% .RunMode='local','background', 'cluster': type of function use +% +% .IndexRange: set the file or frame indices on which the action must be performed +% .FieldTransform: .TransformName: name of the selected transform function +% .TransformPath: path of the selected transform function +% .InputFields: sub structure describing the input fields withfields +% .FieldName: name(s) of the field +% .VelType: velocity type +% .FieldName_1: name of the second field in case of two input series +% .VelType_1: velocity type of the second field in case of two input series +% .Coord_y: name of y coordinate variable +% .Coord_x: name of x coordinate variable +% .ProjObject: %sub structure describing a projection object (read from ancillary GUI set_object) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%======================================================================= +% Copyright 2008-2020, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France +% http://www.legi.grenoble-inp.fr +% Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr +% +% This file is part of the toolbox UVMAT. +% +% UVMAT 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 2 of the license, +% or (at your option) any later version. +% +% UVMAT 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 (see LICENSE.txt) for more details. +%======================================================================= + +function ParamOut=particle_detect(Param) + +%% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed +if isstruct(Param) && isequal(Param.Action.RUN,0) + ParamOut.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default) + ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default) + ParamOut.NbSlice='off'; %nbre of slices ('off' by default) + ParamOut.VelType='one';% menu for selecting the velocity type (options 'off'/'one'/'two', 'off' by default) + ParamOut.FieldName='one';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default) + ParamOut.FieldTransform = 'off';%can use a transform function + ParamOut.TransformPath=fullfile(fileparts(which('uvmat')),'transform_field');% path to transform functions (needed for compilation only) + ParamOut.ProjObject='off';%can use projection object(option 'off'/'on', + ParamOut.Mask='on';%can use mask option (option 'off'/'on', 'off' by default) + ParamOut.OutputDirExt='.detect';%set the output dir extension + ParamOut.OutputFileMode='NbInput';% '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice + %check the input files + ParamOut.CheckOverwriteVisible='on'; % manage the overwrite of existing files (default=1) + first_j=[]; + if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end + PairString=''; + if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end + [i1,i2,j1,j2] = get_file_index(Param.IndexRange.first_i,first_j,PairString); + FirstFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},... + Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2); + if ~exist(FirstFileName,'file') + msgbox_uvmat('WARNING',['the first input file ' FirstFileName ' does not exist']) + end + + prompt = {'threshold(th)';... + 'particle size (sz)' }; + dlg_title = 'get processing parameters'; + num_lines= 2; + def = {'4000';'3'}; + answer = inputdlg(prompt,dlg_title,num_lines,def); + if isempty(answer) + return + end + %check input consistency + ParamOut.ActionInput.th=str2num(answer{1}); + ParamOut.ActionInput.sz=str2num(answer{2}); + return +end + +%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%% +ParamOut=[]; %default output +RUNHandle=[]; +WaitbarHandle=[]; +%% read input parameters from an xml file if input is a file name (batch mode) +checkrun=1; +if ischar(Param) + Param=xml2struct(Param);% read Param as input file (batch case) + checkrun=0; +else + hseries=findobj(allchild(0),'Tag','series'); + RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series + WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series +end + +%% define the directory for result file (with path=RootPath{1}) +OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files + +if ~isfield(Param,'InputFields') + Param.InputFields.FieldName=''; +end + +%% root input file(s) name, type and index series +RootPath=Param.InputTable{1,1}; +RootFile=Param.InputTable{1,3}; +SubDir=Param.InputTable{1,2}; +NomType=Param.InputTable{1,4}; +FileExt=Param.InputTable{1,5}; + +hdisp=disp_uvmat('WAITING...','checking the file series',checkrun); +% gives the series of input file names and indices set by the input parameters: +%[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param); +% filecell{iview,fileindex}: +% iview: line in the table corresponding to a given file series +% fileindex: file index with i and j reshaped as a 1D array +% i1_series(iview,ref_j,ref_i)... are the corresponding arrays of indices i1,i2,j1,j2, depending on the input line iview and the two reference indices ref_i,ref_j +% i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices +if ~isempty(hdisp),delete(hdisp),end;%end the waiting display + + +%% determine the file type on each line from the first input file + + +%%%%%%%%%%%% END STANDARD PART %%%%%%%%%%%% +% EDIT FROM HERE + +RootFileOut=RootFile; + + +%% MAIN LOOP ON FIELDS +%%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%% +% for i_slice=1:NbSlice +% index_slice=i_slice:NbSlice:NbField;% select file indices of the slice +% NbFiles=0; +% nbmissing=0; + +%%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%% +tstart=tic; %used to record the computing time +CheckOverwrite=1;%default +if isfield(Param,'CheckOverwrite') + CheckOverwrite=Param.CheckOverwrite; +end +%%%%%% INPUT %%%%%% +th=Param.ActionInput.th +sz=Param.ActionInput.sz +% th=4000;%threshold on image intensity +% sz=3; %size of particles + +NbImage=Param.IndexRange.last_i-Param.IndexRange.first_i+1; +incr_i=Param.IndexRange.incr_i; +NbBlock=floor(NbImage/incr_i); + +%% MAIN LOOP +for index1=1:NbBlock + OutputFile=fullfile_uvmat(RootPath,OutputDir,RootFileOut,'.mat','_1-2',(index1-1)*incr_i+1,index1*incr_i) + if ~CheckOverwrite && exist(OutputFile,'file') + disp(['existing output file ' OutputFile ' already exists, skip to next field']) + end + for index=(index1-1)*incr_i+1:index1*incr_i + index + + %% reading input file(s) + ImgName=fullfile_uvmat(RootPath,SubDir,RootFile,FileExt,NomType,index) + Im=imread(ImgName); + [Ny,Nx]=size(Im); + % + % + % frame = kframe+idfile*framesperfiles; + % disp(frame); + % ImgName = sprintf(join(['%s/%s_cam%d_' format],''),folderin, ManipName, CamNum, frame); + % fprintf("%s \n",ImgName); + %Im = (cast(imread(ImgName),'like',Background) - Background).*mask; + %% normalizeimage + Im = double(Im); + + % Imfiltered=filter2(transfert_coef,Im); + % + % Imfiltered(:,1:windowsize)=Imfiltered(:,windowsize)*ones(1,windowsize); + % Imfiltered(:,end-windowsize+1:end)=Imfiltered(:,end-windowsize+1)*ones(1,windowsize); + % Imfiltered(1:windowsize,:)=ones(windowsize,1)*Imfiltered(windowsize,:); + % Imfiltered(end-windowsize+1:end,:)=ones(windowsize,1)*Imfiltered(end-windowsize,:); + % Im=Im./(2*Imfiltered); + %% + + out=pkfnd(Im,th,sz); % Provides intensity maxima positionsth,sz,Test,BackgroundType,format + npar = size(out,1); + + %% We keep only spots with a gaussian shape + cnt = 0; + x = []; + y = []; + for j = 1:npar + Nwidth = 1; + if (out(j,2)-Nwidth >0)&&(out(j,1)-Nwidth>0)&&(out(j,2)+Nwidth1 + savefile=OutputFile; + save(savefile,"CC",'-v7.3') + m = matfile(savefile,'Writable',true); + m.nframes = incr_i; + else + figure("NumberTitle","Off","Name",['RAW picture,' SubDir]) + imshow(imread(ImgName),[0,5000]) + colormap gray + % figure("NumberTitle","Off","Name",sprintf("%s, cam %d",BackgroundType,CamNum)) + % imshow(BackgroundMin,[0,5000]) + % colormap gray + % figure("NumberTitle","Off","Name",sprintf("RAW picture - Background, cam %d, frame %d",CamNum,kframe)) + % imshow(Im,[0,th]) + % colormap gray + colorbar + + %% Tracé de l'histogramme des intensités pour définir le seuil + fig = figure('NumberTitle','Off','Name','Intensity histogram'); + histogram(Im,1000) + xlabel("Intensity") + ylabel("Number") + set(gca, 'XScale', 'log') + set(gca, 'YScale', 'log') + + Nx = size(Im,2); + Ny = size(Im,1); + + out=pkfnd(Im,th,sz); % Provides intensity maxima positions + npar = size(out,1); + + %% We keep only spots with a gaussian shape + cnt = 0; + x = []; + y = []; + for j = 1:npar + Nwidth = 1; + if (out(j,2)-Nwidth >0)&&(out(j,1)-Nwidth>0)&&(out(j,2)+Nwidth test mode, false-> classic mode (optional) +%%% BackgroundType (optional) : determine which background is substracted to pictures. By defaut is equal to BackgroundMean, +%%% format (optional) : picture names. By defaut it is '%05d.tif'. +%%% The beginning of picture names has to be %ManipName_cam%CamNum_%format +%-------------------------------------------------------------------------------- +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +close all + +framesperfiles = 1000; + +%% Test if Test exist or not +if ~exist('Test','var') + Test=false; +end + +%% Test if BackgroundType exist or not +if ~exist('BackgroundType','var') + BackgroundType="BackgroundMean"; +end + +% By defaut format='%05d.tif' +if ~exist('format','var') + format='%05d.tif'; +end + +if ~exist('withmask','var') + withmask=false; +end + + + +%% Definition of folders +fprintf(ManipName); +folderin = sprintf("%sDATA/%s/cam%d",session.input_path,ManipName,CamNum) +folderout = sprintf("%sProcessed_DATA/%s",session.output_path,ManipName) +BackgroundFile = sprintf("%s/Background_cam%d.mat",folderout,CamNum); + +%% Find centers +if exist(folderout,'dir')==0 + mkdir(char(folderout)); +end + +% if exist(strcat(folderout,'/Parallel/Matching'),'dir')==0 +% mkdir(char(strcat(folderout,'/Parallel/Matching'))); +% end + +load(BackgroundFile,'BackgroundMin','BackgroundMax','BackgroundMean') +%% Choice of background type +if BackgroundType=="BackgroundMean" + Background=BackgroundMean; +elseif BackgroundType=="BackgroundMax" + Background=BackgroundMax; +elseif BackgroundType=="BackgroundMin" + Background=BackgroundMin; +end +Nx = size(Background,2); +Ny = size(Background,1); + +if withmask + ImgName = sprintf(join(['%s/%s_cam%d_' format],''),folderin, ManipName, CamNum, 1); + Im = imread(ImgName); + imshow(Im,[0,10000]) + BW = roipoly; + mask = cast(BW,'like',Background); +else + mask = cast(ones(Ny,Nx),'like',Background); +end + +%% param for normailze image +windowsize=round(min(Ny,Nx)/20); +transfert_coef = ones(windowsize)/windowsize^2; +%% + +close all +if ~Test + idfirstfile = fix(firstFrame/framesperfiles); + idlastfile = fix((nframes-1)/framesperfiles); + for idfile = idfirstfile:idlastfile + firstFrame = rem(firstFrame,framesperfiles); + lastframe = min(framesperfiles,nframes-idfile*framesperfiles); + for kframe=firstFrame:lastframe + frame = kframe+idfile*framesperfiles; + disp(frame); + ImgName = sprintf(join(['%s/%s_cam%d_' format],''),folderin, ManipName, CamNum, frame); + fprintf("%s \n",ImgName); + Im = (cast(imread(ImgName),'like',Background) - Background).*mask; + %% normalizeimage + Im = double(Im); + + Imfiltered=filter2(transfert_coef,Im); + + Imfiltered(:,1:windowsize)=Imfiltered(:,windowsize)*ones(1,windowsize); + Imfiltered(:,end-windowsize+1:end)=Imfiltered(:,end-windowsize+1)*ones(1,windowsize); + Imfiltered(1:windowsize,:)=ones(windowsize,1)*Imfiltered(windowsize,:); + Imfiltered(end-windowsize+1:end,:)=ones(windowsize,1)*Imfiltered(end-windowsize,:); + Im=Im./(2*Imfiltered); + %% + + + out=pkfnd(Im,th,sz); % Provides intensity maxima positions + npar = size(out,1); + + %% We keep only spots with a gaussian shape + cnt = 0; + x = []; + y = []; + for j = 1:npar + Nwidth = 1; + if (out(j,2)-Nwidth >0)&&(out(j,1)-Nwidth>0)&&(out(j,2)+Nwidth0)&&(out(j,1)-Nwidth>0)&&(out(j,2)+Nwidth th); +n=length(ind); + +if n==0 + out=[];[i,j,ind]=find(im > th); + fprintf('nothing above threshold'); + return; +end +mx=[]; +%convert index from find to row and column +rc=[i,j]; % j corresponds to x axis and i to y axis +% rc=[j,i]; % j corresponds to x axis and i to y axis +for ii=1:n + r=rc(ii,1); + c=rc(ii,2); + %check each pixel above threshold to see if it's brighter than it's neighbors + % THERE'S GOT TO BE A FASTER WAY OF DOING THIS. I'M CHECKING SOME MULTIPLE TIMES, + % BUT THIS DOESN'T SEEM THAT SLOW COMPARED TO THE OTHER ROUTINES, ANYWAY. + if r>1 && r1 && c=im(r-1,c-1) && im(r,c)>=im(r,c-1) && im(r,c)>=im(r+1,c-1) && ... + im(r,c)>=im(r-1,c) && im(r,c)>=im(r+1,c) && ... + im(r,c)>=im(r-1,c+1) && im(r,c)>=im(r,c+1) && im(r,c)>=im(r+1,c+1) + mx=[mx,[r,c]']; %#ok + %tst(ind(i))=im(ind(i)); + end + end +end +%out=tst; +mx=mx'; + +[npks,crap]=size(mx); + +%if size is specified, then get ride of pks within size of boundary +if nargin==3 && npks>0 + %throw out all pks within sz of boundary; + ind=find(mx(:,1)>sz & mx(:,1)<(nr-sz) & mx(:,2)>sz & mx(:,2)<(nc-sz)); + mx=mx(ind,:); +end + +%prevent from finding peaks within size of each other +[npks,crap]=size(mx); +if npks > 1 + %CREATE AN IMAGE WITH ONLY PEAKS + nmx=npks; + tmp=0.*im; + for i=1:nmx + tmp(mx(i,1),mx(i,2))=im(mx(i,1),mx(i,2)); + end + %LOOK IN NEIGHBORHOOD AROUND EACH PEAK, PICK THE BRIGHTEST + for i=1:nmx + roi=tmp( (mx(i,1)-floor(sz/2)):(mx(i,1)+(floor(sz/2)+1)),(mx(i,2)-floor(sz/2)):(mx(i,2)+(floor(sz/2)+1))) ; + [mv,indi]=max(roi); + [mv,indj]=max(mv); + tmp( (mx(i,1)-floor(sz/2)):(mx(i,1)+(floor(sz/2)+1)),(mx(i,2)-floor(sz/2)):(mx(i,2)+(floor(sz/2)+1)))=0; + tmp(mx(i,1)-floor(sz/2)+indi(indj)-1,mx(i,2)-floor(sz/2)+indj-1)=mv; + end + ind=find(tmp>0); + mx=[mod(ind,nr),floor(ind/nr)+1]; +end + +if size(mx)==[0,0] + out=[]; +else + out(:,2)=mx(:,1); + out(:,1)=mx(:,2); +end