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
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 =
        if (''.join(self.image.getbands()) == 'RGB'):
            return np.rot90(np.array(self.image).T[self.NAXIS], 3)
            return np.array(self.image)

    def close(self) -> None:

    def __exit__(self, *args) -> None:

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
        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')

        ax2 = fig.add_subplot(121)
        ax2.imshow(s.blockAverage(SUB), interpolation='none', cmap='gray')
        ax2.set_title('Block averaged')
And the result is as follows.

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.