source: trunk/src/filter_tps.m @ 501

Last change on this file since 501 was 494, checked in by sommeria, 9 years ago

various bugs corrected after testing in Windows OS. Introduction
of filter tps

File size: 6.5 KB
Line 
1%'filter_tps': find the thin plate spline coefficients for interpolation-smoothing
2%------------------------------------------------------------------------
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)
4%
5% OUTPUT:
6% 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
8% FF: false flags
9% 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
11% U_tps,V_tps: weight of the tps for each subdomain
12% to get the interpolated field values, use the function calc_field.m
13%
14% INPUT:
15% coord=[X Y]: matrix whose first column is the x coordinates of the initial data, the second column the y coordiantes
16% U,V: set of velocity components of the initial data
17% Rho: smoothing parameter
18% Threshold: max diff accepted between smoothed and initial data
19% Subdomain: estimated number of data points in each subdomain
20
21function [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
23warning off
24nbvec=size(Coord,1);
25W_tps=[];%default
26W_smooth=[];
27NbCoord=size(Coord,2);
28NbSubDomain=nbvec/SubDomain;
29MinCoord=min(Coord,[],1);
30MaxCoord=max(Coord,[],1);
31Range=MaxCoord-MinCoord;
32AspectRatio=Range(2)/Range(1);
33NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1);
34NbSubDomainY=max(floor(sqrt(NbSubDomain*AspectRatio)),1);
35NbSubDomain=NbSubDomainX*NbSubDomainY;
36Siz(1)=Range(1)/NbSubDomainX;%width of subdomains
37Siz(2)=Range(2)/NbSubDomainY;%height of subdomains
38CentreX=linspace(MinCoord(1)+Siz(1)/2,MaxCoord(1)-Siz(1)/2,NbSubDomainX);
39CentreY=linspace(MinCoord(2)+Siz(2)/2,MaxCoord(2)-Siz(2)/2,NbSubDomainY);
40[CentreX,CentreY]=meshgrid(CentreX,CentreY);
41CentreY=reshape(CentreY,1,[]);% Y positions of subdomain centres
42CentreX=reshape(CentreX,1,[]);% X positions of subdomain centres
43rho=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)
44U_tps_sub=zeros(nbvec,NbSubDomain);%default spline
45V_tps_sub=zeros(nbvec,NbSubDomain);%default spline
46Indices_tps=zeros(nbvec,NbSubDomain);%default indices
47U_smooth=zeros(nbvec,1);
48V_smooth=zeros(nbvec,1);
49nb_select=zeros(nbvec,1);
50FF=zeros(nbvec,1);
51check_empty=zeros(1,NbSubDomain);
52SubRange=zeros(NbCoord,2,NbSubDomain);%initialise the positions of subdomains
53% SubRangy=zeros(NbSubDomain,2);
54for isub=1:NbSubDomain
55    SubRange(1,:,isub)=[CentreX(isub)-0.55*Siz(1) CentreX(isub)+0.55*Siz(1)];
56    SubRange(2,:,isub)=[CentreY(isub)-0.55*Siz(2) CentreY(isub)+0.55*Siz(2)];
57    ind_sel_previous=[];
58    ind_sel=0;
59    while numel(ind_sel)>numel(ind_sel_previous) %increase the subdomain during four iterations at most
60        ind_sel_previous=ind_sel;
61        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));
62        % if no vector in the subdomain, skip the subdomain
63        if isempty(ind_sel)
64            check_empty(isub)=1;   
65            U_tps(1,isub)=0;%define U_tps and V_tps by default
66            V_tps(1,isub)=0;
67            break
68            % if too few selected vectors, increase the subrange for next iteration
69        elseif numel(ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous);
70            SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
71            SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
72        else         
73            [U_smooth_sub,U_tps_sub]=tps_coeff(Coord(ind_sel,:),U(ind_sel),rho);
74            [V_smooth_sub,V_tps_sub]=tps_coeff(Coord(ind_sel,:),V(ind_sel),rho);
75            UDiff=U_smooth_sub-U(ind_sel);
76            VDiff=V_smooth_sub-V(ind_sel);
77            NormDiff=UDiff.*UDiff+VDiff.*VDiff;
78            ind_ind_sel=1:numel(ind_sel);%default
79            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
82            end
83            % if no value exceeds threshold, the result is recorded
84            if isequal(numel(ind_ind_sel),numel(ind_sel))
85                U_smooth(ind_sel)=U_smooth(ind_sel)+U_smooth_sub;
86                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;         
91                nb_select(ind_sel)=nb_select(ind_sel)+1;
92                display('good')
93                break
94            % if too few selected vectors, increase the subrange for next iteration
95            elseif numel(ind_ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous);
96                SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
97                SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
98            % else interpolation-smoothing is done again with the selected vectors
99            else
100                [U_smooth_sub,U_tps_sub]=tps_coeff(Coord(ind_sel(ind_ind_sel),:),U(ind_sel(ind_ind_sel)),rho);
101                [V_smooth_sub,V_tps_sub]=tps_coeff(Coord(ind_sel(ind_ind_sel),:),V(ind_sel(ind_ind_sel)),rho);
102                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;
108                nb_select(ind_sel(ind_ind_sel))=nb_select(ind_sel(ind_ind_sel))+1;
109                display('good2')
110                break
111            end
112        end
113    end
114end
115ind_empty=find(check_empty);
116%remove empty subdomains
117if ~isempty(ind_empty)
118    SubRange(:,:,ind_empty)=[];
119    Coord_tps(:,:,ind_empty)=[];
120    U_tps(:,ind_empty)=[];
121    V_tps(:,ind_empty)=[];
122end
123nb_select(nb_select==0)=1;%ones(size(find(nb_select==0)));
124U_smooth=U_smooth./nb_select;
125V_smooth=V_smooth./nb_select;
126fill=zeros(NbCoord+1,NbCoord,size(SubRange,3)); %matrix of zeros to complement the matrix Data.Civ1_Coord_tps (conveninent for file storage)
127Coord_tps=cat(1,Coord_tps,fill);
128
Note: See TracBrowser for help on using the repository browser.