We present a hybrid method for spectral reflectivity recovery, using 3D extrapolation as a supplemental method for 3D interpolation. The proposed 3D extrapolation is an extended version of 3D interpolation based on the barycentric algorithm. It is faster and more accurate than the conventional spectralrecovery techniques of principalcomponent analysis and nonnegative matrix transformation. Four different extrapolation techniques (based on nearest neighbors, circumcenters, incenters, and centroids) are formulated and applied to recover spectral reflectivity. Under the standard conditions of a D65 illuminant and 1964 10° observer, all reflectivity data from 1269 Munsell color chips are successfully reconstructed. The superiority of the proposed method is demonstrated using statistical data to compare coefficients of correlation and determination. The proposed hybrid method can be applied for fast and accurate spectral reflectivity recovery in image processing.
I. INTRODUCTION
Recently there has been considerable interest in spectral reflectivity recovery methods in the field of color science
[1

7]
because of their fundamental importance, as well as their technological significance. In general, the color of an object is determined by the illumination conditions and the reflectivity of the surface material
[8

10]
. The light spectrum reflected by the material is passed through the eye and quantified by tristimulus values
[8

10]
. The spectral behavior of a colored object is known to be the fingerprint of the color, which is independent of the applied light source and the observer (for nonfluorescent materials). Spectral reflectivity recovery is also important for multispectral image processing, color consistency, color reconstruction, and color image processing
[1

3]
.
Recovering the spectral reflectivity from known tristimulus values is an inverse problem for a given illumination condition
[1
,
7
,
8]
. In general, the solution to such an inverse problem is illposed, meaning that an exact and accurate solution is not easy to find. The solution also depends on the noise and initial conditions
[11

13]
. Two major methods for spectral reflectivity recovery are the pseudoinverse technique
[2

6
,
14

17]
and the interpolation technique
[7
,
18

20]
. In the pseudoinverse technique, principalcomponent analysis
[2

6
,
14

16]
and nonnegative matrix transformation
[7
,
17
,
21

23]
are developed with various adaptation techniques. These are partially successful at reflectivity recovery, but finding an appropriate adaptation set requires extensive computation time, and the adaptation algorithm itself suffers from a localization problem.
A threedimensional (3D) interpolation technique is often applied to minimize calculations when approximating mathematically defined complex functions, or to produce intermediate results from sparse data points
[11
,
24]
. Both situations can be directly employed when we convert images from one devicedependent color space to another, or invert a color space to reflectivity data
[1
,
7
,
18
,
20]
. Nearest neighbors with given tristimulus values and reflectivity data are needed to obtain the reflectivity of the target tristimulus values. 3D interpolation with Delaunay tetrahedration extracts the four vertices of a tetrahedron and calculates the reflectivity value of an unknown point inside. The 3D interpolation technique is faster, in terms of computational time, than other methods
[1
,
7
,
18
,
20
,
24]
. However, the interpolation technique can only be applied to find the tristimulus values of points that are inside the color gamut of a training data set
[7
,
18
,
20]
. The tristimulus values of a point outside the color gamut (convex hull) of the training set cannot be interpolated.
In our previous study
[7]
we applied a hybrid method for spectral reflectivity recovery, interpolating tristimulus values when the data are inside the color gamut, and using nonnegative matrix transformation (NMT) when the data are outside. This technique has overcome the color gamut problem of the interpolation method and is faster than previous pseudoinverse techniques. However, even with this fast and effective hybrid method, outofgamut data still must be recovered through the timeconsuming pseudoinverse technique of NMT. In this report, we apply a simple extrapolation technique
[11
,
24

27]
to the spectral recovery process to address the color gamut problem of interpolation. We start with a brief introduction of the 3D Delaunay interpolation method and proposed algorithms for the 3D extrapolation of tristimulus values. Next, four different algorithms for obtaining the vertices of a tetrahedron are presented and used for 3D extrapolation. We analyze the statistics of the extrapolation technique in terms of efficiency, accuracy, and color difference. Finally, we present a comparison between this extrapolation method and various pseudoinverse techniques.
II. THEORETICAL BACKGROUND
 2.1. 3D Interpolation
The 3D interpolation technique has been used in various fields to deal with scattered data
[1
,
7
,
11
,
18

20
,
24]
. It has also been applied to spectral reflectivity recovery, and is regarded as a fast and reliable method
[7
,
18
,
20]
. In this section, the concept of the 3D interpolation technique is revisited and summarized. This technique will be extended to a 3D extrapolation technique, explained in Section 2.2.
The reference color data points are usually given by CIE
XYZ
values, and they are scattered in 3D coordinates that are enclosed by a color gamut. The color data points can be triangulated with Delaunay triangulation (3D tessellation). To calculate the reflectivity value of a target CIE
XYZ
value that is inside the color gamut, a tetrahedron can be defined to which the target point belongs. Assuming that the four vertices of the tetrahedron are given by
P_{1}
,
P_{2}
,
P_{3}
, and
P_{4},
the CIE
XYZ
values of the corresponding points are given as (X
_{1}
, Y
_{1}
, Z
_{1}
), (X
_{2}
, Y
_{2}
, Z
_{2}
), (X
_{3}
, Y
_{3}
, Z
_{3}
), and (X
_{4}
, Y
_{4}
, Z
_{4}
), and the reflectivity values of the colors are given by
R_{1}
,
R_{2}
,
R_{3}
, and
R_{4}
respectively. Any target point
P_{t}
inside this tetrahedron is characterized by the CIE
XYZ
value (X
_{t,}
Y
_{t}
, Z
_{t}
) and the reflectivity value
R_{t}
. Since the reflectivity is an additive quantity in the CIE
XYZ
coordinates in this case of an additive mixing problem, the task is reduced to finding the appropriate mixing coefficients for the four vertices. It is straightforward to show that the target point can then be found by solving the following linear interpolation polynomial:
The coefficients
a
,
b
,
c
, and
d
can be found by solving the following linear system of equations:
In a 3D interpolation technique, the reflectivity value
R_{t}
can be calculated using the following equation:
In Eq. (3) the four coefficients
a
,
b
,
c
, and
d
are the actual barycentric coordinates of the point
P_{t}
in the tetrahedron consisting of four vertices
P_{1}
,
P_{2}
,
P_{3}
, and
P_{4}
, and can be calculated by solving the following matrix equation:
In the 3D interpolation algorithm, the four coefficients of color matching range between 0 and 1. Implementation of the 3D interpolation algorithm requires the following five steps:

1. Calculating CIEXYZdata from given reflectivity data and other conditions.

2. 3D Delaunay triangulation of data points.

3. Finding a tetrahedron that contains target CIEXYZdata.

4. Calculating the barycentric coordinate of target CIEXYZdata using Eq. (4) and the four vertices of the tetrahedron.

5. Calculating the reflectivity value of the target point using Eq. (3).
Basically, color matching in additive mixing is a linear and additive problem, as shown by Eqs. (1)(3). The 3D interpolation of reflectivity involves a simple, linear technique. However, there are two points to be considered. First, spectral recovery is an inverse problem with many possible solutions. The spectral value of the target point is determined by the characteristics of neighboring points. Depending on the reference data set, a completely different reflectivity with the same CIE
XYZ
value could be obtained. Therefore, the reference data set is critical for 3D interpolation. The second point to be considered is the socalled “color gamut” problem of interpolation. An interpolation algorithm (including 3D interpolation) can be applied only when the target point is inside the color gamut (convex hull) of the reference data set. Therefore, it is preferable to have a large color gamut for 3D interpolation.
 2.2. 3D Extrapolation
In this section we develop a simple 3D extrapolation
[11
,
24

27]
that can be applied to spectralrecovery problems. We reconsider the interpolation problem by examining some basic equations, Eqs. (1)(4). The barycentric coordinate equation of Eq. (4) is the main equation of 3D interpolation. Its limiting condition is that the four coefficients,
a
,
b
,
c
, and
d
of color mixing should be between 0 and 1 for the target point inside the tetrahedron. A careful reexamination of Eq. (4) provides a clue for 3D extrapolation: Even though the target point may be outside the tetrahedron, the barycentric equation still can be applied, as long as the four coefficients can be any real number
[25

27]
. For example, if
d
is a negative value,
d
can be substituted by 
d
’ (while
d
’ can be positive constant), and Eq. (1) becomes
The corresponding reflectivity matching equation (Eq. (3)) is then
The meanings of Eqs. (5) and (6) can be interpreted as follows. To obtain the reflectivity of a target color
R_{t}
, the reflectivities of
R_{1}
,
R_{2}
, and
R_{3}
should be added in the amounts
a
,
b
, and
c
, and then the reflectivity of
R_{4}
should be subtracted by an amount
d
. This statement can be extended to any positive/negative combination of the four coefficients
a
,
b
,
c
, and
d
. Physically, this means that a color can be additive and subtractive in terms of reflectivity in the problem of additive color mixing.
It is essential for implementation of this extrapolation algorithm to find an appropriate tetrahedron for the calculation of the barycentric coordinates. One still can utilize the strategy of 3D interpolation. We tested four different algorithms to define tetrahedra for extrapolating color data. The first choice in defining the tetrahedron is finding the four nearest neighbors of the target color point. Then the Euclidian distance between the target point and reference can be calculated, and a neighborsearching algorithm can be applied to define four points from which to construct a tetrahedron.
The next choice in defining the tetrahedron is the utilization of information from the 3D interpolation. In our 3D interpolation we used the Delaunay tessellation algorithm in MATLAB
^{®}
(version 7.11, MathWorks)
[24]
. In mathematics and computational geometry,
[11
,
24]
a Delaunay triangulation for a set P of points in a plane is a triangulation DT(P) such that no point in P is inside the circumcircle of any triangle in DT(P). Delaunay triangulation maximizes the narrowest angle of all angles of the triangles in the triangulation. Therefore, the circumcenters of all tetrahedra have already been calculated by 3D interpolation. We used the circumcenters of partitioned tetrahedra and calculated the Euclidian distances between the circumcenters and the target point to find an appropriate tetrahedron of neighbors. One can also calculate and tabulate incenters and centroids of tetrahedrons by 3D Delaunay tessellation. We applied three methods (circumcenters, incenters, and centroids) with given tessellation information to define a tetrahedron, and the neighborsearching algorithm was applied to find an appropriate tetrahedron. After finding an appropriate tetrahedron, the four coefficients of mixing can be calculated easily by modifying the algorithm for 3D interpolation.
III. EXPERIMENT
The reflectivity spectra of 1269 color chips from the Munsell color matte finish collection book were downloaded from a website
[28]
and used as a reference data set. The original reflectivity data were measured using a spectrophotometer (Lambda 18, Perkin Elmer), the wavelengths of the data ranging from 380 to 800 nm with intervals of 1 nm. In this work the reflectivity data from 400 to 700 nm with intervals of 10 nm were used to evaluate the proposed reconstruction algorithm. A spectral reflectivity recovery algorithm was developed, which includes the process of rebuilding the reflectivity through a hybrid method using 3D interpolation and 3D extrapolation. To calculate a CIE
XYZ
value we assumed the CIE standard illuminant D65 and 1964 standard observer with a 10° field
[1
,
8
,
9]
. To evaluate the algorithm, the target reflectivity spectrum was removed from the 1269 chip set prior to calculation, and the remaining set of 1268 data points was used. All calculations were carried out using MATLAB
^{®}
(version 7.11, MathWorks)
[24]
. Using 1268 reflectivity data, a 3D Delaunay triangulation was performed. The 3D interpolation method was applied to reconstruct the spectral reflectivity if the target was inside the 3D convex hull of the reference data, or the 3D extrapolation method was applied according to the process if the target was outside the 3D convex hull. This procedure was repeated for 1269 reflectivity data items to obtain statistical significance.
To evaluate the proposed algorithm, the mean and maximum rootmeansquare errors (RMSE) between the original and reconstructed reflectivity spectra were evaluated. We have examined the significance of the reconstructed reflectivity data using the coefficient of correlation (COC)
[12
,
29]
as follows:
where
R_{j}
and
R^{r} _{j}
are the
j
th wavelength components of the original and reconstructed spectra respectively. We also calculated the coefficient of determination
[12
,
29]
. The color difference between the target value and reconstructed value was also compared using CIE LAB
ΔE
, which defines Euclidian distance in CIE LAB coordinates
[30
,
31]
. One should note that due to the characteristics of the interpolation and extrapolation, the color difference is nearly zero (~10
^{11}
).
IV. RESULTS AND DISCUSSION
Figure 1(a)
shows the color gamut and 3D Delaunay tetrahedrization results of 1269 Munsell color chips in CIE
XYZ
space. The color recovery procedures are shown in
Figs. 1(b)
and
(c)
. For example, Munsell 7.5Y 7/8 color (color number 357, denoted as a red dot in
Fig. 1(b)
) was used for the reconstruction algorithm. The coordinates for this color are (32.2774, 34.1852, 9.4060) with standard illuminant D65 and a 1964 10° observer
[1
,
8

10
,
30
,
31]
. The four vertices of the tetrahedron that encloses this color are 10R 6/12, 5Y 7/8, 7.5Y 8/10, and 10Y 8/8, as shown by black dots in
Fig. 1(b)
. The reflectivity values of the tetrahedral vertices are shown as black curves in
Fig. 1(c)
. From their CIE
XYZ
coordinate values, the four colormixing coefficients were calculated to be 0.0715, 0.00558, 0.1292, and 0.7445. Then the reconstructed reflectivity value was calculated and is shown as a red diamond in
Fig. 1(c)
. The original reflectivity value of the target color is shown as a red curve in
Fig. 1(c)
.
(a) CIE XYZ color gamut of 1269 Munsell color data with D65 illuminant and 1964 10° observer. (b) 3D interpolation of a target point using 3D Delaunay tessellation and a barycentric algorithm. (c) Recovered reflectivity of the target point (red diamonds) by barycentric transformation of four reflectivity data points. (d) Coefficients of correlation of 1269 Munsell color data.
This recovery procedure was repeated for all 1269 Munsell colors, and coefficients of correlation are depicted in
Fig. 1(d)
. It is extremely difficult to reconstruct some of the reflectivity data because they belong to the border between different color chips, for which interpolation will have considerable uncertainty
[7
,
18
,
20]
. As previously discussed, 115 data points define the convex hull of the 3D CIE
XYZ
data of 1269 Munsell colors. The reflectivity recovery for these 115 data points cannot be performed because of the 3D interpolation’s limitation, and they are marked as having a zero value in
Fig. 1(d)
.
ΔE
is calculated from the difference in CIE LAB coordinates
, and is summarized in
Table 1
with other statistical measures for 3D interpolation.
ΔE
is less than 10
^{−11}
, which is quite small and originates from numerical errors during calculation.
Statistical summary of various reflectivity reconstruction methods applying 1269 Munsell color data
RMSE is the rootmeansquare error, and Coefficient of correlation is a measure of the strength and direction of the linear relationship between the original and reconstructed spectral reflectivity data. ΔEis the average color difference between the target and reconstructed data. The relative time is the time normalized to the 3D reconstruction method. The method of Simone and the adaptive PCA take nearly same amount of time as NMT, while the proposed algorithm takes nearly the same time as 3D interpolation. Part of the statistical data is the same as in reference 7, because we use the same standard algorithm to calculate statistics.
When the target data points are outside the gamut of the reference reflectivity data, fast and accurate 3D interpolation cannot be applied. Various techniques can be applied in such cases: different versions of principalcomponent analysis (PCA) and nonnegative matrix transformation (NMT) have been applied to the spectral recovery of “outside color gamut data”
[2

7
,
14

20]
. In
Table 1
we report the statistical values for various PCA and NMT techniques, as well as for our interpolation technique. The interpolation for insidegamut data can be performed in 2D CIE
xyY
coordinates and 3D CIE
XYZ
coordinates. For the PCA technique, we have applied two different algorithms: one from general adaptive PCA
[5]
and another from the method of Simone
et. al
[4]
. The twostep adaptive PCA method of Simone
et. al
is known to yield better ΔE than the conventional adaptive version of PCA. We also checked the adaptive version of the NMT technique, which has nearly the same statistical values as the adaptive version of PCA.
As described in Section 2, we implemented four different techniques for 3D extrapolation. Our basic algorithm for spectral recovery is shown in
Fig. 2(a)
. The first implementation found the four nearest neighbors of the target data (shown in
Fig. 2(b)
), and we calculated the Euclidian distance between the target and neighboring data
[32]
and selected the four nearest data points for the spectral reconstruction. The next class of implementation involved using various properties of the tetrahedron that were already partitioned by 3D Delaunay tetrahedrization.
Figures 2(c)

(e)
show the centroid, incenter, and circumcenter of the tetrahedron respectively. The centroid of the tetrahedron is the arithmetic mean of its four vertices. The circumcenter and incenter are the centers of the circumsphere and insphere, respectively, which can be easily found using the builtin functions “circumcenters” and “incenters” of MATLAB
[11
,
13
,
24]
. These points can be calculated and tabulated when the reflectivity data are given, so the additional computational time required for 3D extrapolation using these points can be minimized. It is worthwhile to mention the following points: First, the centroid and the incenter of a tetrahedron always lie within the tetrahedron, but the circumcenter may not. Second, a Delaunay tetrahedrization algorithm is based on circumcenters, the importance of which will be demonstrated later.
(a) Details of the hybrid method algorithm using 3D interpolation and 3D extrapolation. (b) Four nearest neighbors of target tristimulus values in the 3D extrapolation algorithm. (c) Centroid, (d) incenter, and (e) circumcenter of a tetrahedron used in 3D extrapolation.
Having developed the 3D extrapolation techniques, we now describe the details and statistics of these techniques. The results of extrapolating nearest neighbors are shown in
Fig. 3
.
Figures 3(a)
and
(b)
show the spectral reflectivity recovery of the Munsell color 7.5YR 8/10 (#225), the circle in
Fig. 3(a)
. This color belongs to the convex hull made by the 1269 Munsell reference colors and cannot be recovered by 3D interpolation. The four nearest neighbors of 7.5YR 8/10 are 10YR 5/8, 10YR 5/6, 2.5Y 8.5/10, and 2.5Y 8.5/8, plotted as a tetrahedron in
Fig. 3(a)
, and their mixing coefficients in the 3D extrapolation are 0.8870, 1.1033, 0.0467, and 0.9435 respectively. Four reflectivity spectra are plotted as solid curves in
Fig. 3(b)
, and the recovered spectrum is plotted as diamonds, while the original spectrum of 7.5YR 8/10 is depicted as a dotted curve. The spectral recovery of 7.5YR 8/10 is excellent.
Figures 3(c)
and
(d)
show the recovery of the Munsell color 2.5RP 5/12 (#1163). The four nearest neighbors are identified as 2.5RP 5/10, 10P 5/10, 10P 5/12, and 2.5RP 5/8, and they are marked as a tetrahedron; their mixing coefficients are 3.1520, 1.9981, −2.0704, and −2.0798 respectively. The recovered spectrum (diamonds in
Fig. 3(d)
) and original spectrum (dotted curve) do not match well. This disagreement could be a consequence of the tetrahedron being too small in the CIE
Y
coordinates, as shown in
Fig. 3(c)
. The absolute values of the mixing coefficients are greater than 1. Disagreement can also be seen in the original reflectance versus recovered reflectance plot shown in
Fig. 3(e)
. Good agreement data (square symbols) from
Figs. 3(a)
and
(b)
are shown as a straight line in this correlation plot, whereas poor agreement data (diamonds) deviate from a straight line. The coefficients of correlation are plotted in
Fig. 3(f)
for the hybrid method with 3D interpolation (gray bar), and 3D extrapolation using the nearest neighbor (black bar).
(a) An example of 3D extrapolation using a nearestneighbor algorithm. The red dot represents the color for which we want to recover reflectivity in the CIE XYZ coordinate system, and the BLACK squares represent 10 nearestneighbor points for which we know the reflectivity. Four nearestneighbor points around the red dot form a tetrahedron (blue). The (b) reflectivity recovery of Munsell color 7.5YR 8/10 (#225) using a nearestneighbor extrapolation algorithm. Red dots indicate the recovered reflectivity of the target color, and black lines the reflectivities of the four nearest neighbors. (c) Tristimulus representation and (d) reflectivity recovery of 3D extrapolation of Munsell color 2.5RP 5/12 (#1163). (e) Recovered reflectivity versus target reflectivity diagram for Munsell colors 7.5YR 8/10 (#225) (red squares) and 2.5RP 5/12 (#1163) (black diamonds). (f) Histogram plot of absolute value of coefficients of correlation in 3D extrapolation using a nearestneighbor (NN) algorithm. The gray bars are the histogram for the 3D interpolation algorithm, and the black bars that of 3D extrapolation using a NN algorithm.
Information from 3D Delaunay tetrahedrization can be used to obtain better recovery results. We used centroids, incenters, and circumcenters of tetrahedra from 3D Delaunay tetrahedrization, and the results are summarized in
Fig. 4
.
Figures 4(a)
and
(b)
depict the spectral reflectivity reconstruction of the Munsell color 2.5RP 5/12 (##1163) (the same color shown in
Figs. 3(c)
and
(d)
) using the circumcenter algorithm. The circumcenters of all tetrahedra were calculated using the builtin MATLAB function, and the Euclidean distances between circumcenters and the target point were computed. The shortestdistance selection criterion for the tetrahedron is shown by solid lines. The tetrahedron with dashed lines is from the nearestneighbor point (
Fig. 3(c)
versus
Fig. 4(a)
). The Munsell notations of the vertices of the solid tetrahedron are 2.5RP 5/8, 2.5RP 5/10, 10P 5/6, and 10P 5/8. Two of them (2.5RP 5/8 and 2.5RP 5/10) are the same points selected by the nearestneighbor algorithm. The mixing coefficients for the four vertices are −0.5370, 1.5980, −0.3030, and 0.2419, smaller than those of the nearestneighbor selection algorithm (3.1520, 1.9981, −2.0704, and −2.0798). The spectralrecovery result is shown in
Fig. 4(b)
. The reflectances for the four vertices are shown as solid curves. The target reflectance and the recovered reflectance are shown as dotted and diamond curves respectively. Compared to those from previous nearestneighbor algorithms (shown in
Fig. 3(d)
), the recovered spectrum matches the target reflectance fairly well.
(a) Example of 3D extrapolation of Munsell color 2.5RP 5/12 (#1163) using a circumcenter algorithm. The red dot represents the color for which we want to recover reflectivity in CIE XYZ coordinates, and the black squares represent 10 nearestneighbor points for which we know the reflectivity. The tetrahedron with dashed lines is from 3D nearestneighbor extrapolation, while that with solid lines is from 3D circumcenter extrapolation. (b) Reflectivity recovery of 3D extrapolation of Munsell color 2.5RP 5/12 (#1163) using 3D circumcenter extrapolation. The recovered reflectivity is shown with black diamonds, and the target reflectivity is shown as a red dashed line. (c) Histogram plot of coefficients of correlation for 3D extrapolation using a circumcenter (CC) algorithm. Gray bars are the histogram for the 3D interpolation algorithm for comparison, and black bars are that for 3D extrapolation using CC algorithm. (d) Mean values with error bars for coefficients of correlation (blue) and determination (red) for six different recovery algorithms. “3DIN” represents 3D interpolation and “Hybrid” represents the hybrid method with 3D interpolation and a nonnegative matrix transformation. NN, CC, IC, and CE represent the hybrid technique of 3D interpolation plus 3D extrapolation using nearest neighbors, circumcenter, incenter, and centroid respectively.
The statistics for circumcenter extrapolation hybridized with 3D interpolation are shown in
Fig. 4(c)
. The coefficient of correlation for 3D interpolation is shown as a gray bar, and the coefficient for the circumcenter extrapolation is shown as a black bar. We note that there is no reflectance with a zero coefficient of correlation, which appears in the 3D interpolation algorithm due to the outofgamut problem. The statistics for the circumcenter extrapolation are better than those for the nearestneighbor algorithm. The mean values and variations of the coefficient of correlation for different methods are summarized in
Fig. 4(d)
using square symbols. “3DIN” represents 3D interpolation, for which the statistics are calculated without outofgamut data (115 data points). “Hybrid” represents the hybrid technique with 3D interpolation and nonnegative matrix transformation (NMT)
[7]
. NN, CC, IC, and CE represent nearestneighbor, circumcenters, incenters, and centroid extrapolation, respectively. Another statistical value, the coefficient of determination, is computed for the different methods and is shown with diamond symbols. Among the four extrapolation algorithms, the CC, IC, and CE methods are better than the NN method because of 3D Delaunay tetrahedrization. The CC, IC, and CE methods show nearly comparable results com pared to the hybrid technique of 3D interpolation plus NMT for spectral recovery. The statistical analysis data are also derived, and summarized in
Table 1
. We note that the outofthegamut problem of 3D interpolation can be overcome by our 3D extrapolation technique; however, the overall statistics for the outofgamut data are not as good as in the case of the 3D interpolation for insidegamut data. CPU times for reconstruction are also summarized in
Table 1
. The relative time is the time normalized to the 3D interpolation method. The method of Simone and adaptive PCA take nearly the same amount of time as NMT, and the proposed algorithm takes nearly the same as 3D interpolation.
We now consider the merits of 3D extrapolation compared to NMT and PCA. First, 3D extrapolation with CC, IC, and CE uses 3D Delaunay tetrahedrization, which is well known for its superior efficiency. In addition, a fast and accurate mathematical algorithm is used with the 3D Delaunay method. Second, the time to compute the 3D extrapolation is minimal compared to the execution time of the adaptive versions of PCA and NMT. Since the CC, IC, and CE of a tetrahedron are computed and stored as a database or lookup table (LUT), the time required for 3D extrapolation is hardly greater than the time required for 3D interpolation, which is at least two orders of magnitude faster than the other reconstruction methods
[4
,
6
,
7]
. Finally, the color gamut of 3D interpolation can be extended by accumulating more reflectivity data for the LUT. Also, the measure of fitting accuracy can be enhanced by the quality and the quantity of the data in the LUT, because most of the errors of 3D interpolation and extrapolation techniques come from color samples near the border between two colors with different character.
V. CONCLUSION
The recovery of spectral data is important because the spectral reflectivity of a color is a fundamental index which is unchangeable, even if the appearance of the color is affected by the ambient environment or the observer's vision. We propose a 3D extrapolation technique with 3D interpolation based on the Delaunay tetrahedrization algorithm for reflectivity reconstruction from the tristimulus values of colors. A previous 3D interpolation technique using a barycentric algorithm is extended to 3D extrapolation using four points among the nearest neighbors. In the proposed procedure, a 3D extrapolation technique is used to find the tristimulus values of color data that exist outside the color gamut. We applied and analyzed four distinct nearestneighbor searching algorithms using nearest points, and the tristimulus values of the CC, IC, and CE of a tetrahedron. Statistical analysis, the coefficient of correlation, and the coefficient of determination show that the results of the CC, IC, and CE methods yield accuracy as high as that of the pseudoinverse technique. However, the extrapolation technique is much faster than the pseudoinverse technique in terms of processing speed. We believe that the proposed 3D extrapolation technique for spectral reflectivity reconstruction could advance color research and widen applications, including color consistency, color reconstruction, and color image processing.
Acknowledgements
Bog G. Kim was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF2013R1 A1A2004496). Soobeen Park was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF2013R1A1A2013152).
Westland S.
,
Ripamonti C.
2004
Computational Color Science Using Matlab
Wiley
New York, USA
Fairman H. S.
,
Brill M. H.
2004
“The principal components of reflectance,”
Color Res. Appl.
29
104 
110
DOI : 10.1002/col.10230
Zuffi S.
,
Santini S.
,
Schettini R.
2008
“From color sensor space to feasible reflectance spectra,”
IEEE Trans. Signal Process.
56
518 
531
DOI : 10.1109/TSP.2007.907838
Bianco S.
2010
“Reflectance spectra recovery from tristimulus values by adaptive estimation with metameric shape correction,”
J. Opt. Soc. Am. A
27
1868 
1877
DOI : 10.1364/JOSAA.27.001868
Babaei V.
,
Amirshahi S. H.
,
Agahian F.
2011
“Using weighted pseudoinverse method for reconstruction of reflectance spectra and analyzing the dataset in terms of normality,”
Color Res. Appl.
36
295 
305
DOI : 10.1002/col.20613
Kandi S. G.
,
Tehran M. A.
2011
“Applying metamer sets to investigate data dependency of principal component analysis method in recovery of spectral data,”
Color Res. Appl.
36
349 
354
DOI : 10.1002/col.20645
Kim B. G.
,
Han J. W.
,
Park S. B.
2012
“Spectral reflectivity recovery from the tristimulus values using a Hybrid method,”
J. Opt. Soc. Am. A
29
2612 
2621
DOI : 10.1364/JOSAA.29.002612
Berns R. S.
2000
Billmeyer and Saltzman's Principles of Color Technology
3rd ed.
Wiley
New York, USA
Schanda J.
2007
Colorimetry: Understanding the CIE System
1st ed.
Wiley
New York, USA
Ohta N.
,
Robertson A.
2006
Colorimetry: Fundamentals and Applications
1st ed.
Wiley
New York, USA
Press W. H.
,
Teukolsky S. A.
,
Vetterling W. T.
,
Flannery B. P.
2007
Numerical Recipes: The Art of Scientific Computing
3rd ed.
Cambridge University Press
New York, USA
Martinez W. L.
,
Martinez A. R.
,
Solka J. L.
2011
Exploratory Data Analysis with Matlab
2nd ed.
Chapman ＆ Hall
Boca Raton, FL, USA
Magrab E. B.
,
Azarm S.
,
Balachandran B.
,
Duncan J. H.
,
Herold K. E.
,
Walsh G. C.
2011
An Engineer’s Guide to Matlab
Pearson
Upper Saddle River, NJ, USA
Agahian F.
,
Amirshahi S. A.
,
Amirshahi S. H.
2008
“Reconstruction of reflectance spectra using weighted principal component analysis,”
Color Res. Appl.
33
360 
371
DOI : 10.1002/col.20431
Harifi T.
,
Amirshahi S. H.
,
Agahian F.
2008
“Recovery of reflectance spectra from colorimetric data using principal component analysis embedded regression technique,”
Opt. Rev.
15
302 
308
DOI : 10.1007/s1004300800491
Farajikhah S.
,
Madanchi F.
,
Amirshahi S. H.
2011
“Nonlinear principal component analysis for compression of spectral data,”
Proc. Conference on Data Mining and Data warehouses
Ljubljana, Slovenia
http://ailab.ijs.si/dunja/SiKDD2011/Papers/Farajikhah_PCA.pdf
Amirshahi S. H.
,
Amirhahi S. A.
2010
“Adaptive nonnegative bases for reconstruction of spectral data from colorimetric information,”
Opt. Review
17
562 
569
DOI : 10.1007/s1004301001019
Kasson J.
,
Plouffe W.
,
Nin S.
1993
“A tetrahedral interpolation technique for color space conversion,”
Proc. SPIE
1909
127 
138
Amidror I.
2002
“Scattered data interpolation methods for electronic imaging systems: A survey,”
J. Electron. Imaging
11
157 
176
DOI : 10.1117/1.1455013
Abed F. M.
,
Amirshahi S. H.
,
Abed M. R. M.
2009
“Reconstruction of reflectance data using an interpolation technique,”
J. Opt. Soc. Am. A
26
613 
624
DOI : 10.1364/JOSAA.26.000613
Lee D. D.
,
Seung H. S.
1999
“Learning the parts of objects by nonnegative matrix factorization,”
Nature
401
788 
791
DOI : 10.1038/44565
Lee D. D.
,
Seung H. S.
2000
“Algorithms for nonnegative matrix factorization,”
MIT Press
Cambridge, MA, USA
Proc. Neural Information Processing Systems
556 
562
Kim J.
,
Park H.
2011
“Fast nonnegative matrix factorization: An activesetlike method and comparisons,”
SIAM J. of Sci. Comput. (SISC)
33
3261 
3281
DOI : 10.1137/110821172
de Berg M.
,
Cheong O.
,
van Kreveld M.
,
Overmars M.
2010
Computational Geometry: Algorithms and Applications
3rd ed
Springer
Berlin, Germany
Wu Q.
,
Xu H.
,
Zou X.
2005
“An effective method for 3D geological modeling with multi source data integration,”
Computers and Geosciences
31
35 
43
DOI : 10.1016/j.cageo.2004.09.005
Bobach T.
,
Farin G.
,
Hansford D.
,
Umlauf G.
2009
“Natural neighbor extrapolation using ghost points,”
ComputerAided Design
41
350 
365
DOI : 10.1016/j.cad.2008.08.007
0000
“Spectral database,”
University of Joensuu Color Group
http://www.uef.fi/fi/spectral/spectraldatabase
0000
We have used MATLAB function of gfit2.m
http://www.mathworks.com/matlabcentral/fileexchange/22020goodnessoffitmodified/content/gfit2.m
Wagberg J.
2007
“OptProp, Matlab toolbox for calculation of color related optical properties,”
http://apachepersonal.miun.se/~magneu/publications/R0777.pdf
Kim B. G.
,
Han J.W.
,
Park S.B.
2013
“Computer simulation for gradual yellowing of aged lens and its application for test devices,”
Journal of the Optical Society of Korea
17
(4)
344 
349
DOI : 10.3807/JOSK.2013.17.4.344
0000
We have used MATLAB function of nearestneighbor.m
http://www.mathworks.com/matlabcentral/fileexchange/12574nearestneighbourm