U.S. patent application number 10/660570 was filed with the patent office on 2004-03-18 for metid if interpolation.
This patent application is currently assigned to Hewlett-Packard Company. Invention is credited to Hemingway, Peter.
Application Number | 20040051892 10/660570 |
Document ID | / |
Family ID | 8241277 |
Filed Date | 2004-03-18 |
United States Patent
Application |
20040051892 |
Kind Code |
A1 |
Hemingway, Peter |
March 18, 2004 |
Metid if interpolation
Abstract
A method of determining a value for a function which is
particularly useful for mapping values from one colour space to
another comprises a series of steps, as follows. The method is
applicable to n-dimensional spaces, but is particularly described
for three dimensions. The first step is to establish a three
dimensional lattice, the function having values at the lattice
points. The next step is to record values of the function for a
subset of the lattice points, the lattice points of the subset
being known value lattice points. These known value lattice points
form a sparse lattice (preferably the sparse lattice points are
regularly spaced along orthogonal axes). A value of the function
for a given lattice point is established by returning a weighted
average of the values of one or more of four known value lattice
points defining a tetrahedron 42 touching or enclosing the given
lattice point. Each of the lattice points which are intermediate
points in a coarse lattice cube 41 is either within, or on the
boundary of, at least one tetrahedron 42 whose vertices are four of
the vertices of the cube.
Inventors: |
Hemingway, Peter; (Bristol,
GB) |
Correspondence
Address: |
Allan M. Lowe
LOWE HAUPTMAN GILMAN & BERNER, LLP
Suite 300
1700 Diagonal Road
Alexandria
VA
22314
US
|
Assignee: |
Hewlett-Packard Company
Palo Alto
CA
|
Family ID: |
8241277 |
Appl. No.: |
10/660570 |
Filed: |
September 12, 2003 |
Related U.S. Patent Documents
|
|
|
|
|
|
Application
Number |
Filing Date |
Patent Number |
|
|
10660570 |
Sep 12, 2003 |
|
|
|
09527860 |
Mar 17, 2000 |
|
|
|
Current U.S.
Class: |
358/1.9 ;
358/518; 358/525; 382/167 |
Current CPC
Class: |
H04N 1/6019
20130101 |
Class at
Publication: |
358/001.9 ;
358/518; 358/525; 382/167 |
International
Class: |
H04N 001/60 |
Foreign Application Data
Date |
Code |
Application Number |
Mar 19, 1999 |
EP |
99302139.3 |
Claims
1. A method of determining a value for a function, comprising:
establishing an n-dimensional lattice, the function having values
at the lattice points, and where n is greater than or equal to two;
recording values for a subset of the lattice points, the lattice
points of the subset being known value lattice points; and
establishing a value for a given lattice point by returning a
weighted average of the values of one or more of (n+1) known value
lattice points defining an n-simplex touching or enclosing the
given lattice point.
2. A method as claimed in claim 1, wherein n=3 and the n-simplex is
a tetrahedron.
3. A method as claimed in claim 2, wherein a weighted average of
all four known value lattice point values is used if the given
lattice point is enclosed by the tetrahedron but is not touched by
a face of the tetrahedron, a weighted average of three of the four
known value lattice point values is used if the given lattice point
is on a face of the tetrahedron bounded by the three of the four
known value lattice points but is not touched by an edge of the
tetrahedron, a weighted average of two of the four known value
lattice point values is used if the given lattice point is on an
edge of the tetrahedron bounded by the two of the four known value
lattice points but is not at a vertex of the tetrahedron, and
wherein a value of one of the known value lattice points is used if
the given lattice point is also the known value lattice point.
4. A method as claimed in claim 3, wherein a given lattice point
close to a known value lattice point is changed to the known value
lattice point before calculation of the weighted average.
5. A method as claimed in claim 3 or claim 4, wherein a given
lattice point close to an edge or a face of the tetrahedron is
changed to a point lying on the edge or the face of the tetrahedron
before calculation of the weighted average.
6. A method as claimed in any of claims 3 to 5, wherein if the
given lattice point is enclosed by the tetrahedron but is not
touched by a face of the tetrahedron, and the tetrahedron has
vertices of known value lattice points with positions A, B, C, D
and values a, b, c, d at the respective vertices, and wherein the
given lattice point has position P and wherein the volume between
four positions is expressed as Vol(position 1 position 2 position 3
position 4), the value p returned is given by:
p=(Vol(ABCP).d+Vol(ABDP).c+Vol(ACDP).b+Vol(BCDP).a)/Vol(ABCD)
7. A method as claimed in any of claims 3 to 6, wherein if the
given lattice point is on a face of the tetrahedron bounded by the
three of the four known value lattice points but is not touched by
an edge of the tetrahedron, the three of the four known value
lattice points being A, B and C with values a, b and c
respectively, the value p returned is given by p=((Area(BCP).
a)+(Area(ACP). b)+(Area(ABP). c)/Area(ABC).
8. A method as claimed in any of claims 3 to 7, wherein if the
given lattice point is on an edge of the tetrahedron bounded by the
two of the four known value lattice points but is not at a vertex
of the tetrahedron, the two of the four known lattice points being
A and B with values a and b, the value p returned is given by
p=((Distance (AP). b)+(Distance (BP). a))/Distance(AB).
9. A method as claimed in any preceding claim, wherein the known
value lattice points form a sparse lattice with known value lattice
points separated from each other by an integer multiple of the
distance between adjacent lattice points.
10. A method as claimed in claim 9, wherein said integer multiple
is an integer power of two.
11. A method as claimed in claim 10, wherein the integer is 4 and
all given lattice points coincide with a value lattice point or lie
between two adjacent value lattice points or lie within a triangle
described by three adjacent value lattice points.
12. A method as claimed in claim 10, wherein the integer is 8 or
more and all given lattice points coincide with a value lattice
point or lie between two adjacent value lattice points or lie
within a triangle described by three adjacent value lattice points
or lie within or lie within a tetrahedron of four adjacent value
lattice points.
13. A method as claimed in claim 12, where the integer is 8.
14. A method as claimed in any of claims 2 to 13, wherein the step
of establishing a value comprises determining a set of four known
value lattice points which form a tetrahedron touching or enclosing
the given lattice point, and providing the weighted average from
the positions of four known value lattice points, the known values
of one or more of the four known value lattice points, and the
position of the given lattice point.
15. A method as claimed in claim 14, wherein the step of providing
the weighted average comprises using the positions as inputs to a
jump table.
16. A method of mapping values in a first colour space to values in
a second colour space, comprising establishing the value in the
second colour space by the method of determining a value for a
function described in any of claims 1 to 15.
17. A computer programmed to determine a value for a function, by:
establishing an n-dimensional lattice, the function having values
at the lattice points, and where n is greater than or equal to two;
recording values for a subset of the lattice points, the lattice
points of the subset being known value lattice points; and
establishing a value for a given lattice point by returning a
weighted average of the values of one or more of (n+1) known value
lattice points defining an n-simplex touching or enclosing the
given lattice point.
18. A computer as claimed in claim 17, wherein n=3 and the
n-simplex is a tetrahedron.
19. A computer as claimed in claim 18 wherein the computer is
programmed such that a weighted average of all four known value
lattice point values is used if the given lattice point is enclosed
by the tetrahedron but is not touched by a face of the tetrahedron,
a weighted average of three of the four known value lattice point
values is used if the given lattice point is on a face of the
tetrahedron bounded by the three of the four known value lattice
points but is not touched by an edge of the tetrahedron, a weighted
average of two of the four known value lattice point values is used
if the given lattice point is on an edge of the tetrahedron bounded
by the two of the four known value lattice points but is not at a
vertex of the tetrahedron, and wherein a value of one of the known
value lattice points is used if the given lattice point is also the
known value lattice point.
20. A program storage medium readable by a computer, tangibly
embodying a program of instructions executable by the computer to
perform method steps for determining a value for a function, said
method steps comprising: establishing an n-dimensional lattice, the
function having values at the lattice points, and where n is
greater than or equal to two; recording values for a subset of the.
lattice points, the lattice points of the subset being known value
lattice points; and establishing a value for a given lattice point
by returning a weighted average of the values of one or more of
(n+1) known value lattice points defining an n-simplex touching or
enclosing the given lattice point.
21. A program storage medium as claimed in claim 20, wherein n=3
and the n-simplex is a tetrahedron.
22. A program storage medium as claimed in claim 21, wherein in the
step of establishing a value, a weighted average of all four known
value lattice point values is used if the given lattice point is
enclosed by the tetrahedron but is not touched by a face of the
tetrahedron, a weighted average of three of the four known value
lattice point values is used if the given lattice point is on a
face of the tetrahedron bounded by the three of the four known
value lattice points but is not touched by an edge of the
tetrahedron, a weighted average of two of the four known value
lattice point values is used if the given lattice point is on an
edge of the tetrahedron bounded by the two of the four known value
lattice points but is not at a vertex of the tetrahedron, and
wherein a value of one of the known value lattice points is used if
the given lattice point is also the known value lattice point.
Description
FIELD OF INVENTION
[0001] The present invention relates to a method of interpolation
between known values to find intermediate values.
BACKGROUND ART
[0002] Interpolation has a number of applications, particularly in
the field of image processing. A particularly important use is in
colour mapping: the taking of values in one colour space and
mapping them on to corresponding values in another colour space.
Colour mapping is important, because different devices may function
best or most naturally in a particular colour space. For example, a
scanner may advantageously use the additive colour space RGB (Red,
Green, Blue), whereas a printer may more advantageously use the
subtractive colour space CMY (Cyan, Magenta, Yellow). There are a
large number of other colour spaces in common use: CieLab, CMYK,
YC.sub.bC.sub.r, LUV and LHS are examples. To allow devices using
different colour spaces to interact, a method is needed to convert
data from the colour space of one device to the colour space of
another.
[0003] Such conversions can generally be carried out with a
mathematical formula (though not a trivial one, as such mappings
will typically be non-linear). The most satisfactory approach to
computation in such cases is to use a lookup table. The difficulty
with using a lookup table is that this may need to be extremely
large--for example, a full-table for conversion from a 24 bit RGB
colour space to a 24 bit CMY colour space would require 48 MByte of
memory, which is clearly unacceptably large. The solution, set out
in UK Patent 1369702 and U.S. Pat. No. 3,893,166, is to store in
the table only a coarse (sparse) lattice of points in the colour
space, and to use linear interpolation to compute the values
between these points. In this approach, only values for every 2N
points are stored in each of the input dimensions (plus the last
point in every dimension). Choice of N is important--if low, then
the table will be large, and if high, interpolation will not
approximate with sufficient accuracy. If N=4, the conversion above
is reduced in size to a more acceptable 14.4 KByte.
[0004] FIG. 1 illustrates the problem. Input colour space
components 1 are provided in three dimensions (a, b, c). For the
values in each of these dimensions, the higher bits
a.sub.u,b.sub.u,c.sub.u (from the highest bit to bit N) are used to
find a coarse lattice point 2, and more particularly the cube 3 in
the coarse lattice containing the point of interest. The lowest
bits a.sub.1,b.sub.1,c.sub.1 are then used to determine the output
components 4 with values x,y,z in the output colour space using
interpolation.
[0005] A number of interpolation techniques have been developed.
The best known is trilinear interpolation, which uses a weighted
average of values at the eight vertices of a cube in the sparse
lattice to obtain the output. Trilinear interpolation is shown in
FIG. 2. Values are known for the coarse lattice points 21. In a
first linear interpolation step, a weighted average is taken
between one of the sets of points and another of the sets of points
(at opposite faces of the cube formed by the coarse lattice point
21). This generates four first linear interpolation point 22
(though alternative choices 23, 24 could also have been made). The
second linear interpolation step involves interpolation between
these four points to give two second linear interpolation point 25
(using the earlier alternative choices for first linear
interpolation points, again alternative second linear interpolation
point 26, 27 could be chosen). Whichever route is followed, the
third linear interpolation step between the two second linear
interpolation points provides the third linear interpolation point
28 and hence the output value.
[0006] Trilinear interpolation, and most other known interpolation
processes, require eight table accesses (essentially, one access
for each vertex). An improvement on this is provided by radial
interpolation, described in British Patent Application No.
9826975.6. Radial interpolation has the advantage that it requires
only five table accesses, rather than eight. Various other
interpolation approaches are known, but none of these provides any
significant reduction in computational complexity.
[0007] It is desirable to try and still further improve
interpolation by reducing further computational complexity, and in
particular by reducing as far as possible the total memory
consumption for processes such as colour mapping.
SUMMARY OF INVENTION
[0008] Accordingly, the invention provides a method of determining
a value for a function, comprising: establishing an n-dimensional
lattice, the function having values at the lattice points, and
where n is greater than or equal to two; recording values for a
subset of the lattice points, the lattice points of the subset
being known value lattice points; and establishing a value for a
given lattice point by returning a weighted average of the values
of one or more of (n+1) known value lattice points defining an
n-simplex touching or enclosing the given lattice point.
[0009] This approach is highly advantageous, as it ensures that the
value for any given lattice point can be found with a maximum of
(n+1) look up operations from the table of given value lattice
points (the sparse lattice). If these look up operations have a
critical effect on the performance of the overall system--which may
be the case, particularly where the values for the sparse lattice
are packed, to minimise storage requirements--then this
minimisation of look up operations can have a real effect on speed,
in particular.
[0010] While this approach is valid for n-dimensional spaces, it is
of most interest here for three-dimensional spaces (as would
typically be used in colour mapping). The approach is particularly
effective where n=3, and the known value lattice points define a
tetrahedron. Four is now the maximum number of look up operations
required.
[0011] Whereas four is a maximum number of look up operations
required, for many cases preferred embodiments of the invention can
reduce this number to three, two, or one. Advantageously, a
weighted average of all four known value lattice point values is
used if the given lattice point is enclosed by the tetrahedron but
is not touched by a face of the tetrahedron, a weighted average of
three of the four known value lattice point values is used if the
given lattice point is on a face of the tetrahedron bounded by the
three of the four known value lattice points but is not touched by
an edge of the tetrahedron, a weighted average of two of the four
known value lattice point values is used if the given lattice point
is on an edge of the tetrahedron bounded by the two of the four
known value lattice points but is not at a vertex of the
tetrahedron, and wherein a value of one of the known value lattice
points is used if the given lattice point is also the known value
lattice point.
[0012] In particular embodiments, the average number of table look
up operations can be reduced further. For example, where a given
lattice point is close to a known value lattice point, the given
lattice point may be changed to the known value lattice point hence
avoiding need for calculation of any weighted average, thus
reducing the number of table look up operations required to one
(from, typically, four) for such cases.
[0013] Advantageously, the known value lattice points form a sparse
lattice with known value lattice points separated from each other
by an integer multiple (preferably an integer power of two) of the
distance between adjacent lattice points. A particularly preferred
value for the integer is 8.
[0014] In a preferred approach, the step of establishing a value
comprises determining a set of four known value lattice points
which form a tetrahedron touching or enclosing the given lattice
point, and providing the weighted average from the positions of
four known value lattice points, the known values of one or more of
the four known value lattice points, and the position of the given
lattice point. Advantageously, the step of providing the weighted
average comprises using the positions as inputs to a jump table, or
similar. (Case statements are used in the specific example
described below, use of a computed goto is a further
alternative).
[0015] The method can be carried out by using a programmed
computer, or can be facilitated by providing a programmed memory or
storage device which can be accessed by a computer.
[0016] The method can advantageously be used for mapping values in
a first colour space to values in a second colour space. However,
this is far from the only possible application for this technique.
Look up tables are frequently used when the result is known, but
cannot be computed from a function (either because the computation
is too complex to perform in real time, or because the function
itself is unknown). For example, a forward function may exist that,
given the result as its parameter, will yield the input value, but
the function cannot be used to find the result because the inverse
function cannot be found. One solution to the lack of an inverse
function is simply to enumerate all the parameters for the forward
function and tabulate the results. The table can then be used to
perform the inverse function, but if the table is too large to
store, it can be stored in a sparse manner and intermediate values
calculated by interpolation--this interpolation can be carried out
by methods in accordance with the invention.
[0017] Use of large look up tables, best stored sparsely, (and
hence suitable for application of the invention) may arise in a
large number of contexts. In automotive control, engine management
systems may well require such look up tables. Medical imaging is a
further likely area--particularly for image reconstruction in
tomography (especially soft field electrical tomography). In
computer animation, such an approach is particularly appropriate
for recovering the parameters that represent the motion of objects
in the animation. In image analysis, this approach is particularly
suitable to image enhancement and image recognition. The skilled
man will appreciate that this general approach may be applied in
many other contexts.
DESCRIPTION OF FIGURES
[0018] A specific embodiment of the invention will now be
described, by way of example, with reference to the accompanying
drawings, in which:
[0019] FIG. 1 shows conversion of values in a first colour space to
values in a second colour space, where values in the second colour
space are stored for a sparse lattice of points;
[0020] FIG. 2 shows trilinear interpolation for an intermediate
point in a sparse lattice;
[0021] FIG. 3a shows a sparse lattice of points, and FIG. 3b shows
the full set of lattice points in a cube bounded by adjacent sparse
lattice points;
[0022] FIG. 4 shows use of tetrahedra with vertices of adjacent
sparse lattice points in a method according to an embodiment of the
invention;
[0023] FIG. 5 shows calculation of a value for a lattice point
located within a tetrahedron according to an embodiment of the
invention;
[0024] FIG. 6 shows calculation of a value for a lattice point
located on a face of a tetrahedron according to an embodiment of
the invention; and
[0025] FIG. 7 shows storage required for mapping three 8-bit inputs
to three output bytes using an embodiment of the invention.
DESCRIPTION OF SPECIFIC EMBODIMENTS
[0026] The general principles of an embodiment of the invention for
translating values in one space to values in another space are
described below for a three dimensional example (the general
n-dimensional case is given at the end of this description): a
ready application of this is, as indicated above, to colour
spaces.
[0027] FIG. 3a shows a coarse lattice of points for which known
values are stored. These points are separated by a number of
intermediate points, for which values are not stored, but for which
approximate values can be derived by interpolation. For
computational convenience, the spacing between coarse lattice
points is advantageously 2.sup.N times the intermediate point
spacing. An elemental cube in the coarse lattice and its associated
intermediate points, for n=3, in FIG. 3b: if the intermediate point
spacing is considered to be 1, the coarse lattice point spacing is
8, and a total of 512 values can be interpolated for one coarse
lattice point.
[0028] As can be seen from FIG. 4, each of the intermediate points
in a coarse lattice cube 41 is either within, or on the boundary
of, at least one tetrahedron 42 whose vertices are four of the
vertices of the cube. Considering a tetrahedron with vertices of
coarse lattice points ABCD and an intermediate point P as shown in
FIG. 5, it is possible to work out whether point P lies inside or
outside the tetrahedron ABCD: if
Vol(ABCP)+Vol(ABDP)+Vol(ACDP)+Vol (BCDP)>Vol(ABCD)
[0029] (where the volume between four positions is expressed as
Vol(position 1 position 2 position 3 position 4)), then P lies
outside this tetrahedron. If P lies on or within the tetrahedron,
then the value p at P can be found from the values a, b, c, d at A,
B, C, D by
p=(Vol(ABCP).d+Vol(ABDP).c+Vol(ACDP).b+Vol(BCDP).a)/Vol(ABCD)
[0030] The volume of a tetrahedron with one vertex at the origin
can be expressed by the following formula:
Vol=.vertline.(x.sub.1(y.sub.2z.sub.3-z.sub.2y.sub.3)-y.sub.1(x.sub.2z.sub-
.3-z.sub.2x.sub.3)+z.sub.1(x.sub.2y.sub.3-y.sub.2x.sub.3))/6.vertline.
[0031] These equations can be used to determine whether a
particular point lies on or inside a tetrahedron, and also the
relative influences of each point (the ratio of the tetrahedral
volumes). Some points within the coarse lattice cube will have
several equivalent interpolation equations. For example, the point
in the middle of the cube lies at the midpoint of three diagonal
lines. If the space within the coarse lattice cube is linear, then
all equations should yield similar results. In practice the results
may differ, and the degree of difference is a measure of the
distortion within a particular coarse lattice cube. This can be
useful in determining whether the lattice is too coarse to achieve
an acceptable approximation.
[0032] This method is effective in increasing speed and reducing
computational complexity, as only four table lookups (one to find
the value for each of the four vertices of the relevant
tetrahedron) are used. A further advantage of this approach is that
for points lying on any of the boundaries of the tetrahedron, even
fewer points need to be used.
[0033] If a lattice point lies on the face of a tetrahedron, as
shown in FIG. 6 (where point P lies on the plane ABC), then the
value p can be found from d by weighted average of areas (rather
than volumes as for the general case). The expression for the
interpolated value p reduces to:
p=(Area(ABP)*c+Area(ACP)*b+Area(BCP)*a)/Area(ABC)
[0034] where the area of a triangle can be found by, for example,
Heron's formula:
Area=(s(s-a)(s-b)(s-c)).sup.1/2
[0035] where s=(a+b+c)/2 and a b and c are the lengths of the three
sides.
[0036] For P to lie on the plane defined by ABC and within the area
ABC:
Area(ABP)+Area(ACP)+Area(BCP)=Area(ABC)
[0037] In the case of an intermediate point lying on the face of a
tetrahedron, it is therefore possible to calculate the value of the
intermediate point using only three table lookups, with a
consequent further increase in speed.
[0038] For a point lying on a line between two of the vertices of
the tetrahedron, the calculation reduces still further. The
expression for the interpolated value p now reduces to:
p=(Length(BP)*a+Length(AP)*b)/Length(AB))
[0039] again, for P to lie on the line AB, then
Length(AB)=Length(AP)+Length(BP)
[0040] The distance between two points (x.sub.1,y.sub.1,z.sub.1)
and (x.sub.2, y.sub.2, z.sub.2) is equal to
((x.sub.1-x.sub.2).sup.2+(Y.sub.1-
-y.sub.2).sup.2+(z.sub.1-z.sub.2).sup.2).sup.1/2. For such points,
only two table lookups are now required, as a weighted average of
lines is all that is required. Clearly, if the "intermediate point"
in fact lies at a vertex, the general expression reduces to, for
example:
p=a
[0041] as Vol (BCDP)=Vol (ABCD), and Vol (ABCP)=Vol (ABDP)=Vol
(ACDP)=Vol (APCP)=0. In this case, only one table lookup is
required.
[0042] Use of this approach provides effective reduction of the
number of table accesses required for interpolation, which is
particularly important to maximise speed if the coarse lattice
point values are stored in a packed array (in which case less
storage is required, but retrieval is more difficult). Table 1
illustrates the frequency with which intermediate points fall on
vertices, lines or faces of tetrahedra in a coarse lattice
cube.
1TABLE 1 Location of intermediate points in relation to coarse
lattice tetrahedra Inter Lattice Sub Cube On On Within In N
Distance 2.sup.N Volume (2.sup.N).sup.3 Point Line Triangle
Tetrahedron 0 1 1 1 0 0 0 1 2 8 1 7 0 0 2 4 64 1 33 30 0 3 8 512 1
85 378 48 4 16 4096 1 189 2322 1584 5 32 32768 1 397 11202 21168 6
64 262144 1 813 48918 212418 7 128 2097152 1 4579 188586
1903986
[0043] The corresponding probabilities of requiring one, two or
three table accesses rather the general (and maximum) number of
four is shown in Table 2.
2TABLE 2 Probable number of table accesses Average N P(Point)
P(Line) P(Triangle) P(Tetrahedron) Table Access Best Case Worst
Case 0 1.0000000 1.0000000 1.0000000 1.0000000 1 0.1250000
0.8750000 1.8750000 1.0000000 2.0000000 2 0.0156250 0.5156250
0.4687500 2.4531250 1.0000000 3.0000000 3 0.0019531 0.1660156
0.7382813 0.0937500 2.9238281 1.0000000 4.0000000 4 0.0002441
0.0461426 0.5668945 0.3867188 3.3400879 1.0000000 4.0000000 5
0.0000305 0.0121155 0.3418579 0.6459961 3.6338196 1.0000000
4.0000000 6 0.0000038 0.0031013 0.1866074 0.8103104 3.8072701
1.0000000 4.0000000 7 0.0000005 0.0021834 0.0899248 0.9078913
3.9057069 1.0000000 4.0000000
[0044] For small N, on average significantly fewer than four table
accesses are required. For greater efficiency, points close to a
lattice point could be snapped to a lattice point, rather than
treated as general intermediate points. This is described further
below.
[0045] For each possible point P it is possible to derive the
weights and coarse lattice points to be used. For N=2 there are 64
points in the subcube, which can be addressed as P[x][y][z], where
x, y and z are between 0 and 3 inclusive (points P which would have
a value of x, y or z=4, although these would in fact be touched by
the coarse lattice cube under consideration, are in fact treated in
adjoining coarse lattice cubes in which the points would have a
corresponding x, y or z value of 0). It is convenient to write a
program to generate the equation coefficients for evaluation of p
in each case, but the efficiency of generating the equations is not
important: it is the efficiency of execution of the equations, once
generated, that is important for the efficiency of the
interpolation process. If the eight vertices of the cube (which are
coarse lattice points) are referred to as W, X, Y, Z, XY, XZ, YZ,
XYZ, then the interpolations can be given by the equations shown in
Table 3 for the N=2 case.
3TABLE 3 Interpolation Equations for N = 2 P[0][0][0] = (4W + 2)/4
P[0][0][1] = (3W + 1Z + 2)/4 P[0][0][2] = (2W + 2Z + 2)/4
P[0][0][3] = (1W + 3Z + 2)/4 P[0][1][0] = (3W + 1Y + 2)/4
P[0][1][1] = (3W + 1YZ + 2)/4 P[0][1][2] = (2W + 1Z + 1YZ + 2)/4
P[0][1][3] = (1Y + 3Z + 2)/4 P[0][2][0] = (2W + 2Y + 2)/4
P[0][2][1] = (2W + 1Y +1YZ + 2)/4 P[0][2][2] = (2Y + 2Z + 2)/4
P[0][2][3] = (1Y + 2Z +1YZ + 2)/4 P[0][3][0] = (1W + 3Y + 2)/4
P[0][3][1] = (3Y + 1Z + 2)/4 P[0][3][2] = (2Y + 1Z + 1YZ + 2)/4
P[0][3][3] = (1W + 3YZ + 2)/4 P[1][0][0] = (3W + 1X + 2)/4
P[1][0][1] = (3W + 1XZ + 2)/4 P[1][0][2] = (2W + 1Z +1XZ + 2)/4
P[1][0][3] = (1X + 3Z + 2)/4 P[1][1][0] = (3W + 1XY + 2)/4
P[1][1][1] = (3W + 1XYZ + 2)/4 P[1][1][2] = (1X + 1Y + 2Z + 2)/4
P[1][1][3] = (1XY + 3Z + 2)/4 P[1][2][0] = (2W + 1Y + 1XY + 2)/4
P[1][2][1] = (1X + 2Y + 1Z + 2)/4 P[1][2][2] = (2Y + 1Z + 1XZ +
2)/4 P[1][2][3] = (1XY + 2Z + 1YZ + 2)/4 P[1][3][0] = (1X + 3Y +
2)/4 P[1][3][1] = (3Y + 1XZ + 2)/4 P[1][3][2] = (2Y + 1XZ + 1YZ +
2)/4 P[1][3][3] = (1X + 3YZ +2)/4 P[2][0][0] = (2W + 2X + 2)/4
P[2][0][1] = (2W + 1X + 1XZ + 2)14 P[2][0][2] = (2X + 2Z + 2)/4
P[2][1][3] = (1X + 2Z + 1XZ + 2)/4 P[2][1][0] = (2W + 1X + 1XY +
2)/4 P[2][1][1] = (2X + 1Y + 1Z + 2)/4 P[2][1][2] = (2X + 1Z + 1YZ
+ 2)/4 P[2][1][3] = (1XY + 2Z +1XZ + 2)/4 P[2][2][0] = (2X + 2Y +
2)/4 P[2][2][1] = (2X + 1Y +1YZ + 2)/4 P[2][2][2] = (2XY + 2Z +
2)/4 P[2][2][3] = (1XY + 2Z + 1XYZ + 2)/4 P[2][3][0] = (1X + 2Y +
1XY + 2)/4 P[2][3][1] = (2Y + 1XY + 1XZ + 2)/4 P[2][3][2] = (2XY +
1Z + 1YZ + 2)/4 P[2][3][3] = (1XY + 1XZ + 2YZ + 2)/4 P[3][0][0] =
(1W + 3X + 2)/4 P[3][0][1] = (3X + 1Z + 2)/4 P[3][0][2] = (2X + 1Z
+ 1XZ + 2)/4 P[3][0][3] = (IW + 3X2 + 2)/4 P[3][1][0] = (3X + 1Y +
2)/4 P[3][1][1] = (3X + 1YZ + 2)/4 P[3][1][2] = (2X + 1XZ + 1YZ +
2)/4 P[3][1][3] = (1Y + 3XZ + 2)/4 P[3][2][0] = (2X + 1Y + 1XY +
2)/4 P[3][2][1] = (2X + 1XY + 1YZ + 2)/4 P[3][2][2] = (2XY + 1Z +
1XZ + 2)/4 P[3][2][3] = (1XY + 2XZ + 1YZ + 2)/4 P[3][3][0] = (1W +
3XY + 2)/4 P[3][3][1] = (3XY + 1Z + 2)/4 P[3][3][2] = (2XY + 1XZ +
1YZ + 2)/4 P[3][3][3] = (1W + 3XYZ + 2)/4
[0046] For many of these intermediate points, more than one
equation exists. For example, P[2][2][2] lies in the centre of the
8 vertices, and so also lies on the midpoint of four diagonals. Any
of the following four equations could be employed: 1 P [ 2 ] [ 2 ]
[ 2 ] = ( 2 XY + 2 Z + 2 ) / 4 = ( 2 X + 2 YZ + 2 ) / 4 = ( 2 XYZ +
2 W + 2 ) / 4 = ( 2 Y + 2 XZ + 2 ) / 4
[0047] These four equations are all equivalent as they all access
only two elements, and the overall weighting is the same for each
equation. However, each equation may yield a different result if
the volume bounded by W, X, Y, Z, XY, XZ, YZ, XYZ is not linear:
the greater the non-linearity, the greater the difference between
the results. In some cases, it may be preferable to use one of such
alternative equations over another: for example, for colour tables
it is typically better to average along the neutral axis.
[0048] An efficient method of implementing this form of
interpolation in software is by use of a case statement. This
involves listing the appropriate equation for each case, together
with a mechanism to allow a jump to the relevant case. Cases start
from zero and are consecutive--the overhead of a case statement is
typically about six instructions.
[0049] Specifically, the high order bits of the input are used to
provide a table offset (essentially, to determine which coarse
lattice cube to use), and the lower order bits to choose which case
to execute. Each case corresponds to an intermediate point, and
need only access the values for the specific coarse lattice points
required by the relevant case equation. Each case statement can be
generated in advance by a program using the tetrahedral volume
ratios indicated earlier. The interpolation coefficients for the
lattice points are known integers in the range (1 . . .
(2.sup.N-1)), with the sum of the coefficients being 2.sup.N. The
final step (which can be performed outside the case statement) is
to round (by adding (2.sup.N)/2) and normalise (by dividing by
2.sup.N). For example, mapping three input bytes to one output byte
the table size will be (2.sup.8-N+1).sup.3, and there will be
(2.sup.N).sup.3 cases.
[0050] Generation of case statements in advance makes it relatively
easy for the number of table accesses to be reduced further by
"snapping" a lattice point on to a nearby known value lattice
point. All points within some predetermined degree of proximity to
a known value lattice point (for example, all lattice points
adjacent to the known value lattice points) can be determined and
case statements for these points provided which simply return the
value at the known value lattice point. Table accesses can be
reduced still further, albeit with more complex calculation
required in establishing the case statements (though with no added
complexity in the interpolation itself), by snapping lattice points
not only to known value lattice points, but also to lines and
surfaces, where the lattice point is near an edge or a face of the
tetrahedron. Again, this involves at the time of constructing the
case statements determining which points lie within a proximity
threshold, and providing a case statement returning the function of
two or three known value lattice points appropriate to the
equivalent point on the line or surface respectively. Use of this
"snapping" approach does reduce the accuracy of interpolation for
"snapped" points, but also very significantly decreases the average
number of look up operations required.
[0051] An example of code for N=4, providing a representative
collection of case statements (the remainder of the case statements
can be calculated according to principles illustrated above without
difficulty, but are omitted here for reasons of space) is provided
below.
4 Example: Code for N-4 #define N 4 #define SIZE 16 #define STRIDE
((1<<(8-N))+1) #define ROUND 8 unsigned char table [(STRIDE)
* (STRIDE) * (STRIDE)]; #define W offset [0] /* x = 0, y = 0, z = 0
*/ #define X offset [STRIDE*STRIDE] /* x = 1, y = 0, z = 0 */
#define XY offset [STRIDE*STRIDE +STRIDE] /* x = 1, y = 1, z = 0 */
#define XYZ offset [STRIDE*STRIDE +STRIDE+1] /* x = 1, y = 1, z = 1
*/ #define XE offset [STRIDE*STRIDE +1] /* x = 1, y = 0, z = 1 */
#define Y offset [STRIDE] /* x = 0, y = 1, z = 0 */ #define YZ
offset [STRIDE +1] /* x = 0, y = 1, z = 1 */ #define E offset [1]
/* x = 0, y = 0, z = 1 */ char interpolate(xin,yin,zin) int
xin,yin,zin; { char p; int swvar = (((xin&(SIZE-1)) <<
(2*N)).vertline. ((yin&(SIZXE-1)) << N) .vertline.
(zin&(SIZE-1))); register char *offset = &table[
((xin>>N)*STRIDE*STRIDE)) +((yin>>N)*STRIDE) +
(zin>>N)]; switch( swvar){ case 0 : p = 16*W; break; case 1 :
p = 15*W + Z; break; case 2 : p = 14*W + 2*Z; break; case 3 : p =
13*W + 3*Z; break; . . . case 4074 : p = 6*XY + 2*XZ + YZ + 7*XYZ;
break; case 4075 : p = 5*XY + 2*XZ + YZ + 8*XYZ; break; case 4076 :
p = 4*XY + 2*XZ + YZ + 9*XYZ; break; case 4077 : p = 2*X + Y +
13*XYZ; break; case 4078 : p = 2*X + YZ + 13*XYZ; break; case 4079
: p = Y + 2*XZ + 13*XYZ; break; case 4080 : P = W + 15*XY;; break;
case 4081 : p = 15*XY + Z;; break; case 4082 : p = 14*XY + XZ + YZ;
break; case 4083 : p = 13*XY + Z + 2*XYZ; break; case 4084 : p =
12*XY + Z + 3*XYZ; break; case 4085 : p = 11*XY + Z + 4*XYZ; break;
case 4086 : p = 10*XY + Z + 5*XYZ; break; case 4087 : p = 9*XY + Z
+ 6*XYZ; break; case 4088 : p = 8*XY + Z + 7*XYZ; break; case 4089
: p = 7*XY + Z + 8*XYZ; break; case 4090 : p = 6*XY + Z + 9*XYZ;
break; case 4091 : p = 5*XY + Z + 10*XYZ; break; case 4092 : p =
4*XY + Z + 11*XYZ; break; case 4093 : p = 3*XY + Z + 12*XYZ; break;
case 4094 : p = 2*XY + Z + 13*XYZ; break; case 4095 : p = W +
15*XYZ; break; } p = (p + ROUND)>>N; return p; };
[0052] As N increases, a greater number of case statements are
required as the distance between coarse lattice points increases
(though the table size itself decreases relatively). The task of,
for example, mapping three input bytes to three output bytes
(similar to the above example, but requiring three assignment
statements per line instead of one assignment statement per line as
in the example above) can be achieved for different values of N,
but with different effect--in particular, with different overall
storage requirements for the code and the table. FIG. 7 shows the
storage required for different values of N. The code size is for a
program compiled using cc with no flags on a Hewlett-Packard 735
workstation running HP-UX version 9 as reported by the Unix size
program (and would differ with different compilers or target
instruction sets). The table size will be 3.(2.sup.8-N+1).sup.3
bytes and the code will have (2.sup.N).sup.3 cases.
[0053] As can be seen from FIG. 7, in this case the smallest total
storage (200 kBytes) occurs for N=3. Increasing N above 3 has in
this case no advantage and several disadvantages: the table will
have fewer entries, and so will probably be less precise; the
average number of lookups will increase, so the interpolation will
be slower; and the large number of switch statements may defeat
some compilers. If more storage is available, then reducing N will
generally speed up the interpolation and increase precision.
Moreover, the relative advantage of techniques according to the
invention over other interpolation techniques increases for small
N, as even fewer memory accesses are required. This is likely to
render this technique particularly useful where memory is cheap,
and use of large tables common.
[0054] As indicated above, although described for a 3-dimensional
lattice, methods in accordance with the invention can be used for
n-dimensional spaces (and may in particular be useful for
dimensions greater than 3). The general geometric figure providing
the basis for interpolation is an n-simplex in an n-dimensional
space (a tetrahedron is a 3-simplex, and a triangle is a
2-simplex). In the case of a point P within an n-simplex, where the
vertices of the n-simplex are known value lattice points V0, V1, V2
. . . Vn with values v0, v1, v2 . . . vn, and the point is an
intermediate lattice point requiring interpolation to find its
value p, it can first be considered that (n+1) n-simplexes exist
which have n vertices in common with the known value n-simplex and
which have P as the other vertex. If the volumes of these (n+1)
n-simplexes add up to the volume of the known value n-simplex, P
lies within the known value n-simplex (or on its boundary). The
value of p can then be found by:
p=(v0(volume(P,V1,V2. . . Vn)+v1(volume(V0,P,V2. . . Vn)+v2(volume
(V0,V1,P . . . Vn). . . +vn(volume(V0,V1,V2. . .
P)))/volume(V0,V1,V2,Vn)
[0055] If any of the n-simplexes which have P as a vertex should
have a zero volume, then P must lie in the (n-1-simplex that does
not include P--for example, if volume(V0,P,V2 . . . Vn)=0, then P
lies in or on the (n-1-simplex with vertices V0,V2 . . . Vn. If
this is the case, then only n look up operations are required
rather than n+1--correspondingly fewer look up operations are
required if further ones of these n-simplex volumes are zero.
[0056] The same approach to implementing interpolation--use of
earlier generated case statements--can be applied in any
n-dimensional case. Similarly, it is quite as possible to snap to a
vertex, an edge, a plane, or even an (n-1--simplex in the
n-dimensional case as it is in the specific three-dimensional
case.
* * * * *