| 2 |
michaesp |
1 |
function out = int_index (field,indices,mdv)
|
|
|
2 |
|
|
|
3 |
% General interpolation routine for data grid
|
|
|
4 |
% Usage: out = int_index(in,index,mdv)
|
|
|
5 |
% Example: out int_index(inp.TH.data,[10.1 20.4 30.3; 23.6 56.7 78.4],mdv)
|
|
|
6 |
|
|
|
7 |
% Check whether input data are correct
|
|
|
8 |
if ( length(size(field)) > 3 )
|
|
|
9 |
error('Error in int_index: input field can be at most three-dimensional');
|
|
|
10 |
end
|
|
|
11 |
if ( length(size(indices)) ~= 2 | ...
|
|
|
12 |
length(size(indices)) == 2 & size(indices,2) ~= length(size(field)) )
|
|
|
13 |
error('Input fields for int_index are inconsistent');
|
|
|
14 |
end
|
|
|
15 |
|
|
|
16 |
% Get the dimensions of the input field and the indices
|
|
|
17 |
if ( length(size(field)) == 3 )
|
|
|
18 |
nx = size(field,3);
|
|
|
19 |
ny = size(field,2);
|
|
|
20 |
nz = size(field,1);
|
|
|
21 |
xpos = indices(:,3);
|
|
|
22 |
ypos = indices(:,2);
|
|
|
23 |
zpos = indices(:,1);
|
|
|
24 |
elseif ( length(size(field)) == 2 )
|
|
|
25 |
nx = size(field,2);
|
|
|
26 |
ny = size(field,1);
|
|
|
27 |
nz = 1;
|
|
|
28 |
xpos = indices(:,2);
|
|
|
29 |
ypos = indices(:,1);
|
|
|
30 |
zpos = ones(size(xpos));
|
|
|
31 |
end
|
|
|
32 |
|
|
|
33 |
% Make the field one-dimensional (needed for vectorised computation)
|
|
|
34 |
field1 = field(:);
|
|
|
35 |
|
|
|
36 |
% Get the indices of the neighbouring data points (surrounding grid box)
|
|
|
37 |
ix000 = max(floor(xpos),1);
|
|
|
38 |
iy000 = max(floor(ypos),1);
|
|
|
39 |
iz000 = max(floor(zpos),1);
|
|
|
40 |
|
|
|
41 |
ix001 = min(1+floor(xpos),nx);
|
|
|
42 |
iy001 = iy000;
|
|
|
43 |
iz001 = iz000;
|
|
|
44 |
|
|
|
45 |
ix010 = ix000;
|
|
|
46 |
iy010 = min(1+floor(ypos),ny);
|
|
|
47 |
iz010 = iz000;
|
|
|
48 |
|
|
|
49 |
ix011 = min(1+floor(xpos),nx);
|
|
|
50 |
iy011 = min(1+floor(ypos),ny);
|
|
|
51 |
iz011 = iz000;
|
|
|
52 |
|
|
|
53 |
ix100 = ix000;
|
|
|
54 |
iy100 = iy000;
|
|
|
55 |
iz100 = min(1+floor(zpos),nz);
|
|
|
56 |
|
|
|
57 |
ix101 = min(1+floor(xpos),nx);
|
|
|
58 |
iy101 = iy000;
|
|
|
59 |
iz101 = min(1+floor(zpos),nz);
|
|
|
60 |
|
|
|
61 |
ix110 = ix000;
|
|
|
62 |
iy110 = min(1+floor(ypos),ny);
|
|
|
63 |
iz110 = min(1+floor(zpos),nz);
|
|
|
64 |
|
|
|
65 |
ix111 = min(1+floor(xpos),nx);
|
|
|
66 |
iy111 = min(1+floor(ypos),ny);
|
|
|
67 |
iz111 = min(1+floor(zpos),nz);
|
|
|
68 |
|
|
|
69 |
% Get the weights for the individual grid points
|
|
|
70 |
fracx = (1-xpos+ix000);
|
|
|
71 |
fracy = (1-ypos+iy000);
|
|
|
72 |
fracz = (1-zpos+iz000);
|
|
|
73 |
|
|
|
74 |
frac000 = fracx .* fracy .* fracz ;
|
|
|
75 |
frac001 = (1-fracx) .* fracy .* fracz ;
|
|
|
76 |
frac010 = fracx .* (1-fracy) .* fracz ;
|
|
|
77 |
frac011 = (1-fracx) .* (1-fracy) .* fracz ;
|
|
|
78 |
frac100 = fracx .* fracy .* (1-fracz);
|
|
|
79 |
frac101 = (1-fracx) .* fracy .* (1-fracz);
|
|
|
80 |
frac110 = fracx .* (1-fracy) .* (1-fracz);
|
|
|
81 |
frac111 = (1-fracx) .* (1-fracy) .* (1-fracz);
|
|
|
82 |
|
|
|
83 |
% Vectorised interpolation
|
|
|
84 |
out = frac000 .* field1( iz000 + (iy000-1) * nz + (ix000-1) * nz * ny ) + ...
|
|
|
85 |
frac001 .* field1( iz001 + (iy001-1) * nz + (ix001-1) * nz * ny ) + ...
|
|
|
86 |
frac010 .* field1( iz010 + (iy010-1) * nz + (ix010-1) * nz * ny ) + ...
|
|
|
87 |
frac011 .* field1( iz011 + (iy011-1) * nz + (ix011-1) * nz * ny ) + ...
|
|
|
88 |
frac100 .* field1( iz100 + (iy100-1) * nz + (ix100-1) * nz * ny ) + ...
|
|
|
89 |
frac101 .* field1( iz101 + (iy101-1) * nz + (ix101-1) * nz * ny ) + ...
|
|
|
90 |
frac110 .* field1( iz110 + (iy110-1) * nz + (ix110-1) * nz * ny ) + ...
|
|
|
91 |
frac111 .* field1( iz111 + (iy111-1) * nz + (ix111-1) * nz * ny );
|
|
|
92 |
|
|
|
93 |
% Handling of missing data values
|
|
|
94 |
mask = ( field1( iz000 + (iy000-1) * nz + (ix000-1) * nz * ny ) == mdv ) | ...
|
|
|
95 |
( field1( iz001 + (iy001-1) * nz + (ix001-1) * nz * ny ) == mdv ) | ...
|
|
|
96 |
( field1( iz010 + (iy010-1) * nz + (ix010-1) * nz * ny ) == mdv ) | ...
|
|
|
97 |
( field1( iz011 + (iy011-1) * nz + (ix011-1) * nz * ny ) == mdv ) | ...
|
|
|
98 |
( field1( iz100 + (iy100-1) * nz + (ix100-1) * nz * ny ) == mdv ) | ...
|
|
|
99 |
( field1( iz101 + (iy101-1) * nz + (ix101-1) * nz * ny ) == mdv ) | ...
|
|
|
100 |
( field1( iz110 + (iy110-1) * nz + (ix110-1) * nz * ny ) == mdv ) | ...
|
|
|
101 |
( field1( iz111 + (iy111-1) * nz + (ix111-1) * nz * ny ) == mdv );
|
|
|
102 |
out(mask) = mdv;
|
|
|
103 |
|
|
|
104 |
|
|
|
105 |
|
|
|
106 |
|
|
|
107 |
|
|
|
108 |
|