Demo iamagnify

Illustrate the interpolation of magnified images

Description

In this demonstration, the interpolation of magnified image is explained using the frequency domain. The interpolation is a low pass filter that can be applied either in the spatial domain, which is the methods known as nearest neighbor or pixel replication and bi-linear interpolation. Or in the frequency domain, with an ideal filter or a butterworth filter.

Script

Reading and ROI selection

The image is read and a 64x64 ROI is selected and displayed

1 fin = iaread('lenina.pgm')
2 adshow(fin)
3 froi = fin[137:137+64,157:157+64]
4 adshow(froi)
5 print froi.shape
(64, 64)

DFT

The DFT of the small image is taken and its spectrum displayed

1 import numpy as np
2 from ia636 import iadftview
3 
4 fd = froi.astype(float)
5 F = np.fft.fft2(fd)
6 adshow(froi)
7 adshow(iadftview(F))

Expansion by 4 without interpolation

The image is expanded by 4, but filling the new pixels with 0

1 fx4 = zeros(4*array(froi.shape))
2 fx4[::4,::4] = froi
3 adshow(froi)
4 adshow(fx4)

DFT of the expansion without interpolation

Using the expansion propertie of the DFT (only valid for the discrete case), the resulting DFT is a periodical replication of the original DFT.

1 fdx4 = fx4.astype(float)
2 Fx4 = np.fft.fft2(fdx4)
3 
4 adshow(fx4)
5 adshow(iadftview(Fx4))

Filtering by mean filtering - nearest neighbor

Filtering the expanded image using an average filter of size 4x4 is equivalent of applying a nearest neighbor interpolator. The zero pixels are replaced by the nearest non-zero pixel. This is equivalent to interpolation by pixel replication.

1 from ia636 import iapconv
2 
3 k = ones((4,4))
4 fx4nn = iapconv(fdx4, k)
5 adshow(fx4)
6 adshow(fx4nn.astype(int32))

Interpretation of the mean filtering in the frequency domain

Filtering by the average filter in space domain is equivalent to filter in the frequency domain by the sync filter.

1 kzero = zeros(fx4.shape)
2 kzero[0:4,0:4] = k
3 K = np.fft.fft2(kzero)
4 adshow(iadftview(K))
5 Fx4nn = K * Fx4
6 adshow(iadftview(Fx4nn))

Filtering by pyramidal kernel, linear interpolation

Filtering by a pyramidal kernel in space domain is equivalent to make a bi-linear interpolation. The zero pixels are replaced by a weighted sum of the neighbor pixels, the weight is inversely proportional to the non-zero pixel distance.

1 klinear = array([1,2,3,4,3,2,1])/4.
2 k2dlinear = dot(reshape(klinear, (7,1)), reshape(klinear, (1,7)))
3 print 'k2dlinear=',k2dlinear
4 fx4li = iapconv(fdx4, k2dlinear)
5 adshow(fx4)
6 adshow(fx4li.astype(int32))
k2dlinear= [[ 0.0625  0.125   0.1875  0.25    0.1875  0.125   0.0625]
 [ 0.125   0.25    0.375   0.5     0.375   0.25    0.125 ]
 [ 0.1875  0.375   0.5625  0.75    0.5625  0.375   0.1875]
 [ 0.25    0.5     0.75    1.      0.75    0.5     0.25  ]
 [ 0.1875  0.375   0.5625  0.75    0.5625  0.375   0.1875]
 [ 0.125   0.25    0.375   0.5     0.375   0.25    0.125 ]
 [ 0.0625  0.125   0.1875  0.25    0.1875  0.125   0.0625]]

Interpretation of the pyramid filtering in the frequency domain

Filtering by the pyramid filter in space domain is equivalent to filter in the frequency domain by the square of the sync filter.

1 klizero = zeros(fx4.shape).astype(float)
2 klizero[0:7,0:7] = k2dlinear
3 Klinear = np.fft.fft2(klizero)
4 adshow(iadftview(Klinear))
5 Fx4li = Klinear * Fx4
6 adshow(iadftview(Fx4li))

Using an ideal filter

Filtering by cutoff period of 8

 1 from ia636 import iabwlp
 2 from ia636 import ianormalize
 3 
 4 H8 = iabwlp(fx4.shape, 8, 10000)
 5 adshow(iadftview(H8))
 6 G8 = Fx4 * H8
 7 adshow(iadftview(G8))
 8 g_ideal = np.fft.ifft2(G8)
 9 print 'Max of imaginary:', g_ideal.imag.max()
10 g_ideal = ianormalize(g_ideal.real, [0,255])
11 adshow(g_ideal)
Max of imaginary: 2.96873329651e-14

Using a Butterworth filter of order 5

Filtering by cutoff period of 8

1 HB8 = iabwlp(fx4.shape, 8, 5)
2 adshow(iadftview(HB8))
3 GB = Fx4 * HB8
4 adshow(iadftview(GB))
5 g_b = np.fft.ifft2(GB)
6 print 'Max of imaginary:', g_b.imag.max()
7 g_b = ianormalize(g_b.real, [0,255])
8 adshow(g_b)
Max of imaginary: 3.19653094491e-14

Display all four for comparison

Top-left: nearest neighbor, Top-right: linear, Bottom-left: ideal, Bottom-right: Butterworth

1 aux1 = np.concatenate((fx4nn[0:256,0:256], fx4li[0:256,0:256]), 1)
2 aux2 = np.concatenate((g_ideal, g_b), 1)
3 adshow(np.concatenate((aux1, aux2)))