the grayscale intensity value from the source image, 
, at 
, and that we are trying to calculate the interpolated grayscale value in the target\/output image, 
, at 
) (Equations based on those provided in [1]):<\/em><\/p>\n- \n
- The original intensity values of the surrounding pixels (4 values) <\/em>\n
- \n
<\/em><\/li>\n
<\/em><\/li>\n
<\/em><\/li>\n
<\/em><\/li>\n<\/ul>\n<\/li>\n
- The partial derivatives along the x-axis for each of these points (4 values) <\/em>\n
- \n
<\/em><\/li>\n<\/ul>\n<\/li>\n
- The partial derivatives along the y-axis for each of these points (4 values) <\/em>\n
- \n
<\/em><\/li>\n<\/ul>\n<\/li>\n
- The cross-derivatives at each of these points (4 values) <\/em>\n
- \n
<\/em><\/li>\n<\/ul>\n<\/li>\n<\/ul>\n
This leaves us with a total of 16 values to solve for in the case of generalized bicubic interpolation, as opposed to the 4 needed with 1D cubic spline interpolation (
and
).<\/em><\/p>\n
End Quote<\/h1>\n
I also mentioned in the previous section that bicubic spline interpolation is a special case of generalized bicubic interpolation. For the former of the two, we simply carried out 1D cubic spline interpolation along each axis and we were done. We aren’t making use of all the data available to us though, such as the cross-derivatives.<\/p>\n
For each of the neighbouring 4 data points, we need to know its intensity value (
), its partial derivatives along both axes (
![\"\"](\"https:\/\/www.giassa.net\/wp-content\/uploads\/2009\/12\/grid-others1.png\" "\\"grid-others\\"")
<\/a>Grid Coordinates<\/p><\/div>\n
We can calculate the partial derivatives easily enough using centered differences, which are basically slope calculations using the intensity values,
from adjacent pixels. For the partial derivative along the horizontal axis, <\/em>
- \n
<\/li>\n<\/ul>\n
Similarly, we calculate the partial derivative for the vertical axis at each neighbouring pixel by:<\/p>\n
- \n
<\/li>\n<\/ul>\n
Finally, for the cross derivatives, we can calculate them by:<\/p>\n
- \n
<\/li>\n
- Therefore,
<\/li>\n<\/ul>\n
If the spacing between pixels is not equal (ie: you are only carrying out interpolation on part of an image rather than the whole thing), you will need to adjust the denominator values in the above equations. Note that even if you get these values wrong, <\/strong>the final result will still be smooth, just not as accurate!<\/strong><\/p>\n
With these values available, we can defined any interpolated value in our output\/target image,
, as a function of the pixels from our input image according to (refer to the grid coordinates in the figure above for visual definitions of A,B,J,K,h, and w):<\/p>\n
- \n
<\/li>\n
<\/li>\n
<\/li>\n
<\/li>\n<\/ul>\n
However, before we can carry on, we need to solve for the coefficients,
. Assuming that
and
are continuous over the interval
, and that we redefine the values at the gridpoints of the source image,
as
, we can express all 16 coefficients of
as follows (equations provided by [1]):<\/p>\n
- \n
- Function values:\n
- \n
<\/li>\n
<\/li>\n
<\/li>\n
<\/li>\n<\/ul>\n<\/li>\n
- X-Derivatives:\n
- \n
<\/li>\n
<\/li>\n
<\/li>\n
<\/li>\n<\/ul>\n<\/li>\n
- Y-Derivatives:\n
- \n
<\/li>\n
<\/li>\n
<\/li>\n
<\/li>\n<\/ul>\n<\/li>\n
- Cross-Derivatives (XY):\n
- \n
<\/li>\n
<\/li>\n
<\/li>\n
<\/li>\n<\/ul>\n<\/li>\n<\/ul>\n
These coefficient values, combined with the previously calculated partial and cross-derivative values, form a set of linear equations which form a matrix equation of the general form
. We will rewrite it in the form
so that we don’t confuse the values for our matrices with those used in the equations above. Let:<\/p>\n
- \n
<\/li>\n
<\/li>\n
\n
![\"\"](\"https:\/\/www.giassa.net\/wp-content\/uploads\/2009\/12\/M-matrix.png\" "\\"M-matrix\\"")
<\/a>Coefficient Matrix for ‘a’-values<\/p><\/div><\/li>\n<\/ul>\n
Note that the
refers to the transpose of these matrices, rather than an exponent (which would only make sense for a square matrix anyways). With these equations, we can use our favorite method (feel free to pick) of solving this set of linear equations to solve for the inverse of our weighting matrix,
. Once this is done, we can save this matrix for use in any future implementations of bicubic interpolation, as there is no need to recalculate it again in the future. By simple matrix inversion (yes, you can do this by hand if you really want to, the matrix algebra is very simple, just time-consuming; it’d be better to do this with MATLAB), we get the inverted coefficient matrix,
.<\/p>\n
- \n
![\"\"](\"https:\/\/www.giassa.net\/wp-content\/uploads\/2009\/12\/M-matrix-inverted.png\" "\\"M-matrix-inverted\\"")
<\/a>Inverted Coefficient Matrix<\/p><\/div>\n
When
, this gives us the values of each of the 16 coefficients in our
matrix. We can now use the previously defined sigma-notation formula to determine the value of any interpolated data point using the formula:<\/p>\n
- \n
The following MATLAB M-function demonstrates this algorithm in action:<\/p>\n
% GBI.M\r\n% Written December 2009, (C) Matthew Giassa\r\n% teo@giassa.net, www.giassa.net\r\n% Returns an upsampled image using bicubic interpolation\r\nfunction output_image = gbi( input_image,x_res,y_res )\r\n% input_image - an image on which to perform bicubic interpolation\r\n% x_res - the new horizontal dimensions (in pixels)\r\n% y_res - the new vertical dimensions (in pixels)\r\n%Define the inverted weighting matrix, M^(-1), no need to recalculate it\r\n%ever again\r\nM_inv = [\r\n 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;\r\n 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0;\r\n -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0;\r\n 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0;\r\n 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0;\r\n 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0;\r\n 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0;\r\n 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0;\r\n -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0;\r\n 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0;\r\n 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1;\r\n -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1;\r\n 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0;\r\n 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0;\r\n -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1;\r\n 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1\r\n ];\r\n%Make a copy of the input image\r\nI = input_image;\r\n%Convert image to grayscale (intensity) values for simplicity (for now)\r\nI = double(rgb2gray(I));\r\n%Determine the dimensions of the source image\r\n%Note that we will have three values - width, height, and the number\r\n%of color vectors, 3\r\n[j k] = size(I);\r\n%Specify the new image dimensions we want for our larger output image\r\nx_new = x_res;\r\ny_new = y_res;\r\n%Determine the ratio of the old dimensions compared to the new dimensions\r\n%Referred to as S1 and S2 in my tutorial\r\nx_scale = x_new.\/(j-1);\r\ny_scale = y_new.\/(k-1);\r\n%Declare and initialize an output image buffer\r\ntemp_image = zeros(x_new,y_new);\r\n%Calculate the horizontal derivatives for each point (assume that at the\r\n%edge pixels of the source image, all derivatives are zero to maximize\r\n%smoothing). Take note - recreating this in C or C++ might not be so easy,\r\n%as you must pay attention to how this code is using matrix-style\r\n%addressing (x-y positions swapped).\r\nIx = double(zeros(j,k));\r\nfor count1 = 1:j\r\n for count2 = 1:k\r\n if( (count2==1) || (count2==k) )\r\n Ix(count1,count2)=0;\r\n else\r\n Ix(count1,count2)=(0.5).*(I(count1,count2+1)-I(count1,count2-1));\r\n end\r\n end\r\nend\r\n%Similarly, calculate the vertical derivatives\r\nIy = double(zeros(j,k));\r\nfor count1 = 1:j\r\n for count2 = 1:k\r\n if( (count1==1) || (count1==j) )\r\n Iy(count1,count2)=0;\r\n else\r\n Iy(count1,count2)=(0.5).*(I(count1+1,count2)-I(count1-1,count2));\r\n end\r\n end\r\nend\r\n%Finally, calculate the cross derivatives\r\nIxy = double(zeros(j,k));\r\nfor count1 = 1:j\r\n for count2 = 1:k\r\n if( (count1==1) || (count1==j) || (count2==1) || (count2==k) )\r\n Ixy(count1,count2)=0;\r\n else\r\n Ixy(count1,count2)=(0.25).*((I(count1+1,count2+1)+I(count1-1,count2-1)) - (I(count1+1,count2-1)+I(count1-1,count2+1)));\r\n end\r\n end\r\nend\r\n%Generate the output image\r\n%==========================================================================\r\n%Please note that this is definitely not the most efficient way to carry\r\n%out this operation, as we are recalculating our beta and alpha-vectors far\r\n%repeatedly when there is no need. I kept the code like this so you can see\r\n%the order in which we would carry out most operations. If you want to make\r\n%a more efficient version which pre-calculates the alpha vector for each\r\n%pixel in the source image, be my guest.\r\n%==========================================================================\r\nfor count1 = 0:x_new-1\r\n for count2 = 0:y_new-1\r\n %Calculate the normalized distance constants, h and w\r\n W = -(((count1.\/x_scale)-floor(count1.\/x_scale))-1);\r\n H = -(((count2.\/y_scale)-floor(count2.\/y_scale))-1);\r\n %Determine the indexes\/address of the 4 neighbouring pixels from\r\n %the source data\/image\r\n I11_index = [1+floor(count1.\/x_scale),1+floor(count2.\/y_scale)];\r\n I21_index = [1+floor(count1.\/x_scale),1+ceil(count2.\/y_scale)];\r\n I12_index = [1+ceil(count1.\/x_scale),1+floor(count2.\/y_scale)];\r\n I22_index = [1+ceil(count1.\/x_scale),1+ceil(count2.\/y_scale)];\r\n %Calculate the four nearest function values\r\n I11 = I(I11_index(1),I11_index(2));\r\n I21 = I(I21_index(1),I21_index(2));\r\n I12 = I(I12_index(1),I12_index(2));\r\n I22 = I(I22_index(1),I22_index(2));\r\n %Calculate the four nearest horizontal derivatives\r\n Ix11 = Ix(I11_index(1),I11_index(2));\r\n Ix21 = Ix(I21_index(1),I21_index(2));\r\n Ix12 = Ix(I12_index(1),I12_index(2));\r\n Ix22 = Ix(I22_index(1),I22_index(2));\r\n %Calculate the four nearest vertical derivatives\r\n Iy11 = Iy(I11_index(1),I11_index(2));\r\n Iy21 = Iy(I21_index(1),I21_index(2));\r\n Iy12 = Iy(I12_index(1),I12_index(2));\r\n Iy22 = Iy(I22_index(1),I22_index(2));\r\n %Calculate the four nearest cross derivatives\r\n Ixy11 = Ixy(I11_index(1),I11_index(2));\r\n Ixy21 = Ixy(I21_index(1),I21_index(2));\r\n Ixy12 = Ixy(I12_index(1),I12_index(2));\r\n Ixy22 = Ixy(I22_index(1),I22_index(2));\r\n %Create our beta-vector\r\n beta = [I11 I21 I12 I22 Ix11 Ix21 Ix12 Ix22 Iy11 Iy21 Iy12 Iy22 Ixy11 Ixy21 Ixy12 Ixy22];\r\n %Calculate our alpha vector (ie: the a-values) using simple matrix\r\n %multiplication\r\n %==================================================================\r\n %If we wanted to make sure this entire program is written entirely\r\n %from scratch, we would implement our own code to multiple an AxB\r\n %matrix and a BxC matrix right here. This is very simple to do, so\r\n %I won't go over it right now.\r\n %==================================================================\r\n alpha = M_inv*beta';\r\n temp_p=0;\r\n for count3 = 1:16\r\n w_temp = floor((count3-1)\/4);\r\n h_temp = mod(count3-1,4);\r\n %disp(sprintf('aij=%d wsub=%d hsub=%d\\n',count3,w_temp,h_temp))\r\n temp_p = temp_p + alpha(count3).*((1-W)^(w_temp)).*((1-H)^(h_temp));\r\n end\r\n temp_image(count1+1,count2+1)=temp_p;\r\n end\r\nend\r\noutput_image = temp_image;\r\n<\/pre>\n<\/p>\n
Included below is a simple resampling of a 4×4 checkerboard pattern to a much larger resolution using the code above.<\/p>\n
![\"\"](\"https:\/\/www.giassa.net\/wp-content\/uploads\/2009\/12\/bci-mine-300x265.png\" "\\"bci-mine\\"")
<\/a>Bicubic Interpolation on 4×4 Checkerboard Pattern<\/p><\/div>\n
Before concluding this section, I should note that there are many other resampling methods that are definitely worth reading up on, such as Lanczos, Hermite spline, Super Eagle, 2xSAI, HQn, and so on. Many of these algorithms are used in commerical products such as image processing software digital cameras, even console emulators. I will now conclude this section with a brief comparison of “bicubic spline interpolation” and “generalized bicubic interpolation”, and will show why the two are NOT<\/strong> the same.<\/p>\n
First, observe the image pairs below. The first pair was generated by using bicubic spline interpolation, which was covered in the previous subsection, while the latter pair was generated with the generalized bicubic interpolation algorithm covered in this subsection (the two original smiley faces used are slightly different, but no matter):<\/p>\n
![\"\"](\"https:\/\/www.giassa.net\/wp-content\/uploads\/2009\/12\/smile-halo1-300x209.png\" "\\"smile-halo\\"")
<\/a>Resampling Using Bicubic Spline Interpolation<\/p><\/div>\n
![\"\"](\"https:\/\/www.giassa.net\/wp-content\/uploads\/2009\/12\/smile-bci-300x202.png\" "\\"smile-bci\\"")
<\/a>Resampling Using Generalized Bicubic Interpolation<\/p><\/div>\n
If we observe the resampling smiley faces, we can see that both still have halos introduced due to our choice of resampling algorithms. One important difference exists between both of the sampled images, however – the former of the two (generated using bicubic spline interpolation) still demonstrates some sharp edges and even corners, while the latter does not. This is due to bicubic spline interpolation being a rough approximation of generalized bicubic interpolation at best, primarily due to it not making any use of the cross-derivatives of the source data.<\/p>\n
References<\/h1>\n
[1] Numerical recipes in C the art of scientific computing<\/em>. Cambridge: Cambridge UP, 1992. Print.<\/p>\n","protected":false},"excerpt":{"rendered":"
If you haven’t read the previous section on bicubic spline interpolation, now would be a good time! As mentioned in the previous section, when we are performing bicubic interpolation, we are generating a smooth surface that interpolates data points on … Continue reading
- Function values:\n
- The cross-derivatives at each of these points (4 values) <\/em>\n
- The partial derivatives along the y-axis for each of these points (4 values) <\/em>\n