Changeset 246 for trunk/src/civ_uvmat.m
- Timestamp:
- Apr 28, 2011, 10:52:31 AM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/civ_uvmat.m
r243 r246 1 1 % To develop.... 2 2 function [Data,errormsg]= civ_uvmat(Param,ncfile) 3 errormsg=''; 3 4 Data.ListGlobalAttribute={'Conventions','Program','CivStage'}; 4 5 Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes 5 6 Data.Program='civ_uvmat'; 6 7 Data.CivStage=0;%default 8 ListVarCiv1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'}; 9 ListVarFix1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF'}; 7 10 8 11 %% Civ1 … … 20 23 shiftx=str2num(par_civ1.shiftx); 21 24 shifty=str2num(par_civ1.shifty); 22 miniy=max(1+isy2 -shifty,1+iby2);25 miniy=max(1+isy2+shifty,1+iby2); 23 26 minix=max(1+isx2-shiftx,1+ibx2); 24 maxiy=min(size(image1,1)-isy2 -shifty,size(image1,1)-iby2);27 maxiy=min(size(image1,1)-isy2+shifty,size(image1,1)-iby2); 25 28 maxix=min(size(image1,2)-isx2-shiftx,size(image1,2)-ibx2); 26 29 [GridX,GridY]=meshgrid(minix:stepx:maxix,miniy:stepy:maxiy); … … 29 32 30 33 % caluclate velocity data (y and v in indices, reverse to y component) 31 [xtable ytable utable vtable ctable F] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx, shifty,PointCoord,str2num(par_civ1.rho), []);34 [xtable ytable utable vtable ctable F] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,-shifty,PointCoord,str2num(par_civ1.rho), []); 32 35 list_param=(fieldnames(par_civ1))'; 33 36 list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'}; … … 67 70 Data.CivStage=1; 68 71 else 69 Data=nc2struct(ncfile)%read existing netcdf file 72 if isfield(Param,'Fix1') 73 Data=nc2struct(ncfile,ListVarCiv1);%read civ1 data in the existing netcdf file 74 else 75 Data=nc2struct(ncfile,ListVarFix1);%read civ1 and fix1 data in the existing netcdf file 76 end 77 if isfield(Data,'Txt') 78 msgbox_uvmat('ERROR',Data.Txt) 79 return 80 end 81 % read Civx data 70 82 if isfield(Data,'absolut_time_T0')%read civx file 71 83 var={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF';'vec_X','vec_Y','vec_U','vec_V','vec_C','vec_F','vec_FixFlag'}; … … 95 107 ParamName=ListFixParam{ilist}; 96 108 ListName=['Fix1_' ParamName]; 97 ['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];']98 109 eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];']) 99 110 eval(['Data.' ListName '=Param.Fix1.' ParamName ';']) … … 107 118 Data.VarDimName=[Data.VarDimName {'nbvec1'}]; 108 119 nbvar=length(Data.ListVarName); 109 Data.VarAttribute{nbvar}.Role='errorflag'; 110 [Data.Civ1_FF]=fix_uvmat(Param.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V);120 Data.VarAttribute{nbvar}.Role='errorflag'; 121 Data.Civ1_FF=fix_uvmat(Param.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V); 111 122 Data.CivStage=2; 112 123 end … … 117 128 Data.Patch1_Threshold=str2double(Param.Patch1.Threshold); 118 129 Data.Patch1_SubDomain=str2double(Param.Patch1.SubDomain); 119 Data.ListVarName=[Data.ListVarName {'Patch1_U','Patch1_V'}]; 120 Data.VarDimName=[Data.VarDimName {'nbvec1','nbvec1'}]; 130 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'}]; 131 Data.VarDimName=[Data.VarDimName {'NbVec1','NbVec1',{'NbSubDomain1','Two'},{'NbSubDomain1','Two'},'NbSubDomain1',... 132 {'NbVec1Sub','NbSubDomain1'},{'NbVec1Sub','NbSubDomain1'},{'Nbtps1','NbSubDomain1'},{'Nbtps1','NbSubDomain1'}}]; 121 133 nbvar=length(Data.ListVarName); 122 134 Data.VarAttribute{nbvar-1}.Role='vector_x'; 123 135 Data.VarAttribute{nbvar}.Role='vector_y'; 124 Data. Patch1_U=zeros(size(Data.Civ1_X));125 Data. Patch1_V=zeros(size(Data.Civ1_X));136 Data.Civ1_U_Diff=zeros(size(Data.Civ1_X)); 137 Data.Civ1_V_Diff=zeros(size(Data.Civ1_X)); 126 138 if isfield(Data,'Civ1_FF') 127 139 ind_good=find(Data.Civ1_FF==0); … … 129 141 ind_good=1:numel(Data.Civ1_X); 130 142 end 131 Data.Civ1_X 132 [Ures, Vres,FFres]=... 143 [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]=... 133 144 patch_uvmat(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); 134 size(Ures) 135 size(Vres) 136 size(FFres) 137 size(ind_good) 138 Data.Patch1_U(ind_good)=Ures; 139 Data.Patch1_V(ind_good)=Vres; 140 Data.Civ1_FF(ind_good)=FFres 141 Data.CivStage=3; 145 Data.Civ1_U_Diff(ind_good)=Data.Civ1_U(ind_good)-Ures; 146 Data.Civ1_V_Diff(ind_good)=Data.Civ1_V(ind_good)-Vres; 147 Data.Civ1_FF(ind_good)=FFres; 148 Data.CivStage=3; 142 149 end 143 150 %% write result 151 % 'TESTcalc' 152 % [DataOut,errormsg]=calc_field('velocity',Data) 153 if exist('ncfile','var') 144 154 errormsg=struct2nc(ncfile,Data); 145 155 end 146 156 147 157 … … 192 202 FF=FF==1 | (U.*U+V.*V)>thresh; 193 203 end 194 204 FF=double(FF); 195 205 % 196 206 % % criterium on velocity values … … 270 280 %------------------------------------------------------------------------ 271 281 % patch function 272 function [ U_patch,V_patch,FF,SubRangx,SubRangy,X_ctrs,Y_ctrs] =patch_uvmat(X,Y,U,V,Rho,Threshold,SubDomain)282 function [SubRangx,SubRangy,nbpoints,FF,U_smooth,V_smooth,X_tps,Y_tps,U_tps,V_tps] =patch_uvmat(X,Y,U,V,Rho,Threshold,SubDomain) 273 283 %subdomain decomposition 274 284 warning off … … 278 288 Y=reshape(Y,[],1); 279 289 nbvec=numel(X); 280 NbSubDomain=ceil(nbvec/SubDomain) 281 MinX=min(X) 282 MinY=min(Y) 283 MaxX=max(X) 284 MaxY=max(Y) 290 NbSubDomain=ceil(nbvec/SubDomain); 291 MinX=min(X); 292 MinY=min(Y); 293 MaxX=max(X); 294 MaxY=max(Y); 285 295 RangX=MaxX-MinX; 286 296 RangY=MaxY-MinY; 287 AspectRatio=RangY/RangX 288 NbSubDomainX= ceil(sqrt(NbSubDomain/AspectRatio))289 NbSubDomainY= ceil(sqrt(NbSubDomain*AspectRatio))290 297 AspectRatio=RangY/RangX; 298 NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1); 299 NbSubDomainY=max(floor(sqrt(NbSubDomain*AspectRatio)),1); 300 NbSubDomain=NbSubDomainX*NbSubDomainY; 291 301 SizX=RangX/NbSubDomainX;%width of subdomains 292 302 SizY=RangY/NbSubDomainY;%height of subdomains 293 303 CentreX=linspace(MinX+SizX/2,MaxX-SizX/2,NbSubDomainX); 294 304 CentreY=linspace(MinY+SizY/2,MaxY-SizY/2,NbSubDomainY); 295 SubIndexX=ceil((X-MinX)/SizX);%subdomain index of vectors 296 SubIndexY=ceil((Y-MinY)/SizY);297 rho=RangX*RangY*Rho;%optimum rho increase as teh area of the subdomain 298 U_tps=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline 299 V_tps=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline300 X_dist=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline301 Y_dist=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline 302 U_smooth=zeros(length(X),NbSubDomainY,NbSubDomainX);303 V_smooth=zeros(length(X),NbSubDomainY,NbSubDomainX); 304 dist_ctre=zeros(length(X),NbSubDomainY,NbSubDomainX);305 [CentreX,CentreY]=meshgrid(CentreX,CentreY); 306 CentreY=reshape(CentreY,1,[]); 307 CentreX=reshape(CentreX,1,[]); 308 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) 309 U_tps_sub=zeros(length(X),NbSubDomain);%default spline 310 V_tps_sub=zeros(length(X),NbSubDomain);%default spline 311 U_smooth=zeros(length(X),1); 312 V_smooth=zeros(length(X),1); 313 314 nb_select=zeros(length(X),1); 305 315 FF=zeros(length(X),1); 306 for isubx=1:NbSubDomainX 307 for isuby=1:NbSubDomainY 308 SubRangx(isuby,isubx,:)=[CentreX(isubx)-SizX/2 CentreX(isubx)+SizX/2]; 309 SubRangy(isuby,isubx,:)=[CentreY(isubx)-SizY/2 CentreY(isubx)+SizY/2]; 310 for iter=1:3 %increase the subdomain during three iterations at most 311 ind_sel=find(X>SubRangx(isuby,isubx,1) & X<SubRangx(isuby,isubx,2) & Y>SubRangy(isuby,isubx,1) & Y<SubRangy(isuby,isubx,2)); 312 size(ind_sel) 313 if numel(ind_sel)<SubDomain/4;% too few selected vectors, increase the subrange for next iteration 314 SubRangx(isuby,isubx,1)=SubRangx(isuby,isubx,1)-SizX/4; 315 SubRangx(isuby,isubx,2)=SubRangx(isuby,isubx,2)+SizX/4; 316 SubRangy(isuby,isubx,1)=SubRangy(isuby,isubx,1)-SizY/4; 317 SubRangy(isuby,isubx,2)=SubRangy(isuby,isubx,2)+SizY/4; 316 test_empty=zeros(1,NbSubDomain); 317 for isub=1:NbSubDomain 318 SubRangx(isub,:)=[CentreX(isub)-SizX/2 CentreX(isub)+SizX/2]; 319 SubRangy(isub,:)=[CentreY(isub)-SizY/2 CentreY(isub)+SizY/2]; 320 ind_sel_previous=[]; 321 ind_sel=0; 322 while numel(ind_sel)>numel(ind_sel_previous) %increase the subdomain during four iterations at most 323 ind_sel_previous=ind_sel; 324 ind_sel=find(X>SubRangx(isub,1) & X<SubRangx(isub,2) & Y>SubRangy(isub,1) & Y<SubRangy(isub,2)); 325 % if no vector in the subdomain, skip the subdomain 326 if isempty(ind_sel) 327 test_empty(isub)=1; 328 U_tps(1,isub)=0;%define U_tps and V_tps by default 329 V_tps(1,isub)=0; 330 break 331 % if too few selected vectors, increase the subrange for next iteration 332 elseif numel(ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous); 333 SubRangx(isub,1)=SubRangx(isub,1)-SizX/4; 334 SubRangx(isub,2)=SubRangx(isub,2)+SizX/4; 335 SubRangy(isub,1)=SubRangy(isub,1)-SizY/4; 336 SubRangy(isub,2)=SubRangy(isub,2)+SizY/4; 337 else 338 [U_smooth_sub,U_tps_sub]=tps_coeff(X(ind_sel),Y(ind_sel),U(ind_sel),rho); 339 [V_smooth_sub,V_tps_sub]=tps_coeff(X(ind_sel),Y(ind_sel),V(ind_sel),rho); 340 UDiff=U_smooth_sub-U(ind_sel); 341 VDiff=V_smooth_sub-V(ind_sel); 342 NormDiff=UDiff.*UDiff+VDiff.*VDiff; 343 FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination 344 ind_ind_sel=find(FF(ind_sel)==0); % select the indices of ind_sel corresponding to the remaining vectors 345 % no value exceeds threshold, the result is recorded 346 if isequal(numel(ind_ind_sel),numel(ind_sel)) 347 U_smooth(ind_sel)=U_smooth(ind_sel)+U_smooth_sub; 348 V_smooth(ind_sel)=V_smooth(ind_sel)+V_smooth_sub; 349 nbpoints(isub)=numel(ind_sel); 350 X_tps(1:nbpoints(isub),isub)=X(ind_sel); 351 Y_tps(1:nbpoints(isub),isub)=Y(ind_sel); 352 U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub; 353 V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub; 354 nb_select(ind_sel)=nb_select(ind_sel)+1; 355 display('good') 356 break 357 % too few selected vectors, increase the subrange for next iteration 358 elseif numel(ind_ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous); 359 SubRangx(isub,1)=SubRangx(isub,1)-SizX/4; 360 SubRangx(isub,2)=SubRangx(isub,2)+SizX/4; 361 SubRangy(isub,1)=SubRangy(isub,1)-SizY/4; 362 SubRangy(isub,2)=SubRangy(isub,2)+SizY/4; 363 % display('fewsmooth') 364 % interpolation-smoothing is done again with the selected vectors 318 365 else 319 [U_smooth_sub,U_tps_sub]=tps_uvmat(X(ind_sel),Y(ind_sel),U(ind_sel),rho); 320 [V_smooth_sub,V_tps_sub]=tps_uvmat(X(ind_sel),Y(ind_sel),V(ind_sel),rho); 321 size(U_smooth_sub) 322 size(U(ind_sel)) 323 UDiff=U_smooth_sub-U(ind_sel); 324 VDiff=V_smooth_sub-V(ind_sel); 325 NormDiff=UDiff.*UDiff+VDiff.*VDiff; 326 FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination 327 ind_ind_sel=find(FF(ind_sel)==0); 328 if isequal(numel(ind_ind_sel),numel(ind_sel)) 329 U_smooth(ind_sel,isuby,isubx)=U_smooth_sub; 330 V_smooth(ind_sel,isuby,isubx)=V_smooth_sub; 331 break 332 elseif numel(ind_ind_sel)<SubDomain/4;% too few selected vectors, increase the subrange for next iteration 333 SubRangx(isuby,isubx,1)=SubRangx(isuby,isubx,1)-SizX/4; 334 SubRangx(isuby,isubx,2)=SubRangx(isuby,isubx,2)+SizX/4; 335 SubRangy(isuby,isubx,1)=SubRangy(isuby,isubx,1)-SizY/4; 336 SubRangy(isuby,isubx,2)=SubRangy(isuby,isubx,2)+SizY/4; 337 else 338 [U_smooth(ind_sel(ind_ind_sel),isuby,isubx),U_tps(ind_sel(ind_ind_sel),isuby,isubx),U_tps3(isuby,isubx)]=tps_uvmat(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),U(ind_sel(ind_ind_sel)),rho); 339 [V_smooth(ind_sel(ind_ind_sel),isuby,isubx),V_tps(ind_sel(ind_ind_sel),isuby,isubx),V_tps3(isuby,isubx)]=tps_uvmat(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),V(ind_sel(ind_ind_sel)),rho); 340 break 341 end 342 end 343 end 344 X_ctrs(isuby,isubx)=mean(X(ind_sel));%gravity centre of selected points 345 Y_ctrs(isuby,isubx)=mean(Y(ind_sel));%positions of tps sources for the subdomain i, 346 dist_ctre(ind_sel,isuby,isubx)=sqrt(abs(((X(ind_sel)-X_ctrs(isuby,isubx)).*(X(ind_sel)-X_ctrs(isuby,isubx)))+((Y(ind_sel)-Y_ctrs(isuby,isubx)).*(Y(ind_sel)-Y_ctrs(isuby,isubx))))); 347 end 348 end 349 U_patch=sum(sum(U_smooth.*dist_ctre,3),2)./sum(sum(dist_ctre,3),2); 350 V_patch=sum(sum(V_smooth.*dist_ctre,3),2)./sum(sum(dist_ctre,3),2); 351 352 %------------------------------------------------------------------------ 353 %fasshauer@iit.edu MATH 590 ? Chapter 19 32 354 % X,Y initial coordiantes 355 % XI vector, YI column vector for the grid of interpolation points 356 function [U_smooth,U_tps,U_tps3]=tps_uvmat(X,Y,U,rho) 357 %------------------------------------------------------------------------ 358 %rho smoothing parameter 359 ep = 1; 360 X=reshape(X,[],1); 361 Y=reshape(Y,[],1); 362 rhs = reshape(U,[],1); 363 % if exist('FF','var') 364 % test_false=isnan(rhs)|FF~=0; 365 % else 366 % test_false=isnan(rhs); 367 % end 368 % X(test_false)=[]; 369 % Y(test_false)=[]; 370 % rhs(test_false)=[]; 371 %randn('state',3); rhs = rhs + 0.03*randn(size(rhs)); 372 rhs = [rhs; zeros(3,1)]; 373 dsites = [X Y];% coordinates of measurement sites 374 ctrs = dsites;%radial base functions are located at the measurement sites 375 DM_data = DistanceMatrix(dsites,ctrs);%2D matrix of distances between spline centres (=initial points) ctrs 376 % if size(XI,1)==1 && size(YI,2)==1 % XI vector, YI column vector 377 % [XI,YI]=meshgrid(XI,YI); 378 % end 379 % [npy,npx]=size(XI); 380 % epoints = [reshape(XI,[],1) reshape(YI,[],1)]; 381 IM_sites = tps(ep,DM_data);%values of thin plate at site points 382 IM = IM_sites + rho*eye(size(IM_sites));% rho=1/(2*omega) , omega given by fasshauer; 383 PM=[ones(size(dsites,1),1) dsites]; 384 IM=[IM PM; [PM' zeros(3,3)]]; 385 %fprintf('Condition number estimate: %e\n',condest(IM)) 386 %DM_eval = DistanceMatrix(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs 387 %EM = tps(ep,DM_eval);%values of thin plate 388 %PM = [ones(size(epoints,1),1) epoints]; 389 %EM = [EM PM]; 390 U_tps=(IM\rhs); 391 PM = [ones(size(dsites,1),1) dsites]; 392 EM = [IM_sites PM]; 393 U_smooth=EM *U_tps; 394 U_tps3=U_tps(end-2:end); 395 U_tps=U_tps(1:end-3); 366 [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); 367 [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); 368 U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+U_smooth_sub; 369 V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+V_smooth_sub; 370 nbpoints(isub)=numel(ind_ind_sel); 371 X_tps(1:nbpoints(isub),isub)=X(ind_sel(ind_ind_sel)); 372 Y_tps(1:nbpoints(isub),isub)=Y(ind_sel(ind_ind_sel)); 373 U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub; 374 V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub; 375 nb_select(ind_sel(ind_ind_sel))=nb_select(ind_sel(ind_ind_sel))+1; 376 display('good2') 377 break 378 end 379 end 380 end 381 end 382 ind_empty=find(test_empty); 383 %remove empty subdomains 384 if ~isempty(ind_empty) 385 SubRangx(ind_empty,:)=[]; 386 SubRangy(ind_empty,:)=[]; 387 X_tps(:,ind_empty)=[]; 388 Y_tps(:,ind_empty)=[]; 389 U_tps(:,ind_empty)=[]; 390 V_tps(:,ind_empty)=[]; 391 end 392 nb_select(nb_select==0)=1;%ones(size(find(nb_select==0))); 393 U_smooth=U_smooth./nb_select; 394 V_smooth=V_smooth./nb_select; 395 396 397 398 396 399 397 400 … … 421 424 422 425 423 % DM: MxN matrix whose i,j position contains the Euclidean 424 % distance between the i-th data site and j-th center 425 function DM = DistanceMatrix(dsites,ctrs) 426 [M,s] = size(dsites); [N,s] = size(ctrs); 427 DM = zeros(M,N); 428 % Accumulate sum of squares of coordinate differences 429 % The ndgrid command produces two MxN matrices: 430 % dr, consisting of N identical columns (each containing 431 % the d-th coordinate of the M data sites) 432 % cc, consisting of M identical rows (each containing 433 % the d-th coordinate of the N centers) 434 for d=1:s 435 [dr,cc] = ndgrid(dsites(:,d),ctrs(:,d)); 436 DM = DM + (dr-cc).^2; 437 end 438 DM = sqrt(DM); 439 440 441 % rbf = tps(e,r) 442 % Defines thin plate spline RBF 443 function rbf = tps(e,r) 444 rbf = zeros(size(r)); 445 nz = find(r~=0); % to deal with singularity at origin 446 rbf(nz) = (e*r(nz)).^2.*log(e*r(nz)); 447 426 427
Note: See TracChangeset
for help on using the changeset viewer.