Abstract:
A fast method of color interpolation of pixels of an image acquired by a color filtered digital sensor uses a very simple cost function that nevertheless produce interpolated images of good quality. The cost function is computationally simpler because it does not require the calculation of powers and square roots. The triangulation algorithm may be executed in far less time, while practically ensuring the same performance. The triangulation algorithm on average may use only two iteration steps. The interpolation process may be followed by an anti-aliasing processing that effectively removes color artifacts.

Description:
FIELD OF THE INVENTION 
   The invention relates in general to processing techniques of image data acquired by a color digital sensor, and, more in particular, to a method of interpolating the color information for each image pixel. 
   BACKGROUND OF THE INVENTION 
   Modern digital cameras use an electronic sensor for acquiring the spatial variations of light intensity. The two main fabrication technologies of sensors are the so called CCD (Charge Coupled Device) and the CMOS (Complementary Metal Oxide Semiconductor). 
   For example, a CCD is provided by an array of photosensitive diodes that converts photons into electrons (electric charge). Each diode, called a photodiode produces an electric charge proportional to the quantity of light reaching it. 
   The single sensors of a digital camera are normally arranged as a two-dimensional array of millions of cells (photodiodes), which convert an infinitesimal portion of the image into electric charge. The stored generated charge is then transformed into a pixel of digital image by an analog-to-digital converter. 
   A monochromatic sensor is made sensitive to a single color component by a color filter. A filter may act on primary colors or on complementary colors. The primary color filters are commonly arranged according to a scheme known in the art by the name of Bayer pattern in which half of the pixels are made sensitive to the green color, a fourth of the pixels to the red color and the last fourth to the blue color. As depicted in  FIG. 1 , green and blue pixels alternate along a row, while green and red pixels alternate on the successive row and this scheme repeats itself for the whole image. 
   The main reason the number of acquired green pixels are doubled compared to the number of acquired red and blue pixels is the greater sensitivity of human vision of the green chromatic component. Therefore, the value of each pixel of a CCD sensor carries information on a single primary color. 
   To reconstruct a good quality color image it is therefore desired to carry out an operation known as color interpolation for generating a triplet of values RGB for each pixel of the image to be produced at an output. Therefore, for example, for estimating the red component of a pixel for which only the green component has been acquired by the filtered light sensitive cell, it is necessary to take into consideration the adjacent pixels carrying information on the red component. 
   Bicubic and bilinear color interpolation are simple and efficient algorithms, but they may introduce large errors in the contour zones of the images. To address this problem several authors have proposed optimization algorithms that are sensitive to the input data. For example, Adams introduced an edge oriented method, and Kimmel, Adams, and Pei-Tam suggested other color correlation methods. 
   Triangulation is another well-known technique of linear interpolation that has been proposed for performing zooming on a digital image. A brief review of triangulation-based linear interpolation approaches is presented herein below. 
   Triangulation: Basic Considerations 
   Basically, the approach is to partition an image into triangles, the vertices of which are the known pixels, according to certain criteria. The color rendering or reconstruction phase follows during which each pixel value of the image to be produced is calculated through a linear interpolation of the vertices of the pertinent triangle. 
   Given a set of distinct points V={(x i , y i )}, what is of interest is a convex hull triangulation of V, that is a set T={T i } i=1 . . . q  of non degenerate triangles that satisfies the following conditions:
     every triangle vertex is an element of V, and every element of V is a triangle vertex;   every edge of a triangle in T contains exactly two points from V;   the intersection of any two different triangles in T is either empty, or is a shared edge or a shared vertex; and   the union of all triangles in T is the convex hull of V.   

   The criteria for choosing one among all the possible triangulations is to minimize a certain cost function C(T) of the entire triangulation T. Dyn-Levin-Rippa [5] suggested, among others, the following overall function cost:
 
 C ( T )=Σ|cost( e )|  (1)
 
for each edge e of T with a chosen cost(e). Once a possible triangulation is chosen, the values of the pixels will be interpolated by the unique 3D linear polynomial passing for the three vertices of the triangles. The 3D liner polynomial will be:
 
 P   i ( x,y )= a   i   x+b   i   y+c   i   (2)
 
for each triangle T i  of T.
 
Delaunay Triangulation
 
   According to the Delaunay triangulation each triangle is such that the circle passing for his three vertices does not contain any other vertex of V, where V={v i } set of distinct vertices. The construction of a Delaunay triangulation is defined through its Voronoi diagram. 
   The Voronoi diagram for V is the partition Vor(V) of the plane into the Voronoi polygons associated to V. The Voronoi polygon associated to v i  is the locus of points in the plane that are closer to v i  than to any other member of V. 
   The Delaunay triangulation can be defined as follow:
     partitioning of the set V in two subsets V 1  and V 2  of the same dimension;   recursively constructing Vor(V 1 ) and Vor(V 2 ); and   merging of Vor(V 1 ) and Vor(V 2 ) to obtain Vor(V).
 
 FIGS. 2 and 3  show the duality of the two problems.
   

   If the vertices of V are equally spaced (as it would be in the case of color interpolation), the resulting triangles will be equilateral triangles. Dyn-Levin-Rippa [9] demonstrated that triangulations involving long, narrow triangles give better results in linear interpolation. 
   Data Dependent Triangulation 
   Several optimization algorithms for performing what may be defined as Data Dependent Triangulation have been developed. For example, Lawson&#39;s algorithm [10] converges to the globally optimal solution if the cost function is based on the Delaunay criterion. Schumaker&#39;s simulated annealing method [11], gives good results, but is relatively slow. The “edge swap with look-ahead” method [8], is relatively faster than other algorithms. 
   The Yu-Morse-Sederberg&#39;s “edge swap with look-ahead” method includes the following steps:
         a) the phase of triangle construction (step  0 ) is basically a Delaunay Triangulation, conducted in a very simple way: all the vertices are joined to form quadrilaterals (for vertices equally spaced they would be squares) and the diagonal (edge) having the lower cost according to the particular cost function chosen in the quadrilateral will determine the initial triangulation, as depicted in  FIG. 4 .  FIG. 5  illustrates a single diagonal swapping inside a considered quadrilateral; and   b) an iterative phase (step  1  . . . n) during which each quadrilateral and its four adjacent triangles are considered, trying all the possible single swapping of the diagonals to minimize the local cost (that is the cost of all the 13 edges involved). The routine ends when all quadrilaterals have been processed.  FIG. 6  illustrates, at bottom, the possible 4 triangulations originated from a starting triangulation (depicted at the top of the figure) by a single edge swap.
 
The process is repeated n times, using the results of the previous iteration step.
 
Cost Functions
       

   A number of cost functions for characterizing a satisfactory triangulation interpolation have been proposed. 
   In particular, Dyn-Levin-Rippa [5] identify four cost functions called NC 1  (Nearly C 1 ), as particularly efficient. The cost for an edge e=[(x 1 ,y 1 ), (x 2 ,y 2 )] in common with two triangles T 1  and T 2 , the interpolating planes of which are respectively P 1 (x,y)=a 1 x+b 1 y+c 1  and P 2 (x,y)=a 2 x+b 2 y+c 2 , may be calculated according to a selected cost function:
     a) ABN (Angle Between Normals) cost function disclosed in Choi, et al. [12]
 
cost( e )=cos(α)  (3)
       where α is the angle between the normals to the interpolating planes P 1  and P 2 .   
       b) JND (Jump in Normal Derivatives) cost function
 
cost( e )=| n   x ( a   1   −a   2 )+ n   y ( b   1   −b   2 )|  (4)
       where (n x , n y ) is a unit vector orthogonal to the projection of the edge e.   
       c) DLP (Deviation from Linear Polynomials) cost function
 
cost( e )=∥ h∥   (5)
       where   
       

   
     
       
         
           
             h 
             = 
             
               [ 
               
                 
                   
                     
                        
                       
                         
                           
                             P 
                             1 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 x 
                                 2 
                               
                               , 
                               
                                 y 
                                 2 
                               
                             
                             ) 
                           
                         
                         - 
                         
                           F 
                           2 
                         
                       
                        
                     
                   
                 
                 
                   
                     
                        
                       
                         
                           
                             P 
                             2 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 x 
                                 1 
                               
                               , 
                               
                                 y 
                                 1 
                               
                             
                             ) 
                           
                         
                         - 
                         
                           F 
                           1 
                         
                       
                        
                     
                   
                 
               
               ] 
             
           
           ; 
         
       
     
       
       
         
            F i =value (x i , y i ) 
         
       
       d) DP (Distance from Planes) cost function
 
cost( e )=∥ g∥   (6)
       where   
     
     
  
   
     
       
         
           
             g 
             = 
             
               [ 
               
                 
                   
                     
                       dist 
                       ⁡ 
                       
                         ( 
                         
                           
                             P 
                             1 
                           
                           , 
                           
                             w 
                             2 
                           
                         
                         ) 
                       
                     
                   
                 
                 
                   
                     
                       dist 
                       ⁡ 
                       
                         ( 
                         
                           
                             P 
                             2 
                           
                           , 
                           
                             w 
                             1 
                           
                         
                         ) 
                       
                     
                   
                 
               
               ] 
             
           
           ; 
         
       
     
       
       
         
            w i =(x i , y i , F i ); 
         
       
     
  
   
     
       
         
           
             dist 
             ⁡ 
             
               ( 
               
                 
                   P 
                   j 
                 
                 , 
                 
                   w 
                   t 
                 
               
               ) 
             
           
           = 
           
             
                
               
                 
                   
                     P 
                     j 
                   
                   ⁡ 
                   
                     ( 
                     
                       
                         x 
                         t 
                       
                       , 
                       
                         y 
                         t 
                       
                     
                     ) 
                   
                 
                 - 
                 
                   F 
                   t 
                 
               
                
             
             
               sqrt 
               ⁡ 
               
                 ( 
                 
                   
                     a 
                     j 
                     2 
                   
                   + 
                   
                     b 
                     j 
                     2 
                   
                   + 
                   1 
                 
                 ) 
               
             
           
         
       
     
   
   Another cost function has been proposed by Yu-Morse-Sederberg [8]:
 
cost( e )=∥□ P   1   ∥*∥□P   2   ∥−□P   1   *□P   2   (7)
 
where
 
∥□ P   i ∥=sqrt( a   i   2   +b   i   2 )
 
 P   i ( x,y )= a   i   x+b   i   y+c   i 
 
□ P   i =( a   i   , b   i ) gradient of  P   i  
 
   DISCUSSION OF THE PRIOR ART 
   A first approach to use pixel level data dependent triangulation [5] even for color interpolation, was proposed by Su-Wills [6], following a previous work [7]. This method of color interpolation matches the edge orientation of the images and correlates the three distinct pixel channels (RGB). 
   Su-Wills Algorithm 
   Su-Wills algorithm [6] is based on the Delaunay Triangulation. As depicted in  FIG. 7 , the construction of triangles is relatively simple. For R and B channels (or patterns), they are constructed as in the case of a zooming by 2. Differently, for the G patterns, they are constructed with a rotation by 45 degrees. In practice, due to the simplicity of triangulation, the triangles do not require an appreciable burden for their definition. 
   To determine which diagonal should be chosen, a simple cost function is used as indicated in  FIG. 8 . As usual, the cost function for the G pattern or channel will be the same, considering a rotation of the triangles by 45 degrees. The interpolated points of the triangles (the black dots in  FIG. 8 ) are easily calculated as the average of the relative vertices. Practically, this known technique does not contemplate any iteration and though relatively simple to implement, does not produce a satisfactory color rendering in the interpolated images. 
   On the other hand, the use of more sophisticated triangulation algorithms as described above and of a cost function chosen among the above-mentioned four NC1 cost functions and the Yu-Morse-Sederberg [8] cost function for color interpolation would be excessively costly in terms of computation complexity, and in some cases, excessively slow. 
   It is generally perceived from these published works that it is not possible to have good interpolated images with a simple cost function and a fast triangulation algorithm, and that it is necessary to employ a sophisticated triangulation algorithm and a relatively complex cost function for obtaining satisfactory performance. 
   SUMMARY OF THE INVENTION 
   Contrary to these generally accepted conclusions, the present applicants have found simpler cost function and faster triangulation algorithms than those disclosed in Yu-Morse-Sederberg [8] that nevertheless produce interpolated images of comparably good quality. Surprisingly, applicants found that the quality of color interpolated images does not seem to degrade with respect to images obtained with the method disclosed in [8] when a computationally simpler cost function is used that does not require the calculation of powers and square roots. 
   Moreover, they found an efficient triangulation algorithm that may be executed in far less time than the triangulation algorithm disclosed in [8], while providing similar performance. The peculiarity of this triangulation algorithm may require on average only two iteration steps, while the triangulation algorithm disclosed in [8] may be completed only after four iteration steps. 
   Optionally, the interpolation process may be followed by anti-aliasing processing that effectively removes color artifacts that may be created during the interpolation process. 
   According to an innovative aspect of this invention in the case of a Bayer pattern, this step is based on the use of only two difference arrays, namely green-blue and green-red, respectively, and the consequent calculation of only two median arrays therefrom. In addition, the invention may be easily implemented with a software computer program. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
       FIG. 1  shows a Bayer Pattern Array according to the prior art; 
       FIG. 2  depicts a sample Voronoi diagram according to the prior art; 
       FIG. 3  illustrates how to obtain a Delaunay triangulation related to the Voronoi diagram according to the prior art; 
       FIG. 4  illustrates a possible construction phase (step  0 ) for vertices equally spaced according to the prior art; 
       FIG. 5  shows a single diagonal swapping inside the considered quadrilateral according to the prior art; 
       FIG. 6  shows all the possible 4 triangulation (bottom) originated from a starting triangulation (up) by doing a single edge swap according to the prior art; 
       FIG. 7  illustrates the triangles “construction” in Su-Wills algorithm: Left: R or B case; Right: G case according to the prior art; 
       FIG. 8  basically illustrates the Su-Wills algorithm according to the prior art; 
       FIG. 9  depicts a color interpolation scheme of according to the invention; 
       FIG. 10  illustrates the process steps performed in the block Data Dependent Triangulation of  FIG. 9 ; 
       FIG. 11  illustrates an exemplary process in accordance with the invention; 
       FIG. 12  shows the Bayer Data splitting according to the invention; 
       FIG. 13  shows a central square composed of two triangles and its four adjacent triangles obtained according to the invention; 
       FIGS. 14   a  and  14   b  illustrate possible edge swappings from the triangulation in  FIG. 13 ; 
       FIGS. 15   a  and  15   b  illustrate two cases that may occur while performing a data dependent triangulation technique (DDT) according to the invention; 
       FIG. 16  illustrates how to merge the three channels of an interpolated image according to the invention; 
       FIG. 17  shows the Freeman&#39;s algorithm according to the prior art; 
       FIG. 18  summarizes the differences between the Freeman&#39;s algorithm (left) and the algorithm of this invention (right). 
   

   DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
   The color interpolation method of this invention will be described by way of an example when considering the interpolation of a Bayer pattern image data as acquired by a digital color sensor. Naturally, the method is equally applicable to any other different scheme of acquisition using different filters, for example, for discriminating among complementary colors instead of primary colors. 
     FIG. 9  shows the flow scheme of the whole data dependent triangulation color interpolation according to the method of this invention, including also the optional anti-aliasing processing step following the interpolation step.  FIG. 10  illustrates in detail the operations performed in the step Data Dependent Triangulation shown in  FIG. 9 , while  FIG. 11  illustrates by way of an example an embodiment of the method of the invention. 
   By referring to the flow chart of  FIG. 9 , starting from the Bayer Pattern Data, three channels arrays are constructed, as shown in  FIG. 12 . The R, G and B pixels in the respective three channels are the vertices of the triangles to be constructed. 
   Improved Cost Function 
   The cost function that has been found to be outstandingly efficient in terms of quality vs. computational complexity is the following
 
cost( e )=(| a   1   |+|b   1 |)*(| a   2   |+|b   2 |)−□ P   1   *□P   2   (8)
 
where
 
 P   i ( x,y )= a   i   x+b   i   y+c   i 
 
and
 
 □P   i =( a   i   , b   i )
 
is the gradient of the interpolating function P i , which depends on the coordinates x, y of the pixel to be interpolated.
 
   It may be immediately noticed that this function is much simpler to calculate than the Yu-Morse-Sederberg cost function because it uses the calculation of absolute values, additions and multiplications. Using this cost function in the method of color interpolation disclosed in [8], the same quality of images is obtained, while at the same time reducing the computational cost of the interpolation algorithm. Another method in which this cost function may be conveniently used is disclosed hereinafter. 
   Improved Triangulation 
   According to another embodiment of the invention, a color interpolation method that uses simple data dependent triangulation has been found. Different from the known method of color interpolation by DDT of Su-Wills [6], the novel method of DDT color interpolation of this invention uses the “edge swap with look-ahead” technique of triangulation. The method of the invention contemplates the following operations:
         a) splitting the image in primary colors or complementary colors;   b) subdividing in blocks of 4×4 pixels or also 2×2 pixels an image in primary colors, generating a Delaunay data triangulation using a cost function; and   c) each square at the center of the block and its four adjacent triangles are considered, testing whether the at most two possible single edge swappings from the triangulation shown in  FIG. 13  of the diagonals, as depicted in  FIGS. 14   a  and  14   b , minimize the cost function, and eventually performing these swappings to minimize the local cost.       

   According to a preferred embodiment, the method of this invention is applied on all three images in primary colors and blocks of 4×4 pixels are considered. Of course, this method may be applied also to two images in primary colors and/or considering blocks of 2×2 pixels. 
   The process ends when all squares have been processed, and this may be done at most in two iteration steps. This is evident in  FIG. 13 , but let us consider the case of  FIG. 15   a , which is possible even if quite rare. 
   There are four triangles adjacent to the central square that form a quadrilateral with a triangle composing the central square. Only in this case, all the four adjacent triangles should be considered, but at most two swappings are possible. In fact, for each of such adjacent triangles there is another adjacent triangle that form a quadrilateral with the same triangle of the central square. 
   Therefore, if a diagonal of one of such quadrilaterals is swapped, the diagonal of the other quadrilateral involving the same triangle of the central square cannot be swapped anymore. If this were possible, there would be diagonals crossing each other, which is forbidden in any triangulation technique because each unknown pixel value is calculated by a function of only two known pixels values. 
   Only in the case of a triangulation as illustrated in  FIG. 15   a , which rarely occurs, the method of this invention would not produce substantial advantages of reduced computational complexity with respect to the method disclosed in [8]. 
   In case of a triangulation as illustrated in  FIG. 15   b , which also rarely occurs, there is not any edge swapping possible. Therefore, the method of this invention does not contemplate iteration and is even faster than the prior art method [8]. 
   The cases in which the central square will have only one, two or three adjacent triangles that form a quadrilateral with a triangle that compose it, are more probable than the case of  FIG. 15   a . On average, according to the method of this invention, the central square will have only two adjacent triangles that form a quadrilateral with a triangle that compose it, while in the algorithm of Yu, Morse and Sederberg, all four adjacent triangles should be considered. 
   According to a preferred embodiment of this method, the cost function to be used in steps b) and c) is given by eq. (8). Simulations indicate that more than one iteration step may not produce any visual perceptive improvement on the color interpolated image. Therefore, according to a preferred embodiment of this invention, only one iteration step may be performed. 
   Merging of Channels 
   The channels R and B are shifted compared to the G channel and then to the original image. For this reason a smart cropping of the borders may be necessary to avoid mixing of colors. 
     FIG. 16  shows how the merging may be done: only the central pixels are considered, while the peripheral ones are removed. In this implementation four rows and four columns are removed, so the interpolated image will be smaller than that of the original (input) Bayer Pattern. 
   Anti-Aliasing Algorithm 
   To remove color artifacts, due to the color interpolation step, an anti-aliasing algorithm is necessary. An effective algorithm is the Freeman&#39;s algorithm [13], which is summarized in the  FIG. 17 . 
   In particular for the RGB channels of a Bayer pattern input image, the Freeman&#39;s algorithm uses the following formulas: 
   for green (G) pixels in the original Bayer Pattern:
 
 Ĝ=G; {circumflex over (R)}=G−ν   GR   ; {circumflex over (B)}=G−ν   GB ;
 
for red (R) pixels in the original Bayer Pattern:
 
 {circumflex over (R)}=R; Ĝ=R+ν   GR   ; {circumflex over (B)}=R−ν   RB ;
 
for blue (B) pixels in the original Bayer Pattern:
 
 {circumflex over (B)}=B; Ĝ=B+ν   GB   ; {circumflex over (R)}=B+ν   RB ;
 
wherein ({circumflex over (R)}, {circumflex over (B)}, Ĝ) are the output RGB pixels and ν xy  is the median filter (the selection window of which is usually a 3×3 or 5×5 window centered on the pixel to be filtered) applied on the image difference (X−Y). Freeman&#39;s approach works fairly well to remove color artifacts, but it is slow due to the median filter elaboration of the three difference arrays.
 
   According to an effective approach, only the two difference arrays (G−B) and (G−R) may be used, and only the two median arrays ν GB  and ν GR  are calculated. The difference array (R−B) and the related median array ν RB  is not calculated, thus elaboration time is drastically reduced. 
   In detail, the color artifact removal after the color interpolation is performed through the steps of:
     for each pixel to be filtered, selecting a block of pixels from each component color image with a selection window of pre-established dimensions centered on the pixel to be filtered;   calculating a first median difference between the values of the selected green pixels and the selected red pixels v GR  and a second median difference between the values of the selected green pixels and the selected blue pixels v GB ; and   if the pixel to be filtered corresponds to:   a green pixel in the detected image, keeping its original green component and choosing as the value of a new red component thereof {circumflex over (R)} the difference between the green value and the first median difference v GR  and as the value of a new blue component thereof {circumflex over (B)} the difference between the green value and the second median difference v GB ,   a red pixel in the detected image, keeping its original red component and choosing as the value of a new green component thereof Ĝ the sum between the red value and the first median difference v GR  and as the value of a new blue component thereof {circumflex over (B)} the difference between the new green component Ĝ and the second median difference v GB , and   a blue pixel in the detected image, keeping its original blue component and choosing as the value of a new green component thereof Ĝ the sum between the blue value and the second median difference v GB  and as the value of a new red component thereof {circumflex over (R)} the difference between the new green component Ĝ and the first median difference v GR .   

     FIG. 18  summarizes this innovative aspect and highlights the difference between the just described algorithm and the Freeman&#39;s algorithm. Considering that the G patterns are double the number of R or B patterns, the output of the proposed algorithm is ⅚ equal to the output of Freeman&#39;s algorithm. In fact, if the repetitive pattern 
             [         GR           BG         ]               
with final 4*3=12 RGB samples is considered, it is possible to note that G patterns in the two outputs are equal and that the R or B patterns are ⅔ equal, so 6/12+ 2/12+ 2/12=⅚. Simulations revealed that the quality of the final image in the proposed algorithm is practically the same as that obtained with the Freeman&#39;s algorithm.
 
   REFERENCES 
   The following references were cited in foregoing description of the invention and the entire contents of each are hereby incorporated by reference.
     [1] J. E. Adams, Jr. “Interactions between Color Plane Interpolation and Other Image Processing Functions in Electronic Photography”  Proc. of SPIE  Vol. 2416, pp. 144-151, 1995   [2] R. Kimmel, “Demosaicing: Image reconstruction from colour CCD samples,”  IEEE Trans. Image Processing  Vol. 8, pp. 1221-1228, September 1999   [3] J. E. Adams, Jr. “Design of Practical Colour Filter Array Interpolation Algorithms for Digital Cameras”  Proc. Of SPIE  Vol. 3028, pp. 117-125, 1997   [4] S. C. Pei and I. K. Tam, “Effective colour interpolation in CCD colour filter array using signal correlation”  Proc. IEEE Int. Conf. Image Processing  Vol. 3, 2000, pp. 488-491   [5] N. Dyn, D. Levin, S. Rippa, “Data Dependent Triangulations for Piecewise Linear Interpolation” IMA Journal of Numerical Analysis, Vol. 10, pp. 137-154, 1990.   [6] D. Su, P. Willis, “Demosaicing of Colour Images Using Pixel Level Data-Dependent Triangulation”, Theory and Practice of Computer Graphics 2003, University of Birmingham (UK), June 2003.   [7] D. Su, P. Willis, “Image Interpolation by Pixel Level Data-Dependent Triangulation”, Computer Graphics Forum, October 2002.   [8] X. Yu, B. Morse, T. W. Sederberg, “Image Reconstruction Using Data-Dependent Triangulation”, IEEE Transactions on Computer Graphics and Applications, pp. 62-68, May/June 2001.   [9] N. Dyn, D. Levin, S. Rippa, “Long, thin triangles can be good for linear interpolation”, SIAM Journal of Numerical Analysis, 1992.   [10] C. L. Lawson, “Software for c1 Surface Interpolations,”  Mathematical Software III,  J. R. Rice, ed., Academic Press, New York, July 1977, pp. 161-194.   [11] L. Schumaker, “Computing Optimal Triangulations using Simulated Annealing,”  Computer Aided Geometric Design , vol. 10, nos. 3-4, 1993, pp. 329-345.   [12] B. Choi et al., “Triangulation of Scattered Data in 3DSpace,”  Computer Aided Design , vol. 20, no. 5, 1988, pp. 239-248.   [13] W. T. Freeman “Median filter for reconstructing missing color samples” U.S. Pat. No. 4,724,395 (1988).