[1088] | 1 | %'merge_proj': concatene several fields from series, can project them on a regular grid in phys coordinates |
---|
| 2 | %------------------------------------------------------------------------ |
---|
| 3 | % function ParamOut=merge_proj(Param) |
---|
| 4 | %------------------------------------------------------------------------ |
---|
| 5 | %%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 6 | % |
---|
| 7 | %OUTPUT |
---|
| 8 | % ParamOut: sets options in the GUI series.fig needed for the function |
---|
| 9 | % |
---|
| 10 | %INPUT: |
---|
| 11 | % In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series. |
---|
| 12 | % In batch mode, Param is the name of the corresponding xml file containing the same information |
---|
| 13 | % when Param.Action.RUN=0 (as activated when the current Action is selected |
---|
| 14 | % in series), the function ouput paramOut set the activation of the needed GUI elements |
---|
| 15 | % |
---|
| 16 | % Param contains the elements:(use the menu bar command 'export/GUI config' in series to |
---|
| 17 | % see the current structure Param) |
---|
| 18 | % .InputTable: cell of input file names, (several lines for multiple input) |
---|
| 19 | % each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension} |
---|
| 20 | % .OutputSubDir: name of the subdirectory for data outputs |
---|
| 21 | % .OutputDirExt: directory extension for data outputs |
---|
| 22 | % .Action: .ActionName: name of the current activated function |
---|
| 23 | % .ActionPath: path of the current activated function |
---|
| 24 | % .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled Matlab fct |
---|
| 25 | % .RUN =0 for GUI input, =1 for function activation |
---|
| 26 | % .RunMode='local','background', 'cluster': type of function use |
---|
| 27 | % |
---|
| 28 | % .IndexRange: set the file or frame indices on which the action must be performed |
---|
| 29 | % .FieldTransform: .TransformName: name of the selected transform function |
---|
| 30 | % .TransformPath: path of the selected transform function |
---|
| 31 | % .InputFields: sub structure describing the input fields withfields |
---|
| 32 | % .FieldName: name(s) of the field |
---|
| 33 | % .VelType: velocity type |
---|
| 34 | % .FieldName_1: name of the second field in case of two input series |
---|
| 35 | % .VelType_1: velocity type of the second field in case of two input series |
---|
| 36 | % .Coord_y: name of y coordinate variable |
---|
| 37 | % .Coord_x: name of x coordinate variable |
---|
| 38 | % .ProjObject: %sub structure describing a projection object (read from ancillary GUI set_object) |
---|
| 39 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 40 | |
---|
| 41 | %======================================================================= |
---|
[1093] | 42 | % Copyright 2008-2021, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France |
---|
[1088] | 43 | % http://www.legi.grenoble-inp.fr |
---|
| 44 | % Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr |
---|
| 45 | % |
---|
| 46 | % This file is part of the toolbox UVMAT. |
---|
| 47 | % |
---|
| 48 | % UVMAT is free software; you can redistribute it and/or modify |
---|
| 49 | % it under the terms of the GNU General Public License as published |
---|
| 50 | % by the Free Software Foundation; either version 2 of the license, |
---|
| 51 | % or (at your option) any later version. |
---|
| 52 | % |
---|
| 53 | % UVMAT is distributed in the hope that it will be useful, |
---|
| 54 | % but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
| 55 | % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
| 56 | % GNU General Public License (see LICENSE.txt) for more details. |
---|
| 57 | %======================================================================= |
---|
| 58 | |
---|
| 59 | function ParamOut=particle_detect(Param) |
---|
| 60 | |
---|
| 61 | %% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed |
---|
| 62 | if isstruct(Param) && isequal(Param.Action.RUN,0) |
---|
| 63 | ParamOut.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default) |
---|
| 64 | ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default) |
---|
| 65 | ParamOut.NbSlice='off'; %nbre of slices ('off' by default) |
---|
| 66 | ParamOut.VelType='one';% menu for selecting the velocity type (options 'off'/'one'/'two', 'off' by default) |
---|
| 67 | ParamOut.FieldName='one';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default) |
---|
| 68 | ParamOut.FieldTransform = 'off';%can use a transform function |
---|
| 69 | ParamOut.TransformPath=fullfile(fileparts(which('uvmat')),'transform_field');% path to transform functions (needed for compilation only) |
---|
| 70 | ParamOut.ProjObject='off';%can use projection object(option 'off'/'on', |
---|
| 71 | ParamOut.Mask='on';%can use mask option (option 'off'/'on', 'off' by default) |
---|
| 72 | ParamOut.OutputDirExt='.detect';%set the output dir extension |
---|
| 73 | 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 |
---|
| 74 | %check the input files |
---|
| 75 | ParamOut.CheckOverwriteVisible='on'; % manage the overwrite of existing files (default=1) |
---|
| 76 | first_j=[]; |
---|
| 77 | if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end |
---|
| 78 | PairString=''; |
---|
| 79 | if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end |
---|
| 80 | [i1,i2,j1,j2] = get_file_index(Param.IndexRange.first_i,first_j,PairString); |
---|
| 81 | FirstFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},... |
---|
| 82 | Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2); |
---|
| 83 | if ~exist(FirstFileName,'file') |
---|
| 84 | msgbox_uvmat('WARNING',['the first input file ' FirstFileName ' does not exist']) |
---|
| 85 | end |
---|
| 86 | |
---|
| 87 | prompt = {'threshold(th)';... |
---|
| 88 | 'particle size (sz)' }; |
---|
| 89 | dlg_title = 'get processing parameters'; |
---|
| 90 | num_lines= 2; |
---|
| 91 | def = {'4000';'3'}; |
---|
| 92 | answer = inputdlg(prompt,dlg_title,num_lines,def); |
---|
| 93 | if isempty(answer) |
---|
| 94 | return |
---|
| 95 | end |
---|
| 96 | %check input consistency |
---|
| 97 | ParamOut.ActionInput.th=str2num(answer{1}); |
---|
| 98 | ParamOut.ActionInput.sz=str2num(answer{2}); |
---|
| 99 | return |
---|
| 100 | end |
---|
| 101 | |
---|
| 102 | %%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%% |
---|
| 103 | ParamOut=[]; %default output |
---|
| 104 | RUNHandle=[]; |
---|
| 105 | WaitbarHandle=[]; |
---|
| 106 | %% read input parameters from an xml file if input is a file name (batch mode) |
---|
| 107 | checkrun=1; |
---|
| 108 | if ischar(Param) |
---|
| 109 | Param=xml2struct(Param);% read Param as input file (batch case) |
---|
| 110 | checkrun=0; |
---|
| 111 | else |
---|
| 112 | hseries=findobj(allchild(0),'Tag','series'); |
---|
| 113 | RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series |
---|
| 114 | WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series |
---|
| 115 | end |
---|
| 116 | |
---|
| 117 | %% define the directory for result file (with path=RootPath{1}) |
---|
| 118 | OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files |
---|
| 119 | |
---|
| 120 | if ~isfield(Param,'InputFields') |
---|
| 121 | Param.InputFields.FieldName=''; |
---|
| 122 | end |
---|
| 123 | |
---|
| 124 | %% root input file(s) name, type and index series |
---|
| 125 | RootPath=Param.InputTable{1,1}; |
---|
| 126 | RootFile=Param.InputTable{1,3}; |
---|
| 127 | SubDir=Param.InputTable{1,2}; |
---|
| 128 | NomType=Param.InputTable{1,4}; |
---|
| 129 | FileExt=Param.InputTable{1,5}; |
---|
| 130 | |
---|
| 131 | hdisp=disp_uvmat('WAITING...','checking the file series',checkrun); |
---|
| 132 | % gives the series of input file names and indices set by the input parameters: |
---|
| 133 | %[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param); |
---|
| 134 | % filecell{iview,fileindex}: |
---|
| 135 | % iview: line in the table corresponding to a given file series |
---|
| 136 | % fileindex: file index with i and j reshaped as a 1D array |
---|
| 137 | % 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 |
---|
| 138 | % i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices |
---|
| 139 | if ~isempty(hdisp),delete(hdisp),end;%end the waiting display |
---|
| 140 | |
---|
| 141 | |
---|
| 142 | %% determine the file type on each line from the first input file |
---|
| 143 | |
---|
| 144 | |
---|
| 145 | %%%%%%%%%%%% END STANDARD PART %%%%%%%%%%%% |
---|
| 146 | % EDIT FROM HERE |
---|
| 147 | |
---|
| 148 | RootFileOut=RootFile; |
---|
| 149 | |
---|
| 150 | |
---|
| 151 | %% MAIN LOOP ON FIELDS |
---|
| 152 | %%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%% |
---|
| 153 | % for i_slice=1:NbSlice |
---|
| 154 | % index_slice=i_slice:NbSlice:NbField;% select file indices of the slice |
---|
| 155 | % NbFiles=0; |
---|
| 156 | % nbmissing=0; |
---|
| 157 | |
---|
| 158 | %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%% |
---|
| 159 | tstart=tic; %used to record the computing time |
---|
| 160 | CheckOverwrite=1;%default |
---|
| 161 | if isfield(Param,'CheckOverwrite') |
---|
| 162 | CheckOverwrite=Param.CheckOverwrite; |
---|
| 163 | end |
---|
| 164 | %%%%%% INPUT %%%%%% |
---|
| 165 | th=Param.ActionInput.th |
---|
| 166 | sz=Param.ActionInput.sz |
---|
| 167 | % th=4000;%threshold on image intensity |
---|
| 168 | % sz=3; %size of particles |
---|
| 169 | |
---|
| 170 | NbImage=Param.IndexRange.last_i-Param.IndexRange.first_i+1; |
---|
| 171 | incr_i=Param.IndexRange.incr_i; |
---|
| 172 | NbBlock=floor(NbImage/incr_i); |
---|
| 173 | |
---|
| 174 | %% MAIN LOOP |
---|
| 175 | for index1=1:NbBlock |
---|
| 176 | OutputFile=fullfile_uvmat(RootPath,OutputDir,RootFileOut,'.mat','_1-2',(index1-1)*incr_i+1,index1*incr_i) |
---|
| 177 | if ~CheckOverwrite && exist(OutputFile,'file') |
---|
| 178 | disp(['existing output file ' OutputFile ' already exists, skip to next field']) |
---|
| 179 | end |
---|
| 180 | for index=(index1-1)*incr_i+1:index1*incr_i |
---|
| 181 | index |
---|
| 182 | |
---|
| 183 | %% reading input file(s) |
---|
| 184 | ImgName=fullfile_uvmat(RootPath,SubDir,RootFile,FileExt,NomType,index) |
---|
| 185 | Im=imread(ImgName); |
---|
| 186 | [Ny,Nx]=size(Im); |
---|
| 187 | % |
---|
| 188 | % |
---|
| 189 | % frame = kframe+idfile*framesperfiles; |
---|
| 190 | % disp(frame); |
---|
| 191 | % ImgName = sprintf(join(['%s/%s_cam%d_' format],''),folderin, ManipName, CamNum, frame); |
---|
| 192 | % fprintf("%s \n",ImgName); |
---|
| 193 | %Im = (cast(imread(ImgName),'like',Background) - Background).*mask; |
---|
| 194 | %% normalizeimage |
---|
| 195 | Im = double(Im); |
---|
| 196 | |
---|
| 197 | % Imfiltered=filter2(transfert_coef,Im); |
---|
| 198 | % |
---|
| 199 | % Imfiltered(:,1:windowsize)=Imfiltered(:,windowsize)*ones(1,windowsize); |
---|
| 200 | % Imfiltered(:,end-windowsize+1:end)=Imfiltered(:,end-windowsize+1)*ones(1,windowsize); |
---|
| 201 | % Imfiltered(1:windowsize,:)=ones(windowsize,1)*Imfiltered(windowsize,:); |
---|
| 202 | % Imfiltered(end-windowsize+1:end,:)=ones(windowsize,1)*Imfiltered(end-windowsize,:); |
---|
| 203 | % Im=Im./(2*Imfiltered); |
---|
| 204 | %% |
---|
| 205 | |
---|
| 206 | out=pkfnd(Im,th,sz); % Provides intensity maxima positionsth,sz,Test,BackgroundType,format |
---|
| 207 | npar = size(out,1); |
---|
| 208 | |
---|
| 209 | %% We keep only spots with a gaussian shape |
---|
| 210 | cnt = 0; |
---|
| 211 | x = []; |
---|
| 212 | y = []; |
---|
| 213 | for j = 1:npar |
---|
| 214 | Nwidth = 1; |
---|
| 215 | if (out(j,2)-Nwidth >0)&&(out(j,1)-Nwidth>0)&&(out(j,2)+Nwidth<Ny)&&(out(j,1)+Nwidth<Nx) |
---|
| 216 | cnt = cnt+1; |
---|
| 217 | |
---|
| 218 | Ip = double(Im(out(j,2)-Nwidth:out(j,2)+Nwidth,out(j,1)-Nwidth:out(j,1)+Nwidth)); |
---|
| 219 | |
---|
| 220 | x(end+1) = out(j, 1) + 0.5*log(Ip(2,3)/Ip(2,1))/(log((Ip(2,2)*Ip(2,2))/(Ip(2,1)*Ip(2,3)))); |
---|
| 221 | y(end+1) = out(j, 2) + 0.5*log(Ip(3,2)/Ip(1,2))/(log((Ip(2,2)*Ip(2,2))/(Ip(1,2)*Ip(3,2)))); |
---|
| 222 | end |
---|
| 223 | end |
---|
| 224 | |
---|
| 225 | CC(index).X=x; |
---|
| 226 | CC(index).Y=y; |
---|
| 227 | end |
---|
| 228 | |
---|
| 229 | %% Centers saving into a .mat file |
---|
| 230 | if Param.IndexRange.last_i>1 |
---|
| 231 | savefile=OutputFile; |
---|
| 232 | save(savefile,"CC",'-v7.3') |
---|
| 233 | m = matfile(savefile,'Writable',true); |
---|
| 234 | m.nframes = incr_i; |
---|
| 235 | else |
---|
| 236 | figure("NumberTitle","Off","Name",['RAW picture,' SubDir]) |
---|
| 237 | imshow(imread(ImgName),[0,5000]) |
---|
| 238 | colormap gray |
---|
| 239 | % figure("NumberTitle","Off","Name",sprintf("%s, cam %d",BackgroundType,CamNum)) |
---|
| 240 | % imshow(BackgroundMin,[0,5000]) |
---|
| 241 | % colormap gray |
---|
| 242 | % figure("NumberTitle","Off","Name",sprintf("RAW picture - Background, cam %d, frame %d",CamNum,kframe)) |
---|
| 243 | % imshow(Im,[0,th]) |
---|
| 244 | % colormap gray |
---|
| 245 | colorbar |
---|
| 246 | |
---|
| 247 | %% TracÃ© de l'histogramme des intensitÃ©s pour dÃ©finir le seuil |
---|
| 248 | fig = figure('NumberTitle','Off','Name','Intensity histogram'); |
---|
| 249 | histogram(Im,1000) |
---|
| 250 | xlabel("Intensity") |
---|
| 251 | ylabel("Number") |
---|
| 252 | set(gca, 'XScale', 'log') |
---|
| 253 | set(gca, 'YScale', 'log') |
---|
| 254 | |
---|
| 255 | Nx = size(Im,2); |
---|
| 256 | Ny = size(Im,1); |
---|
| 257 | |
---|
| 258 | out=pkfnd(Im,th,sz); % Provides intensity maxima positions |
---|
| 259 | npar = size(out,1); |
---|
| 260 | |
---|
| 261 | %% We keep only spots with a gaussian shape |
---|
| 262 | cnt = 0; |
---|
| 263 | x = []; |
---|
| 264 | y = []; |
---|
| 265 | for j = 1:npar |
---|
| 266 | Nwidth = 1; |
---|
| 267 | if (out(j,2)-Nwidth >0)&&(out(j,1)-Nwidth>0)&&(out(j,2)+Nwidth<Ny)&&(out(j,1)+Nwidth<Nx) |
---|
| 268 | cnt = cnt+1; |
---|
| 269 | |
---|
| 270 | Ip = double(Im(out(j,2)-Nwidth:out(j,2)+Nwidth,out(j,1)-Nwidth:out(j,1)+Nwidth)); |
---|
| 271 | |
---|
| 272 | x(end+1) = out(j, 1) + 0.5*log(Ip(2,3)/Ip(2,1))/(log((Ip(2,2)*Ip(2,2))/(Ip(2,1)*Ip(2,3)))); |
---|
| 273 | y(end+1) = out(j, 2) + 0.5*log(Ip(3,2)/Ip(1,2))/(log((Ip(2,2)*Ip(2,2))/(Ip(1,2)*Ip(3,2)))); |
---|
| 274 | end |
---|
| 275 | end |
---|
| 276 | CC(1).X=x; |
---|
| 277 | CC(1).Y=y; |
---|
| 278 | |
---|
| 279 | fprintf("%d treated \n",1) |
---|
| 280 | |
---|
| 281 | %% Let's plot picture and detected points on a graph !!! Be careful the vertical axis is reversed compared to reality !!! |
---|
| 282 | figure('NumberTitle','Off','Name',sprintf("frame %d, %d detected points",1,numel(x))) |
---|
| 283 | imshow(Im,[0,th]) |
---|
| 284 | colormap gray |
---|
| 285 | |
---|
| 286 | hold on |
---|
| 287 | plot(flip(x),flip(y),'r+') |
---|
| 288 | end |
---|
| 289 | end |
---|
| 290 | |
---|
| 291 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 292 | function CC = CenterFinding2D(session,ManipName,CamNum,firstFrame,nframes,th,sz,Test,BackgroundType,format,withmask) |
---|
| 293 | %%% Detect particles position in picture and provides their positions in |
---|
| 294 | %%% px. |
---|
| 295 | %-------------------------------------------------------------------------------- |
---|
| 296 | %%% Parameters : |
---|
| 297 | %%% session : Path to the achitecture root (2 fields: session.input_path |
---|
| 298 | % and session.output_path) |
---|
| 299 | %%% ManipName : Name of the folder experiment |
---|
| 300 | %%% NumCam : number of the camera studied |
---|
| 301 | %%% nframes : total number of pictures |
---|
| 302 | %%% th : threshold |
---|
| 303 | %%% sz : typical size of the particles |
---|
| 304 | %%% Test : true-> test mode, false-> classic mode (optional) |
---|
| 305 | %%% BackgroundType (optional) : determine which background is substracted to pictures. By defaut is equal to BackgroundMean, |
---|
| 306 | %%% format (optional) : picture names. By defaut it is '%05d.tif'. |
---|
| 307 | %%% The beginning of picture names has to be %ManipName_cam%CamNum_%format |
---|
| 308 | %-------------------------------------------------------------------------------- |
---|
| 309 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 310 | close all |
---|
| 311 | |
---|
| 312 | framesperfiles = 1000; |
---|
| 313 | |
---|
| 314 | %% Test if Test exist or not |
---|
| 315 | if ~exist('Test','var') |
---|
| 316 | Test=false; |
---|
| 317 | end |
---|
| 318 | |
---|
| 319 | %% Test if BackgroundType exist or not |
---|
| 320 | if ~exist('BackgroundType','var') |
---|
| 321 | BackgroundType="BackgroundMean"; |
---|
| 322 | end |
---|
| 323 | |
---|
| 324 | % By defaut format='%05d.tif' |
---|
| 325 | if ~exist('format','var') |
---|
| 326 | format='%05d.tif'; |
---|
| 327 | end |
---|
| 328 | |
---|
| 329 | if ~exist('withmask','var') |
---|
| 330 | withmask=false; |
---|
| 331 | end |
---|
| 332 | |
---|
| 333 | |
---|
| 334 | |
---|
| 335 | %% Definition of folders |
---|
| 336 | fprintf(ManipName); |
---|
| 337 | folderin = sprintf("%sDATA/%s/cam%d",session.input_path,ManipName,CamNum) |
---|
| 338 | folderout = sprintf("%sProcessed_DATA/%s",session.output_path,ManipName) |
---|
| 339 | BackgroundFile = sprintf("%s/Background_cam%d.mat",folderout,CamNum); |
---|
| 340 | |
---|
| 341 | %% Find centers |
---|
| 342 | if exist(folderout,'dir')==0 |
---|
| 343 | mkdir(char(folderout)); |
---|
| 344 | end |
---|
| 345 | |
---|
| 346 | % if exist(strcat(folderout,'/Parallel/Matching'),'dir')==0 |
---|
| 347 | % mkdir(char(strcat(folderout,'/Parallel/Matching'))); |
---|
| 348 | % end |
---|
| 349 | |
---|
| 350 | load(BackgroundFile,'BackgroundMin','BackgroundMax','BackgroundMean') |
---|
| 351 | %% Choice of background type |
---|
| 352 | if BackgroundType=="BackgroundMean" |
---|
| 353 | Background=BackgroundMean; |
---|
| 354 | elseif BackgroundType=="BackgroundMax" |
---|
| 355 | Background=BackgroundMax; |
---|
| 356 | elseif BackgroundType=="BackgroundMin" |
---|
| 357 | Background=BackgroundMin; |
---|
| 358 | end |
---|
| 359 | Nx = size(Background,2); |
---|
| 360 | Ny = size(Background,1); |
---|
| 361 | |
---|
| 362 | if withmask |
---|
| 363 | ImgName = sprintf(join(['%s/%s_cam%d_' format],''),folderin, ManipName, CamNum, 1); |
---|
| 364 | Im = imread(ImgName); |
---|
| 365 | imshow(Im,[0,10000]) |
---|
| 366 | BW = roipoly; |
---|
| 367 | mask = cast(BW,'like',Background); |
---|
| 368 | else |
---|
| 369 | mask = cast(ones(Ny,Nx),'like',Background); |
---|
| 370 | end |
---|
| 371 | |
---|
| 372 | %% param for normailze image |
---|
| 373 | windowsize=round(min(Ny,Nx)/20); |
---|
| 374 | transfert_coef = ones(windowsize)/windowsize^2; |
---|
| 375 | %% |
---|
| 376 | |
---|
| 377 | close all |
---|
| 378 | if ~Test |
---|
| 379 | idfirstfile = fix(firstFrame/framesperfiles); |
---|
| 380 | idlastfile = fix((nframes-1)/framesperfiles); |
---|
| 381 | for idfile = idfirstfile:idlastfile |
---|
| 382 | firstFrame = rem(firstFrame,framesperfiles); |
---|
| 383 | lastframe = min(framesperfiles,nframes-idfile*framesperfiles); |
---|
| 384 | for kframe=firstFrame:lastframe |
---|
| 385 | frame = kframe+idfile*framesperfiles; |
---|
| 386 | disp(frame); |
---|
| 387 | ImgName = sprintf(join(['%s/%s_cam%d_' format],''),folderin, ManipName, CamNum, frame); |
---|
| 388 | fprintf("%s \n",ImgName); |
---|
| 389 | Im = (cast(imread(ImgName),'like',Background) - Background).*mask; |
---|
| 390 | %% normalizeimage |
---|
| 391 | Im = double(Im); |
---|
| 392 | |
---|
| 393 | Imfiltered=filter2(transfert_coef,Im); |
---|
| 394 | |
---|
| 395 | Imfiltered(:,1:windowsize)=Imfiltered(:,windowsize)*ones(1,windowsize); |
---|
| 396 | Imfiltered(:,end-windowsize+1:end)=Imfiltered(:,end-windowsize+1)*ones(1,windowsize); |
---|
| 397 | Imfiltered(1:windowsize,:)=ones(windowsize,1)*Imfiltered(windowsize,:); |
---|
| 398 | Imfiltered(end-windowsize+1:end,:)=ones(windowsize,1)*Imfiltered(end-windowsize,:); |
---|
| 399 | Im=Im./(2*Imfiltered); |
---|
| 400 | %% |
---|
| 401 | |
---|
| 402 | |
---|
| 403 | out=pkfnd(Im,th,sz); % Provides intensity maxima positions |
---|
| 404 | npar = size(out,1); |
---|
| 405 | |
---|
| 406 | %% We keep only spots with a gaussian shape |
---|
| 407 | cnt = 0; |
---|
| 408 | x = []; |
---|
| 409 | y = []; |
---|
| 410 | for j = 1:npar |
---|
| 411 | Nwidth = 1; |
---|
| 412 | if (out(j,2)-Nwidth >0)&&(out(j,1)-Nwidth>0)&&(out(j,2)+Nwidth<Ny)&&(out(j,1)+Nwidth<Nx) |
---|
| 413 | cnt = cnt+1; |
---|
| 414 | |
---|
| 415 | Ip = double(Im(out(j,2)-Nwidth:out(j,2)+Nwidth,out(j,1)-Nwidth:out(j,1)+Nwidth)); |
---|
| 416 | |
---|
| 417 | x(end+1) = out(j, 1) + 0.5*log(Ip(2,3)/Ip(2,1))/(log((Ip(2,2)*Ip(2,2))/(Ip(2,1)*Ip(2,3)))); |
---|
| 418 | y(end+1) = out(j, 2) + 0.5*log(Ip(3,2)/Ip(1,2))/(log((Ip(2,2)*Ip(2,2))/(Ip(1,2)*Ip(3,2)))); |
---|
| 419 | end |
---|
| 420 | end |
---|
| 421 | |
---|
| 422 | CC(kframe).X=x; |
---|
| 423 | CC(kframe).Y=y; |
---|
| 424 | end |
---|
| 425 | %% Centers saving into a .mat file |
---|
| 426 | firstFramefile = framesperfiles*idfile+1; |
---|
| 427 | lastFramefile = framesperfiles*(idfile+1); |
---|
| 428 | savefile = sprintf(['%s/Parallel/Matching/centers_cam%d_',format(1:end-4),'-',format(1:end-4),'.mat'],folderout,CamNum,firstFramefile,lastFramefile); |
---|
| 429 | if exist(savefile,'file') |
---|
| 430 | m = matfile(savefile,'Writable',true); |
---|
| 431 | m.CC(1,firstFrame:lastframe) = CC(firstFrame:lastframe); |
---|
| 432 | m.nframes = framesperfiles; |
---|
| 433 | else |
---|
| 434 | save(savefile,"CC","nframes",'-v7.3') |
---|
| 435 | m = matfile(savefile,'Writable',true); |
---|
| 436 | m.nframes = framesperfiles; |
---|
| 437 | end |
---|
| 438 | firstFrame = 1; |
---|
| 439 | end |
---|
| 440 | else |
---|
| 441 | kframe=1 |
---|
| 442 | ImgName = sprintf(join(['%s/%s_cam%d_' format],''),folderin, ManipName, CamNum, kframe); |
---|
| 443 | fprintf("%s \n",ImgName); |
---|
| 444 | Im = (cast(imread(ImgName),'like',Background) - Background).*mask; |
---|
| 445 | %% normalizeimage |
---|
| 446 | Im = double(Im); |
---|
| 447 | |
---|
| 448 | Imfiltered=filter2(transfert_coef,Im); |
---|
| 449 | |
---|
| 450 | Imfiltered(:,1:windowsize)=Imfiltered(:,windowsize)*ones(1,windowsize); |
---|
| 451 | Imfiltered(:,end-windowsize+1:end)=Imfiltered(:,end-windowsize+1)*ones(1,windowsize); |
---|
| 452 | Imfiltered(1:windowsize,:)=ones(windowsize,1)*Imfiltered(windowsize,:); |
---|
| 453 | Imfiltered(end-windowsize+1:end,:)=ones(windowsize,1)*Imfiltered(end-windowsize,:); |
---|
| 454 | Im=Im./(2*Imfiltered); |
---|
| 455 | %% |
---|
| 456 | figure("NumberTitle","Off","Name",sprintf("RAW picture, cam %d, frame %d",CamNum,kframe)) |
---|
| 457 | imshow(imread(ImgName),[0,5000]) |
---|
| 458 | colormap gray |
---|
| 459 | figure("NumberTitle","Off","Name",sprintf("%s, cam %d",BackgroundType,CamNum)) |
---|
| 460 | imshow(BackgroundMin,[0,5000]) |
---|
| 461 | colormap gray |
---|
| 462 | figure("NumberTitle","Off","Name",sprintf("RAW picture - Background, cam %d, frame %d",CamNum,kframe)) |
---|
| 463 | imshow(Im,[0,th]) |
---|
| 464 | colormap gray |
---|
| 465 | colorbar |
---|
| 466 | |
---|
| 467 | % if exist('Erosion','var') |
---|
| 468 | % se = strel('disk',1); |
---|
| 469 | % Imerode = imerode(Im,se); |
---|
| 470 | % Imdilate = imdilate(Imerode,se); |
---|
| 471 | % figure() |
---|
| 472 | % imagesc(Imdilate) |
---|
| 473 | % axis image |
---|
| 474 | % colormap gray |
---|
| 475 | % |
---|
| 476 | % Image = Im; |
---|
| 477 | % Im = Imdilate; |
---|
| 478 | % end |
---|
| 479 | |
---|
| 480 | %% TracÃ© de l'histogramme des intensitÃ©s pour dÃ©finir le seuil |
---|
| 481 | fig = figure('NumberTitle','Off','Name','Intensity histogram'); |
---|
| 482 | histogram(Im,1000) |
---|
| 483 | xlabel("Intensity") |
---|
| 484 | ylabel("Number") |
---|
| 485 | set(gca, 'XScale', 'log') |
---|
| 486 | set(gca, 'YScale', 'log') |
---|
| 487 | |
---|
| 488 | Nx = size(Im,2); |
---|
| 489 | Ny = size(Im,1); |
---|
| 490 | |
---|
| 491 | out=pkfnd(Im,th,sz); % Provides intensity maxima positions |
---|
| 492 | npar = size(out,1); |
---|
| 493 | |
---|
| 494 | %% We keep only spots with a gaussian shape |
---|
| 495 | cnt = 0; |
---|
| 496 | x = []; |
---|
| 497 | y = []; |
---|
| 498 | for j = 1:npar |
---|
| 499 | Nwidth = 1; |
---|
| 500 | if (out(j,2)-Nwidth >0)&&(out(j,1)-Nwidth>0)&&(out(j,2)+Nwidth<Ny)&&(out(j,1)+Nwidth<Nx) |
---|
| 501 | cnt = cnt+1; |
---|
| 502 | |
---|
| 503 | Ip = double(Im(out(j,2)-Nwidth:out(j,2)+Nwidth,out(j,1)-Nwidth:out(j,1)+Nwidth)); |
---|
| 504 | |
---|
| 505 | x(end+1) = out(j, 1) + 0.5*log(Ip(2,3)/Ip(2,1))/(log((Ip(2,2)*Ip(2,2))/(Ip(2,1)*Ip(2,3)))); |
---|
| 506 | y(end+1) = out(j, 2) + 0.5*log(Ip(3,2)/Ip(1,2))/(log((Ip(2,2)*Ip(2,2))/(Ip(1,2)*Ip(3,2)))); |
---|
| 507 | end |
---|
| 508 | end |
---|
| 509 | CC(kframe).X=x; |
---|
| 510 | CC(kframe).Y=y; |
---|
| 511 | |
---|
| 512 | fprintf("%d treated \n",kframe) |
---|
| 513 | |
---|
| 514 | %% Let's plot picture and detected points on a graph !!! Be careful the vertical axis is reversed compared to reality !!! |
---|
| 515 | figure('NumberTitle','Off','Name',sprintf("frame %d, %d detected points",kframe,numel(x))) |
---|
| 516 | imshow(Im,[0,th]) |
---|
| 517 | colormap gray |
---|
| 518 | |
---|
| 519 | hold on |
---|
| 520 | plot(flip(x),flip(y),'r+') |
---|
| 521 | end |
---|
| 522 | |
---|
| 523 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 524 | function out=pkfnd(im,th,sz) |
---|
| 525 | % finds local maxima in an image to pixel level accuracy. |
---|
| 526 | % this provides a rough guess of particle |
---|
| 527 | % centers to be used by cntrd.m. Inspired by the lmx subroutine of Grier |
---|
| 528 | % and Crocker's feature.pro |
---|
| 529 | %-------------------------------------------------------------------------------- |
---|
| 530 | % INPUTS: |
---|
| 531 | % im : image to process, particle should be bright spots on dark background with little noise |
---|
| 532 | % ofen an bandpass filtered brightfield image (fbps.m, fflt.m or bpass.m) or a nice |
---|
| 533 | % fluorescent image |
---|
| 534 | % th : the minimum brightness of a pixel that might be local maxima. |
---|
| 535 | % (NOTE: Make it big and the code runs faster |
---|
| 536 | % but you might miss some particles. Make it small and you'll get |
---|
| 537 | % everything and it'll be slow.) |
---|
| 538 | % sz : if your data's noisy, (e.g. a single particle has multiple local |
---|
| 539 | % maxima), then set this optional keyword to a value slightly larger than the diameter of your blob. if |
---|
| 540 | % multiple peaks are found withing a radius of sz/2 then the code will keep |
---|
| 541 | % only the brightest. Also gets rid of all peaks within sz of boundary |
---|
| 542 | % OUTPUT: a N x 2 array containing, [row,column] coordinates of local maxima |
---|
| 543 | % out(:,1) are the x-coordinates of the maxima |
---|
| 544 | % out(:,2) are the y-coordinates of the maxima |
---|
| 545 | %-------------------------------------------------------------------------------- |
---|
| 546 | %CREATED: Eric R. Dufresne, Yale University, Feb 4 2005 |
---|
| 547 | %MODIFIED: ERD, 5/2005, got rid of ind2rc.m to reduce overhead on tip by |
---|
| 548 | % Dan Blair; added sz keyword |
---|
| 549 | % ERD, 6/2005: modified to work with one and zero peaks, removed automatic |
---|
| 550 | % normalization of image |
---|
| 551 | % ERD, 6/2005: due to popular demand, altered output to give x and y |
---|
| 552 | % instead of row and column |
---|
| 553 | % ERD, 8/24/2005: pkfnd now exits politely if there's nothing above |
---|
| 554 | % threshold instead of crashing rudely |
---|
| 555 | % ERD, 6/14/2006: now exits politely if no maxima found |
---|
| 556 | % ERD, 10/5/2006: fixed bug that threw away particles with maxima |
---|
| 557 | % consisting of more than two adjacent points |
---|
| 558 | |
---|
| 559 | |
---|
| 560 | |
---|
| 561 | %find all the pixels above threshold |
---|
| 562 | %im=im./max(max(im)); |
---|
| 563 | [nr,nc] = size(im); |
---|
| 564 | [i,j,ind]=find(im > th); |
---|
| 565 | n=length(ind); |
---|
| 566 | |
---|
| 567 | if n==0 |
---|
| 568 | out=[];[i,j,ind]=find(im > th); |
---|
| 569 | fprintf('nothing above threshold'); |
---|
| 570 | return; |
---|
| 571 | end |
---|
| 572 | mx=[]; |
---|
| 573 | %convert index from find to row and column |
---|
| 574 | rc=[i,j]; % j corresponds to x axis and i to y axis |
---|
| 575 | % rc=[j,i]; % j corresponds to x axis and i to y axis |
---|
| 576 | for ii=1:n |
---|
| 577 | r=rc(ii,1); |
---|
| 578 | c=rc(ii,2); |
---|
| 579 | %check each pixel above threshold to see if it's brighter than it's neighbors |
---|
| 580 | % THERE'S GOT TO BE A FASTER WAY OF DOING THIS. I'M CHECKING SOME MULTIPLE TIMES, |
---|
| 581 | % BUT THIS DOESN'T SEEM THAT SLOW COMPARED TO THE OTHER ROUTINES, ANYWAY. |
---|
| 582 | if r>1 && r<nr && c>1 && c<nc |
---|
| 583 | if im(r,c)>=im(r-1,c-1) && im(r,c)>=im(r,c-1) && im(r,c)>=im(r+1,c-1) && ... |
---|
| 584 | im(r,c)>=im(r-1,c) && im(r,c)>=im(r+1,c) && ... |
---|
| 585 | im(r,c)>=im(r-1,c+1) && im(r,c)>=im(r,c+1) && im(r,c)>=im(r+1,c+1) |
---|
| 586 | mx=[mx,[r,c]']; %#ok<AGROW> |
---|
| 587 | %tst(ind(i))=im(ind(i)); |
---|
| 588 | end |
---|
| 589 | end |
---|
| 590 | end |
---|
| 591 | %out=tst; |
---|
| 592 | mx=mx'; |
---|
| 593 | |
---|
| 594 | [npks,crap]=size(mx); |
---|
| 595 | |
---|
| 596 | %if size is specified, then get ride of pks within size of boundary |
---|
| 597 | if nargin==3 && npks>0 |
---|
| 598 | %throw out all pks within sz of boundary; |
---|
| 599 | ind=find(mx(:,1)>sz & mx(:,1)<(nr-sz) & mx(:,2)>sz & mx(:,2)<(nc-sz)); |
---|
| 600 | mx=mx(ind,:); |
---|
| 601 | end |
---|
| 602 | |
---|
| 603 | %prevent from finding peaks within size of each other |
---|
| 604 | [npks,crap]=size(mx); |
---|
| 605 | if npks > 1 |
---|
| 606 | %CREATE AN IMAGE WITH ONLY PEAKS |
---|
| 607 | nmx=npks; |
---|
| 608 | tmp=0.*im; |
---|
| 609 | for i=1:nmx |
---|
| 610 | tmp(mx(i,1),mx(i,2))=im(mx(i,1),mx(i,2)); |
---|
| 611 | end |
---|
| 612 | %LOOK IN NEIGHBORHOOD AROUND EACH PEAK, PICK THE BRIGHTEST |
---|
| 613 | for i=1:nmx |
---|
| 614 | 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))) ; |
---|
| 615 | [mv,indi]=max(roi); |
---|
| 616 | [mv,indj]=max(mv); |
---|
| 617 | 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; |
---|
| 618 | tmp(mx(i,1)-floor(sz/2)+indi(indj)-1,mx(i,2)-floor(sz/2)+indj-1)=mv; |
---|
| 619 | end |
---|
| 620 | ind=find(tmp>0); |
---|
| 621 | mx=[mod(ind,nr),floor(ind/nr)+1]; |
---|
| 622 | end |
---|
| 623 | |
---|
| 624 | if size(mx)==[0,0] |
---|
| 625 | out=[]; |
---|
| 626 | else |
---|
| 627 | out(:,2)=mx(:,1); |
---|
| 628 | out(:,1)=mx(:,2); |
---|
| 629 | end |
---|