condat_shift.py

# Implements the 1D fractional-pixel shift operator from
# Condat, L. (2008). FULLY REVERSIBLE IMAGE ROTATION BY 1-D FILTERING. Signal Processing

import numpy as np
from scipy.misc import comb


def make_b(N, tau):

    b = np.zeros(N)

    for k in range(1, N + 1):

        # Note; we use zero-based indexing for k whereas the paper uses one-based indexing
        C = comb(N, k)
        numerator = tau - np.arange(0, N + 1)
        denominator = tau - np.arange(0, N + 1) - k
        b[k - 1] = ((-1) ** k) * C * np.prod(numerator / denominator)

    return b


def condat_shift(s, N, tau):

    b = make_b(N, tau)
    T = s.size
    for i in range(T - 1, -1, -1):
        for k in range(1, N + 1):
            s[i] += b[k - 1] * (s[(i - k) % T] - s[(i + k) % T])

    return s


s = np.zeros(100)
s[20:50] = 1.
condat_shift(s, 100, 0.5)

import matplotlib.pyplot as plt
plt.plot(s)