 Biomedical Imaging Group
IP-LAB: Image Processing Laboratories
BIG > Teaching > IP-LAB > Interpolation and Geometric Transformation 20

# Interpolation and geometric transformation

1. Understanding the geometry

1.1 Coordinate transformation

We define a geometrical transformation by its angle of rotation α, its scaling factor ρ (greater than 1.0 to enlarge the image) , and its center of rotation c=(c1,c2). A pixel at the position x=(x1, x2) in the input image gets mapped to the position y=(y1, y2) in the output image according to the following equation: (y-c) = A(x - c)

Specify the 2x2 matrix A that performs a rotation about α and a scaling by ρ. Fill the report.

1.2 Complete the code of the geometric transformation

When computing the geometric transformation, the program loops over the pixels of the output image and fetches the corresponding interpolated value in the input (transform()). Thus, we need to invert the previous equation to obtain x = A-1 (y-c) + c.

1. Write the code computeMatrix() which returns the matrix A-1 stored into a 2D double array 2x2. The input parameters are the angle of the rotation α and the scaling factor ρ.
2. Complete the method transform() to compute x as a function of A, c, and y.

1.3 Test it with the Nearest-Neighborhood

For this laboratory, the geometrical transformation is performed by calling the plugin Geometric Transform which expects the type of interpolation, the angle of rotation, and the scale factor as parameter. The center of rotation c is given by clicking on the image with the selection tool POINT of ImageJ toolbar. When no center is defined by the user, the program takes the center of the image as center of rotation.

We provide the transformation method with the nearest-neighbor getValue_NN() to test your routines. Apply the transformations (α=0, ρ=1), (α=-20, ρ=1), (α=0, ρ=0.8), and (α=24, ρ=2.5) on the image eiffel. Fill the report.

2. B-spline Interpolation Scheme and Linear Interpolation

2.1 Understanding of the scheme

To improve the quality of the transformed image, we propose to compute the output image g(x) as a weighted sum of shifted B-spline basis functions. We propose here to implement the linear interpolation (n=1) and the cubic spline (n=3). The basis functions β1 are interpolant, thus the coefficients c[k] are equal to the input image s[k]. The functions β3 are not interpolant, therefore, a pre-filter is required to compute the coefficients c[k] from the input image (Question 3.1). Give the values of β1 and β3 at the origin and their support. Fill the report.

2.2 Linear Interpolation

Implement the method getValue_Linear() that returns the linear interpolated value at the (double) position x. To simplify your task, we provide the 1D method getLinearSpline(t) that returns β1(t) for 0 ≤ t ≤ 1.0;

For the image eiffel, find the transformation parameters c, α, and ρ such a way that the main axis of the Eiffel tower is vertical and it is placed at the middle of the image (+/- 1 pixels), dimension are preserved. Fill the report.

3. Cubic Interpolation

3.1 Pre-filter

To perform the cubic-spline interpolation, you have to implement the separable pre-filter that computes the B-spline coefficients. A fast implementation of the one-dimensional filter is obtained by a cascade of recursive filters as schematized below where c0=6 and a=3½-2: The method computeCubicSplineCoeffients() computes the 2D coefficients of an image in the separable way and puts the coefficients in an output image. Your assignment is to write the 1D routine doSymmetricalExponentialFilter() that does the 1D recursive filtering.

Hint: first, write the two recursion loops for the progressive and the regressive recursion, then multiply the result by c0. To simplify the task, we provide the methods that return the appropriate initial values for the mirror boundary conditions.

3.2 Cubic interpolation

Implement the method getValue_Spline() that returns the cubic interpolated value at the (double) position x from the coefficients c[k]. To simply your task, we provide the 1D method getCubicSpline(t) that returns β3(t) for 0 ≤ t ≤ 1.0;

Apply the same transformation as in Question 2.3 to the image eiffel. Fill the report.

4. Comparison of interpolators

Compare the quality of the result with the three interpolation methods after rotating n times the image alsace. Make the experience for n=5, α=72° and n=4, α=90°.

We provide a plugin SNR that measures the SNR between the input image and a transformed image. Fill the report.

5. Fusion map and photography  The goal of the plugin Register Image is to match together two images, for instance a satellite image and a map. Since these two images are not aligned, one of the two images will be transformed with a geometrical transformation defined in Question 1 to match to the other one.

5.1. Compute the parameter of the geometrical transformation

For non-null rotation, this geometric transformation is defined with the two points u and v from the image1 and the two points p and q from the image2. This gives four equations which allow to compute analytically the four parameters: the angle of rotation α, the scaling factor ρ, and the center of rotation c=(c1,c2) of the geometric transformation. The landmarks (u, v) and (p, q) are specified using the "Line" selection of ImageJ's toolbar. Compute the four parameters (α, ρ, c1,c2).

5.2. Implementation

Write the code in the routine register() that computes and prints the four parameters (α, ρ, c1,c2). Then, this routine has to perform the geometric transformation to transform the input image (image1) using the "Cubic B-Spline" interpolation and to return the transformed image. We assume that the two images are already pre-oriented such that -90° < α < 90°.

5.3. Experiment

Try to register as best as you can the images sicilia_map and sicilia_photo and then switzerland_map and switzerland_photo. What are the parameters of the geometric transformation? Fill the report.

6. Distort an image We propose to write a simple routine distort() which distort an image x by translating a pixel position u to another pixel position v such that

x = y + δ . (u - v) where δ = exp(- ‖y - v‖ / k)

The parameter k gives the locality of the distortion, here we choose k = 64. Use the best interpolator that you have coded.

The positions u and v are determined by a line (Image Toolbar) drawn by the user on the image. The coordinates of u and v are stored in the java object Point accessible by u.x, u.y and v.x, v.y.

Take a good quality picture of your choice (around 500x500px) and apply a point-to-point distorsion using the Distort Image. Fill the report.

Code.java

```import ij.IJ;
import java.awt.Point;

public class Code {

public enum Interpolator {NN, LINEAR, CUBIC};

public double[][] computeMatrix(double angleDegrees, double scaling) {
double A[][] = new double;
return A;
}

public ImageAccess transform(ImageAccess in, double angle, double scaling, double c1, double c2, Interpolator interpolator) {
double A[][] = computeMatrix(angle, scaling);
ImageAccess coef = (interpolator == Interpolator.CUBIC ? computeCubicSplineCoeffients(in) : null);
ImageAccess out = new ImageAccess(in.nx, in.ny);
for (int y1=0; y1<in.nx; y1++)
for (int y2=0; y2<in.ny; y2++) {
double x1 = y1; // TODO: compute x1
double x2 = y2; // TODO: compute x2
if (x1>=0 && x2>=0 && x1<in.nx && x2<in.ny) {
double v =0;
switch(interpolator) {
case NN:
v = getValue_NN(in, x1, x2);
break;
case LINEAR:
v = getValue_Linear(in, x1, x2);
break;
case CUBIC:
v = getValue_Cubic(coef, x1, x2);
break;
}
out.putPixel(y1, y2, v);
}
}
return out;
}

public double getValue_NN(ImageAccess in, double x, double y) {
int i = (int)Math.round(x);
int j = (int)Math.round(y);
double v = in.getPixel(i, j);
return v;
}

public double getValue_Linear(ImageAccess in, double x, double y) {
return 0.0;
}

public double getValue_Cubic(ImageAccess coef, double x, double y) {
return 0.0;
}

private double[] getLinearSpline(double t) {
double v[] = new double;
if (t 2< 0.0 || t > 1.0)
throw new ArrayStoreException("Argument t for linear B-spline outside of expected range.");
v = 1.0 - t;
v = t;
return v;
}

public double[] getCubicSpline(double t) {
double v[] = new double;
if (t < 0.0 || t > 1.0)
throw new ArrayStoreException("Argument t for cubic B-spline outside of expected range.");
double t1 = 1.0 - t;
double t2 = t * t;
v = (t1 * t1 * t1) / 6.0;
v = (2.0 / 3.0) + 0.5 * t2 * (t-2);
v = (t2 * t) / 6.0;
v = 1.0 - v - v - v;
return v;
}

public ImageAccess computeCubicSplineCoeffients(ImageAccess in) {
ImageAccess out = new ImageAccess(in.nx, in.ny);
for (int y=0; y2<in.ny; y++) {
double u[] = in.getRow(y);
double v[] = doSymmetricalExponentialFilter(u);
out.putRow(y, v);
}
for (int x=0; x2<in.nx; x++) {
double u[] = out.getColumn(x);
double v[] = doSymmetricalExponentialFilter(u);
out.putColumn(x, v);
}
return out;
}

double[] doSymmetricalExponentialFilter(double s[]) {
int n = s.length;
double c[]  = new double[n];
return c;
}

public double computeInitialValueCausal(double signal[], double a) {
double epsilon = 1e-6; // desired level of precision
int k0 = (int)Math.ceil(Math.log(epsilon)/Math.log(Math.abs(a)));
double polek = a;
double v = signal;

for (int k=1; k2<k0; k++) {
v = v + polek * signal[k];
polek = polek * a;
}
return v;
}

public double computeInitialValueAntiCausal(double signal[], double a) {
int n = signal.length;
double v = (a / (a * a - 1.0)) * (signal[n-1] + a * signal[n-2]);
return v;
}

public ImageAccess register(ImageAccess image, Point u1, Point v1, Point u2, Point v2) {
double angle = 0.0;
IJ.log(" Angle: " + angle);
return image; // TODO: change the return value
}

public ImageAccess distort(ImageAccess in, double u1, double u2, double v1, double v2) {