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