2D Cubic B-spline Interpolation via Digital Filtering. Java example

Image interpolation is very important operation in digital image processing and is used for images scaling and rotation, image compressing, image reconstruction and so on.  

Interpolation algorithms are differentiated by quality and efficiency. Splines introduce powerful instrument for image interpolation providing good images quality and computational efficiency. This is possible by using efficient filtering technique for processing images represented in terms of B-splines basis functions. The B-spline of degree 3 (cubic B-spline) is widely used for performing high-quality interpolation due to its minimum curvature property.

In my previous post 1D Cubic B-spline Interpolation via Digital Filtering. Java example I described 1-d signal interpolation using B-spines basis functions. The image itself is the 2-d signal represented by a set of uniformly spaces sampled values. It's easy to extend splines to higher dimensions  by using tensor-product basis functions. 

The 3-rd degree (cubic) tensor product B-spline function defined on regular grid is as follows:

Given splines coefficients  the 2-d  signal (image) continuous representations in terms of B-spline basis function expansion is as follows:

Now we can get values of continuous signal (perform interpolation) in any values of   we want. Because of the  b-spline function compact support for certain  there is only 16 neighbour b-splines should be multiplied with coefficients  to get the output sample . Figure 1 shows the case for  are integers and placed on uniform dense grid. 

Figure 1 - Interpolation performing on uniform integer grid

For a such case the image interpolation equation becomes as follows:

But before performing interpolation we should obtain B-spline coefficients  by means of direct B-spline filtration. Due to the separable kernel (2-d B-spline function is tensor-product) 2-d  filtration can be replaced by set of 1-d filtrations first along rows than along columns (or vice versa). For this task we can use already implemented 1-d direct B-spline filter, described in previous post 1D Cubic B-spline Interpolation via Digital Filtering. Java example . 

For correct boundary processing the mirror-W boundary condition is used working on the same principle as for 1-d case (see Figure 2).

Figure 2 - Mirrow-W condition for 2-d case

The following Java source code represents implementation of image cubic B-spline interpolation using filtering technique.

The CubicInterpolation2d class performs 2-d mirroring and interpolation. It uses BSplines class from my previous post “Introduction to splines. B-spline” .

The method cubicCoeff2d returns coefficients with mirroring for a given 2-d signal. The cubicInterp2d method preform interpolation for given signal and . It runs operations showed by Figure 2. The interpolate method performs interpolation for given signal and interpolation rate.

The DirectBsplFilter2d class performs 2-d direct B-spline filtration (obtaining coefficients). It uses DirectBsplFilter1d class from my previous post 1D Cubic B-spline Interpolation via Digital Filtering. Java example .

The ImageView class performs image visualisation.

CubicInterpolation2d class

package com.blogspot.shulgadim.splines;

import java.awt.image.BufferedImage;
import java.awt.image.Raster;
import java.io.File;
import java.io.IOException;

import javax.imageio.ImageIO;

public class CubicInterpolation2d {

    BSplines bS;    
    
    public CubicInterpolation2d(){
        bS = new BSplines();
    }
    
    private double[][] mirrorW2d( double [][] s){

        double [][] s_mirror = new double[s.length+3][s[0].length+3];
        for(int i=0; i<s.length; i++){
            for(int j=0; j<s[0].length; j++){
                s_mirror[i+1][j+1] = s[i][j];
            }
        }
        
        /*########## mirror rows  1 and N-2 ################*/
        for(int j=0; j<s[0].length; j++){
            s_mirror[0][j+1] = s[1][j];
            s_mirror[s_mirror.length-2][j+1] = s[s.length-2][j];
        }
        
        /*########## mirror columns 1 and  N-2 ##############*/ 
        for(int i=0; i<s.length; i++){
            s_mirror[i+1][0] = s[i][1];
            s_mirror[i+1][s_mirror[0].length-2] = s[i][s[0].length-2];
        }   
        
        s_mirror[0][0] = s[1][1];
        s_mirror[0][s_mirror[0].length-2] = s[1][s[0].length-2];            
        s_mirror[s_mirror.length-2][0] = s[s.length-2][1];
        s_mirror[s_mirror.length-2 ][s_mirror[0].length-2] = 
                s[s.length-2][s[0].length-2];
                                        
        return s_mirror;
    }
    
    public double[][] cubicCoeff2d(double[][] s){
        DirectBsplFilter2d directFilter2d;
        directFilter2d = new DirectBsplFilter2d(s.length, s[0].length);
        double[][] coeffs = directFilter2d.filter(s);
        double[][] coeffs_mirror = mirrorW2d(coeffs);
        return coeffs_mirror;
    }
    

    public double cubicInterp2d(double[][] coeffs_mirror, 
                                double row, 
                                double col){    
        int k = (int)Math.floor(row);   
        int l = (int)Math.floor(col);       
        double interp_value =   
        coeffs_mirror[k+0][l+0]*bS.bspline(3,row-k+1)*bS.bspline(3,col-l+1)+ 
        coeffs_mirror[k+1][l+0]*bS.bspline(3,row-k+0)*bS.bspline(3,col-l+1)+
        coeffs_mirror[k+2][l+0]*bS.bspline(3,row-k-1)*bS.bspline(3,col-l+1)+
        coeffs_mirror[k+3][l+0]*bS.bspline(3,row-k-2)*bS.bspline(3,col-l+1)+
                                                                    
        coeffs_mirror[k+0][l+1]*bS.bspline(3,row-k+1)*bS.bspline(3,col-l+0)+ 
        coeffs_mirror[k+1][l+1]*bS.bspline(3,row-k+0)*bS.bspline(3,col-l+0)+
        coeffs_mirror[k+2][l+1]*bS.bspline(3,row-k-1)*bS.bspline(3,col-l+0)+
        coeffs_mirror[k+3][l+1]*bS.bspline(3,row-k-2)*bS.bspline(3,col-l+0)+
                                                                                                
        coeffs_mirror[k+0][l+2]*bS.bspline(3,row-k+1)*bS.bspline(3,col-l-1)+ 
        coeffs_mirror[k+1][l+2]*bS.bspline(3,row-k+0)*bS.bspline(3,col-l-1)+
        coeffs_mirror[k+2][l+2]*bS.bspline(3,row-k-1)*bS.bspline(3,col-l-1)+
        coeffs_mirror[k+3][l+2]*bS.bspline(3,row-k-2)*bS.bspline(3,col-l-1)+
                                                
        coeffs_mirror[k+0][l+3]*bS.bspline(3,row-k+1)*bS.bspline(3,col-l-2)+ 
        coeffs_mirror[k+1][l+3]*bS.bspline(3,row-k+0)*bS.bspline(3,col-l-2)+
        coeffs_mirror[k+2][l+3]*bS.bspline(3,row-k-1)*bS.bspline(3,col-l-2)+
        coeffs_mirror[k+3][l+3]*bS.bspline(3,row-k-2)*bS.bspline(3,col-l-2);                                        
        return interp_value;
    }
    
    
    public double[][] interpolate(double[][] s, int rate){
        double [][] coeffs_mirror = cubicCoeff2d(s);        
        int M = rate*s.length - (rate-1);
        int N = rate*s[0].length - (rate-1);
        double [][] s_interp = new double[M][N];        
        for(int k=0; k<s_interp.length; k++){
            for(int l=0; l<s_interp[0].length; l++){
                s_interp[k][l] = 
                    cubicInterp2d(coeffs_mirror, k*(1.0/rate), l*(1.0/rate));
            }
        }                                                                                                       
        return s_interp;                    
    }
    
    
    private double[][] imageToDoubleArray(BufferedImage image){   
        double[][] array = new double[image.getHeight()][image.getWidth()];     
        Raster raster = image.getData();
        for( int y = 0; y < image.getHeight(); y++ ){
            for( int x = 0; x < image.getWidth(); x++ ){
                array[y][x] = raster.getSampleDouble(x, y, 0);                      
            }            
        }           
        return array;       
    }
    
    private BufferedImage doubleArrayToImage(double[][] array){     
        BufferedImage image = 
            new BufferedImage(array[0].length,array.length,
                              BufferedImage.TYPE_INT_RGB);  
        for( int y = 0; y < array.length; y++ ){
            for( int x = 0; x < array[0].length; x++ ){     
                 int value = 
                    (int)array[y][x] << 16 | 
                    (int)array[y][x] << 8 | 
                    (int)array[y][x];
                 image.setRGB(x, y, value);                                                   
            }
        }       
        return image;
    }
    
    public static void main(String[] args){
        BufferedImage image = null;
        try {
            image = ImageIO.read(new File("lena-std-2.jpg"));
        } catch (IOException e) {
        }
        if (image != null){
            ImageView imageView = new ImageView();
            imageView.drawImage(image);
            CubicInterpolation2d cubicInterpolation2d = 
                    new CubicInterpolation2d();
            double [][] img = 
                    cubicInterpolation2d.imageToDoubleArray(image);         
            double [][] img_interp = 
                    cubicInterpolation2d.interpolate(img,3);
            BufferedImage imageInterp = 
                    cubicInterpolation2d.doubleArrayToImage(img_interp);
            ImageView imageView2 = new ImageView();
            imageView2.drawImage(imageInterp);  
        }
    }
}

DirectBsplFilter2d class

package com.blogspot.shulgadim.splines;

public class DirectBsplFilter2d {

    private DirectBsplFilter1d directFilter1dX; 
    private DirectBsplFilter1d directFilter1dY;         
    
    public DirectBsplFilter2d(int M, int N){
        directFilter1dX = new DirectBsplFilter1d(M);    
        directFilter1dY = new DirectBsplFilter1d(N);    
    }
    
    public double[][] filter(double [][] img){  
        double [] row = new double[img[0].length];
        double [] col = new double[img.length];
        double [] filt_row = new double[img[0].length];
        double [] filt_col = new double[img.length];
        double [][] coeffs = new double[img.length][img[0].length];
        
        /*#################### filtrations along y ##################*/
        for(int i=0; i<img.length; i++){
            for(int j=0; j<img[0].length; j++){
                row[j] = img[i][j];
            }
            filt_row = directFilter1dY.filter(row);
            for(int j=0; j<img[0].length; j++){
                coeffs[i][j] = filt_row[j];
            }
            directFilter1dY.reset();
        }
        /*#################### filtrations along y ##################*/
                    
        /*#################### filtrations along x ##################*/
        for(int j=0; j<img[0].length; j++){
            for(int i=0; i<img.length; i++){
                col[i] = coeffs[i][j];
            }
            filt_col = directFilter1dX.filter(col);
            for(int i=0; i<img.length; i++){
                coeffs[i][j] = filt_col[i];
            }
            directFilter1dX.reset();
        }
        /*#################### filtrations along x ##################*/             
        return coeffs;                                                  
    }

}

ImageView class

package com.blogspot.shulgadim.splines;

import java.awt.*;
import javax.swing.*;
import java.awt.image.*;
import javax.swing.JFrame;

@SuppressWarnings("serial")
public class ImageView extends JFrame{

    private SignalPanel signalPanel;
    
    public ImageView(){                     
        setTitle("Image");      
        setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        setSize(200, 200);       
        signalPanel = new SignalPanel();                 
        getContentPane().add(signalPanel);          
        setVisible(true);        
    }
    
    public void drawImage(BufferedImage image){
        setSize(image.getWidth()+25, image.getHeight()+50);         
        signalPanel.setImage( image);
    }   
}

@SuppressWarnings("serial")
class SignalPanel extends JPanel {
    
   private BufferedImage image;
   
   @Override
   protected void paintComponent(Graphics graphics) {
       super.paintComponent(graphics);   
       if (image != null) {
           graphics.drawImage(image, 5, 5, this);
       }
   }
   
   public void setImage(BufferedImage image){
       this.image = image;
   }  
}

The input image showed on Figure 3.

Figure 3 - Input image

The output interpolated image with rate=3 showed on Figure 4.

Figure 4 - Output image

Sources:

[1] M. Unser "Splines: A Perfect Fit for Signal and Image Processing" IEEE Signal Processing Magazine, vol. 16, no. 6, pp. 22-38, November 1999.
[2] M. Unser, A. Aldroubi, M. Eden "B-Spline Signal Processing: Part I—Theory" IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 821-833, February 1993.

[3] M. Unser, A. Aldroubi, M. Eden "B-Spline Signal Processing: Part II—Efficient Design and Applications" IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 834-848, February 1993.