source: trunk/src/set_subdomains.m @ 643

Last change on this file since 643 was 587, checked in by sommeria, 11 years ago

tps_coeff_field introduced and several bugs corrected

File size: 3.3 KB
Line 
1%'set_subdomains': sort a set of points defined by scattered coordinates in subdomains, as needed for tps interpolation
2%------------------------------------------------------------------------
3% [SubRange,NbSites,IndSelSubDomain] =set_subdomains(Coord,SubDomainNbPoint)
4%
5% OUTPUT:
6% SubRange(NbCoord,NbSubdomain,2): range (min, max) of the coordinates x and y respectively, for each subdomain
7% NbSites(NbSubdomain): number of source points for each subdomain
8% IndSelSubDomain(SubDomainNbPointMax,NbSubdomain): set of indices of the input point array
9% selected in each subdomain, =0 beyond NbSites points
10%
11% INPUT:
12% coord=[X Y]: matrix whose first column is the x coordinates of the input data points, the second column the y coordinates
13% SubdomainNbPoint: estimated number of data points whished for each subdomain
14
15function [SubRange,NbSites,IndSelSubDomain] =set_subdomains(Coord,SubDomainNbPoint)
16
17%% adjust subdomain decomposition
18NbVec=size(Coord,1);
19NbCoord=size(Coord,2);
20MinCoord=min(Coord,[],1);
21MaxCoord=max(Coord,[],1);
22Range=MaxCoord-MinCoord;
23AspectRatio=Range(2)/Range(1);
24NbSubDomain=NbVec/SubDomainNbPoint;
25NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1);
26NbSubDomainY=max(floor(sqrt(NbSubDomain*AspectRatio)),1);
27NbSubDomain=NbSubDomainX*NbSubDomainY;
28Siz(1)=Range(1)/NbSubDomainX;%width of subdomains
29Siz(2)=Range(2)/NbSubDomainY;%height of subdomains
30CentreX=linspace(MinCoord(1)+Siz(1)/2,MaxCoord(1)-Siz(1)/2,NbSubDomainX);
31CentreY=linspace(MinCoord(2)+Siz(2)/2,MaxCoord(2)-Siz(2)/2,NbSubDomainY);
32[CentreX,CentreY]=meshgrid(CentreX,CentreY);
33CentreY=reshape(CentreY,1,[]);% Y positions of subdomain centres
34CentreX=reshape(CentreX,1,[]);% X positions of subdomain centres
35
36%% default output
37SubRange=zeros(NbCoord,2,NbSubDomain);%initialise the positions of subdomains
38NbSites=zeros(NbSubDomain);%number of interpolated values per subdomain, =0 by default
39%Coord_tps=zeros(NbVec,NbCoord,NbSubDomain);% default positions of the tps source= initial positions of the good vectors sorted by subdomain
40check_empty=zeros(1,NbSubDomain);
41
42%% adjust subdomains
43for isub=1:NbSubDomain
44    SubRange(1,:,isub)=[CentreX(isub)-0.55*Siz(1) CentreX(isub)+0.55*Siz(1)];
45    SubRange(2,:,isub)=[CentreY(isub)-0.55*Siz(2) CentreY(isub)+0.55*Siz(2)];
46    IndSel_previous=[];
47    IndSel=0;
48    while numel(IndSel)>numel(IndSel_previous) %increase the subdomain if it contains less than SubDomainNbPoint/4 sources
49        IndSel_previous=IndSel;
50        IndSel=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));
51        % if no source in the subdomain, skip the subdomain
52        if isempty(IndSel)
53            check_empty(isub)=1;   
54            break
55        % if too few selected sources, increase the subrange for next iteration
56        elseif numel(IndSel)<SubDomainNbPoint/4 && ~isequal( IndSel,IndSel_previous);
57            SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
58            SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
59        else
60            break
61        end
62    end
63    NbSites(isub)=numel(IndSel);
64    IndSelSubDomain(:,isub)=IndSel;
65end
66
67
68%% remove empty subdomains
69ind_empty=find(check_empty);
70if ~isempty(ind_empty)
71    SubRange(:,:,ind_empty)=[];
72    NbSites(ind_empty)=[];
73    IndSelSubDomain(:,ind_empty)=[];
74end
75
76
Note: See TracBrowser for help on using the repository browser.