Changeset 382 for trunk/src/civ_matlab.m
 Timestamp:
 Feb 6, 2012, 11:46:39 PM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/civ_matlab.m
r372 r382 3 3 % civ: PIV function itself 4 4 % fix: removes false vectors after detection by various criteria 5 % patch: make interpolationsmoothing5 % filter_tps: make interpolationsmoothing 6 6 % 7 7 % function [Data,errormsg,result_conv]= civ_uvmat(Param,ncfile) … … 169 169 Data.Patch1_Threshold=Param.Patch1.MaxDiff; 170 170 Data.Patch1_SubDomain=Param.Patch1.SubdomainSize; 171 Data.ListVarName=[Data.ListVarName {'Civ1_U_Diff','Civ1_V_Diff','Civ1_ X_SubRange','Civ1_Y_SubRange','Civ1_NbSites','Civ1_X_tps','Civ1_Y_tps','Civ1_U_tps','Civ1_V_tps'}];172 Data.VarDimName=[Data.VarDimName {'NbVec1','NbVec1',{'Nb SubDomain1','Two'},{'NbSubDomain1','Two'},'NbSubDomain1',...173 {'NbVec1Sub','Nb SubDomain1'},{'NbVec1Sub','NbSubDomain1'},{'Nbtps1','NbSubDomain1'},{'Nbtps1','NbSubDomain1'}}];171 Data.ListVarName=[Data.ListVarName {'Civ1_U_Diff','Civ1_V_Diff','Civ1_SubRange','Civ1_NbSites','Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps'}]; 172 Data.VarDimName=[Data.VarDimName {'NbVec1','NbVec1',{'NbCoord','Two','NbSubDomain1'},{'NbSubDomain1'},... 173 {'NbVec1Sub','NbCoord','NbSubDomain1'},{'Nbtps1','NbSubDomain1'},{'Nbtps1','NbSubDomain1'}}]; 174 174 nbvar=length(Data.ListVarName); 175 175 Data.VarAttribute{nbvar1}.Role='vector_x'; … … 182 182 ind_good=1:numel(Data.Civ1_X); 183 183 end 184 [Data.Civ1_X_SubRange,Data.Civ1_Y_SubRange,Data.Civ1_NbSites,FFres,Ures, Vres,Data.Civ1_X_tps,Data.Civ1_Y_tps,Data.Civ1_U_tps,Data.Civ1_V_tps]=... 185 patch(Data.Civ1_X(ind_good)',Data.Civ1_Y(ind_good)',Data.Civ1_U(ind_good)',Data.Civ1_V(ind_good)',Data.Patch1_Rho,Data.Patch1_Threshold,Data.Patch1_SubDomain); 184 % [SubRange,NbPoints,Coord_tps,U_tps,V_tps,W_tps,U_smooth,V_smooth,W_smooth,FF] 185 [Data.Civ1_SubRange,Data.Civ1_NbSites,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=... 186 filter_tps([Data.Civ1_X(ind_good) Data.Civ1_Y(ind_good)],Data.Civ1_U(ind_good),Data.Civ1_V(ind_good),[],Data.Patch1_SubDomain,Data.Patch1_Threshold,Data.Patch1_Rho); 186 187 Data.Civ1_U_Diff(ind_good)=Data.Civ1_U(ind_good)Ures; 187 188 Data.Civ1_V_Diff(ind_good)=Data.Civ1_V(ind_good)Vres; … … 335 336 Data.Patch2_Threshold=Param.Patch2.MaxDiff; 336 337 Data.Patch2_SubDomain=Param.Patch2.SubdomainSize; 337 Data.ListVarName=[Data.ListVarName {'Civ2_U_Diff','Civ2_V_Diff','Civ2_X_SubRange','Civ2_Y_SubRange','Civ2_ NbSites','Civ2_X_tps','Civ2_Y_tps','Civ2_U_tps','Civ2_V_tps'}];338 Data.VarDimName=[Data.VarDimName {'NbVec2','NbVec2',{'NbSubDomain2','Two'},{'NbSubDomain2','Two'}, 'NbSubDomain2',...339 {'NbVec2Sub','NbSubDomain2'},{'NbVec2Sub','NbSubDomain2'},{'Nbtps2','NbSubDomain2'},{'Nbtps2','NbSubDomain2'} }];338 Data.ListVarName=[Data.ListVarName {'Civ2_U_Diff','Civ2_V_Diff','Civ2_X_SubRange','Civ2_Y_SubRange','Civ2_X_tps','Civ2_Y_tps','Civ2_U_tps','Civ2_V_tps','Civ2_Indices_tps'}]; 339 Data.VarDimName=[Data.VarDimName {'NbVec2','NbVec2',{'NbSubDomain2','Two'},{'NbSubDomain2','Two'},... 340 {'NbVec2Sub','NbSubDomain2'},{'NbVec2Sub','NbSubDomain2'},{'Nbtps2','NbSubDomain2'},{'Nbtps2','NbSubDomain2'},{'NbVec2Sub','NbSubDomain2'}}]; 340 341 nbvar=length(Data.ListVarName); 341 342 Data.VarAttribute{nbvar1}.Role='vector_x'; … … 348 349 ind_good=1:numel(Data.Civ2_X); 349 350 end 350 [Data.Civ2_ X_SubRange,Data.Civ2_Y_SubRange,Data.Civ2_NbSites,FFres,Ures, Vres,Data.Civ2_X_tps,Data.Civ2_Y_tps,Data.Civ2_U_tps,Data.Civ2_V_tps]=...351 patch(Data.Civ2_X(ind_good)',Data.Civ2_Y(ind_good)',Data.Civ2_U(ind_good)',Data.Civ2_V(ind_good)',Data.Patch2_Rho,Data.Patch2_Threshold,Data.Patch2_SubDomain);351 [Data.Civ2_SubRange,Data.Civ2_NbSites,Data.Civ2_Coord_tps,Data.Civ2_U_tps,Data.Civ2_V_tps,tild,Ures, Vres,tild,FFres]=... 352 filter_tps([Data.Civ2_X(ind_good) Data.Civ2_Y(ind_good)],Data.Civ2_U(ind_good),Data.Civ2_V(ind_good),Data.Patch2_SubDomain,Data.Patch2_Rho,Data.Patch2_Threshold); 352 353 Data.Civ2_U_Diff(ind_good)=Data.Civ2_U(ind_good)Ures; 353 354 Data.Civ2_V_Diff(ind_good)=Data.Civ2_V(ind_good)Vres; … … 697 698 698 699 699 % 700 % patch function 701 % OUTPUT: 702 % SubRangx,SubRangy(NbSubdomain,2): range (min, max) of the coordiantes x and y respectively, for each subdomain 703 % nbpoints(NbSubdomain): number of source points for each subdomain 704 % FF: false flags 705 % U_smooth, V_smooth: filtered velocity components at the positions of the initial data 706 % X_tps,Y_tps,U_tps,V_tps: positions and weight of the tps for each subdomain 707 % 708 % INPUT: 709 % X, Y: set of coordinates of the initial data 710 % U,V: set of velocity components of the initial data 711 % Rho: smoothing parameter 712 % Threshold: max diff accepted between smoothed and initial data 713 % Subdomain: estimated number of data points in each subdomain 714 715 function [SubRangx,SubRangy,nbpoints,FF,U_smooth,V_smooth,X_tps,Y_tps,U_tps,V_tps] =patch(X,Y,U,V,Rho,Threshold,SubDomain) 716 %subdomain decomposition 717 warning off 718 U=reshape(U,[],1); 719 V=reshape(V,[],1); 720 X=reshape(X,[],1); 721 Y=reshape(Y,[],1); 722 nbvec=numel(X); 723 NbSubDomain=ceil(nbvec/SubDomain); 724 MinX=min(X); 725 MinY=min(Y); 726 MaxX=max(X); 727 MaxY=max(Y); 728 RangX=MaxXMinX; 729 RangY=MaxYMinY; 730 AspectRatio=RangY/RangX; 731 NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1); 732 NbSubDomainY=max(floor(sqrt(NbSubDomain*AspectRatio)),1); 733 NbSubDomain=NbSubDomainX*NbSubDomainY; 734 SizX=RangX/NbSubDomainX;%width of subdomains 735 SizY=RangY/NbSubDomainY;%height of subdomains 736 CentreX=linspace(MinX+SizX/2,MaxXSizX/2,NbSubDomainX); 737 CentreY=linspace(MinY+SizY/2,MaxYSizY/2,NbSubDomainY); 738 [CentreX,CentreY]=meshgrid(CentreX,CentreY); 739 CentreY=reshape(CentreY,1,[]);% Y positions of subdomain centres 740 CentreX=reshape(CentreX,1,[]);% X positions of subdomain centres 741 rho=SizX*SizY*Rho/1000000;%optimum rho increase as the area of the subdomain (division by 10^6 to reach good values with the default GUI input) 742 U_tps_sub=zeros(length(X),NbSubDomain);%default spline 743 V_tps_sub=zeros(length(X),NbSubDomain);%default spline 744 U_smooth=zeros(length(X),1); 745 V_smooth=zeros(length(X),1); 746 747 nb_select=zeros(length(X),1); 748 FF=zeros(length(X),1); 749 check_empty=zeros(1,NbSubDomain); 750 SubRangx=zeros(NbSubDomain,2);%initialise the positions of subdomains 751 SubRangy=zeros(NbSubDomain,2); 752 for isub=1:NbSubDomain 753 SubRangx(isub,:)=[CentreX(isub)0.55*SizX CentreX(isub)+0.55*SizX]; 754 SubRangy(isub,:)=[CentreY(isub)0.55*SizY CentreY(isub)+0.55*SizY]; 755 ind_sel_previous=[]; 756 ind_sel=0; 757 while numel(ind_sel)>numel(ind_sel_previous) %increase the subdomain during four iterations at most 758 ind_sel_previous=ind_sel; 759 ind_sel=find(X>=SubRangx(isub,1) & X<=SubRangx(isub,2) & Y>=SubRangy(isub,1) & Y<=SubRangy(isub,2)); 760 % if no vector in the subdomain, skip the subdomain 761 if isempty(ind_sel) 762 check_empty(isub)=1; 763 U_tps(1,isub)=0;%define U_tps and V_tps by default 764 V_tps(1,isub)=0; 765 break 766 % if too few selected vectors, increase the subrange for next iteration 767 elseif numel(ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous); 768 SubRangx(isub,1)=SubRangx(isub,1)SizX/4; 769 SubRangx(isub,2)=SubRangx(isub,2)+SizX/4; 770 SubRangy(isub,1)=SubRangy(isub,1)SizY/4; 771 SubRangy(isub,2)=SubRangy(isub,2)+SizY/4; 772 else 773 774 [U_smooth_sub,U_tps_sub]=tps_coeff([X(ind_sel) Y(ind_sel)],U(ind_sel),rho); 775 [V_smooth_sub,V_tps_sub]=tps_coeff([X(ind_sel) Y(ind_sel)],V(ind_sel),rho); 776 UDiff=U_smooth_subU(ind_sel); 777 VDiff=V_smooth_subV(ind_sel); 778 NormDiff=UDiff.*UDiff+VDiff.*VDiff; 779 FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination 780 ind_ind_sel=find(FF(ind_sel)==0); % select the indices of ind_sel corresponding to the remaining vectors 781 % no value exceeds threshold, the result is recorded 782 if isequal(numel(ind_ind_sel),numel(ind_sel)) 783 U_smooth(ind_sel)=U_smooth(ind_sel)+U_smooth_sub; 784 V_smooth(ind_sel)=V_smooth(ind_sel)+V_smooth_sub; 785 nbpoints(isub)=numel(ind_sel); 786 X_tps(1:nbpoints(isub),isub)=X(ind_sel); 787 Y_tps(1:nbpoints(isub),isub)=Y(ind_sel); 788 U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub; 789 V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub; 790 nb_select(ind_sel)=nb_select(ind_sel)+1; 791 display('good') 792 break 793 % too few selected vectors, increase the subrange for next iteration 794 elseif numel(ind_ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous); 795 SubRangx(isub,1)=SubRangx(isub,1)SizX/4; 796 SubRangx(isub,2)=SubRangx(isub,2)+SizX/4; 797 SubRangy(isub,1)=SubRangy(isub,1)SizY/4; 798 SubRangy(isub,2)=SubRangy(isub,2)+SizY/4; 799 % display('fewsmooth') 800 % interpolationsmoothing is done again with the selected vectors 801 else 802 [U_smooth_sub,U_tps_sub]=tps_coeff([X(ind_sel(ind_ind_sel)) Y(ind_sel(ind_ind_sel))],U(ind_sel(ind_ind_sel)),rho); 803 [V_smooth_sub,V_tps_sub]=tps_coeff([X(ind_sel(ind_ind_sel)) Y(ind_sel(ind_ind_sel))],V(ind_sel(ind_ind_sel)),rho); 804 U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+U_smooth_sub; 805 V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+V_smooth_sub; 806 nbpoints(isub)=numel(ind_ind_sel); 807 X_tps(1:nbpoints(isub),isub)=X(ind_sel(ind_ind_sel)); 808 Y_tps(1:nbpoints(isub),isub)=Y(ind_sel(ind_ind_sel)); 809 U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub; 810 V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub; 811 nb_select(ind_sel(ind_ind_sel))=nb_select(ind_sel(ind_ind_sel))+1; 812 display('good2') 813 break 814 end 815 end 816 end 817 end 818 ind_empty=find(check_empty); 819 %remove empty subdomains 820 if ~isempty(ind_empty) 821 SubRangx(ind_empty,:)=[]; 822 SubRangy(ind_empty,:)=[]; 823 X_tps(:,ind_empty)=[]; 824 Y_tps(:,ind_empty)=[]; 825 U_tps(:,ind_empty)=[]; 826 V_tps(:,ind_empty)=[]; 827 end 828 nb_select(nb_select==0)=1;%ones(size(find(nb_select==0))); 829 U_smooth=U_smooth./nb_select; 830 V_smooth=V_smooth./nb_select; 831 832 833 834 835 700 701 702 703
Note: See TracChangeset
for help on using the changeset viewer.