source: trunk/src/filter_tps.m @ 503

Last change on this file since 503 was 494, checked in by sommeria, 12 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.