Changeset 1158


Ignore:
Timestamp:
Jul 12, 2024, 8:46:38 AM (3 months ago)
Author:
sommeria
Message:

civ_3D corrected, automatic vector position modified in civ_series for better symmetry

Location:
trunk/src
Files:
1 added
1 deleted
6 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/mouse_motion.m

    r1143 r1158  
    165165                            end
    166166                            %display the field values
     167                                if isfield (CellInfo{icell},'XName')
     168                                    XName=CellInfo{icell}.XName;
     169                                    text_displ_2=[XName '=' num2str(Field.(XName)(ivec),4)];
     170                                end
     171                                if isfield (CellInfo{icell},'YName')
     172                                    YName=CellInfo{icell}.YName;
     173                                    text_displ_2=[text_displ_2 ',' YName '=' num2str(Field.(YName)(ivec),4)];
     174                                end
    167175                            for ivar=1:numel(CellInfo{icell}.VarIndex)
    168176                                VarName=Field.ListVarName{CellInfo{icell}.VarIndex(ivar)};
    169177                                VarVal=Field.(VarName)(ivec);
    170178                                var_text=[VarName '=' num2str(VarVal,4) ','];
     179 
    171180                                if isequal(ivar,CellInfo{icell}.CoordIndex(end))||isequal(ivar,CellInfo{icell}.CoordIndex(end-1))||isequal(ivar,CellInfo{icell}.CoordIndex(1))
    172181                                    text_displ_1=[text_displ_1 var_text];
  • trunk/src/plot_field.m

    r1152 r1158  
    13311331    n=size(xc);
    13321332    xN=NaN*ones(size(xc));
    1333     matx=[xc(:)-uc(:)/2 xc(:)+uc(:)/2 xN(:)]';
     1333    matx=[xc(:) xc(:)+uc(:) xN(:)]';
    13341334    matx=reshape(matx,1,3*n(2));
    1335     maty=[yc(:)-vc(:)/2 yc(:)+vc(:)/2 xN(:)]';
     1335    maty=[yc(:) yc(:)+vc(:) xN(:)]';
    13361336    maty=reshape(maty,1,3*n(2));
    13371337   
     
    13391339    arrowplus=rot*[uc;vc];
    13401340    arrowmoins=rot'*[uc;vc];
    1341     x1=xc+uc/2-arrowplus(1,:);
    1342     x2=xc+uc/2;
    1343     x3=xc+uc/2-arrowmoins(1,:);
    1344     y1=yc+vc/2-arrowplus(2,:);
    1345     y2=yc+vc/2;
    1346     y3=yc+vc/2-arrowmoins(2,:);
     1341    % x1=xc+uc/2-arrowplus(1,:);
     1342    % x2=xc+uc/2;
     1343    % x3=xc+uc/2-arrowmoins(1,:);
     1344    % y1=yc+vc/2-arrowplus(2,:);
     1345    % y2=yc+vc/2;
     1346    % y3=yc+vc/2-arrowmoins(2,:);
     1347    x1=xc+uc-arrowplus(1,:);
     1348    x2=xc+uc;
     1349    x3=xc+uc-arrowmoins(1,:);
     1350    y1=yc+vc-arrowplus(2,:);
     1351    y2=yc+vc;
     1352    y3=yc+vc-arrowmoins(2,:);
    13471353    matxar=[x1(:) x2(:) x3(:) xN(:)]';
    13481354    matxar=reshape(matxar,1,4*n(2));
  • trunk/src/proj_field.m

    r1135 r1158  
    13601360                        if ~(check3D && ivar==CellInfo{icell}.CoordIndex(1))
    13611361                            ListVarName=[ListVarName VarName];
    1362                             VarDimName=[VarDimName DimCell];
     1362                            VarDimName=[VarDimName {DimCell}];
    13631363                            nbvar=nbvar+1;
    13641364                            if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
  • trunk/src/series/civ_3D.m

    r1157 r1158  
    180180Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param];
    181181% set the list of variables
    182 Data.ListVarName={'Coord_z','Coord_y','Coord_x','Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};%  cell array containing the names of the fields to record
    183 Data.VarDimName={'npz','npy','npx',{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},...
     182Data.ListVarName={'Coord_z','Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};%  cell array containing the names of the fields to record
     183Data.VarDimName={'npz',{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},...
    184184    {'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'}};
    185 Data.VarAttribute{1}.Role='coord_x';
    186 Data.VarAttribute{2}.Role='coord_y';
    187 Data.VarAttribute{3}.Role='vector_x';
    188 Data.VarAttribute{4}.Role='vector_y';
    189 Data.VarAttribute{5}.Role='vector_z';
    190 Data.VarAttribute{6}.Role='ancillary';
    191 Data.VarAttribute{7}.Role='errorflag';
     185Data.VarAttribute{1}.Role='coord_z';
     186Data.VarAttribute{2}.Role='coord_x';
     187ata.VarAttribute{3}.Role='coord_y';
     188Data.VarAttribute{4}.Role='vector_x';
     189Data.VarAttribute{5}.Role='vector_y';
     190Data.VarAttribute{6}.Role='vector_z';
     191Data.VarAttribute{7}.Role='ancillary';
     192Data.VarAttribute{8}.Role='errorflag';
    192193Data.Coord_z=j1_series_Civ1;
    193194% path for output nc file
     
    297298        par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx);
    298299       
    299         z_index=1;%first vertical block centered at image index z_index=par_civ1.Dz
    300         for iz=1:par_civ1.Dz+SearchRange_z
     300        z_index=1;%first vertical block centered at image index z_index=SearchRange_z+1
     301        for iz=1:2*SearchRange_z+1
    301302            j_image_index=j1_series_Civ1(iz,1)% j index of the current image
    302303            ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);%
     
    304305            ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j_image_index);
    305306            B= read_image(ImageName_B,FileType_B);
    306             par_civ1.ImageA(iz+par_civ1.Dz-1,:,:) = A;
    307             par_civ1.ImageB(iz+par_civ1.Dz-1,:,:) = B;
     307            par_civ1.ImageA(iz,:,:) = A;
     308            par_civ1.ImageB(iz,:,:) = B;
    308309        end
    309310        % calculate velocity data at the first position z, corresponding to j_index=par_civ1.Dz
     
    338339        Data.Civ1_C(Data.Civ1_C==Inf)=0;
    339340        [npz,npy,npx]=size(Data.Civ1_X);
    340         Data.Coord_z=1:npz;
     341        Data.Coord_z=SearchRange_z+1:par_civ1.Dz:NbSlice-1;
    341342        Data.Coord_y=1:npy;
    342343        Data.Coord_x=1:npx;
     
    788789
    789790%% prepare measurement grid
    790 
    791 minix=floor(par_civ.Dx/2)-0.5;
    792 maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
    793 miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
    794 maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
    795 [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
     791nbinterv_x=floor((par_civ.ImageWidth-1)/par_civ.Dx);
     792gridlength_x=nbinterv_x*par_civ.Dx;
     793minix=ceil((par_civ.ImageWidth-gridlength_x)/2);
     794nbinterv_y=floor((par_civ.ImageHeight-1)/par_civ.Dy);
     795gridlength_y=nbinterv_y*par_civ.Dy;
     796miniy=ceil((par_civ.ImageHeight-gridlength_y)/2);
     797[GridX,GridY]=meshgrid(minix:par_civ.Dx:par_civ.ImageWidth-1,miniy:par_civ.Dy:par_civ.ImageHeight-1);
     798% minix:par_civ.Dx:par_civ.ImageWidth-1
     799% miniy:par_civ.Dy:par_civ.ImageHeight-1
     800% minix=floor(par_civ.Dx/2)-0.5;
     801% maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
     802% miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
     803% maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
     804% [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
    796805par_civ.Grid(:,:,1)=GridX;
    797806par_civ.Grid(:,:,2)=GridY;% increases with array index,
     
    931940                    end
    932941
    933                
    934                    
    935                     for kz=1:par_civ.SearchBoxSize(3)
     942                    npxycorr=par_civ.SearchBoxSize(1:2)-par_civ.CorrBoxSize(1:2)+1;
     943                    result_conv=zeros([par_civ.SearchBoxSize(3) npxycorr]);%initialise the conv product
     944                    max_xy=zeros(par_civ.SearchBoxSize(3),1);%initialise the max correlation vs z
     945                    xk=ones(par_civ.SearchBoxSize(3),1);%initialise the x displacement corresponding to the max corr
     946                    yk=ones(par_civ.SearchBoxSize(3),1);%initialise the y displacement corresponding to the max corr
     947                 
     948                    for kz=kref:kref%par_civ.SearchBoxSize(3)
    936949                        subima2=squeeze(image2_crop(kz,:,:));
    937950                        subima1=squeeze(image1_crop(kref,:,:));
    938951                        correl_xy=conv2(subima2,flip(flip(subima1,2),1),'valid');
    939                           result_conv(kz,:,:)= correl_xy;
     952                        result_conv(kz,:,:)= 255*correl_xy;
    940953                        max_xy(kz)=max(max(correl_xy));
    941                     [xk(kz),yk(kz)]=find(correl_xy==max_xy(kz),1);
    942                
     954                        % max_xy(kz)=max(max(correl_xy));
     955                         [yk(kz),xk(kz)]=find(correl_xy==max_xy(kz),1);
    943956                    end
    944                     [corrmax,zmax]=max(max_xy);
    945                
    946                     x=xk(zmax);
    947                     y=yk(zmax);
    948                     result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
    949                     subimage2_crop=squeeze(image2_crop(zmax,y:y+2*iby2/mesh,x:x+2*ibx2/mesh));%subimage of image 2 corresponding to the optimum displacement of first image
    950                     sum_square=sum(sum(squeeze(image1_crop(zmax,:,:).*image1_crop(zmax,:,:))));
     957                    [corrmax,dz]=max(max_xy);
     958                    dx=xk(dz);
     959                    dy=yk(dz);
     960                 
     961                    %
     962                    % result_conv=255*result_conv/corrmax; %normalize, peak=always 255
     963                    % for kz=1:par_civ.SearchBoxSize(3)
     964                    %     [dy,dx] = find(result_conv(kz,:,:)==255,1)
     965                    %     if ~isempty(dy)
     966                    %         dz=kz;
     967                    %         break
     968                    %     end
     969                    % end
     970                    subimage2_crop=squeeze(image2_crop(dz,dy:dy+2*iby2/mesh,dx:dx+2*ibx2/mesh));%subimage of image 2 corresponding to the optimum displacement of first image
     971                    sum_square=sum(sum(squeeze(image1_crop(dz,:,:).*image1_crop(dz,:,:))));
    951972                    sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2
    952973                    sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation
    953                     if ~isempty(y) && ~isempty(x)
    954        
    955                             if par_civ.CorrSmooth==1
    956                                 [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (squeeze(result_conv(zmax,:,:)),x,y);%TODO: improve by max optimisation along z
    957                             elseif par_civ.CorrSmooth==2
    958                                 [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (squeeze(result_conv(zmax,:,:)),x,y);
    959                             else
    960                                 [vector,FF(ivec_y,ivec_x)] = quadr_fit(squeeze(result_conv(zmax,:,:)),x,y);
    961                             end
    962                             utable(ivec_y,ivec_x)=vector(1)*mesh+shiftx(ivec_y,ivec_x);
    963                             vtable(ivec_y,ivec_x)=vector(2)*mesh+shifty(ivec_y,ivec_x);
    964                             xtable(ivec_y,ivec_x)=iref+utable(ivec_y,ivec_x)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
    965                             ytable(ivec_y,ivec_x)=jref+vtable(ivec_y,ivec_x)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
    966                             iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector
    967                             jref=round(ytable(ivec_y,ivec_x)+0.5);
    968                             wtable(ivec_y,ivec_x)=zmax-kref;
    969                             % eliminate vectors located in the mask
    970                             if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
    971                                 utable(ivec_y,ivec_x)=0;
    972                                 vtable(ivec_y,ivec_x)=0;
    973                                 FF(ivec_y,ivec_x)=1;
    974                             end
    975                             ctable(ivec_y,ivec_x)=corrmax/sum_square;% correlation value
    976 
     974                    if ~isempty(dz)&& ~isempty(dy) && ~isempty(dx)
     975                        if par_civ.CorrSmooth==1
     976                            [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (result_conv,dx,dy,dz);
     977                        elseif par_civ.CorrSmooth==2
     978                            [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (result_conv,dx,dy,dz);
     979                        else
     980                            [vector,FF(ivec_y,ivec_x)] = quadr_fit(result_conv,dx,dy,dz);
     981                        end
     982                        utable(ivec_y,ivec_x)=vector(1)+shiftx(ivec_y,ivec_x);
     983                        vtable(ivec_y,ivec_x)=vector(2)+shifty(ivec_y,ivec_x);
     984                        wtable(ivec_y,ivec_x)=vector(3);
     985                        xtable(ivec_y,ivec_x)=iref+utable(ivec_y,ivec_x)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
     986                        ytable(ivec_y,ivec_x)=jref+vtable(ivec_y,ivec_x)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
     987                        iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector
     988                        jref=round(ytable(ivec_y,ivec_x)+0.5);
     989                 
     990                        % eliminate vectors located in the mask
     991                        if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
     992                            utable(ivec_y,ivec_x)=0;
     993                            vtable(ivec_y,ivec_x)=0;
     994                            FF(ivec_y,ivec_x)=1;
     995                        end
     996                        ctable(ivec_y,ivec_x)=corrmax/sum_square;% correlation value
    977997                    else
    978998                        FF(ivec_y,ivec_x)=true;
     
    9931013% x,y: position of the maximum correlation at integer values
    9941014
    995 function [vector,FF] = SUBPIXGAUSS (result_conv,x,y)
     1015function [vector,FF] = SUBPIXGAUSS (result_conv,x,y,z)
    9961016%------------------------------------------------------------------------
    997 vector=[0 0]; %default
    998 [npy,npx]=size(result_conv);
     1017
     1018[npz,npy,npx]=size(result_conv);
    9991019result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1  (=0 by discretisation, to avoid divergence in the log)
    10001020%the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology
    10011021%http://urapiv.wordpress.com
    10021022FF=false;
    1003     peaky = y;
    1004     if y < npy && y > 1
    1005         f0 = log(result_conv(y,x));
    1006         f1 = log(result_conv(y-1,x));
    1007         f2 = log(result_conv(y+1,x));
    1008         peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
    1009     else
    1010         FF=true; % error flag for vector truncated by the limited search box in y
    1011     end
    1012     peakx=x;
    1013     if x < npx-1 && x > 1
    1014         f0 = log(result_conv(y,x));
    1015         f1 = log(result_conv(y,x-1));
    1016         f2 = log(result_conv(y,x+1));
    1017         peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
    1018     else
    1019         FF=true; % warning flag for vector truncated by the limited search box in x
    1020     end
    1021     vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
     1023peakz=z;peaky=y;peakx=x;
     1024if z < npz && z > 1 && y < npy && y > 1 && x < npx && x > 1
     1025    f0 = log(result_conv(z,y,x));
     1026    f1 = log(result_conv(z-1,y,x));
     1027    f2 = log(result_conv(z+1,y,x));
     1028    peakz = peakz+ (f1-f2)/(2*f1-4*f0+2*f2);
     1029
     1030    f1 = log(result_conv(z,y-1,x));
     1031    f2 = log(result_conv(z,y+1,x));
     1032    peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
     1033
     1034    f1 = log(result_conv(z,y,x-1));
     1035    f2 = log(result_conv(z,y,x+1));
     1036    peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
     1037else
     1038    FF=true;
     1039end
     1040
     1041vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1 peakz-floor(npz/2)-1];
    10221042
    10231043%------------------------------------------------------------------------
  • trunk/src/series/civ_series.m

    r1157 r1158  
    962962    % else par_civ.Grid is already an array, no action here
    963963else% automatic grid in x, y image coordinates
    964     minix=floor(par_civ.Dx/2)-0.5;
    965     maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
    966     miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
    967     maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
    968     [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
     964nbinterv_x=floor((par_civ.ImageWidth-1)/par_civ.Dx);
     965gridlength_x=nbinterv_x*par_civ.Dx;
     966minix=ceil((par_civ.ImageWidth-gridlength_x)/2);
     967nbinterv_y=floor((par_civ.ImageHeight-1)/par_civ.Dy);
     968gridlength_y=nbinterv_y*par_civ.Dy;
     969miniy=ceil((par_civ.ImageHeight-gridlength_y)/2);
     970[GridX,GridY]=meshgrid(minix:par_civ.Dx:par_civ.ImageWidth-1,miniy:par_civ.Dy:par_civ.ImageHeight-1);
    969971    par_civ.Grid(:,1)=reshape(GridX,[],1);
    970972    par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index
     973    %
     974    %
     975    % minix=floor(par_civ.Dx/2)-0.5;
     976    % maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
     977    % miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
     978    % maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
     979    % [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
     980    % par_civ.Grid(:,1)=reshape(GridX,[],1);
     981    % par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index
    971982end
    972983nbvec=size(par_civ.Grid,1);
     
    10071018par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images
    10081019par_civ.ImageB=sum(double(par_civ.ImageB),3);
    1009 [npy_ima npx_ima]=size(par_civ.ImageA);
     1020[npy_ima,npx_ima]=size(par_civ.ImageA);
    10101021if ~isequal(size(par_civ.ImageB),[npy_ima npx_ima])
    10111022    errormsg='image pair with unequal size';
     
    11161127                sum_square=sum(sum(image1_crop.*image1_crop));
    11171128                %reference: Oliver Pust, PIV: Direct Cross-Correlation
    1118                 result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid');
     1129                result_conv= conv2(image2_crop,flip(flip(image1_crop,2),1),'valid');
    11191130                corrmax= max(max(result_conv));
    11201131                result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
     
    11761187%the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology
    11771188%http://urapiv.wordpress.com
    1178 peaky = y;
    1179 if y < npy && y > 1
     1189peaky = y;peakx=x;
     1190if y < npy && y > 1 && x < npx-1 && x > 1
    11801191    f0 = log(result_conv(y,x));
    11811192    f1 = log(result_conv(y-1,x));
    11821193    f2 = log(result_conv(y+1,x));
    11831194    peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
    1184 else
    1185     F=1; % warning flag for vector truncated by the limited search box
    1186 end
    1187 peakx=x;
    1188 if x < npx-1 && x > 1
    1189     f0 = log(result_conv(y,x));
    11901195    f1 = log(result_conv(y,x-1));
    11911196    f2 = log(result_conv(y,x+1));
     
    11941199    F=1; % warning flag for vector truncated by the limited search box
    11951200end
     1201
    11961202vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
    11971203
  • trunk/src/series/filter_time.m

    r1157 r1158  
    195195    % DataOut.Wfilter=squeeze(mean(Wblock(end-sumindex:end,:,:),1,'omitnan'));
    196196    %updating output the netcdf file
    197     [i1,i2,j1,j2] = get_file_index(i1_series{1}(index),[],PairString{1});
     197    [i1,i2,j1,j2] = get_file_index(i1_series{1}(index),[],PairString);
    198198    OutputFile=fullfile_uvmat(OutputPath,OutputDir,Param.InputTable{1,3},'.nc',NomType{1},i1,i2,j1,j2);
    199199    errormsg=struct2nc(OutputFile, DataOut);
Note: See TracChangeset for help on using the changeset viewer.