Changeset 1154


Ignore:
Timestamp:
Jul 7, 2024, 11:22:00 PM (3 months ago)
Author:
sommeria
Message:

various improvements

Location:
trunk/src
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/filter_tps.m

    r1152 r1154  
    66% SubRange(NbCoord,2,NbSubdomain): range (min, max) of the coordinates x and y respectively, for each subdomain
    77% NbCentre(NbSubdomain): number of source points for each subdomain
    8 % FF: false flags preserved from the input, or equal to 20 for vectors excluded by the difference with the smoothed field
     8% FF: false flags preserved from the input, or equal to true for vectors excluded by the difference with the smoothed field
    99% U_smooth, V_smooth: filtered velocity components at the positions of the initial data
    1010% Coord_tps(NbCentre,NbCoord,NbSubdomain): positions of the tps centres
     
    6161
    6262%% smoothing parameter: CHANGED 03 May 2024 TO GET RESULTS INDEPENDENT OF SUBDOMAINSIZE
    63 %smoothing=Siz(1)*Siz(2)*FieldSmooth/1000%optimum smoothing increase as the area of the subdomain (division by 1000 to reach good values with the default GUI input)
     63%smoothing=Siz(1)*Siz(2)*FieldSmooth/1000%old calculation before 03 May < r1129
    6464NbVecSub=NbVec/NbSubDomain;% refined estimation of the nbre of vectors per subdomain
    6565smoothing=sqrt(Siz(1)*Siz(2)/NbVecSub)*FieldSmooth;%optimum smoothing increase as the typical mesh size =sqrt(SizX*SizY/NbVecSub)^1/2
     66Threshold=Threshold*Threshold;% take the square of the threshold to work with the modulus squared (not done before r1154)
    6667
    6768%% default output
     
    7576V_smooth=zeros(NbVec,1);% smoothed velocity V at the initial positions
    7677W_smooth=[];%default (2 component case)
    77 FF=zeros(NbVec,1);
     78FF=false(NbVec,1);%false flag=0 (false) by default
    7879nb_select=zeros(NbVec,1);
    7980check_empty=zeros(1,NbSubDomain);
     
    8889    while numel(ind_sel)>numel(ind_sel_previous)
    8990        ind_sel_previous=ind_sel;% record the set of selected vector indices for next iteration
    90         ind_sel= find(FF==0 & Coord(:,1)>=SubRange(1,1,isub) & Coord(:,1)<=SubRange(1,2,isub) & Coord(:,2)>=SubRange(2,1,isub) & Coord(:,2)<=SubRange(2,2,isub));% indices of vectors in the subdomain #isub
     91        ind_sel= find(~FF & Coord(:,1)>=SubRange(1,1,isub) & Coord(:,1)<=SubRange(1,2,isub) & Coord(:,2)>=SubRange(2,1,isub) & Coord(:,2)<=SubRange(2,2,isub));% indices of vectors in the subdomain #isub
     92        %ind_sel= find( Coord(:,1)>=SubRange(1,1,isub) & Coord(:,1)<=SubRange(1,2,isub) & Coord(:,2)>=SubRange(2,1,isub) & Coord(:,2)<=SubRange(2,2,isub));% indices of vectors in the subdomain #isub
    9193        % if no vector in the subdomain  #isub, skip the subdomain
    9294        if isempty(ind_sel)
     
    9799            SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
    98100            SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
    99         % subdomain includes enough vectors, perform tps interpolation
     101        % if subdomain includes enough vectors, perform tps interpolation
    100102        else
    101103            [U_smooth_sub,U_tps_sub]=tps_coeff(Coord(ind_sel,:),U(ind_sel),smoothing);
     
    106108            ind_ind_sel=1:numel(ind_sel);%default
    107109            if exist('Threshold','var')&&~isempty(Threshold)
    108                 FF(ind_sel)=4*(NormDiff>Threshold);%put FF value to 4 to identify the criterium of elimmination
    109                 ind_ind_sel=find(FF(ind_sel)==0); % select the indices of remaining vectors in the subset of ind_sel vectors
     110                FF(ind_sel)=(NormDiff>Threshold);%put FF value to 1 to identify the criterium of elimmination
     111                ind_ind_sel=find(~FF(ind_sel)); % select the indices of remaining vectors in the subset of ind_sel vectors
    110112            end
    111113            % if no value exceeds threshold, the result is recorded
     
    116118                y_dist=(Coord(ind_sel,2)-CentreY(isub))/y_width;% relative ydistance to the retangle centre
    117119                weight=cos(x_dist).*cos(y_dist);%weighting fct =1 at the rectangle center and 0 at edge
     120                %weight=1;% case for r1129 and before
    118121                U_smooth(ind_sel)=U_smooth(ind_sel)+weight.*U_smooth_sub;
    119122                V_smooth(ind_sel)=V_smooth(ind_sel)+weight.*V_smooth_sub;
     
    138141                y_dist=(Coord(ind_sel(ind_ind_sel),2)-CentreY(isub))/y_width;% relative ydistance to the retangle centre
    139142                weight=cos(x_dist).*cos(y_dist);%weighting fct =1 at the rectangle center and 0 at edge
     143                %weight=1;
    140144                U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+weight.*U_smooth_sub;
    141145                V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+weight.*V_smooth_sub;
     
    166170U_smooth=U_smooth./nb_select;% take the average at the intersection of several subdomains
    167171V_smooth=V_smooth./nb_select;
    168 U_smooth(FF==4)=U(FF==4);% set to the initial values the eliminated vectors (flagged as false)
    169 V_smooth(FF==4)=V(FF==4);
     172
     173%eliminate the vectors with diff>threshold not yet eliminated
     174if exist('Threshold','var')&&~isempty(Threshold)
     175UDiff=U_smooth-U;% difference between interpolated U component and initial value
     176VDiff=V_smooth-V;% difference between interpolated V component and initial value
     177NormDiff=UDiff.*UDiff+VDiff.*VDiff;% Square of difference norm
     178FF(NormDiff>Threshold)=true;%put FF value to 1 to identify the criterium of elimmination
     179end
     180
     181U_smooth(FF)=U(FF);% set to the initial values the eliminated vectors (flagged as false)
     182V_smooth(FF)=V(FF);
    170183fill=zeros(NbCoord+1,NbCoord,size(SubRange,3)); %matrix of zeros to complement the matrix Data.Civ1_Coord_tps (conveninent for file storage)
    171184Coord_tps=cat(1,Coord_tps,fill);
  • trunk/src/keyboard_callback.m

    r1127 r1154  
    2626    switch xx
    2727        case {29,28,30,31}    %arrows for displacement
     28            hhh=get(hObject,'CurrentObject');
    2829            AxeData=get(cur_axes,'UserData');
    2930            if isfield(AxeData,'ZoomAxes')&&ishandle(AxeData.ZoomAxes)
     
    3132                axes(cur_axes)
    3233            end
    33             if ~isempty(cur_axes)
     34            if ~isempty(cur_axes) && ~strcmp(get(hhh,'Type'),'uicontrol')
    3435                xlimit=get(cur_axes,'XLim');
    3536                ylimit=get(cur_axes,'Ylim');
  • trunk/src/series.m

    r1152 r1154  
    572572        PairString=get(handles.PairString,'Data');
    573573        if numel(PairString)>=iview
    574             checkpair=strfind(PairString{iview},'j=');
    575             if checkpair
    576                 j1=str2double(PairString{iview}(4));
    577                 j2=str2double(PairString{iview}(6));
    578             end
     574                r=regexp(PairString{iview},'(?<num1>\d+)-(?<num2>\d+)' ,'names');
     575                if ~isempty(r)
     576                j1=str2double(r.num1);
     577                j2=str2double(r.num2);
     578                end
    579579        end
    580580        InputFile=fullfile_uvmat('','',InputTable{iview,3},InputTable{iview,5},InputTable{iview,4},i1,[],j1,j2);
     
    14001400    errormsg='input field name(s) not defined, select add_field...';
    14011401    return
     1402end
     1403[status,result]=system(['svn info ' Param.Action.ActionPath]);
     1404if status==0
     1405    t=regexp(result,'R.vision\s*:\s*(?<rev>\d+)','names');%detect 'revision' or 'Revision' in the text
     1406    if ~isempty(t)
     1407        Param.UvmatRevision=t.rev; %version nbre of the current package
     1408    end
    14021409end
    14031410
     
    23052312        return
    23062313    end
    2307     [tild,ActionName,ActionExt]=fileparts(FileName);
     2314    [~,ActionName,ActionExt]=fileparts(FileName);
    23082315
    23092316    % insert the choice in the menu ActionName
     
    23972404%% Activate the Action fct to adapt the configuration of the GUI series and bring specific parameters in SeriesData
    23982405Param=read_GUI_series(handles); % read the parameters from the GUI series
    2399 Param.Action.RUN=0;
    2400 Param.SeriesData=SeriesData;
     2406Param.Action.RUN=0;% indicate that we are in the mode of parameter input, not program run
     2407Param.SeriesData=SeriesData;% info stored in 'UserData' of the fig 'series'
    24012408ParamOut=h_fun(Param); % run the selected Action function to get the relevant input
    24022409
     
    24652472        ind_var=get(handles.FieldName,'Value'); % indices of previously selected variables
    24662473        for ilist=1:numel(ind_var)
    2467             if isempty(find(strcmp(FieldList{ind_var(ilist)},ListVarName)))
     2474            if isempty(find(strcmp(FieldList{ind_var(ilist)},ListVarName), 1))
    24682475                FieldList={}; % previous choice not consistent with new input field
    24692476                set(handles.FieldName,'Value',1)
     
    24712478            end
    24722479        end
    2473         if ~isempty(FieldList)iview_netcdf
    2474             if isempty(find(strcmp(get(handles.Coord_x,'String'),ListVarName)))||...
    2475                     isempty(find(strcmp(get(handles.Coord_y,'String'),ListVarName)))
     2480        if ~isempty(FieldList)
     2481            if isempty(find(strcmp(get(handles.Coord_x,'String'),ListVarName), 1))||...
     2482                    isempty(find(strcmp(get(handles.Coord_y,'String'),ListVarName), 1))
    24762483                FieldList={};
    24772484                set(handles.Coord_x,'String','')
     
    24792486            end
    24802487            Coord_z=get(handles.Coord_z,'String');
    2481             if ~isempty(Coord_z) && isempty(find(strcmp(Coord_z,ListVarName)))REFRESH
     2488            if ~isempty(Coord_z) && isempty(find(strcmp(Coord_z,ListVarName), 1))
    24822489                FieldList={};
    24832490                set(handles.Coord_z,'String','')
     
    24892496    end
    24902497    set(handles_coord,'Visible','on')
    2491     if isempty(find(strcmp('add_field...',FieldList)))
     2498    if isempty(find(strcmp('add_field...',FieldList), 1))
    24922499        FieldList=[FieldList;{'add_field...'}];%add 'add_field...' to the menu FieldName if it is not already
    24932500    end
     
    24962503        set(handles.Field_text_1,'Visible','on')
    24972504        if CheckPivData_1==0        % not civ input made
    2498             FieldList_1={'add_field...'}
     2505            FieldList_1={'add_field...'};
    24992506            ListVarName=SeriesData.FileInfo{iview_netcdf(2)}.ListVarName;
    25002507            ind_var=get(handles.FieldName,'Value'); % indices of previously selected variables
    25012508            for ilist=1:numel(ind_var)
    2502                 if isempty(find(strcmp(FieldList{ind_var(ilist)},ListVarName)))
     2509                if isempty(find(strcmp(FieldList{ind_var(ilist)},ListVarName), 1))
    25032510                    %FieldList_1={}; % previous choice not consistent with new input field
    25042511                    set(handles.FieldName_1,'Value',1)
     
    25072514            end
    25082515            warn_coord=0;
    2509             if isempty(find(strcmp(get(handles.Coord_x,'String'),ListVarName)))||...
    2510                     isempty(find(strcmp(get(handles.Coord_y,'String'),ListVarName)))
     2516            if isempty(find(strcmp(get(handles.Coord_x,'String'),ListVarName), 1))||...
     2517                    isempty(find(strcmp(get(handles.Coord_y,'String'),ListVarName), 1))
    25112518                warn_coord=1;
    25122519            end
    2513             if ~isempty(Coord_z) && isempty(find(strcmp(Coord_z,ListVarName)))
     2520            if ~isempty(Coord_z) && isempty(find(strcmp(Coord_z,ListVarName), 1))
    25142521                FieldList_1={'add_field...'};
    25152522                warn_coord=1;
     
    25472554%% Check whether alphabetical sorting of input Subdir is allowed by the Action fct  (for multiples series entries)
    25482555if isfield(ParamOut,'AllowInputSort')&&isequal(ParamOut.AllowInputSort,'on')&& size(Param.InputTable,1)>1
    2549     [tild,iview]=sort(Param.InputTable(:,2)); % subdirectories sorted in alphabetical order
     2556    [~,iview]=sort(Param.InputTable(:,2)); % subdirectories sorted in alphabetical order
    25502557    set(handles.InputTable,'Data',Param.InputTable(iview,:));
    25512558    MinIndex_i=get(handles.MinIndex_i,'Data');
     
    26322639
    26332640%% NbSlice visibility
    2634 % if isfield(ParamOut,'OutputFileMode')&& strcmp(ParamOut.OutputFileMode,'NbSlice')
    2635 %     ParamOut.NbSlice='on';
    2636 % end
    26372641if isfield(ParamOut,'NbSlice') && (strcmp(ParamOut.NbSlice,'on')||isnumeric(ParamOut.NbSlice))
    26382642    set(handles.num_NbSlice,'Visible','on')
     
    26902694end
    26912695set(handles.CheckMask,'Visible',MaskVisible);
     2696set(handles.MaskTable,'Visible',MaskVisible);
     2697
    26922698%% Setting of expected iteration time
    26932699if isfield(ParamOut,'CPUTime')
     
    26982704InputTable=get(handles.InputTable,'Data');
    26992705[OutputPath,Device,DeviceExt]=fileparts(InputTable{1,1});
    2700 [OutputPath,Experiment,ExperimentExt]=fileparts(OutputPath);
     2706[~,Experiment,ExperimentExt]=fileparts(OutputPath);
    27012707set(handles.Device,'String',[Device DeviceExt])
    27022708set(handles.Device,'Visible','on')
     
    27102716
    27112717%% definition of the subdirectory containing the output files
    2712 
    27132718if  ~(isfield(SeriesData,'ActionName') && strcmp(ActionName,SeriesData.ActionName))
    27142719    OutputDirExt='.series'; % default
     
    27522757set(handles.OutputDirExt,'Visible',OutputDirVisible)
    27532758set(handles.OutputSubDir,'Visible',OutputDirVisible)
    2754 % set(handles.OutputDir_title,'Visible',OutputDirVisible)
    27552759SeriesData.ActionName=ActionName; % record ActionName for next use
    27562760
     
    27762780%% definition of an additional parameter set, determined by an ancillary GUI
    27772781if isfield(ParamOut,'ActionInput')
    2778 %     set(handles.ActionInput,'Visible','on')
    27792782    ParamOut.ActionInput.Program=ActionName; % record the program in ActionInput
    27802783    SeriesData.ActionInput=ParamOut.ActionInput;
    27812784else
    2782 %     set(handles.ActionInput,'Visible','off')
    27832785    if isfield(SeriesData,'ActionInput')
    27842786        SeriesData=rmfield(SeriesData,'ActionInput');
  • trunk/src/series/civ_input.m

    r1153 r1154  
    141141        iview_image=2;%line # for the input images
    142142    otherwise
    143         if ~strcmp(FieldType,'image')
    144         msgbox_uvmat('ERROR','civ_series needs images, scalar fields in netcdf format, or civ data as input')
    145         return
    146         end
     143        % if ~strcmp(FileType,'image')
     144        % msgbox_uvmat('ERROR','civ_series needs images, scalar fields in netcdf format, or civ data as input')
     145        % return
     146        % end
    147147end
    148148       
  • trunk/src/series/civ_series.m

    r1153 r1154  
    127127    end
    128128    if isfield(Param,'InputTable')
    129         [tild,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
     129        [~,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
    130130        iview_A=0;% series index (iview) for the first image series
    131131        iview_B=0;% series index (iview) for the second image series (only non zero for option 'shift' comparing two image series )
     
    195195                j2_series_Civ2=j2_series_Civ1;
    196196                NomTypeNc='_1';
    197             case 'PIV volume'
    198                 % TODO, TODO
    199197        end
    200198        %determine frame indices for input with movie or other multiframe input file
     
    252250Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
    253251Data.Program='civ_series';
     252if isfield(Param,'UvmatRevision')
     253    Data.Program=[Data.Program ', uvmat r' Param.UvmatRevision];
     254end
    254255Data.CivStage=0;%default
    255256
     
    509510            % caluclate velocity data (y and v in indices, reverse to y component)
    510511            tstart_civ1=tic;
    511             [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1);
     512            [xtable, ytable, utable, vtable, ctable, FF, result_conv, errormsg] = civ (par_civ1);
    512513            if ~isempty(errormsg)
    513514                disp_uvmat('ERROR',errormsg,checkrun)
     
    519520            Data.Civ1_V=reshape(-vtable,[],1);
    520521            Data.Civ1_C=reshape(ctable,[],1);
    521             Data.Civ1_FF=reshape(F,[],1);
     522            Data.Civ1_FF=reshape(FF,[],1);
    522523            time_civ1=toc(tstart_civ1);
    523524        end
     
    601602
    602603        % perform Patch calculation using the UVMAT fct 'filter_tps'
    603        
    604604        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=...
    605605            filter_tps([Data.Civ1_X(ind_good) Data.Civ1_Y(ind_good)],Data.Civ1_U(ind_good),Data.Civ1_V(ind_good),[],Data.Patch1_SubDomainSize,Data.Patch1_FieldSmooth,Data.Patch1_MaxDiff);
    606606        Data.Civ1_U_smooth(ind_good)=Ures;% take the interpolated (smoothed) velocity values for good vectors, keep civ1 data for the other
    607607        Data.Civ1_V_smooth(ind_good)=Vres;
    608         Data.Civ1_FF(ind_good)=uint8(FFres);
     608        Data.Civ1_FF(ind_good)=uint8(4*FFres);
    609609        time_patch1=toc(tstart_patch1);
    610610        disp('patch1 performed')
     
    770770        % calculate velocity data (y and v in image indices, reverse to y component)
    771771       
    772         [xtable, ytable, utable, vtable, ctable, F,result_conv,errormsg] = civ (par_civ2);
     772        [xtable, ytable, utable, vtable, ctable, FF,result_conv,errormsg] = civ (par_civ2);
    773773       
    774774        list_param=(fieldnames(Param.ActionInput.Civ2))';
     
    810810        Data.Civ2_V=reshape(-vtable,[],1);
    811811        Data.Civ2_C=reshape(ctable,[],1);
    812         Data.Civ2_FF=reshape(F,[],1);
     812        Data.Civ2_FF=reshape(FF,[],1);
    813813        disp('civ2 performed')
    814814        time_civ2=toc(tstart_civ2);
     
    821821                return
    822822            end
    823 %         elseif isfield(Param,'Civ2_X')% use Civ2 data as input in Param (test mode)
    824 %             Data.ListGlobalAttribute={};
    825 %             Data.ListVarName={};
    826 %             Data.VarDimName={};
    827 %             Data.Civ2_X=Param.Civ2_X;
    828 %             Data.Civ2_Y=Param.Civ2_Y;
    829 %             Data.Civ2_U=Param.Civ2_U;
    830 %             Data.Civ2_V=Param.Civ2_V;
    831 %             Data.Civ2_FF=Param.Civ2_FF;
    832823        end
    833824    end
     
    886877        Data.Civ2_U_smooth(ind_good)=Ures;
    887878        Data.Civ2_V_smooth(ind_good)=Vres;
    888         Data.Civ2_FF(ind_good)=FFres;
     879        Data.Civ2_FF(ind_good)=uint8(4*FFres);
    889880        Data.CivStage=Data.CivStage+1;
    890881        time_patch2=toc(tstart_patch2);
     
    921912% OUTPUT:
    922913% xtable: set of x coordinates
    923 % ytable: set of y coordiantes
     914% ytable: set of y coordinates
    924915% utable: set of u displacements (along x)
    925916% vtable: set of v displacements (along y)
     
    12781269% FF=2, for too small correlation
    12791270% FF=3, for velocity outside bounds
    1280 % FF=4 for exclusion by difference with the smoothed field, not set in this function
     1271% FF=4 for exclusion by difference with the smoothed field, set by call to function filter_tps
    12811272
    12821273if isfield (Param,'MinCorr')
  • trunk/src/uvmat.m

    r1153 r1154  
    20162016%------------------------------------------------------------------------
    20172017% open the GUI 'series'
    2018 function MenuSeries_Callback(hObject, eventdata, handles)
     2018function MenuRun1_Callback(hObject, eventdata, handles)
    20192019%------------------------------------------------------------------------
    20202020Param=read_GUI(handles.uvmat);
     
    20232023
    20242024% --------------------------------------------------------------------
    2025 function MenuPIV_Callback(hObject, eventdata, handles)
     2025function MenuRun2_Callback(hObject, eventdata, handles)
    20262026Param=read_GUI(handles.uvmat);
    20272027Param.HiddenData=get(handles.uvmat,'UserData');
     
    20332033series('ActionName_Callback',hObject,eventdata,hhseries); %file input with xml reading  in uvmat, show the image in phys coordinates
    20342034
    2035 %------------------------------------------------------------------------
    2036 % -- open the GUI civ.fig for PIV
    2037 function MenuCIVx_Callback(hObject, eventdata, handles)
    2038 %------------------------------------------------------------------------
    2039  [RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
    2040  FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];
    2041 civ(FileName);% interface de civ(not in the uvmat file)
     2035% %------------------------------------------------------------------------
     2036% % -- open the GUI civ.fig for PIV
     2037% function MenuCIVx_Callback(hObject, eventdata, handles)
     2038% %------------------------------------------------------------------------
     2039% [RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
     2040% FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];
     2041% civ(FileName);% interface de civ(not in the uvmat file)
    20422042
    20432043%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
    26032603        set(handles.ColorScalar,'String',ColorList)
    26042604        set(handles.Vectors,'Visible','on')
    2605         %set(handles.Coord_x,'Value',1);
    26062605        set(handles.Coord_x,'String','X');
    26072606        set(handles.Coord_y,'String','Y');
     2607        set(handles.MenuRun3,'Label','test_filter_tps')
    26082608    case {'netcdf','mat'}
    26092609        set(handles_Fields,'Value',1)
     
    61086108
    61096109
    6110 
    61116110% --------------------------------------------------------------------
    61126111% --- Executes on button press in CheckTable.
     
    61216120
    61226121
    6123 
    6124 
    6125 
    6126 
    6127 
     6122% --------------------------------------------------------------------
     6123function MenuRun3_Callback(hObject, eventdata, handles)
     6124Param=read_GUI(handles.uvmat);
     6125Param.HiddenData=get(handles.uvmat,'UserData');
     6126hseries=series(Param);
     6127hhseries=guidata(hseries);
     6128ActionMenu=get(hhseries.ActionName,'String');
     6129index_action=find(strcmp('test_filter_tps',ActionMenu));
     6130set(hhseries.ActionName,'Value',index_action);
     6131series('ActionName_Callback',hObject,eventdata,hhseries); %file input with xml reading  in uvmat, show the image in phys coordinates
Note: See TracChangeset for help on using the changeset viewer.