- Timestamp:
- Jul 16, 2024, 9:57:54 AM (3 months ago)
- Location:
- trunk/src
- Files:
-
- 1 added
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/script_delete_rdvision.m
r1159 r1160 85 85 if exist(ExtractFolder,'dir') && exist(PngFolder,'dir') 86 86 filename_seq=fullfile(ExtractFolder,'im.seq'); 87 NumberOfFrames=0; 87 88 try 88 s=ini2struct(filename_seq); 89 FileInfo=s.sequenceSettings; 90 if isfield(s.sequenceSettings,'numberoffiles') 91 NumberOfFrames=str2double(s.sequenceSettings.numberoffiles); 92 else 93 status='bad seq file'; 94 end 89 s=ini2struct(filename_seq); 90 FileInfo=s.sequenceSettings; 91 92 if isfield(s.sequenceSettings,'numberoffiles') 93 NumberOfFrames=str2double(s.sequenceSettings.numberoffiles); 94 else 95 status='bad seq file'; 96 end 95 97 catch ME 96 98 disp(['error in ' filename_seq]) 97 99 end 98 DirPng=dir(PngFolder); 100 DirPng=dir(PngFolder); 101 99 102 if numel(DirPng)==NumberOfFrames+2 100 Checkdelete=1; 103 PngCells=struct2cell(DirPng); 104 mm=cell2mat(PngCells(4,:));% check the sizes of extracted images 105 sizemax=max(mm); 106 sizemin=min(mm(3:end)); 107 disp(['max size(Mbytes)=' num2str(sizemax/1000000)]); 108 if min(mm(3:end))<0.9*sizemax 109 status=['WARNING' 'min size(Mbytes)=' num2str(sizemin/1000000)]; 110 else 111 Checkdelete=1;% approve deletion of the source bin files 112 end 101 113 else 102 114 status=['extraction not finished,' num2str(numel(DirPng)-2) ' images extracted']; … … 141 153 end 142 154 end 155 'END' -
trunk/src/series/civ_3D.m
r1158 r1160 16 16 % if absent, the fct looks for input data in Param.ActionInput (test mode) 17 17 % Param.OutputSubDir: sets the folder name of output file(s, 18 % if absent no file is produced, result in the output structure Data (test mode)19 18 % Param.ActionInput: substructure with the parameters provided by the GUI civ_input 20 19 % .Civ1: parameters for civ1cc … … 185 184 Data.VarAttribute{1}.Role='coord_z'; 186 185 Data.VarAttribute{2}.Role='coord_x'; 187 ata.VarAttribute{3}.Role='coord_y';186 Data.VarAttribute{3}.Role='coord_y'; 188 187 Data.VarAttribute{4}.Role='vector_x'; 189 188 Data.VarAttribute{5}.Role='vector_y'; … … 310 309 % calculate velocity data at the first position z, corresponding to j_index=par_civ1.Dz 311 310 [Data.Civ1_X(1,:,:),Data.Civ1_Y(1,:,:), Data.Civ1_U(1,:,:), Data.Civ1_V(1,:,:),Data.Civ1_W(1,:,:),... 312 Data.Civ1_C(1,:,:), Data.Civ1_FF(1,:,:), result_conv, errormsg] = civ3D (par_civ1);311 Data.Civ1_C(1,:,:), Data.Civ1_FF(1,:,:), ~, errormsg] = civ3D (par_civ1); 313 312 if ~isempty(errormsg) 314 313 disp_uvmat('ERROR',errormsg,checkrun) 315 314 return 316 315 end 317 for z_index=2:floor((NbSlice-1)/par_civ1.Dz)% loop on slices 316 %for z_index=2:floor((NbSlice-1)/par_civ1.Dz)% loop on slices 317 for z_index=2:floor((NbSlice-SearchRange_z)/par_civ1.Dz)% loop on slices 318 318 par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz);%shift the indices in the image block upward by par_civ1.Dz 319 319 par_civ1.ImageB=circshift(par_civ1.ImageA,-par_civ1.Dz); … … 336 336 end 337 337 end 338 Data.Civ1_ C=double(Data.Civ1_C);339 Data.Civ1_ C(Data.Civ1_C==Inf)=0;338 Data.Civ1_V=-Data.Civ1_V;%reverse v 339 Data.Civ1_Y=npy-Data.Civ1_Y+1;%reverse y 340 340 [npz,npy,npx]=size(Data.Civ1_X); 341 341 Data.Coord_z=SearchRange_z+1:par_civ1.Dz:NbSlice-1; 342 Data.Coord_y=1:npy;343 Data.Coord_x=1:npx;342 % Data.Coord_y=flip(1:npy); 343 % Data.Coord_x=1:npx; 344 344 % case of mask TO ADAPT 345 345 if par_civ1.CheckMask&&~isempty(par_civ1.Mask) … … 445 445 % perform Patch calculation using the UVMAT fct 'filter_tps' 446 446 447 [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps, tild,Ures, Vres,tild,FFres]=...447 [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,~,Ures, Vres,~,FFres]=... 448 448 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); 449 449 Data.Civ1_U_smooth(ind_good)=Ures;% take the interpolated (smoothed) velocity values for good vectors, keep civ1 data for the other … … 659 659 if exist('ncfile','var') 660 660 CivFile=ncfile; 661 [Data, tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and detect_false1 data in the existing netcdf file661 [Data,~,~,errormsg]=nc2struct(CivFile);%read civ1 and detect_false1 data in the existing netcdf file 662 662 if ~isempty(errormsg) 663 663 disp_uvmat('ERROR',errormsg,checkrun) 664 664 return 665 665 end 666 % elseif isfield(Param,'Civ2_X')% use Civ2 data as input in Param (test mode)667 % Data.ListGlobalAttribute={};668 % Data.ListVarName={};669 % Data.VarDimName={};670 % Data.Civ2_X=Param.Civ2_X;671 % Data.Civ2_Y=Param.Civ2_Y;672 % Data.Civ2_U=Param.Civ2_U;673 % Data.Civ2_V=Param.Civ2_V;674 % Data.Civ2_FF=Param.Civ2_FF;675 666 end 676 667 end … … 787 778 788 779 function [xtable,ytable,utable,vtable,wtable,ctable,FF,result_conv,errormsg] = civ3D (par_civ) 780 %% check image sizes 781 [npz,npy_ima,npx_ima]=size(par_civ.ImageA);% npz = 782 if ~isequal(size(par_civ.ImageB),[npz npy_ima npx_ima]) 783 errormsg='image pair with unequal size'; 784 return 785 end 789 786 790 787 %% prepare measurement grid 791 nbinterv_x=floor(( par_civ.ImageWidth-1)/par_civ.Dx);788 nbinterv_x=floor((npx_ima-1)/par_civ.Dx); 792 789 gridlength_x=nbinterv_x*par_civ.Dx; 793 minix=ceil(( par_civ.ImageWidth-gridlength_x)/2);794 nbinterv_y=floor(( par_civ.ImageHeight-1)/par_civ.Dy);790 minix=ceil((npx_ima-gridlength_x)/2); 791 nbinterv_y=floor((npy_ima-1)/par_civ.Dy); 795 792 gridlength_y=nbinterv_y*par_civ.Dy; 796 miniy=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); 793 miniy=ceil((npy_ima-gridlength_y)/2); 794 [GridX,GridY]=meshgrid(minix:par_civ.Dx:npx_ima-1,miniy:par_civ.Dy:npy_ima-1); 805 795 par_civ.Grid(:,:,1)=GridX; 806 796 par_civ.Grid(:,:,2)=GridY;% increases with array index, 807 797 [nbvec_y,nbvec_x,~]=size(par_civ.Grid); 808 % 809 % 810 % minix=floor(par_civ.Dx/2)-0.5; 811 % maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx); 812 % miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy 813 % maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy); 814 % [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy); 815 % par_civ.Grid(:,1)=reshape(GridX,[],1); 816 % par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index 798 817 799 818 800 %% prepare correlation and search boxes … … 826 808 shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases 827 809 if numel(shiftx)==1% case of a unique shift for the whole field( civ1) 828 shiftx=shiftx*ones(nbvec_y,nbvec_x ,1);829 shifty=shifty*ones(nbvec_y,nbvec_x ,1);810 shiftx=shiftx*ones(nbvec_y,nbvec_x); 811 shifty=shifty*ones(nbvec_y,nbvec_x); 830 812 end 831 813 832 814 %% Array initialisation and default output if par_civ.CorrSmooth=0 (just the grid calculated, no civ computation) 833 815 xtable=round(par_civ.Grid(:,:,1)+0.5)-0.5; 834 ytable=round(par_civ.ImageHeight-par_civ.Grid(:,:,2)+0.5)-0.5;% y index corresponding to the position in image coordiantes 835 utable=shiftx;%zeros(nbvec,1); 836 vtable=shifty;%zeros(nbvec,1); 816 %ytable=round(npy_ima-par_civ.Grid(:,:,2)+0.5)-0.5;% y index corresponding to the position in image coordiantes 817 ytable=round(par_civ.Grid(:,:,2)+0.5)-0.5; 818 utable=shiftx; 819 vtable=shifty; 837 820 wtable=zeros(size(utable)); 838 ctable=zeros(nbvec_y,nbvec_x ,1);839 FF= zeros(nbvec_y,nbvec_x,1);821 ctable=zeros(nbvec_y,nbvec_x); 822 FF=false(nbvec_y,nbvec_x); 840 823 result_conv=[]; 841 824 errormsg=''; … … 844 827 check_MinIma=isfield(par_civ,'MinIma');% test for image luminosity threshold 845 828 check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma); 846 847 [npz,npy_ima,npx_ima]=size(par_civ.ImageA);848 if ~isequal(size(par_civ.ImageB),[npz npy_ima npx_ima])849 errormsg='image pair with unequal size';850 return851 end852 829 853 830 %% Apply mask … … 860 837 checkmask=0; 861 838 MinA=min(min(min(par_civ.ImageA))); 862 %MinB=min(min(par_civ.ImageB));863 %check_undefined=false(size(par_civ.ImageA));864 839 if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask) 865 840 checkmask=1; … … 868 843 return 869 844 end 870 check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 ); 845 check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 );% logical of the size of the image block =true for masked pixels 871 846 end 872 847 … … 875 850 sum_square=1;% default 876 851 mesh=1;% default 877 CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;878 if CheckDeformation879 mesh=0.25;%mesh in pixels for subpixel image interpolation (x 4 in each direction)880 par_civ.CorrSmooth=2;% use SUBPIX2DGAUSS (take into account more points near the max)881 end852 % CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1; 853 % if CheckDeformation 854 % mesh=0.25;%mesh in pixels for subpixel image interpolation (x 4 in each direction) 855 % par_civ.CorrSmooth=2;% use SUBPIX2DGAUSS (take into account more points near the max) 856 % end 882 857 883 858 if par_civ.CorrSmooth~=0 % par_civ.CorrSmooth=0 implies no civ computation (just input image and grid points given) … … 885 860 for ivec_y=1:nbvec_y 886 861 iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5);% xindex on the image A for the middle of the correlation box 887 jref=round( par_civ.ImageHeight-par_civ.Grid(ivec_y,ivec_x,2)+0.5);% j index for the middle of the correlation box in the image A862 jref=round(npy_ima-par_civ.Grid(ivec_y,ivec_x,2)+0.5);% j index for the middle of the correlation box in the image A 888 863 subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage 889 864 subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage … … 892 867 image1_crop=MinA*ones(npz,numel(subrange1_y),numel(subrange1_x));% default value=min of image A 893 868 image2_crop=MinA*ones(npz,numel(subrange2_y),numel(subrange2_x));% default value=min of image A 894 check1_x=subrange1_x>=1 & subrange1_x<= par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1895 check1_y=subrange1_y>=1 & subrange1_y<= par_civ.ImageHeight;896 check2_x=subrange2_x>=1 & subrange2_x<= par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2897 check2_y=subrange2_y>=1 & subrange2_y<= par_civ.ImageHeight;869 check1_x=subrange1_x>=1 & subrange1_x<=npx_ima;% check which points in the subimage 1 are contained in the initial image 1 870 check1_y=subrange1_y>=1 & subrange1_y<=npy_ima; 871 check2_x=subrange2_x>=1 & subrange2_x<=npx_ima;% check which points in the subimage 2 are contained in the initial image 2 872 check2_y=subrange2_y>=1 & subrange2_y<=npy_ima; 898 873 image1_crop(:,check1_y,check1_x)=par_civ.ImageA(:,subrange1_y(check1_y),subrange1_x(check1_x));%extract a subimage (correlation box) from image A 899 874 image2_crop(:,check2_y,check2_x)=par_civ.ImageB(:,subrange2_y(check2_y),subrange2_x(check2_x));%extract a larger subimage (search box) from image B 900 if checkmask 901 mask1_crop= ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask902 mask2_crop= ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask903 mask1_crop( check1_y,check1_x)=check_undefined(subrange1_y(check1_y),subrange1_x(check1_x));%extract a mask subimage (correlation box) from image A904 mask2_crop( check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from image B905 sizemask=sum(sum(mask1_crop ))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image875 if checkmask% case of mask images 876 mask1_crop=true(npz,2*iby2+1,2*ibx2+1);% default value=1 for mask 877 mask2_crop=true(npz,2*isy2+1,2*isx2+1);% default value=1 for mask 878 mask1_crop(npz,check1_y,check1_x)=check_undefined(npz,subrange1_y(check1_y),subrange1_x(check1_x));%extract a mask subimage (correlation box) from image A 879 mask2_crop(npz,check2_y,check2_x)=check_undefined(npz,subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from image B 880 sizemask=sum(sum(mask1_crop,2),3)/(npz*(2*ibx2+1)*(2*iby2+1));%size of the masked part relative to the correlation sub-image 906 881 if sizemask > 1/2% eliminate point if more than half of the correlation box is masked 907 882 FF(ivec_y,ivec_x)=1; % … … 909 884 vtable(ivec_y,ivec_x)=NaN; 910 885 else 911 FF(ivec_y,ivec_x)=0;912 886 image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1) 913 887 image2_crop=image2_crop.*~mask2_crop; 914 image1_mean=mean(mean(image1_crop ))/(1-sizemask);915 image2_mean=mean(mean(image2_crop ))/(1-sizemask);888 image1_mean=mean(mean(image1_crop,2),3)/(1-sizemask); 889 image2_mean=mean(mean(image2_crop,2),3)/(1-sizemask); 916 890 end 917 891 else 918 image1_mean=mean(mean(image1_crop ));919 image2_mean=mean(mean(image2_crop ));920 end 921 %threshold on image minimum892 image1_mean=mean(mean(image1_crop,2),3); 893 image2_mean=mean(mean(image2_crop,2),3); 894 end 895 922 896 if FF(ivec_y,ivec_x)==0 923 897 if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma) 924 FF(ivec_y,ivec_x)=1; 925 %threshold on image maximum898 FF(ivec_y,ivec_x)=1; %threshold on image minimum 899 926 900 elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma) 927 FF(ivec_y,ivec_x)=1; 901 FF(ivec_y,ivec_x)=1; %threshold on image maximum 928 902 end 929 903 if FF(ivec_y,ivec_x)==1 930 utable(ivec_y,ivec_x)=NaN; 904 utable(ivec_y,ivec_x)=NaN; 931 905 vtable(ivec_y,ivec_x)=NaN; 932 906 else 933 % mask934 if checkmask935 image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts936 image2_crop=(image2_crop-image2_mean).*~mask2_crop;937 else938 939 940 end907 % case of mask 908 % if checkmask 909 % image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts 910 % image2_crop=(image2_crop-image2_mean).*~mask2_crop;% TO MODIFY !!!!!!!!!!!!!!!!!!!!! 911 % else 912 % % image1_crop=(image1_crop-image1_mean); 913 % % image2_crop=(image2_crop-image2_mean); 914 % end 941 915 942 916 npxycorr=par_civ.SearchBoxSize(1:2)-par_civ.CorrBoxSize(1:2)+1; 943 917 result_conv=zeros([par_civ.SearchBoxSize(3) npxycorr]);%initialise the conv product 944 918 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) 949 subima2=squeeze(image2_crop(kz,:,:)); 950 subima1=squeeze(image1_crop(kref,:,:)); 919 xk=ones(npz,1);%initialise the x displacement corresponding to the max corr 920 yk=ones(npz,1);%initialise the y displacement corresponding to the max corr 921 subima1=squeeze(image1_crop(kref,:,:))-image1_mean(kref); 922 for kz=1:npz 923 subima2=squeeze(image2_crop(kz,:,:))-image2_mean(kz); 951 924 correl_xy=conv2(subima2,flip(flip(subima1,2),1),'valid'); 952 result_conv(kz,:,:)= 255*correl_xy;925 result_conv(kz,:,:)=correl_xy; 953 926 max_xy(kz)=max(max(correl_xy)); 954 % max_xy(kz)=max(max(correl_xy)); 955 [yk(kz),xk(kz)]=find(correl_xy==max_xy(kz),1); 927 [yk(kz),xk(kz)]=find(correl_xy>=0.999*max_xy(kz),1);%find the indices of the maximum, use 0.999 to avoid rounding errors 956 928 end 957 929 [corrmax,dz]=max(max_xy); 930 result_conv=result_conv*255/corrmax;% normalise with a max =255 958 931 dx=xk(dz); 959 932 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,:,:)))); 933 subimage2_crop=squeeze(image2_crop(dz,dy:dy+2*iby2/mesh,dx:dx+2*ibx2/mesh))-image2_mean(dz);%subimage of image 2 corresponding to the optimum displacement of first image 934 sum_square=sum(sum(subima1.*subima1)); 972 935 sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2 973 936 sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation 937 % if ivec_x==3 && ivec_y==4 938 % 'test' 939 % end 974 940 if ~isempty(dz)&& ~isempty(dy) && ~isempty(dx) 975 941 if par_civ.CorrSmooth==1 … … 987 953 iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector 988 954 jref=round(ytable(ivec_y,ivec_x)+0.5); 989 955 % if utable(ivec_y,ivec_x)==0 && vtable(ivec_y,ivec_x)==0 956 % 'test' 957 % end 990 958 % eliminate vectors located in the mask 991 959 if checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100)) … … 1017 985 1018 986 [npz,npy,npx]=size(result_conv); 1019 result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1 ( =0 by discretisation,to avoid divergence in the log)987 result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1 (to avoid divergence in the log) 1020 988 %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology 1021 989 %http://urapiv.wordpress.com 1022 990 FF=false; 1023 991 peakz=z;peaky=y;peakx=x; 1024 if z < npz && z > 1 && y < npy && y > 1 && x < npx && x > 11025 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);1037 else1038 FF=true;1039 end992 % if z < npz && z > 1 && y < npy && y > 1 && x < npx && x > 1 993 % f0 = log(result_conv(z,y,x)); 994 % f1 = log(result_conv(z-1,y,x)); 995 % f2 = log(result_conv(z+1,y,x)); 996 % peakz = peakz+ (f1-f2)/(2*f1-4*f0+2*f2); 997 % 998 % f1 = log(result_conv(z,y-1,x)); 999 % f2 = log(result_conv(z,y+1,x)); 1000 % peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2); 1001 % 1002 % f1 = log(result_conv(z,y,x-1)); 1003 % f2 = log(result_conv(z,y,x+1)); 1004 % peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2); 1005 % else 1006 % FF=true; 1007 % end 1040 1008 1041 1009 vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1 peakz-floor(npz/2)-1]; -
trunk/src/series/civ_series.m
r1158 r1160 952 952 % .DVDY 953 953 954 function [xtable,ytable,utable,vtable,ctable,F ,result_conv,errormsg] = civ (par_civ)954 function [xtable,ytable,utable,vtable,ctable,FF,result_conv,errormsg] = civ (par_civ) 955 955 956 956 %% prepare measurement grid … … 972 972 par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index 973 973 % 974 %975 974 % minix=floor(par_civ.Dx/2)-0.5; 976 975 % maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx); … … 1001 1000 vtable=shifty;%zeros(nbvec,1); 1002 1001 ctable=zeros(nbvec,1); 1003 F =zeros(nbvec,1);1002 FF=zeros(nbvec,1); 1004 1003 result_conv=[]; 1005 1004 errormsg=''; … … 1058 1057 iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box 1059 1058 jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% j index for the middle of the correlation box in the image A 1060 F (ivec)=0;1059 FF(ivec)=0; 1061 1060 subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage 1062 1061 subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage … … 1078 1077 sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image 1079 1078 if sizemask > 1/2% eliminate point if more than half of the correlation box is masked 1080 F (ivec)=1; %1079 FF(ivec)=1; % 1081 1080 utable(ivec)=NaN; 1082 1081 vtable(ivec)=NaN; … … 1092 1091 end 1093 1092 %threshold on image minimum 1094 if F (ivec)~=11093 if FF(ivec)==0 1095 1094 if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma) 1096 F (ivec)=1;1095 FF(ivec)=1; 1097 1096 %threshold on image maximum 1098 1097 elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma) 1099 F (ivec)=1;1100 end 1101 if F (ivec)==11098 FF(ivec)=1; 1099 end 1100 if FF(ivec)==1 1102 1101 utable(ivec)=NaN; 1103 1102 vtable(ivec)=NaN; … … 1127 1126 sum_square=sum(sum(image1_crop.*image1_crop)); 1128 1127 %reference: Oliver Pust, PIV: Direct Cross-Correlation 1128 %%%%%% correlation calculation 1129 1129 result_conv= conv2(image2_crop,flip(flip(image1_crop,2),1),'valid'); 1130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1130 1131 corrmax= max(max(result_conv)); 1131 1132 result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 … … 1138 1139 try 1139 1140 if par_civ.CorrSmooth==1 1140 [vector,F (ivec)] = SUBPIXGAUSS (result_conv,x,y);1141 [vector,FF(ivec)] = SUBPIXGAUSS (result_conv,x,y); 1141 1142 elseif par_civ.CorrSmooth==2 1142 [vector,F (ivec)] = SUBPIX2DGAUSS (result_conv,x,y);1143 [vector,FF(ivec)] = SUBPIX2DGAUSS (result_conv,x,y); 1143 1144 else 1144 [vector,F (ivec)] = quadr_fit(result_conv,x,y);1145 [vector,FF(ivec)] = quadr_fit(result_conv,x,y); 1145 1146 end 1146 1147 utable(ivec)=vector(1)*mesh+shiftx(ivec); … … 1154 1155 utable(ivec)=0; 1155 1156 vtable(ivec)=0; 1156 F (ivec)=1;1157 FF(ivec)=1; 1157 1158 end 1158 1159 ctable(ivec)=corrmax/sum_square;% correlation value 1159 1160 catch ME 1160 F (ivec)=1;1161 FF(ivec)=1; 1161 1162 disp(ME.message) 1162 1163 end 1163 1164 else 1164 F (ivec)=1;1165 FF(ivec)=1; 1165 1166 end 1166 1167 end … … 1188 1189 %http://urapiv.wordpress.com 1189 1190 peaky = y;peakx=x; 1190 if y < npy && y > 1 && x < npx-1 && x > 11191 f0 = log(result_conv(y,x));1192 f1 = log(result_conv(y-1,x));1193 f2 = log(result_conv(y+1,x));1194 peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);1195 f1 = log(result_conv(y,x-1));1196 f2 = log(result_conv(y,x+1));1197 peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);1198 else1199 F=1; % warning flag for vector truncated by the limited search box1200 end1191 % if y < npy && y > 1 && x < npx-1 && x > 1 1192 % f0 = log(result_conv(y,x)); 1193 % f1 = log(result_conv(y-1,x)); 1194 % f2 = log(result_conv(y+1,x)); 1195 % peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2); 1196 % f1 = log(result_conv(y,x-1)); 1197 % f2 = log(result_conv(y,x+1)); 1198 % peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2); 1199 % else 1200 % F=1; % warning flag for vector truncated by the limited search box 1201 % end 1201 1202 1202 1203 vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
Note: See TracChangeset
for help on using the changeset viewer.