### 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 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).
#! /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.imshow(s.subSample(SUB), interpolation='none', cmap='gray')
ax1.set_title('Subsampled')