### Subsampling image arrays in Python

Recently, as I was writing a class that performs Dokkum's algorithm for cosmic ray detection, I came across a nifty math trick. One of the steps is to

By

\begin{align}A=\begin{bmatrix}a&b&c\\d&e&f\\g&h&i\end{bmatrix},\end{align}

our goal (if we wish to subsample only once) is to make it look like

\begin{align}B=\begin{bmatrix}a&a&b&b&c&c\\a&a&b&b&c&c\\d&d&e&e&f&f\\d&d&e&e&f&f\\g&g&h&h&i&i\\g&g&h&h&i&i\end{bmatrix}.\end{align} This is where the Kronecker product comes in. It was an odd step at first because I usually don't think of whole image arrays (first) as matrices. Thus, letting \(I_n\) be the \(n\times n\) ones matrix, we have

\begin{align}A\otimes_{K} I_{2}&=\begin{bmatrix}a\cdot I_{2}&b\cdot I_{2}\\c\cdot I_{2}&d\cdot I_{2}\end{bmatrix}\\&=B.\end{align}This is, of course, exactly what we were looking for.

To do the opposite, i.e. `block average`, we partition the matrix into blocks

\begin{align}B=\begin{bmatrix}B_{11}&B_{12}&B_{13}\\B_{21}&B_{22}&B_{23}\\B_{31}&B_{32}&B_{33}\end{bmatrix},\end{align} and reshape the partitioned matrix into a vector

\begin{align}\vec{B}:=\begin{bmatrix}B_{11}&\cdots&B_{1n^2}\end{bmatrix}^{\text{T}}.\end{align} Now, multiplying \(\langle \vec{B}, \begin{bmatrix}I_{n}\end{bmatrix}/n^2\rangle_{F}\) (and waving the details on the fact that the former does not have the same dimensions as the latter -- computers are far less picky), yields \(A\). Surely there's a faster way, but for moderately large images it doesn't take more than a few seconds. This all gets abstracted away to loops (and the reshape isn't entirely necessary).

The image on the right looks the same because it's the same relative size as the original, but each pixel has been replaced with 4 of the same size.

*subsample*the image. At first thought, one may interpret that word as meaning that the image will be scaled down.*Au contraire.*By

*subsample*the author actually means to subsample an individual pixel to a square array of pixels before convolution, each the same count as the original. So if we had an array that looked like\begin{align}A=\begin{bmatrix}a&b&c\\d&e&f\\g&h&i\end{bmatrix},\end{align}

our goal (if we wish to subsample only once) is to make it look like

\begin{align}B=\begin{bmatrix}a&a&b&b&c&c\\a&a&b&b&c&c\\d&d&e&e&f&f\\d&d&e&e&f&f\\g&g&h&h&i&i\\g&g&h&h&i&i\end{bmatrix}.\end{align} This is where the Kronecker product comes in. It was an odd step at first because I usually don't think of whole image arrays (first) as matrices. Thus, letting \(I_n\) be the \(n\times n\) ones matrix, we have

\begin{align}A\otimes_{K} I_{2}&=\begin{bmatrix}a\cdot I_{2}&b\cdot I_{2}\\c\cdot I_{2}&d\cdot I_{2}\end{bmatrix}\\&=B.\end{align}This is, of course, exactly what we were looking for.

To do the opposite, i.e. `block average`, we partition the matrix into blocks

\begin{align}B=\begin{bmatrix}B_{11}&B_{12}&B_{13}\\B_{21}&B_{22}&B_{23}\\B_{31}&B_{32}&B_{33}\end{bmatrix},\end{align} and reshape the partitioned matrix into a vector

\begin{align}\vec{B}:=\begin{bmatrix}B_{11}&\cdots&B_{1n^2}\end{bmatrix}^{\text{T}}.\end{align} Now, multiplying \(\langle \vec{B}, \begin{bmatrix}I_{n}\end{bmatrix}/n^2\rangle_{F}\) (and waving the details on the fact that the former does not have the same dimensions as the latter -- computers are far less picky), yields \(A\). Surely there's a faster way, but for moderately large images it doesn't take more than a few seconds. This all gets abstracted away to loops (and the reshape isn't entirely necessary).

And the result is as follows.#! /usr/bin/env python3.4 # -*- coding: utf-8 -*- """ Testing the outer product on an image array """ import numpy as np from PIL import Image IMAGE = '/home/brandon/Pictures/Abominable.jpg' SUB = 4 class OpenImage: """ Image opening context manager """ def __init__(self, imageName: str ='', NAXIS: int =0) -> None: self.imageName = imageName self.NAXIS = NAXIS def __enter__(self) -> np.ndarray: self.image = Image.open(self.imageName) if (''.join(self.image.getbands()) == 'RGB'): return np.rot90(np.array(self.image).T[self.NAXIS], 3) else: return np.array(self.image) def close(self) -> None: self.image.close() def __exit__(self, *args) -> None: self.close() class Sampling: """ Sampling and subsampling routines """ SUBBED = False CACHE = [] def __init__(self, image: np.ndarray) -> None: self.image = image self.__X, self.__Y, *_ = image.shape def subSample(self, pix: int =3) -> np.ndarray: """ Subsample array in the context of Dokkum's cosmic ray removal algorithm """ id = np.ones((pix, pix)) self.image = np.kron(self.image, id) # update internal state self.SUBBED = True self.CACHE.append(pix) return self.image def blockAverage(self, pix: int =3) -> np.ndarray: """ Block average array after it's been subsampled """ assert self.SUBBED, 'Must subSample before blockAverage' # if the former is true, then the following should be as well assert not all(map(lambda x: x % pix, self.image.shape)),\ 'image dimensions not cleanly divisible by pix (?)' assert pix == self.CACHE[-1], 'Must use last subsample {0}, received {1}'\ .format(self.CACHE[-1], pix) # now block average---takes a few seconds avrgd = np.empty((self.__X, self.__Y)) for i in range(self.__X): for j in range(self.__Y): Xs = i * pix, i * pix + pix Ys = j * pix, j * pix + pix avrgd[i][j] = np.mean(self.image[slice(*Xs), slice(*Ys)]) self.image = avrgd self.SUBBED = False return self.image if __name__ == '__main__': import matplotlib.pyplot as plt with OpenImage(IMAGE) as im: s = Sampling(im) fig = plt.figure() ax1 = fig.add_subplot(122) ax1.imshow(s.subSample(SUB), interpolation='none', cmap='gray') ax1.set_title('Subsampled') ax2 = fig.add_subplot(121) ax2.imshow(s.blockAverage(SUB), interpolation='none', cmap='gray') ax2.set_title('Block averaged') plt.show()