Changeset 856 for trunk/src/series/civ_series.m
- Timestamp:
- Jan 26, 2015, 12:37:56 AM (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/series/civ_series.m
r855 r856 58 58 Data.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default) 59 59 Data.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default) 60 if isfield(Data,'ActionInput') && strcmp(Data.ActionInput.PairIndices.ListPairMode,'pair j1-j2')60 if isfield(Data,'ActionInput') && isfield(Data.ActionInput,'PairIndices')&& strcmp(Data.ActionInput.PairIndices.ListPairMode,'pair j1-j2') 61 61 Data.Desable_j_index='on';% hide the j index in series (set by the pair choice in civ_input) 62 62 end … … 174 174 j1_series_Civ2=j1_series_Civ1; 175 175 j2_series_Civ2=j2_series_Civ1; 176 NomTypeNc= NomType;176 NomTypeNc=Param.InputTable{1,4}; 177 177 case 'PIV volume' 178 178 % TODO, TODO … … 329 329 330 330 %%%%% MAIN LOOP %%%%%% 331 maskname='';% initiate the mask name 331 332 for ifield=1:NbField 332 333 if ~isempty(RUNHandle)% update the waitbar in interactive mode with GUI series (checkrun=1) … … 360 361 if isfield (Param.ActionInput,'Civ1') 361 362 par_civ1=Param.ActionInput.Civ1; 362 if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1 )363 if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B) 363 364 try 364 365 ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(ifield),[],j1_series_Civ1(ifield)); … … 422 423 Data.VarAttribute{4}.Role='vector_y'; 423 424 Data.VarAttribute{5}.Role='warnflag'; 424 425 if par_civ1.CheckMask&&~isempty(par_civ1.Mask) 426 if strcmp(maskname,par_civ1.Mask)% mask exist, not already read in civ1 427 par_civ1.Mask=mask; %use mask already opened 428 else 429 par_civ1.Mask=imread(par_civ1.Mask);%update the mask, an store it for future use 430 mask=par_civ1.Mask; 431 maskname=par_civ1.Mask; 432 end 433 end 425 434 if strcmp(Param.ActionInput.ListCompareMode, 'PIV volume') 426 435 Data.ListVarName=[Data.ListVarName 'Civ1_Z']; … … 429 438 for ivol=1:NbSlice 430 439 % caluclate velocity data (y and v in indices, reverse to y component) 431 [xtable ytable utable vtable ctable F result_converrormsg] = civ (par_civ1);440 [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1); 432 441 if ~isempty(errormsg) 433 442 disp_uvmat('ERROR',errormsg,checkrun) … … 444 453 else %usual PIV 445 454 % caluclate velocity data (y and v in indices, reverse to y component) 446 [xtable ytable utable vtable ctable F result_converrormsg] = civ (par_civ1);455 [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1); 447 456 if ~isempty(errormsg) 448 457 disp_uvmat('ERROR',errormsg,checkrun) … … 473 482 [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file 474 483 end 475 else 484 elseif isfield(Param,'Civ1_X') 476 485 Data.ListGlobalAttribute={}; 477 486 Data.ListVarName={}; … … 572 581 if isfield (Param.ActionInput,'Civ2') 573 582 par_civ2=Param.ActionInput.Civ2; 574 par_civ2.ImageA=[]; 575 par_civ2.ImageB=[]; 576 i1=i1_series_Civ2(ifield); 577 i2=i1; 578 if ~isempty(i2_series_Civ2) 579 i2=i2_series_Civ2(ifield); 580 end 581 j1=1; 582 if ~isempty(j1_series_Civ2) 583 j1=j1_series_Civ2(ifield); 584 end 585 j2=j1; 586 if ~isempty(j2_series_Civ2) 587 j2=j2_series_Civ2(ifield); 588 end 589 ImageName_A_Civ2=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1,[],j1); 590 591 if strcmp(ImageName_A_Civ2,ImageName_A) && isequal(FrameIndex_A_Civ1(ifield),FrameIndex_A_Civ2(ifield)) 592 par_civ2.ImageA=par_civ1.ImageA; 593 else 594 [par_civ2.ImageA,VideoObject_A] = read_image(ImageName_A_Civ2,FileType_A,VideoObject_A,FrameIndex_A_Civ2(ifield)); 595 end 596 ImageName_B_Civ2=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2,[],j2); 597 if strcmp(ImageName_B_Civ2,ImageName_B) && isequal(FrameIndex_B_Civ1(ifield),FrameIndex_B_Civ2) 598 par_civ2.ImageB=par_civ1.ImageB; 599 else 600 [par_civ2.ImageB,VideoObject_B] = read_image(ImageName_B_Civ2,FileType_B,VideoObject_B,FrameIndex_B_Civ2(ifield)); 601 end 602 603 ncfile=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i1,i2,j1,j2); 604 par_civ2.ImageWidth=FileInfo_A.Width; 605 par_civ2.ImageHeight=FileInfo_A.Height; 606 607 if isfield(par_civ2,'Grid')% grid points set as input file 608 if ischar(par_civ2.Grid)%read the grid file if the input is a file name 609 par_civ2.Grid=dlmread(par_civ2.Grid); 610 par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file) 611 end 612 else% automatic grid 613 minix=floor(par_civ2.Dx/2)-0.5; 614 maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx); 615 miniy=floor(par_civ2.Dy/2)-0.5; 616 maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy); 617 [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy); 618 par_civ2.Grid(:,1)=reshape(GridX,[],1); 619 par_civ2.Grid(:,2)=reshape(GridY,[],1); 620 end 621 Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data 622 Shifty=zeros(size(par_civ2.Grid,1),1); 623 nbval=zeros(size(par_civ2.Grid,1),1); 583 if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B) 584 par_civ2.ImageA=[]; 585 par_civ2.ImageB=[]; 586 i1=i1_series_Civ2(ifield); 587 i2=i1; 588 if ~isempty(i2_series_Civ2) 589 i2=i2_series_Civ2(ifield); 590 end 591 j1=1; 592 if ~isempty(j1_series_Civ2) 593 j1=j1_series_Civ2(ifield); 594 end 595 j2=j1; 596 if ~isempty(j2_series_Civ2) 597 j2=j2_series_Civ2(ifield); 598 end 599 ImageName_A_Civ2=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1,[],j1); 600 601 if strcmp(ImageName_A_Civ2,ImageName_A) && isequal(FrameIndex_A_Civ1(ifield),FrameIndex_A_Civ2(ifield)) 602 par_civ2.ImageA=par_civ1.ImageA; 603 else 604 [par_civ2.ImageA,VideoObject_A] = read_image(ImageName_A_Civ2,FileType_A,VideoObject_A,FrameIndex_A_Civ2(ifield)); 605 end 606 ImageName_B_Civ2=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2,[],j2); 607 if strcmp(ImageName_B_Civ2,ImageName_B) && isequal(FrameIndex_B_Civ1(ifield),FrameIndex_B_Civ2) 608 par_civ2.ImageB=par_civ1.ImageB; 609 else 610 [par_civ2.ImageB,VideoObject_B] = read_image(ImageName_B_Civ2,FileType_B,VideoObject_B,FrameIndex_B_Civ2(ifield)); 611 end 612 if strcmp(Param.ActionInput.ListCompareMode,'PIV') 613 ncfile=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i1,i2,j1,j2); 614 else % displacement 615 ncfile=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i2,[],j2); 616 end 617 par_civ2.ImageWidth=FileInfo_A.Width; 618 par_civ2.ImageHeight=FileInfo_A.Height; 619 if isfield(par_civ2,'Grid')% grid points set as input file 620 if ischar(par_civ2.Grid)%read the grid file if the input is a file name 621 par_civ2.Grid=dlmread(par_civ2.Grid); 622 par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file) 623 end 624 else% automatic grid 625 minix=floor(par_civ2.Dx/2)-0.5; 626 maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx); 627 miniy=floor(par_civ2.Dy/2)-0.5; 628 maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy); 629 [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy); 630 par_civ2.Grid(:,1)=reshape(GridX,[],1); 631 par_civ2.Grid(:,2)=reshape(GridY,[],1); 632 end 633 end 624 634 if par_civ2.CheckDeformation 625 635 DUDX=zeros(size(par_civ2.Grid,1),1); … … 629 639 end 630 640 631 % get the guess from patch1 or patch2 (case 'CheckCiv3') 632 if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from patch2 633 SubRange= Data.Civ2_SubRange; 634 NbCentres=Data.Civ2_NbCentres; 635 Coord_tps=Data.Civ2_Coord_tps; 636 U_tps=Data.Civ2_U_tps; 637 V_tps=Data.Civ2_V_tps; 638 CivStage=Data.CivStage; 639 Civ1_Dt=Data.Civ2_Dt; 640 Data=[];%reinitialise the result structure Data 641 Data.ListGlobalAttribute={'Conventions','Program','CivStage'}; 642 Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes 643 Data.Program='civ_series'; 644 Data.CivStage=CivStage; 645 Data.ListVarName={}; 646 Data.VarDimName={}; 647 else % get the guess from patch1 648 SubRange= Data.Civ1_SubRange; 649 NbCentres=Data.Civ1_NbCentres; 650 Coord_tps=Data.Civ1_Coord_tps; 651 U_tps=Data.Civ1_U_tps; 652 V_tps=Data.Civ1_V_tps; 653 Civ1_Dt=Data.Civ1_Dt; 654 end 655 NbSubDomain=size(SubRange,3); 641 % get the guess from patch1 or patch2 (case 'CheckCiv3') 642 if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B) 643 if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from patch2 644 SubRange= Data.Civ2_SubRange; 645 NbCentres=Data.Civ2_NbCentres; 646 Coord_tps=Data.Civ2_Coord_tps; 647 U_tps=Data.Civ2_U_tps; 648 V_tps=Data.Civ2_V_tps; 649 CivStage=Data.CivStage; 650 Civ1_Dt=Data.Civ2_Dt; 651 Data=[];%reinitialise the result structure Data 652 Data.ListGlobalAttribute={'Conventions','Program','CivStage'}; 653 Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes 654 Data.Program='civ_series'; 655 Data.CivStage=CivStage; 656 Data.ListVarName={}; 657 Data.VarDimName={}; 658 else % get the guess from patch1 659 SubRange= Data.Civ1_SubRange; 660 NbCentres=Data.Civ1_NbCentres; 661 Coord_tps=Data.Civ1_Coord_tps; 662 U_tps=Data.Civ1_U_tps; 663 V_tps=Data.Civ1_V_tps; 664 Civ1_Dt=Data.Civ1_Dt; 665 end 666 else 667 SubRange= par_civ2.Civ1_SubRange; 668 NbCentres=par_civ2.Civ1_NbCentres; 669 Coord_tps=par_civ2.Civ1_Coord_tps; 670 U_tps=par_civ2.Civ1_U_tps; 671 V_tps=par_civ2.Civ1_V_tps; 672 Civ1_Dt=par_civ2.Civ1_Dt; 673 Data.ListVarName={}; 674 Data.VarDimName={}; 675 end 676 Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data 677 Shifty=zeros(size(par_civ2.Grid,1),1); 678 nbval=zeros(size(par_civ2.Grid,1),1);% nbre of interpolated values at each grid point (from the different patch subdomains) 679 NbSubDomain=size(SubRange,3); 656 680 for isub=1:NbSubDomain% for each sub-domain of Patch1 657 681 nbvec_sub=NbCentres(isub);% nbre of Civ vectors in the subdomain … … 660 684 epoints = par_civ2.Grid(ind_sel,:);% coordinates of interpolation sites 661 685 ctrs=Coord_tps(1:nbvec_sub,:,isub) ;%(=initial points) ctrs 662 nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eac ninterpolation point (in case of subdomain overlap)663 EM = tps_eval(epoints,ctrs); 686 nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for each interpolation point (in case of subdomain overlap) 687 EM = tps_eval(epoints,ctrs); 664 688 Shiftx(ind_sel)=Shiftx(ind_sel)+EM*U_tps(1:nbvec_sub+3,isub); 665 689 Shifty(ind_sel)=Shifty(ind_sel)+EM*V_tps(1:nbvec_sub+3,isub); … … 672 696 end 673 697 end 674 mask=''; 675 if par_civ2.CheckMask&&~isempty(par_civ2.Mask)&& ~strcmp(maskname,par_civ2.Mask)% mask exist, not already read in civ1 676 mask=imread(par_civ2.Mask); 698 if par_civ2.CheckMask&&~isempty(par_civ2.Mask) 699 if strcmp(maskname,par_civ2.Mask)% mask exist, not already read in civ1 700 par_civ2.Mask=mask; %use mask already opened 701 else 702 par_civ2.Mask=imread(par_civ2.Mask);%update the mask, and store it for future use 703 mask=par_civ2.Mask; 704 maskname=par_civ2.Mask; 705 end 677 706 end 678 707 ibx2=ceil(par_civ2.CorrBoxSize(1)/2); … … 680 709 par_civ2.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess 681 710 par_civ2.SearchBoxSize(2)=2*iby2+9; 682 Civ2_Dt=time(i2+1,j2+1)-time(i1+1,j1+1); 711 if strcmp(Param.ActionInput.ListCompareMode,'displacement') 712 Civ1_Dt=1; 713 Civ2_Dt=1; 714 elseif exist('time','var') 715 Civ2_Dt=time(i2+1,j2+1)-time(i1+1,j1+1); 716 else 717 Civ2_Dt=1;%TO CHECK, TEST 718 end 683 719 par_civ2.SearchBoxShift=(Civ2_Dt/Civ1_Dt)*[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)]; 684 720 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];% grid taken at the extrapolated origin of the displacement vectors … … 690 726 end 691 727 % calculate velocity data (y and v in indices, reverse to y component) 692 [xtable ytable utable vtable ctable F] = civ (par_civ2); 728 [xtable, ytable, utable, vtable, ctable, F,result_conv,errormsg] = civ (par_civ2); 729 693 730 list_param=(fieldnames(Param.ActionInput.Civ2))'; 694 731 Civ2_param=regexprep(list_param,'^.+','Civ2_$0');% insert 'Civ2_' before each string in list_param 695 732 Civ2_param=[{'Civ2_ImageA','Civ2_ImageB','Civ2_Time','Civ2_Dt'} Civ2_param]; %insert the names of the two input images 696 %indicate the values of all the global attributes in the output data 733 %indicate the values of all the global attributes in the output data 734 if exist('ImageName_A','var') 697 735 Data.Civ2_ImageA=ImageName_A; 698 736 Data.Civ2_ImageB=ImageName_B; 699 737 Data.Civ2_Time=(time(i2+1,j2+1)+time(i1+1,j1+1))/2; 700 738 Data.Civ2_Dt=Civ2_Dt; 701 % Data.Civ2_Time=1; 702 % Data.Civ2_Dt=1; 739 end 740 % Data.Civ2_Time=1; 741 % Data.Civ2_Dt=1; 703 742 for ilist=1:length(list_param) 704 743 Data.(Civ2_param{4+ilist})=Param.ActionInput.Civ2.(list_param{ilist}); … … 725 764 %% Fix2 726 765 if isfield (Param.ActionInput,'Fix2') 727 766 list_param=fieldnames(Param.ActionInput.Fix2)'; 728 767 Fix2_param=regexprep(list_param,'^.+','Fix2_$0');% insert 'Fix1_' before each string in ListFixParam 729 768 %indicate the values of all the global attributes in the output data … … 824 863 %INPUT: 825 864 % par_civ: structure of input parameters, with fields: 826 % .CorrBoxSize 827 % .SearchBoxSize 828 % .SearchBoxShift 829 % .ImageHeight 830 % .ImageWidth 831 % .Dx, Dy 832 % .Grid 833 % .Mask 834 % .MinIma 865 % .ImageA: first image for correlation (matrix) 866 % .ImageB: second image for correlation(matrix) 867 % .CorrBoxSize: 1,2 vector giving the size of the correlation box in x and y 868 % .SearchBoxSize: 1,2 vector giving the size of the search box in x and y 869 % .SearchBoxShift: 1,2 vector or 2 column matrix (for civ2) giving the shift of the search box in x and y 870 % .CorrSmooth: =1 or 2 determines the choice of the sub-pixel determination of the correlation max 871 % .ImageWidth: nb of pixels of the image in x 872 % .Dx, Dy: mesh for the PIV calculation 873 % .Grid: grid giving the PIV calculation points (alternative to .Dx .Dy) 874 % .Mask: name of a mask file or mask image matrix itself 875 % .MinIma: thresholds for image luminosity 835 876 % .MaxIma 836 % .image1:first image (matrix) 837 % image2: second image (matrix) 838 % ibx2,iby2: half size of the correlation box along x and y, in px (size=(2*iby2+1,2*ibx2+1) 839 % isx2,isy2: half size of the search box along x and y, in px (size=(2*isy2+1,2*isx2+1) 840 % shiftx, shifty: shift of the search box (in pixel index, yshift reversed) 841 % step: mesh of the measurement points (in px) 842 % subpixfinder=1 or 2 controls the curve fitting of the image correlation 843 % mask: =[] for no mask 844 % roi: 4 element vector defining a region of interest: x position, y position, width, height, (in image indices), for the whole image, roi=[]; 877 % .CheckDeformation=1 for subpixel interpolation and image deformation (linear transform) 878 % .DUDX: matrix of deformation obtained from patch at each grid point 879 % .DUDY 880 % .DVDX: 881 % .DVDY 882 845 883 function [xtable,ytable,utable,vtable,ctable,F,result_conv,errormsg] = civ (par_civ) 846 884 847 885 %% prepare measurement grid 848 886 if isfield(par_civ,'Grid')% grid points set as input 849 if ischar(par_civ.Grid)%read the drid file if the input is a file name887 if ischar(par_civ.Grid)%read the grid file if the input is a file name 850 888 par_civ.Grid=dlmread(par_civ.Grid); 851 889 par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file) 852 890 end 891 % else par_civ.Grid is already an array 853 892 else% automatic grid 854 893 minix=floor(par_civ.Dx/2)-0.5; … … 889 928 return % get the grid only, no civ calculation 890 929 elseif ischar(par_civ.Mask) 891 par_civ.Mask=imread(par_civ.Mask); 930 par_civ.Mask=imread(par_civ.Mask);% read the mask if not allready done 892 931 end 893 932 end … … 948 987 sum_square=1;% default 949 988 mesh=1;% default 950 CheckDecimal=isfield(par_civ,'CheckDecimal')&& par_civ.CheckDecimal==1; 951 if CheckDecimal 952 mesh=0.2;%mesh in pixels for subpixel image interpolation 953 CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1; 954 end 955 % vector=[0 0];%default 956 for ivec=1:nbvec 957 iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box 958 jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% yindex on the image B for the middle of the correlation box 959 %if ~(checkmask && par_civ.Mask(jref,iref)<=20) %velocity not set to zero by the black mask 960 % if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||... 961 % jref+shifty(ivec)-isy2<1||jref+shifty(ivec)+isy2>par_civ.ImageHeight|| iref+shiftx(ivec)-isx2<1 || iref+shiftx(ivec)+isx2>par_civ.ImageWidth % we are outside the image 962 % F(ivec)=3; 963 % else 964 F(ivec)=0; 965 subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage 966 subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage 967 subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage 968 subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage 969 image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A 970 image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A 971 check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1 972 check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight; 973 check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2 974 check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight; 975 976 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 977 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 978 image1_mean=mean(mean(image1_crop)); 979 image2_mean=mean(mean(image2_crop)); 980 %threshold on image minimum 981 if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma) 982 F(ivec)=3; 983 end 984 %threshold on image maximum 985 if check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma) 986 F(ivec)=3; 987 end 988 % end 989 if F(ivec)~=3 990 image1_crop=image1_crop-image1_mean;%substract the mean 991 image2_crop=image2_crop-image2_mean; 992 if CheckDecimal 993 xi=(1:mesh:size(image1_crop,2)); 994 yi=(1:mesh:size(image1_crop,1))'; 989 CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1; 990 if CheckDeformation 991 mesh=0.25;%mesh in pixels for subpixel image interpolation (x 4 in each direction) 992 end 993 if par_civ.CorrSmooth~=0 % par_civ.CorrSmooth=0 implies no civ computation (just input image and grid points viven) 994 for ivec=1:nbvec 995 iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box 996 jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% yindex on the image B for the middle of the correlation box 997 F(ivec)=0; 998 subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage 999 subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage 1000 subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage 1001 subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage 1002 image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A 1003 image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A 1004 check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1 1005 check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight; 1006 check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2 1007 check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight; 1008 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 1009 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 1010 image1_mean=mean(mean(image1_crop)); 1011 image2_mean=mean(mean(image2_crop)); 1012 %threshold on image minimum 1013 if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma) 1014 F(ivec)=3; 1015 end 1016 %threshold on image maximum 1017 if check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma) 1018 F(ivec)=3; 1019 end 1020 1021 if F(ivec)~=3 1022 image1_crop=image1_crop-image1_mean;%substract the mean 1023 image2_crop=image2_crop-image2_mean; 995 1024 if CheckDeformation 1025 xi=(1:mesh:size(image1_crop,2)); 1026 yi=(1:mesh:size(image1_crop,1))'; 996 1027 [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2)); 997 1028 XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); 998 1029 YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2); 999 1030 image1_crop=interp2(image1_crop,XIant,YIant); 1000 else 1001 image1_crop=interp2(image1_crop,xi,yi); 1002 end 1003 xi=(1:mesh:size(image2_crop,2)); 1004 yi=(1:mesh:size(image2_crop,1))'; 1005 image2_crop=interp2(image2_crop,xi,yi); 1006 end 1007 sum_square=sum(sum(image1_crop.*image1_crop)); 1008 %reference: Oliver Pust, PIV: Direct Cross-Correlation 1009 result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid'); 1010 corrmax= max(max(result_conv)); 1011 result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 1012 %Find the correlation max, at 255 1013 [y,x] = find(result_conv==255,1); 1014 if ~isempty(y) && ~isempty(x) 1015 try 1016 if par_civ.CorrSmooth==1 1017 [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y); 1018 elseif par_civ.CorrSmooth==2 1019 [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y); 1020 end 1021 % if ~isfield(par_civ,'CheckDeformation') 1022 utable(ivec)=vector(1)*mesh+shiftx(ivec); 1023 vtable(ivec)=vector(2)*mesh+shifty(ivec); 1024 % else 1025 % utable(ivec)=shiftx(ivec);% TEST TEST !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1026 % vtable(ivec)=shifty(ivec); 1027 % end 1028 xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2) 1029 ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge) 1030 iref=round(xtable(ivec));% image index for the middle of the vector 1031 jref=round(ytable(ivec)); 1032 if checkmask && par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100 1033 utable(ivec)=0; 1034 vtable(ivec)=0; 1031 xi=(1:mesh:size(image2_crop,2)); 1032 yi=(1:mesh:size(image2_crop,1))'; 1033 image2_crop=interp2(image2_crop,xi,yi); 1034 end 1035 sum_square=sum(sum(image1_crop.*image1_crop)); 1036 %reference: Oliver Pust, PIV: Direct Cross-Correlation 1037 result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid'); 1038 corrmax= max(max(result_conv)); 1039 result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 1040 %Find the correlation max, at 255 1041 [y,x] = find(result_conv==255,1); 1042 if ~isempty(y) && ~isempty(x) 1043 try 1044 if par_civ.CorrSmooth==1 1045 [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y); 1046 elseif par_civ.CorrSmooth==2 1047 [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y); 1048 end 1049 utable(ivec)=vector(1)*mesh+shiftx(ivec); 1050 vtable(ivec)=vector(2)*mesh+shifty(ivec); 1051 xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2) 1052 ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge) 1053 iref=round(xtable(ivec));% image index for the middle of the vector 1054 jref=round(ytable(ivec)); 1055 if checkmask && par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100 1056 utable(ivec)=0; 1057 vtable(ivec)=0; 1058 F(ivec)=3; 1059 end 1060 ctable(ivec)=corrmax/sum_square;% correlation value 1061 catch ME 1035 1062 F(ivec)=3; 1036 1063 end 1037 ctable(ivec)=corrmax/sum_square;% correlation value 1038 catch ME 1064 else 1039 1065 F(ivec)=3; 1040 1066 end 1041 else1042 F(ivec)=3;1043 1067 end 1044 1068 end
Note: See TracChangeset
for help on using the changeset viewer.