source: trunk/src/set_subdomains.m @ 903

Last change on this file since 903 was 809, checked in by g7moreau, 10 years ago
  • Add license
File size: 4.2 KB
Line 
1%'set_subdomains': sort a set of points defined by scattered coordinates in subdomains, as needed for tps interpolation
2%------------------------------------------------------------------------
3% [SubRange,NbCentre,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% NbCentre(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 NbCentre 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
15%=======================================================================
16% Copyright 2008-2014, LEGI UMR 5519 / CNRS UJF G-INP, Grenoble, France
17%   http://www.legi.grenoble-inp.fr
18%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
19%
20%     This file is part of the toolbox UVMAT.
21%
22%     UVMAT is free software; you can redistribute it and/or modify
23%     it under the terms of the GNU General Public License as published
24%     by the Free Software Foundation; either version 2 of the license,
25%     or (at your option) any later version.
26%
27%     UVMAT is distributed in the hope that it will be useful,
28%     but WITHOUT ANY WARRANTY; without even the implied warranty of
29%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
30%     GNU General Public License (see LICENSE.txt) for more details.
31%=======================================================================
32
33function [SubRange,NbCentre,IndSelSubDomain] =set_subdomains(Coord,SubDomainNbPoint)
34
35%% adjust subdomain decomposition
36NbVec=size(Coord,1);
37NbCoord=size(Coord,2);
38MinCoord=min(Coord,[],1);
39MaxCoord=max(Coord,[],1);
40Range=MaxCoord-MinCoord;
41AspectRatio=Range(2)/Range(1);
42NbSubDomain=NbVec/SubDomainNbPoint;
43NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1);
44NbSubDomainY=max(floor(sqrt(NbSubDomain*AspectRatio)),1);
45NbSubDomain=NbSubDomainX*NbSubDomainY;
46Siz(1)=Range(1)/NbSubDomainX;%width of subdomains
47Siz(2)=Range(2)/NbSubDomainY;%height of subdomains
48CentreX=linspace(MinCoord(1)+Siz(1)/2,MaxCoord(1)-Siz(1)/2,NbSubDomainX);
49CentreY=linspace(MinCoord(2)+Siz(2)/2,MaxCoord(2)-Siz(2)/2,NbSubDomainY);
50[CentreX,CentreY]=meshgrid(CentreX,CentreY);
51CentreY=reshape(CentreY,1,[]);% Y positions of subdomain centres
52CentreX=reshape(CentreX,1,[]);% X positions of subdomain centres
53
54%% default output
55SubRange=zeros(NbCoord,2,NbSubDomain);%initialise the positions of subdomains
56NbCentre=zeros(1,NbSubDomain);%number of interpolated values per subdomain, =0 by default
57%Coord_tps=zeros(NbVec,NbCoord,NbSubDomain);% default positions of the tps source= initial positions of the good vectors sorted by subdomain
58check_empty=zeros(1,NbSubDomain);
59
60%% adjust subdomains
61for isub=1:NbSubDomain
62    SubRange(1,:,isub)=[CentreX(isub)-0.55*Siz(1) CentreX(isub)+0.55*Siz(1)];
63    SubRange(2,:,isub)=[CentreY(isub)-0.55*Siz(2) CentreY(isub)+0.55*Siz(2)];
64    IndSel_previous=[];
65    IndSel=0;
66    while numel(IndSel)>numel(IndSel_previous) %increase the subdomain if it contains less than SubDomainNbPoint/4 sources
67        IndSel_previous=IndSel;
68        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));
69        % if no source in the subdomain, skip the subdomain
70        if isempty(IndSel)
71            check_empty(isub)=1;   
72            break
73        % if too few selected sources, increase the subrange for next iteration
74        elseif numel(IndSel)<SubDomainNbPoint/4 && ~isequal( IndSel,IndSel_previous);
75            SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
76            SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
77        else
78            break
79        end
80    end
81    NbCentre(isub)=numel(IndSel);
82    IndSelSubDomain(1:numel(IndSel),isub)=IndSel;
83end
84
85
86%% remove empty subdomains
87ind_empty=find(check_empty);
88if ~isempty(ind_empty)
89    SubRange(:,:,ind_empty)=[];
90    NbCentre(ind_empty)=[];
91    IndSelSubDomain(:,ind_empty)=[];
92end
93
94
Note: See TracBrowser for help on using the repository browser.