Changeset 243 for trunk/src


Ignore:
Timestamp:
Apr 19, 2011, 5:41:35 PM (13 years ago)
Author:
sommeria
Message:

bug repair in uvmat

Location:
trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/civ_uvmat.m

    r242 r243  
    129129        ind_good=1:numel(Data.Civ1_X);
    130130    end
     131    Data.Civ1_X
    131132    [Ures, Vres,FFres]=...
    132133                            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);
     
    269270%------------------------------------------------------------------------
    270271% patch function
    271 function [U_smooth,V_smooth,FF] =patch_uvmat(X,Y,U,V,Rho,Threshold,SubDomain)
     272function [U_patch,V_patch,FF,SubRangx,SubRangy,X_ctrs,Y_ctrs] =patch_uvmat(X,Y,U,V,Rho,Threshold,SubDomain)
    272273%subdomain decomposition
    273274warning off
     275U=reshape(U,[],1);
     276V=reshape(V,[],1);
     277X=reshape(X,[],1);
     278Y=reshape(Y,[],1);
    274279nbvec=numel(X);
    275 NbSubDomain=ceil(nbvec/SubDomain);
    276 MinX=min(X);
    277 MinY=min(Y);
    278 MaxX=max(X);
    279 MaxY=max(Y);
     280NbSubDomain=ceil(nbvec/SubDomain)
     281MinX=min(X)
     282MinY=min(Y)
     283MaxX=max(X)
     284MaxY=max(Y)
    280285RangX=MaxX-MinX;
    281286RangY=MaxY-MinY;
    282 AspectRatio=RangY/RangX;
     287AspectRatio=RangY/RangX
    283288NbSubDomainX=ceil(sqrt(NbSubDomain/AspectRatio))
    284289NbSubDomainY=ceil(sqrt(NbSubDomain*AspectRatio))
     
    291296SubIndexY=ceil((Y-MinY)/SizY);
    292297rho=RangX*RangY*Rho;%optimum rho increase as teh area of the subdomain
    293 U_tps=zeros(length(X)+3,1);%default spline
    294 V_tps=zeros(length(X)+3,1);%default spline
     298U_tps=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
     299V_tps=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
     300X_dist=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
     301Y_dist=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
     302U_smooth=zeros(length(X),NbSubDomainY,NbSubDomainX);
     303V_smooth=zeros(length(X),NbSubDomainY,NbSubDomainX);
     304dist_ctre=zeros(length(X),NbSubDomainY,NbSubDomainX);
     305FF=zeros(length(X),1);
    295306for isubx=1:NbSubDomainX
    296307    for isuby=1:NbSubDomainY
    297         SubRangx(isuby,isubx,:)=[CentreX(isubx)-SizX/2 CentreX(isubx)+SizX/2]
    298         SubRangy(isuby,isubx,:)=[CentreY(isubx)-SizY/2 CentreY(isubx)+SizY/2]
     308        SubRangx(isuby,isubx,:)=[CentreX(isubx)-SizX/2 CentreX(isubx)+SizX/2];
     309        SubRangy(isuby,isubx,:)=[CentreY(isubx)-SizY/2 CentreY(isubx)+SizY/2];
    299310        for iter=1:3 %increase the subdomain during three iterations at most
    300311            ind_sel=find(X>SubRangx(isuby,isubx,1) & X<SubRangx(isuby,isubx,2) & Y>SubRangy(isuby,isubx,1) & Y<SubRangy(isuby,isubx,2)); 
    301             [U_smooth(ind_sel)]=tps_uvmat(X(ind_sel),Y(ind_sel),U(ind_sel),rho);
    302             [V_smooth(ind_sel)]=tps_uvmat(X(ind_sel),Y(ind_sel),V(ind_sel),rho);
    303             iter
    304             numel(ind_sel)
     312            size(ind_sel)
    305313            if numel(ind_sel)<SubDomain/4;% too few selected vectors, increase the subrange for next iteration
    306314                SubRangx(isuby,isubx,1)=SubRangx(isuby,isubx,1)-SizX/4;
     
    309317                SubRangy(isuby,isubx,2)=SubRangy(isuby,isubx,2)+SizY/4;
    310318            else
    311                 UDiff=U_smooth(ind_sel)-U(ind_sel);
    312                 VDiff=U_smooth(ind_sel)-U(ind_sel);
     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);
    313325                NormDiff=UDiff.*UDiff+VDiff.*VDiff;
    314326                FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination
    315                 ind_ind_false=find(FF(ind_sel)~=0);
    316                 ind_ind_sel=find(FF(ind_sel)==0);
    317                 U_smooth(ind_sel(ind_ind_false))=zeros(numel(ind_ind_false),1);
    318                 V_smooth(ind_sel(ind_ind_false))=zeros(numel(ind_ind_false),1);         
    319                 [U_smooth(ind_sel(ind_ind_sel))]=tps_uvmat(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),U(ind_sel(ind_ind_sel)),rho);
    320                 [V_smooth(ind_sel(ind_ind_sel))]=tps_uvmat(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),V(ind_sel(ind_ind_sel)),rho);
    321                  if numel(ind_ind_sel)<SubDomain/4;% too few selected vectors, increase the subrange for next iteration
     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
    322333                    SubRangx(isuby,isubx,1)=SubRangx(isuby,isubx,1)-SizX/4;
    323334                    SubRangx(isuby,isubx,2)=SubRangx(isuby,isubx,2)+SizX/4;
    324335                    SubRangy(isuby,isubx,1)=SubRangy(isuby,isubx,1)-SizY/4;
    325336                    SubRangy(isuby,isubx,2)=SubRangy(isuby,isubx,2)+SizY/4;
    326                  else
    327                      break
    328                  end
     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
    329342            end       
    330         end
    331         %intrpolate to new points XI, YI
    332 %         X_ctrs{isuby,isubx}=X(ind_sel(ind_ind_sel));%positions of tps sources for the subdomain i,j
    333 %         Y_ctrs{isuby,isubx}=Y(ind_sel(ind_ind_sel));     
    334     end
    335 end
    336 
     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
     348end
     349U_patch=sum(sum(U_smooth.*dist_ctre,3),2)./sum(sum(dist_ctre,3),2);
     350V_patch=sum(sum(V_smooth.*dist_ctre,3),2)./sum(sum(dist_ctre,3),2);
    337351
    338352%------------------------------------------------------------------------
     
    340354% X,Y initial coordiantes
    341355% XI vector, YI column vector for the grid of interpolation points
    342 function [U_smooth,U_tps]=tps_uvmat(X,Y,U,rho)
     356function [U_smooth,U_tps,U_tps3]=tps_uvmat(X,Y,U,rho)
    343357%------------------------------------------------------------------------
    344358%rho smoothing parameter
     
    378392EM = [IM_sites PM];
    379393U_smooth=EM *U_tps;
     394U_tps3=U_tps(end-2:end);
     395U_tps=U_tps(1:end-3);
    380396
    381397
  • trunk/src/pivlab.m

    r238 r243  
    136136    if y <= npy-1 && y >= 1
    137137        f0 = log(result_conv(y,x));
    138         f1 = log(result_conv(y-1,x));
    139         f2 = log(result_conv(y+1,x));
     138        f1 = real(log(result_conv(y-1,x)));
     139        f2 = real(log(result_conv(y+1,x)));
    140140        peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
    141141    else
     
    145145    if x <= npx-1 && x >= 1
    146146        f0 = log(result_conv(y,x));
    147         f1 = log(result_conv(y,x-1));
    148         f2 = log(result_conv(y,x+1));
     147        f1 = real(log(result_conv(y,x-1)));
     148        f2 = real(log(result_conv(y,x+1)));
    149149        peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
    150150    else
  • trunk/src/uvmat.m

    r242 r243  
    24272427   UvData.Field=Field{1};
    24282428end
    2429 UvData.Field
    2430 max(max(UvData.Field.A))
    2431 min(min(UvData.Field.A))
     2429
    24322430%% get bounds and mesh (needed for mouse action and to open set_object)
    24332431test_x=0;
Note: See TracChangeset for help on using the changeset viewer.