Spectral Reflectivity Recovery from Tristimulus Values Using 3D Extrapolation with 3D Interpolation
Spectral Reflectivity Recovery from Tristimulus Values Using 3D Extrapolation with 3D Interpolation
Journal of the Optical Society of Korea. 2014. Oct, 18(5): 507-516
Copyright © 2014, Journal of the Optical Society of Korea
  • Received : June 05, 2014
  • Accepted : September 05, 2014
  • Published : October 25, 2014
Export by style
Cited by
About the Authors
Bog G. Kim
Department of Physics, Pusan National University, Busan 609-735, Korea
John S. Werner
Department of Ophthalmology and Vision Science, University of California, Davis, CA 95817, USA
Michael Siminovitch
Design Department, California Lighting Technology Center, University of California, Davis, CA 95618, USA
Kostantinos Papamichael
Design Department, California Lighting Technology Center, University of California, Davis, CA 95618, USA
Jeongwon Han
Department of Housing and Interior design, Pusan National University, Busan 609-735, Korea
Soobeen Park
Department of Housing and Interior design, Pusan National University, Busan 609-735, Korea

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 spectral-recovery techniques of principal-component analysis and nonnegative matrix transformation. Four different extrapolation techniques (based on nearest neighbors, circumcenters, in-centers, 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.
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 ill-posed, 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 pseudo-inverse technique [2 - 6 , 14 - 17] and the interpolation technique [7 , 18 - 20] . In the pseudo-inverse technique, principal-component 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 three-dimensional (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 device-dependent 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 pseudo-inverse techniques. However, even with this fast and effective hybrid method, out-of-gamut data still must be recovered through the time-consuming 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 pseudo-inverse techniques.
- 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 P1 , P2 , P3 , and P4, 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 R1 , R2 , R3 , and R4 respectively. Any target point Pt inside this tetrahedron is characterized by the CIE XYZ value (X t, Y t , Z t ) and the reflectivity value Rt . 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:
PPT Slide
Lager Image
The coefficients a , b , c , and d can be found by solving the following linear system of equations:
PPT Slide
Lager Image
In a 3D interpolation technique, the reflectivity value Rt can be calculated using the following equation:
PPT Slide
Lager Image
In Eq. (3) the four coefficients a , b , c , and d are the actual barycentric coordinates of the point Pt in the tetrahedron consisting of four vertices P1 , P2 , P3 , and P4 , and can be calculated by solving the following matrix equation:
PPT Slide
Lager Image
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 so-called “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 spectral-recovery 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
PPT Slide
Lager Image
The corresponding reflectivity matching equation (Eq. (3)) is then
PPT Slide
Lager Image
The meanings of Eqs. (5) and (6) can be interpreted as follows. To obtain the reflectivity of a target color Rt , the reflectivities of R1 , R2 , and R3 should be added in the amounts a , b , and c , and then the reflectivity of R4 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 neighbor-searching 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 in-centers and centroids of tetrahedrons by 3D Delaunay tessellation. We applied three methods (circumcenters, in-centers, and centroids) with given tessellation information to define a tetrahedron, and the neighbor-searching 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.
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 root-mean-square 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:
PPT Slide
Lager Image
where Rj and Rr 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 ).
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 color-mixing 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) .
PPT Slide
Lager Image
(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
PPT Slide
Lager Image
, 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
PPT Slide
Lager Image
RMSE is the root-mean-square 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 principal-component 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 inside-gamut 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 two-step 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, in-center, and circumcenter of the tetrahedron respectively. The centroid of the tetrahedron is the arithmetic mean of its four vertices. The circumcenter and in-center are the centers of the circumsphere and in-sphere, respectively, which can be easily found using the built-in functions “circumcenters” and “in-centers” 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 in-center 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.
PPT Slide
Lager Image
(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) in-center, 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).
PPT Slide
Lager Image
(a) An example of 3D extrapolation using a nearest-neighbor 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 nearest-neighbor points for which we know the reflectivity. Four nearest-neighbor points around the red dot form a tetrahedron (blue). The (b) reflectivity recovery of Munsell color 7.5YR 8/10 (#225) using a nearest-neighbor 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 nearest-neighbor (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, in-centers, 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 built-in MATLAB function, and the Euclidean distances between circumcenters and the target point were computed. The shortest-distance selection criterion for the tetrahedron is shown by solid lines. The tetrahedron with dashed lines is from the nearest-neighbor 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 nearest-neighbor 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 spectral-recovery 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 nearest-neighbor algorithms (shown in Fig. 3(d) ), the recovered spectrum matches the target reflectance fairly well.
PPT Slide
Lager Image
(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 nearest-neighbor points for which we know the reflectivity. The tetrahedron with dashed lines is from 3D nearest-neighbor 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, in-center, 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 out-of-gamut problem. The statistics for the circumcenter extrapolation are better than those for the nearest-neighbor 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 out-of-gamut data (115 data points). “Hybrid” represents the hybrid technique with 3D interpolation and nonnegative matrix transformation (NMT) [7] . NN, CC, IC, and CE represent nearest-neighbor, circumcenters, in-centers, 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 out-of-the-gamut problem of 3D interpolation can be overcome by our 3D extrapolation technique; however, the overall statistics for the out-of-gamut data are not as good as in the case of the 3D interpolation for inside-gamut 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 look-up 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.
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 nearest-neighbor 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 pseudo-inverse 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.
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 (NRF-2013R1 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 (NRF-2013R1A1A2013152).
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 pseudo-inverse 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/s10043-008-0049-1
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
Amirshahi S. H. , Amirhahi S. A. 2010 “Adaptive non-negative bases for reconstruction of spectral data from colorimetric information,” Opt. Review 17 562 - 569    DOI : 10.1007/s10043-010-0101-9
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 non-negative matrix factorization,” Nature 401 788 - 791    DOI : 10.1038/44565
Lee D. D. , Seung H. S. 2000 “Algorithms for non-negative matrix factorization,” MIT Press Cambridge, MA, USA Proc. Neural Information Processing Systems 556 - 562
Kim J. , Park H. 2011 “Fast nonnegative matrix factorization: An active-set-like 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
Hiyoshi H. 2008 “Stable computation of natural neighbor interpolation,” Int. J. Comput. Geom. Appl. 18 321 -    DOI : 10.1142/S0218195908002659
Bobach T. , Farin G. , Hansford D. , Umlauf G. 2009 “Natural neighbor extrapolation using ghost points,” Computer-Aided Design 41 350 - 365    DOI : 10.1016/j.cad.2008.08.007
0000 “Spectral database,” University of Joensuu Color Group
0000 We have used MATLAB function of gfit2.m
Wagberg J. 2007 “OptProp, Matlab toolbox for calculation of color related optical properties,”
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