Changeset 326 for trunk/src/civ_matlab.m
- Timestamp:
- Dec 8, 2011, 8:28:38 AM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/civ_matlab.m
r321 r326 68 68 % if exist('ncfile','var')% TEST for use interactively with mouse_motion (no file created) 69 69 Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'};% cell array containing the names of the fields to record 70 Data.VarDimName={' nbvec1','nbvec1','nbvec1','nbvec1','nbvec1','nbvec1'};70 Data.VarDimName={'NbVec1','NbVec1','NbVec1','NbVec1','NbVec1','NbVec1'}; 71 71 Data.VarAttribute{1}.Role='coord_x'; 72 72 Data.VarAttribute{2}.Role='coord_y'; … … 124 124 else 125 125 Data.ListVarName=[Data.ListVarName {'Civ1_FF'}]; 126 Data.VarDimName=[Data.VarDimName {' nbvec1'}];126 Data.VarDimName=[Data.VarDimName {'NbVec1'}]; 127 127 nbvar=length(Data.ListVarName); 128 128 Data.VarAttribute{nbvar}.Role='errorflag'; … … 133 133 %% Patch1 134 134 if isfield (Param,'Patch1') 135 if check_civx 136 errormsg='Civ Matlab input needed for patch'; 137 return 138 end 135 139 check_patch1=1; 136 140 Param.Patch1 … … 233 237 Data.Civ1_Dt=str2double(par_civ1.Dt); 234 238 Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'};% cell array containing the names of the fields to record 235 Data.VarDimName={' nbvec1','nbvec1','nbvec1','nbvec1','nbvec1','nbvec1'};239 Data.VarDimName={'NbVec1','NbVec1','NbVec1','NbVec1','NbVec1','NbVec1'}; 236 240 Data.VarAttribute{1}.Role='coord_x'; 237 241 Data.VarAttribute{2}.Role='coord_y'; … … 596 600 if isfield(Param,FlagName{iflag}) && Param.(FlagName{iflag}) 597 601 FF=(FF==1| F==FlagVal(iflag)); 598 % if isfield (Param,'WarnFlags')599 % for iflag=1:numel(Param.WarnFlags)600 % FF=(FF==1| F==Param.WarnFlags(iflag));601 % end602 % end603 602 end 604 603 end … … 618 617 return 619 618 620 %621 % if isfield (Param,'LowerBoundVel')&& ~isequal(Param.LowerBoundVel,0)622 % thresh=Param.LowerBoundVel*Param.LowerBoundVel;623 % FF=FF==1 | (U.*U+V.*V)<thresh;624 % end625 % if isfield (Param,'UpperBoundVel')&& ~isequal(Param.UpperBoundVel,0)626 % thresh=Param.UpperBoundVel*Param.UpperBoundVel;627 % FF=FF==1 | (U.*U+V.*V)>thresh;628 % end629 % if isfield(Param,'MaskName')630 % M=imread(Param.MaskName);631 % nxy=size(M);632 % M=reshape(M,1,[]);633 % rangx0=[0.5 nxy(2)-0.5];634 % rangy0=[0.5 nxy(1)-0.5];635 % vec_x1=X-U/2;%beginning points636 % vec_x2=X+U/2;%end points of vectors637 % vec_y1=Y-V/2;%beginning points638 % vec_y2=Y+V/2;%end points of vectors639 % indx=1+round((nxy(2)-1)*(vec_x1-rangx0(1))/(rangx0(2)-rangx0(1)));% image index x at abcissa vec_x1640 % indy=1+round((nxy(1)-1)*(vec_y1-rangy0(1))/(rangy0(2)-rangy0(1)));% image index y at ordinate vec_y1641 % check_in=~(indx < 1 |indy < 1 | indx > nxy(2) |indy > nxy(1)); %=0 out of the mask image, 1 inside642 % indx=indx.*check_in+(1-check_in); %replace indx by 1 out of the mask range643 % indy=indy.*check_in+(1-check_in); %replace indy by 1 out of the mask range644 % ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));%determine the indices in the image reshaped in a Matlab vector645 % Mvalues=M(ICOMB);646 % flag7b=((20 < Mvalues) & (Mvalues < 200))| ~check_in';647 % indx=1+round((nxy(2)-1)*(vec_x2-rangx0(1))/(rangx0(2)-rangx0(1)));% image index x at abcissa vec_x2648 % indy=1+round((nxy(1)-1)*(vec_y2-rangy0(1))/(rangy0(2)-rangy0(1)));% image index y at ordinate vec_y2649 % check_in=~(indx < 1 |indy < 1 | indx > nxy(2) |indy > nxy(1)); %=0 out of the mask image, 1 inside650 % indx=indx.*check_in+(1-check_in); %replace indx by 1 out of the mask range651 % indy=indy.*check_in+(1-check_in); %replace indy by 1 out of the mask range652 % ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));%determine the indices in the image reshaped in a Matlab vector653 % Mvalues=M(ICOMB);654 % flag7e=((Mvalues > 20) & (Mvalues < 200))| ~check_in';655 % FF=FF==1 |(flag7b|flag7e)';656 % end657 % % flag7=0;658 % % end659 %660 619 661 620 FF=double(FF); 662 %663 % % criterium on velocity values664 % delta_u=Field.U;%default without ref file665 % delta_v=Field.V;666 % if exist('fileref','var') && ~isempty(fileref)667 % if ~exist(fileref,'file')668 % error='reference file not found in RUN_FIX.m';669 % display(error);670 % return671 % end672 % FieldRef=read_civxdata(fileref,[],fieldref);673 % if isfield(FieldRef,'FF')674 % index_true=find(FieldRef.FF==0);675 % FieldRef.X=FieldRef.X(index_true);676 % FieldRef.Y=FieldRef.Y(index_true);677 % FieldRef.U=FieldRef.U(index_true);678 % FieldRef.V=FieldRef.V(index_true);679 % end680 % if ~isfield(FieldRef,'X') || ~isfield(FieldRef,'Y') || ~isfield(FieldRef,'U') || ~isfield(FieldRef,'V')681 % error='reference file is not a velocity field in RUN_FIX.m '; %bad input file682 % return683 % end684 % if length(FieldRef.X)<=1685 % errordlg('reference field with one vector or less in RUN_FIX.m')686 % return687 % end688 % vec_U_ref=griddata_uvmat(FieldRef.X,FieldRef.Y,FieldRef.U,Field.X,Field.Y); %interpolate vectors in the ref field689 % vec_V_ref=griddata_uvmat(FieldRef.X,FieldRef.Y,FieldRef.V,Field.X,Field.Y); %interpolate vectors in the ref field to the positions of the main field690 % delta_u=Field.U-vec_U_ref;%take the difference with the interpolated ref field691 % delta_v=Field.V-vec_V_ref;692 % end693 % thresh_vel_x=thresh_vel;694 % thresh_vel_y=thresh_vel;695 % if isequal(inf_sup,1)696 % flag5=abs(delta_u)<thresh_vel_x & abs(delta_v)<thresh_vel_y &(flag1~=1)&(flag2~=1)&(flag3~=1)&(flag4~=1);697 % elseif isequal(inf_sup,2)698 % flag5=(abs(delta_u)>thresh_vel_x | abs(delta_v)>thresh_vel_y) &(flag1~=1)&(flag2~=1)&(flag3~=1)&(flag4~=1);699 % end700 %701 % % flag7 introduce a grey mask, matrix M702 703 % flagmagenta=flag1|flag2|flag3|flag4|flag5|flag7;704 % fixflag_unit=Field.FF-10*floor(Field.FF/10); %unity term of fix_flag705 621 706 622 … … 708 624 %------------------------------------------------------------------------ 709 625 % patch function 626 % OUTPUT: 627 % SubRangx,SubRangy(NbSubdomain,2): range (min, max) of the coordiantes x and y respectively, for each subdomain 628 % nbpoints(NbSubdomain): number of source points for each subdomain 629 % FF: false flags 630 % U_smooth, V_smooth: filtered velocity components at the positions of the initial data 631 % X_tps,Y_tps,U_tps,V_tps: positions and weight of the tps for each subdomain 632 % 633 % INPUT: 634 % X, Y: set of coordinates of the initial data 635 % U,V: set of velocity components of the initial data 636 % Rho: smoothing parameter 637 % Threshold: max diff accepted between smoothed and initial data 638 % Subdomain: estimated number of data points in each subdomain 639 710 640 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) 711 641 %subdomain decomposition … … 744 674 check_empty=zeros(1,NbSubDomain); 745 675 for isub=1:NbSubDomain 746 SubRangx(isub,:)=[CentreX(isub)- SizX/2 CentreX(isub)+SizX/2];747 SubRangy(isub,:)=[CentreY(isub)- SizY/2 CentreY(isub)+SizY/2];676 SubRangx(isub,:)=[CentreX(isub)-0.55*SizX CentreX(isub)+0.55*SizX]; 677 SubRangy(isub,:)=[CentreY(isub)-0.55*SizY CentreY(isub)+0.55*SizY]; 748 678 ind_sel_previous=[]; 749 679 ind_sel=0; 750 680 while numel(ind_sel)>numel(ind_sel_previous) %increase the subdomain during four iterations at most 751 681 ind_sel_previous=ind_sel; 752 ind_sel=find(X> SubRangx(isub,1) & X<SubRangx(isub,2) & Y>SubRangy(isub,1) & Y<SubRangy(isub,2));682 ind_sel=find(X>=SubRangx(isub,1) & X<=SubRangx(isub,2) & Y>=SubRangy(isub,1) & Y<=SubRangy(isub,2)); 753 683 % if no vector in the subdomain, skip the subdomain 754 684 if isempty(ind_sel)
Note: See TracChangeset
for help on using the changeset viewer.