Index: trunk/src/civ.m =================================================================== --- trunk/src/civ.m (revision 1167) +++ trunk/src/civ.m (revision 1168) @@ -257,5 +257,5 @@ end ytable=npy_ima-ytable+1;%reverse from j index to image coordinate y -result_conv=result_conv*corrmax/sum_square;% keep the last (not normalised) correlation matrix for output +result_conv=result_conv/sum_square;% keep the last correlation matrix for output Index: trunk/src/filter2_uvmat.m =================================================================== --- trunk/src/filter2_uvmat.m (revision 1168) +++ trunk/src/filter2_uvmat.m (revision 1168) @@ -0,0 +1,19 @@ +% Apply a sliding average with ny*nx points excluding NaN and outerbound values +function Afilt=filter2_omitnan(npx,npy,A) +Afilt=A;%default +npy2=ceil(npy/2); +npx2=ceil(npx/2); +for iline=1:size(A,2) + for icolumn=1:size(A,1) + lowboundy=max(icolumn-npy2,1); + trunc_y1=lowboundy-icolumn+npy2;%nbre of elements removed low + highboundy=min(icolumn+npy2,size(A,1)); + trunc_y2=icolumn+npy2-highboundy;%nbre of elements removed high + lowboundx=max(iline-npx2,1); + trunc_x1=lowboundx-iline+npx2;%nbre of elements removed low + highboundx=min(iline+npx2,size(A,2)); + trunc_x2=iline+npx2-highboundx;%nbre of elements removed high + Asub=A(lowboundy:highboundy,lowboundx:highboundx);%nbre of elements removed high + Afilt(icolumn,iline)=mean(Asub,'all','omitnan'); + end +end Index: trunk/src/mouse_motion.m =================================================================== --- trunk/src/mouse_motion.m (revision 1167) +++ trunk/src/mouse_motion.m (revision 1168) @@ -276,14 +276,15 @@ 'LineStyle','-','Tag','PIV_box_marker'); rectangle('Curvature',[0 0],... - 'Position',[xround-isx2+shiftx yround-isy2+shifty 2*isx2 2*isy2],'EdgeColor','m',... + 'Position',[xround-isx2-ibx2+shiftx yround-isy2-iby2+shifty 2*isx2+2*ibx2 2*isy2+2*iby2],'EdgeColor','m',... 'LineStyle','- -','Tag','PIV_search_marker'); else set(hhh,'Position',[xround-ibx2 yround-iby2 2*ibx2 2*iby2]) - set(hhhh,'Position',[xround-isx2+shiftx yround-isy2+shifty 2*isx2 2*isy2]) + set(hhhh,'Position',[xround-isx2-ibx2+shiftx yround-isy2-iby2+shifty 2*isx2+2*ibx2 2*isy2+2*iby2]) end % perform the PIV calculation as the single point xround yround - Param.CheckFix1=0; - Param.CheckPatch1=0;%desactivate all calculations except Civ2 or Civ1 + Param.ActionInput.CheckCiv1=1; +% Param.CheckFix1=0; +% Param.CheckPatch1=0;%desactivate all calculations except Civ2 or Civ1 Param.Action.RUN=1; Param.ActionInput.ListCompareMode='PIV'; @@ -299,8 +300,12 @@ text_displ_5=errormsg; else - rangx(1)=-(isx2-ibx2)+shiftx; - rangx(2)=isx2-ibx2+shiftx; - rangy(1)=-(isy2-iby2)-shifty; - rangy(2)=(isy2-iby2)-shifty; +% rangx(1)=-(isx2-ibx2)+shiftx; +% rangx(2)=isx2-ibx2+shiftx; +% rangy(1)=-(isy2-iby2)-shifty; +% rangy(2)=(isy2-iby2)-shifty; + rangx(1)=-isx2+shiftx; + rangx(2)=isx2+shiftx; + rangy(1)=-isy2-shifty; + rangy(2)=isy2-shifty; hcorr=[]; if isfield(AxeData,'CurrentCorrImage') Index: trunk/src/series.m =================================================================== --- trunk/src/series.m (revision 1167) +++ trunk/src/series.m (revision 1168) @@ -3224,8 +3224,14 @@ FirstFileName=fullfile_uvmat(Param.InputTable{iview,1},Param.InputTable{iview,2},Param.InputTable{iview,3},... Param.InputTable{iview,5},Param.InputTable{iview,4},i1,i2,j1,j2); - Data=nc2struct(FirstFileName); + + [FileInfo,VideoObject]=get_file_info(FirstFileName); + + +% Data=nc2struct(FirstFileName); if isfield(Data,'Civ2_Mask') + if exist(Data.Civ2_Mask,'file') MaskTable{iview,1}=Data.Civ2_Mask; checkmask=1; + end end if ~checkmask Index: trunk/src/series/civ_series.m =================================================================== --- trunk/src/series/civ_series.m (revision 1167) +++ trunk/src/series/civ_series.m (revision 1168) @@ -284,5 +284,6 @@ end end - if CheckInputFile +% if CheckInputFile + if CheckOutputFile OutputPath=fullfile(Param.OutputPath,Param.Experiment,Param.Device); if iview_A==0 % no nc file has been entered @@ -488,5 +489,5 @@ % caluclate velocity data - tstart_civ1=tic; + % tstart_civ1=tic; [Data.Civ1_X,Data.Civ1_Y,Data.Civ1_U,Data.Civ1_V,Data.Civ1_C,Data.Civ1_FF, result_conv, errormsg] = civ (par_civ1); if ~isempty(errormsg) Index: trunk/src/transform_field/FFT2_detrend.m =================================================================== --- trunk/src/transform_field/FFT2_detrend.m (revision 1167) +++ trunk/src/transform_field/FFT2_detrend.m (revision 1168) @@ -1,3 +1,3 @@ -% 'FFT2_detrend': calculate the 2D spectrum of the input scalar after removing the linear trend +% 'FFT2_detrend': calculate the 2D spectrum of the input scalar or image after removing the linear trend %------------------------------------------------------------------------ %%%% Use the general syntax for transform fields with a single input %%%% @@ -34,5 +34,5 @@ end %%%%%%%%%%%%%%%%%%%%%%%%% -[CellInfo,NbDim,errormsg]=find_field_cells(DataIn); +[CellInfo,NbDim,errormsg]=find_field_cells(DataIn);% detect the input fields if ~isempty(errormsg) DataOut.Txt=errormsg; @@ -42,12 +42,19 @@ DataOut.VarDimName={}; for ilist=1:numel(CellInfo) - if NbDim(ilist)==2 && numel(CellInfo{ilist}.CoordIndex)==2 % field with structured coordinates + if NbDim(ilist)==2 && numel(CellInfo{ilist}.CoordIndex)==2 % field with structured 2D coordinates + %process scalar + ivar=CellInfo{ilist}.VarIndex_scalar(1); + VarName=DataIn.ListVarName{ivar};% name o the scalar field + z=double(DataIn.(VarName));% transform integer input to double real + z(isnan(z))=0;% set NaN values to 0 + [npy,npx]=size(z); + %process coordinates CoordName=DataIn.ListVarName(CellInfo{ilist}.CoordIndex); - x1 = DataIn.(CoordName{2}); y1 = DataIn.(CoordName{1}); - [x y] = meshgrid(x1,y1); - coeff(1,1) = sum(sum(x.^2)); coeff(1,2) = sum(sum(x.*y)); coeff(1,3) = sum(sum(x)); - coeff(2,1) = sum(sum(x.*y)); coeff(2,2) = sum(sum(y.^2)); coeff(2,3) = sum(sum(y)); - coeff(3,1) = sum(sum(x)); coeff(3,2) = sum(sum(y)); coeff(3,3) = length(x1)*length(y1); + x1 = linspace(DataIn.(CoordName{2})(1),DataIn.(CoordName{2})(end),npx); + y1 =linspace(DataIn.(CoordName{1})(1),DataIn.(CoordName{1})(end),npy); + [x,y] = meshgrid(x1,y1);% 2D grid of coordinates + + % prepare the grid of wave vectors delta_x = x1(2) - x1(1); delta_y = y1(2) - y1(1); Nx = length(x1); Ny = length(y1); @@ -55,32 +62,34 @@ ssx = find(Nxa=0); ky1 = ky1(ss); - [kx ky] = meshgrid(kx1,ky1); + [kx,ky] = meshgrid(kx1,ky1); DataOut.(CoordName{2}) = kx1; DataOut.(CoordName{1}) = ky1; if isfield(DataIn,'CoordUnit') - DataOut.CoordUnit=[DataIn.CoordUnit '^{-1}']; + DataOut.CoordUnit=[DataIn.CoordUnit '^{-1}'];%unit of wave vectors DataOut.ListGlobalAttribute={'CoordUnit'}; end - %process scalar - ivar=CellInfo{ilist}.VarIndex_scalar(1); - VarName=DataIn.ListVarName{ivar}; - z=DataIn.(VarName); - z(isnan(z))=0;% set to 0 NaN values + + % get the coeff to calulate the linear trend + coeff(1,1) = sum(sum(x.^2)); coeff(1,2) = sum(sum(x.*y)); coeff(1,3) = sum(sum(x)); + coeff(2,1) = sum(sum(x.*y)); coeff(2,2) = sum(sum(y.^2)); coeff(2,3) = sum(sum(y)); + coeff(3,1) = sum(sum(x)); coeff(3,2) = sum(sum(y)); coeff(3,3) = length(x1)*length(y1); rhs(1) = sum(sum(x.*z)); rhs(2) = sum(sum(y.*z)); rhs(3) = sum(sum(z)); - lin_coeff = inv(coeff)*rhs'; + % lin_coeff = inv(coeff)*rhs'; + lin_coeff = coeff\rhs'; lin_trend = lin_coeff(1)*x + lin_coeff(2)*y + lin_coeff(3); - z2 = z - lin_trend; - spec2 = abs(fft2(z2)).^2; + z2 = z - lin_trend;% substract the linear trend to the input field + + spec2 = abs(fft2(z2)).^2;% get spectrum as squared of fft modulus spec2 = spec2(Iy,Ix); spec2 = spec2(ss,:); - DataOut.(VarName) = log10(spec2); + DataOut.(VarName) = log10(spec2);% take the log10 of spectrum %DataOut.(VarName) = spec2; - spec_sum=sum(sum(spec2)); - kx_mean=sum(sum(spec2.*kx))/spec_sum; - ky_mean=sum(sum(spec2.*ky))/spec_sum; - theta=atand(ky_mean/kx_mean); - lambda=2*pi/(sqrt(kx_mean*kx_mean+ky_mean*ky_mean)); +% spec_sum=sum(sum(spec2)); +% kx_mean=sum(sum(spec2.*kx))/spec_sum; +% ky_mean=sum(sum(spec2.*ky))/spec_sum; +% theta=atand(ky_mean/kx_mean); +% lambda=2*pi/(sqrt(kx_mean*kx_mean+ky_mean*ky_mean)); DataOut.ListVarName=[CoordName {VarName}];%list of variables Index: trunk/src/transform_field/filter2_omitnan.m =================================================================== --- trunk/src/transform_field/filter2_omitnan.m (revision 1168) +++ trunk/src/transform_field/filter2_omitnan.m (revision 1168) @@ -0,0 +1,92 @@ +% 'ima_filter': low-pass filter of an image or other 2D fields defined on a regular grid +% the size of the filtering window in x and y is interactivement defined + +%------------------------------------------------------------------------ +%%%% Use the general syntax for transform fields with a single input and parameters %%%% +% OUTPUT: +% DataOut: output field structure +% +%INPUT: +% DataIn: input field structure +% +% Param: matlab structure whose field Param.TransformInput contains the filter parameters +% DataIn_1: variables possibly introduced as a second input field +%----------------------------------- + +%======================================================================= +% Copyright 2008-2024, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France +% http://www.legi.grenoble-inp.fr +% Joel.Sommeria - Joel.Sommeria (A) univ-grenoble-alpes.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 DataOut=ima_filter(DataIn,Param,DataIn_1) + +%% request input parameters +if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0) + prompt = {'npx';'npy'}; + dlg_title = 'get the filter size in x and y'; + num_lines= 2; + def = { '21';'21'}; + if isfield(Param,'TransformInput')&&isfield(Param.TransformInput,'FilterBoxSize_x')&&... + isfield(Param.TransformInput,'FilterBoxSize_y') + def={num2str(Param.TransformInput.FilterBoxSize_x);num2str(Param.TransformInput.FilterBoxSize_y)}; + end + answer = inputdlg(prompt,dlg_title,num_lines,def); + DataOut.TransformInput.FilterBoxSize_x=str2num(answer{1}); %size of the filtering window + DataOut.TransformInput.FilterBoxSize_y=str2num(answer{2}); %size of the filtering window + return +end + +DataOut=DataIn; %default + +%definition of the cos shape matrix filter +ix=1/2-Param.TransformInput.FilterBoxSize_x/2:-1/2+Param.TransformInput.FilterBoxSize_x/2;% +iy=1/2-Param.TransformInput.FilterBoxSize_y/2:-1/2+Param.TransformInput.FilterBoxSize_y/2;% +%del=np/3; +%fct=exp(-(ix/del).^2); +fct2_x=cos(ix/((Param.TransformInput.FilterBoxSize_x-1)/2)*pi/2); +fct2_y=cos(iy/((Param.TransformInput.FilterBoxSize_y-1)/2)*pi/2); +%Mfiltre=(ones(5,5)/5^2); +Mfiltre=fct2_y'*fct2_x; +Mfiltre=Mfiltre/(sum(sum(Mfiltre)));%normalize filter + +[CellInfo,NbDim,errormsg]=find_field_cells(DataIn) +for icell=1:numel(CellInfo) + if isfield(CellInfo{icell},'CoordType')&& strcmp(CellInfo{icell}.CoordType,'grid') + for ivar=1:numel(CellInfo{icell}.VarIndex) + VarName=DataIn.ListVarName{CellInfo{icell}.VarIndex(ivar)}; + Atype=class(DataIn.(VarName));% detect integer 8 or 16 bits + + % DataOut.(VarName)=filter2(Mfiltre,DataIn.(VarName)); + DataOut.(VarName)=filter2_uvmat(Param.TransformInput.FilterBoxSize_x,Param.TransformInput.FilterBoxSize_y,DataIn.(VarName)); + DataOut.(VarName)=feval(Atype,DataOut.(VarName));%transform to the initial image format + end + end +end +if exist('DataIn_1','var') + [CellInfo,NbDim,errormsg]=find_field_cells(DataIn_1); + for icell=1:numel(CellInfo) + if isfield(CellInfo{icell},'CoordType')&& strcmp(CellInfo{icell}.CoordType,'grid') + for ivar=1:numel(CellInfo{icell}.VarIndex) + VarName=DataIn_1.ListVarName{CellInfo{icell}.VarIndex(ivar)}; + Atype=class(DataIn_1.(VarName));% detect integer 8 or 16 bits + DataOut.(VarName)=filter2_uvmat(Param.TransformInput.FilterBoxSize_x,Param.TransformInput.FilterBoxSize_y,DataIn_1.(VarName)); + DataOut.(VarName)=feval(Atype,DataOut.(VarName));%transform to the initial image format + end + end + end +end + + Index: trunk/src/transform_field/ima_filter.m =================================================================== --- trunk/src/transform_field/ima_filter.m (revision 1167) +++ trunk/src/transform_field/ima_filter.m (revision 1168) @@ -39,5 +39,5 @@ dlg_title = 'get the filter size in x and y'; num_lines= 2; - def = { '20';'20'}; + def = { '21';'21'}; if isfield(Param,'TransformInput')&&isfield(Param.TransformInput,'FilterBoxSize_x')&&... isfield(Param.TransformInput,'FilterBoxSize_y') @@ -74,4 +74,5 @@ else DataOut.(VarName)=filter2(Mfiltre,DataIn.(VarName)); + %DataOut.(VarName)=filter2_uvmat(Param.TransformInput.FilterBoxSize_x,Param.TransformInput.FilterBoxSize_y,DataIn.(VarName)); DataOut.(VarName)=feval(Atype,DataOut.(VarName));%transform to the initial image format end Index: trunk/src/uvmat.m =================================================================== --- trunk/src/uvmat.m (revision 1167) +++ trunk/src/uvmat.m (revision 1168) @@ -2294,5 +2294,9 @@ set(handles.j2,'String',num2stra(j2,NomType)); if isfield(FileInfo,'MaskFile') + if exist(FileInfo.MaskFile) set(handles.CheckMask,'Value',1) + else + set(handles.CheckMask,'Value',0) + end Mask.File=FileInfo.MaskFile; if isfield(FileInfo,'MaskNbSlice') @@ -2932,5 +2936,5 @@ Mask.File=filemask; [RootPath,SubDir,RootFile,i1_series,i2,j1,j2,NomType]=find_file_series(MaskPath,[FileName FileExt]); - Mask.NbSlice=1;%default + Mask.NbSlice=[];%default Mask.VolumeScan=0;% TO UPDATE *** if strcmp(NomType,'_1') @@ -2944,4 +2948,7 @@ set(handles.CheckMask,'UserData',Mask); errormsg=update_mask(handles); + if ~isempty(errormsg) + msgbox_uvmat(['ERROR','error in displaying mask, ' errormsg]); + end testmask=1; end @@ -2969,5 +2976,5 @@ MaskInfo=get(handles.CheckMask,'UserData'); if isfield(MaskInfo,'File') - if isfield(MaskInfo,'NbSlice') + if isfield(MaskInfo,'NbSlice')&& ~isempty(MaskInfo.NbSlice) if isfield(MaskInfo,'VolumeScan') && MaskInfo.VolumeScan MaskIndex_i=str2num(get(handles.j1,'String'));