Index: /trunk/src/get_field.m =================================================================== --- /trunk/src/get_field.m (revision 875) +++ /trunk/src/get_field.m (revision 876) @@ -268,27 +268,4 @@ drawnow uiwait(handles.get_field); - -%------------------------------------------------------------------------ -% --- Executes when user attempts to close get_field. -%------------------------------------------------------------------------ -%------------------------------------------------------------------------ - - -% Use UIRESUME instead of delete because the OutputFcn needs -% to get the updated handles structure. - -% function get_field_CloseRequestFcn(hObject, eventdata) -% handles.output=[]; -% guidata(hObject, handles);% Update handles structure -% %delete(handles.get_field); -% uiresume(handles.get_field); -%drawnow -% if isequal(get(handles.get_field, 'waitstatus'), 'waiting') -% % The GUI is still in UIWAIT, us UIRESUME -% uiresume(handles.get_field); -% else -% % The GUI is no longer waiting, just close it -% delete(handles.get_field); -% end % ----------------------------------------------------------------------- @@ -775,8 +752,8 @@ end if numel(var_component)<2 - if numel(test_coord)<2 + if numel(find(test_coord))<2 ListCoord={''}; else - if numel(test_coord)>=3 + if numel(find(test_coord))>=3 set(handles.Coord_x,'Value',3) set(handles.Coord_y,'Value',2) @@ -1068,4 +1045,5 @@ set(handles.vector_z,'Visible',status) set(handles.W_title,'Visible',status) +if strcmp(status,'on') Field=get(handles.get_field,'UserData'); if Field.MaxDim>=3% for 3D fields, propose to use the third variable as time @@ -1077,5 +1055,5 @@ end end - +end %------------------------------------------------------------------------ % --- Executes on button press in OK. Index: /trunk/src/series.m =================================================================== --- /trunk/src/series.m (revision 875) +++ /trunk/src/series.m (revision 876) @@ -2419,4 +2419,5 @@ GetFieldData=get_field(FirstFileName,ParamIn); FieldList={}; + if isfield(GetFieldData,'FieldOption')% if a field has been selected switch GetFieldData.FieldOption case 'vectors' @@ -2482,4 +2483,5 @@ set(handles.Coord_x,'Visible','on') set(handles.Coord_y,'Visible','on') + end else msgbox_uvmat('ERROR',[FirstFileName ' does not exist']) Index: /trunk/src/series/civ_input.m =================================================================== --- /trunk/src/series/civ_input.m (revision 875) +++ /trunk/src/series/civ_input.m (revision 876) @@ -1860,8 +1860,17 @@ end if isfield(Param,'OutputSubDir') - Param=rmfield(Param,'OutputSubDir'); %remove output file option from civ_series + Param=rmfield(Param,'OutputSubDir'); %remove output file option from civ_series end ParamPatch1=Param.ActionInput.Patch1; %store the patch1 parameters Param.ActionInput=rmfield(Param.ActionInput,'Patch1');% does not execute Patch + Param.IndexRange.first_i=str2num(get(handles.ref_i,'String')); + Param.IndexRange.last_i=Param.IndexRange.first_i; + if strcmp(get(handles.ref_j,'Visible'),'on') + Param.IndexRange.first_j=str2num(get(handles.ref_j,'String')); + Param.IndexRange.last_j=Param.IndexRange.first_j; + else + Param.IndexRange.first_j=1; + Param.IndexRange.last_j=1; + end [Data,errormsg]=civ_series(Param);% get the civ1+fix1 results bckcolor=get(handles.civ_input,'Color'); @@ -1938,4 +1947,13 @@ set(handles.Fix1,'BackgroundColor',[1 1 0]) set(handles.Patch1,'BackgroundColor',[1 1 0]) + Param.IndexRange.first_i=str2num(get(handles.ref_i,'String')); + Param.IndexRange.last_i=Param.IndexRange.first_i; + if strcmp(get(handles.ref_j,'Visible'),'on') + Param.IndexRange.first_j=str2num(get(handles.ref_j,'String')); + Param.IndexRange.last_j=Param.IndexRange.first_j; + else + Param.IndexRange.first_j=1; + Param.IndexRange.last_j=1; + end [Data,errormsg]=civ_series(Param);% get the civ1+fix1 results @@ -2031,4 +2049,13 @@ ParamPatch2=Param.ActionInput.Patch2; %store the patch1 parameters Param.ActionInput=rmfield(Param.ActionInput,'Patch2');% does not execute Patch + Param.IndexRange.first_i=str2num(get(handles.ref_i,'String')); + Param.IndexRange.last_i=Param.IndexRange.first_i; + if strcmp(get(handles.ref_j,'Visible'),'on') + Param.IndexRange.first_j=str2num(get(handles.ref_j,'String')); + Param.IndexRange.last_j=Param.IndexRange.first_j; + else + Param.IndexRange.first_j=1; + Param.IndexRange.last_j=1; + end [Data,errormsg]=civ_series(Param);% get the civ1+fix1 results bckcolor=get(handles.civ_input,'Color'); Index: /trunk/src/series/civ_series.m =================================================================== --- /trunk/src/series/civ_series.m (revision 875) +++ /trunk/src/series/civ_series.m (revision 876) @@ -691,8 +691,4 @@ end end -% ibx2=ceil(par_civ2.CorrBoxSize(1)/2); -% iby2=ceil(par_civ2.CorrBoxSize(2)/2); - % par_civ2.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess - %par_civ2.SearchBoxSize(2)=2*iby2+9; if CheckInputFile % else Dt given by par_civ2 if strcmp(Param.ActionInput.ListCompareMode,'displacement') @@ -707,8 +703,8 @@ par_civ2.Grid=[par_civ2.Grid(nbval>=1,1)-par_civ2.SearchBoxShift(:,1)/2 par_civ2.Grid(nbval>=1,2)-par_civ2.SearchBoxShift(:,2)/2]; if par_civ2.CheckDeformation - par_civ2.DUDX=DUDX./nbval; - par_civ2.DUDY=DUDY./nbval; - par_civ2.DVDX=DVDX./nbval; - par_civ2.DVDY=DVDY./nbval; + par_civ2.DUDX=DUDX(nbval>=1)./nbval(nbval>=1); + par_civ2.DUDY=DUDY(nbval>=1)./nbval(nbval>=1); + par_civ2.DVDX=DVDX(nbval>=1)./nbval(nbval>=1); + par_civ2.DVDY=DVDY(nbval>=1)./nbval(nbval>=1); end @@ -910,8 +906,8 @@ %% prepare correlation and search boxes -ibx2=ceil(par_civ.CorrBoxSize(1)/2); -iby2=ceil(par_civ.CorrBoxSize(2)/2); -isx2=ceil(par_civ.SearchBoxSize(1)/2); -isy2=ceil(par_civ.SearchBoxSize(2)/2); +ibx2=floor(par_civ.CorrBoxSize(1)/2); +iby2=floor(par_civ.CorrBoxSize(2)/2); +isx2=floor(par_civ.SearchBoxSize(1)/2); +isy2=floor(par_civ.SearchBoxSize(2)/2); shiftx=round(par_civ.SearchBoxShift(:,1));%use the input shift estimate, rounded to the next integer value shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases @@ -968,6 +964,6 @@ end check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 ); - par_civ.ImageA(check_undefined)=0;% put image A to zero (i.e. the min image value) in the undefined area - par_civ.ImageB(check_undefined)=0;% put image B to zero (i.e. the min image value) in the undefined area +% par_civ.ImageA(check_undefined)=0;% put image A to zero (i.e. the min image value) in the undefined area +% par_civ.ImageB(check_undefined)=0;% put image B to zero (i.e. the min image value) in the undefined area end @@ -993,5 +989,5 @@ image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask - mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=min for mask + mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1 check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight; @@ -1006,4 +1002,6 @@ F(ivec)=3; % else + image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1) + image2_crop=image2_crop.*~mask2_crop; image1_mean=mean(mean(image1_crop))/(1-sizemask); image2_mean=mean(mean(image2_crop))/(1-sizemask); @@ -1026,10 +1024,14 @@ XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2); + image1_crop_1=image1_crop; image1_crop=interp2(image1_crop,XIant,YIant); + image1_crop_2=image1_crop; image1_crop(isnan(image1_crop))=0; + image1_crop_3=image1_crop; xi=(1:mesh:size(image2_crop,2)); yi=(1:mesh:size(image2_crop,1))'; image2_crop=interp2(image2_crop,xi,yi); image2_crop(isnan(image2_crop))=0; + par_civ.CorrSmooth=3; %% end @@ -1041,4 +1043,7 @@ %Find the correlation max, at 255 [y,x] = find(result_conv==255,1); + subimage2_crop=image2_crop(y:y+2*iby2/mesh,x:x+2*ibx2/mesh);%subimage of image 2 corresponding to the optimum displacement of first image + sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2 + sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation if ~isempty(y) && ~isempty(x) try @@ -1047,4 +1052,6 @@ elseif par_civ.CorrSmooth==2 [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y); + else + [vector,F(ivec)] = quadr_fit(result_conv,x,y); end utable(ivec)=vector(1)*mesh+shiftx(ivec); @@ -1052,6 +1059,7 @@ xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2) ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge) - iref=round(xtable(ivec));% image index for the middle of the vector + iref=round(xtable(ivec));% nearest image index for the middle of the vector jref=round(ytable(ivec)); + % eliminate vectors located in the mask if checkmask && par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100 utable(ivec)=0; @@ -1150,4 +1158,42 @@ end vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1]; + +%------------------------------------------------------------------------ +% --- Find the maximum of the correlation function after quadratic interpolation +function [vector,F] = quadr_fit(result_conv,x,y) +[npy,npx]=size(result_conv); +if x<4 || y<4 || npx-x<4 ||npy-y <4 + F=-2; + vector=[x y]; +else + F=0; + x_ind=x-4:x+4; + y_ind=y-4:y+4; + x_vec=0.25*(x_ind-x); + y_vec=0.25*(y_ind-y); + [X,Y]=meshgrid(x_vec,y_vec); + coord=[reshape(X,[],1) reshape(Y,[],1)]; + result_conv=reshape(result_conv(y_ind,x_ind),[],1); + + + % n=numel(X); + % x=[X Y]; + % X=X-0.5; + % Y=Y+0.5; + % y = (X.*X+2*Y.*Y+X.*Y+6) + 0.1*rand(n,1); + p = polyfitn(coord,result_conv,2); + A(1,1)=2*p.Coefficients(1); + A(1,2)=p.Coefficients(2); + A(2,1)=p.Coefficients(2); + A(2,2)=2*p.Coefficients(4); + vector=[x y]'-A\[p.Coefficients(3) p.Coefficients(5)]'; + vector=vector'-[floor(npx/2) floor(npy/2)]-1 ; + % zg = polyvaln(p,coord); + % figure + % surf(x_vec,y_vec,reshape(zg,9,9)) + % hold on + % plot3(X,Y,reshape(result_conv,9,9),'o') + % hold off +end %'RUN_FIX': function for fixing velocity fields: Index: /trunk/src/series/stereo_civ.m =================================================================== --- /trunk/src/series/stereo_civ.m (revision 875) +++ /trunk/src/series/stereo_civ.m (revision 876) @@ -1,3 +1,3 @@ -%'civ_series': PIV function activated by the general GUI series +%'stereo_civ': determination of topography by image correlation of two stereo views % --- call the sub-functions: % civ: PIV function itself Index: /trunk/src/series/stereo_input.m =================================================================== --- /trunk/src/series/stereo_input.m (revision 875) +++ /trunk/src/series/stereo_input.m (revision 876) @@ -1,3 +1,3 @@ -%'stereo_input': function associated with the GUI 'stereo_input.fig' to set the input parameters for civ_series +%'stereo_input': function associated with the GUI 'stereo_input.fig' to set the input parameters for stereo_civ %------------------------------------------------------------------------ % function ParamOut = stereo_input(Param) Index: /trunk/src/sub_field.m =================================================================== --- /trunk/src/sub_field.m (revision 875) +++ /trunk/src/sub_field.m (revision 876) @@ -43,4 +43,7 @@ SubData=Field; return +end +if ~isfield(Field_1,'VarAttribute') + Field_1.VarAttribute={}; end @@ -107,4 +110,5 @@ end + %% rename the dimensions of the second field if identical to those of the first for ilist=1:numel(Field_1.VarDimName) @@ -127,5 +131,5 @@ ind_remove=false(size(Field_1.ListVarName)); % loop on the variables of the second field Field_1 -for ilist=1:numel(Field_1.ListVarName) +for ilist=1:numel(Field_1.VarAttribute) % case of variable with a single dimension if ~isempty(Field_1.VarAttribute{ilist}) && ~isempty(regexp(Field_1.VarAttribute{ilist}.Role,'^coord'))% if variable with Role coord... is found. @@ -173,4 +177,9 @@ ListVarNameNew=ListVarNameNew(~check_remove); % %list of renaimed varaibles corresponding to ListVarNameSub VarDimNameSub=Field_1.VarDimName(~check_remove); +if numel(Field_1.VarAttribute)