A units of red {u1-uA}, B units of green {v1-vB} and C units of blue {w1-wC} are sampled and measurements are performed, and an interpolation equation is found from the A×B×C units of data in three-dimensions (L*, a*, b*) that are obtained. Firstly, A units of point sequence data in 3×B×C dimensions are least squares approximated for the parameters u, and Mu control points in 3×B×C dimensions are found. Then, B units of point sequence data in 3×C×Mu dimensions are least squares approximated for the parameters v, and Mv control points in 3×C×Mu dimensions are found. Next, C units of point sequence data in 3×Mu×Mv dimensions are least squares approximated for the parameters w, and Mw control points in 3×Mu×Mv dimensions are found. Then, output values (L*,a*,b*) corresponding to arbitrary input values (Red,Green,Blue) can be obtained by an interpolation equation according to Mu×Mv×Mw final control points in three dimensions and basis functions used in the least squares approximations.