Changeset 581 for trunk/src/filter_tps.m


Ignore:
Timestamp:
Mar 12, 2013, 12:52:13 PM (8 years ago)
Author:
sommeria
Message:

clean the transform field functions

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/filter_tps.m

    r494 r581  
    11%'filter_tps': find the thin plate spline coefficients for interpolation-smoothing
    22%------------------------------------------------------------------------
    3 % [SubRange,NbSites,Coord_tps,U_tps,V_tps,W_tps,U_smooth,V_smooth,W_smooth,FF] =filter_tps(Coord,U,V,W,SubDomain,Rho,Threshold)
     3% [SubRange,NbCentres,Coord_tps,U_tps,V_tps,W_tps,U_smooth,V_smooth,W_smooth,FF] =filter_tps(Coord,U,V,W,SubDomain,Rho,Threshold)
    44%
    55% OUTPUT:
    66% SubRange(NbCoord,NbSubdomain,2): range (min, max) of the coordiantes x and y respectively, for each subdomain
    7 % NbSites(NbSubdomain): number of source points for each subdomain
     7% NbCentres(NbSubdomain): number of source points for each subdomain
    88% FF: false flags
    99% U_smooth, V_smooth: filtered velocity components at the positions of the initial data
    10 % Coord_tps(NbSites,NbCoord,NbSubdomain): positions of the tps centres
     10% Coord_tps(NbCentres,NbCoord,NbSubdomain): positions of the tps centres
    1111% U_tps,V_tps: weight of the tps for each subdomain
    1212% to get the interpolated field values, use the function calc_field.m
     
    1919% Subdomain: estimated number of data points in each subdomain
    2020
    21 function [SubRange,NbSites,Coord_tps,U_tps,V_tps,W_tps,U_smooth,V_smooth,W_smooth,FF] =filter_tps(Coord,U,V,W,SubDomain,Rho,Threshold)
    22 %subdomain decomposition
     21function [SubRange,NbCentres,Coord_tps,U_tps,V_tps,W_tps,U_smooth,V_smooth,W_smooth,FF] =filter_tps(Coord,U,V,W,SubDomain,Rho,Threshold)
     22
     23%% adjust subdomain decomposition
    2324warning off
    24 nbvec=size(Coord,1);
    25 W_tps=[];%default
    26 W_smooth=[];
     25NbVec=size(Coord,1);
    2726NbCoord=size(Coord,2);
    28 NbSubDomain=nbvec/SubDomain;
    29 MinCoord=min(Coord,[],1);
    30 MaxCoord=max(Coord,[],1);
     27MinCoord=min(Coord,[],1);%lower coordinate bounds
     28MaxCoord=max(Coord,[],1);%upper coordinate bounds
    3129Range=MaxCoord-MinCoord;
    3230AspectRatio=Range(2)/Range(1);
     31NbSubDomain=NbVec/SubDomain;
    3332NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1);
    3433NbSubDomainY=max(floor(sqrt(NbSubDomain*AspectRatio)),1);
     
    4140CentreY=reshape(CentreY,1,[]);% Y positions of subdomain centres
    4241CentreX=reshape(CentreX,1,[]);% X positions of subdomain centres
    43 rho=Siz(1)*Siz(2)*Rho/1000;%optimum rho increase as the area of the subdomain (division by 10^6 to reach good values with the default GUI input)
    44 U_tps_sub=zeros(nbvec,NbSubDomain);%default spline
    45 V_tps_sub=zeros(nbvec,NbSubDomain);%default spline
    46 Indices_tps=zeros(nbvec,NbSubDomain);%default indices
    47 U_smooth=zeros(nbvec,1);
    48 V_smooth=zeros(nbvec,1);
    49 nb_select=zeros(nbvec,1);
    50 FF=zeros(nbvec,1);
     42
     43%% smoothing parameter
     44rho=Siz(1)*Siz(2)*Rho/1000;%optimum rho increase as the area of the subdomain (division by 1000 to reach good values with the default GUI input)
     45
     46%% default output
     47SubRange=zeros(NbCoord,2,NbSubDomain);%initialise the positions of subdomains
     48NbCentres=zeros(NbSubDomain);%number of interpolated values per subdomain, =0 by default
     49Coord_tps=zeros(NbVec,NbCoord,NbSubDomain);% default positions of the tps source= initial positions of the good vectors sorted by subdomain
     50U_tps=zeros(NbVec,NbSubDomain);%default spline
     51V_tps=zeros(NbVec,NbSubDomain);%default spline
     52W_tps=[];%default (2 component case)
     53U_smooth=zeros(NbVec,1); % smoothed velocity U at the initial positions
     54V_smooth=zeros(NbVec,1);% smoothed velocity V at the initial positions
     55W_smooth=[];%default (2 component case)
     56FF=zeros(NbVec,1);
     57nb_select=zeros(NbVec,1);
    5158check_empty=zeros(1,NbSubDomain);
    52 SubRange=zeros(NbCoord,2,NbSubDomain);%initialise the positions of subdomains
    53 % SubRangy=zeros(NbSubDomain,2);
     59
     60%% calculate tps coeff in each subdomain
    5461for isub=1:NbSubDomain
    5562    SubRange(1,:,isub)=[CentreX(isub)-0.55*Siz(1) CentreX(isub)+0.55*Siz(1)];
     
    5764    ind_sel_previous=[];
    5865    ind_sel=0;
    59     while numel(ind_sel)>numel(ind_sel_previous) %increase the subdomain during four iterations at most
     66    %increase iteratively the subdomain if it contains less than SubDomainNbVec/4 source vectors
     67    while numel(ind_sel)>numel(ind_sel_previous)
    6068        ind_sel_previous=ind_sel;
    6169        ind_sel=find(Coord(:,1)>=SubRange(1,1,isub) & Coord(:,1)<=SubRange(1,2,isub) & Coord(:,2)>=SubRange(2,1,isub) & Coord(:,2)<=SubRange(2,2,isub));
    6270        % if no vector in the subdomain, skip the subdomain
    6371        if isempty(ind_sel)
    64             check_empty(isub)=1;   
     72            check_empty(isub)=1;
    6573            U_tps(1,isub)=0;%define U_tps and V_tps by default
    6674            V_tps(1,isub)=0;
     
    7078            SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
    7179            SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
    72         else         
     80        else
    7381            [U_smooth_sub,U_tps_sub]=tps_coeff(Coord(ind_sel,:),U(ind_sel),rho);
    7482            [V_smooth_sub,V_tps_sub]=tps_coeff(Coord(ind_sel,:),V(ind_sel),rho);
     
    7886            ind_ind_sel=1:numel(ind_sel);%default
    7987            if exist('Threshold','var')
    80             FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination
    81             ind_ind_sel=find(FF(ind_sel)==0); % select the indices of ind_sel corresponding to the remaining vectors
     88                FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination
     89                ind_ind_sel=find(FF(ind_sel)==0); % select the indices of ind_sel corresponding to the remaining vectors
    8290            end
    8391            % if no value exceeds threshold, the result is recorded
     
    8593                U_smooth(ind_sel)=U_smooth(ind_sel)+U_smooth_sub;
    8694                V_smooth(ind_sel)=V_smooth(ind_sel)+V_smooth_sub;
    87                 NbSites(isub)=numel(ind_sel);
    88                 Coord_tps(1:NbSites(isub),:,isub)=Coord(ind_sel,:);
    89                 U_tps(1:NbSites(isub)+3,isub)=U_tps_sub;
    90                 V_tps(1:NbSites(isub)+3,isub)=V_tps_sub;         
     95                NbCentres(isub)=numel(ind_sel);
     96                Coord_tps(1:NbCentres(isub),:,isub)=Coord(ind_sel,:);
     97                U_tps(1:NbCentres(isub)+3,isub)=U_tps_sub;
     98                V_tps(1:NbCentres(isub)+3,isub)=V_tps_sub;
    9199                nb_select(ind_sel)=nb_select(ind_sel)+1;
    92100                display('good')
    93101                break
    94             % if too few selected vectors, increase the subrange for next iteration
     102                % if too few selected vectors, increase the subrange for next iteration
    95103            elseif numel(ind_ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous);
    96104                SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
    97105                SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
    98             % else interpolation-smoothing is done again with the selected vectors
     106                % else interpolation-smoothing is done again with the selected vectors
    99107            else
    100108                [U_smooth_sub,U_tps_sub]=tps_coeff(Coord(ind_sel(ind_ind_sel),:),U(ind_sel(ind_ind_sel)),rho);
    101109                [V_smooth_sub,V_tps_sub]=tps_coeff(Coord(ind_sel(ind_ind_sel),:),V(ind_sel(ind_ind_sel)),rho);
    102110                U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+U_smooth_sub;
    103                 V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+V_smooth_sub;   
    104                 NbSites(isub)=numel(ind_ind_sel);
    105                 Coord_tps(1:NbSites(isub),:,isub)=Coord(ind_sel(ind_ind_sel),:);
    106                 U_tps(1:NbSites(isub)+3,isub)=U_tps_sub;
    107                 V_tps(1:NbSites(isub)+3,isub)=V_tps_sub;
     111                V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+V_smooth_sub;
     112                NbCentres(isub)=numel(ind_ind_sel);
     113                Coord_tps(1:NbCentres(isub),:,isub)=Coord(ind_sel(ind_ind_sel),:);
     114                U_tps(1:NbCentres(isub)+3,isub)=U_tps_sub;
     115                V_tps(1:NbCentres(isub)+3,isub)=V_tps_sub;
    108116                nb_select(ind_sel(ind_ind_sel))=nb_select(ind_sel(ind_ind_sel))+1;
    109117                display('good2')
     
    113121    end
    114122end
     123
     124%% remove empty subdomains
    115125ind_empty=find(check_empty);
    116 %remove empty subdomains
    117126if ~isempty(ind_empty)
    118127    SubRange(:,:,ind_empty)=[];
     
    121130    V_tps(:,ind_empty)=[];
    122131end
    123 nb_select(nb_select==0)=1;%ones(size(find(nb_select==0)));
    124 U_smooth=U_smooth./nb_select;
     132
     133%% final adjustments
     134nb_select(nb_select==0)=1;
     135U_smooth=U_smooth./nb_select;% take the average at the intersection of several subdomains
    125136V_smooth=V_smooth./nb_select;
    126137fill=zeros(NbCoord+1,NbCoord,size(SubRange,3)); %matrix of zeros to complement the matrix Data.Civ1_Coord_tps (conveninent for file storage)
Note: See TracChangeset for help on using the changeset viewer.