numpy/scipy, loop over subarrays

Question

Lately I've been doing a lot of processing on 8x8 blocks of image-data. Standard approach has been to use nested for-loops to extract the blocks, e.g.

for y in xrange(0,height,8):
    for x in xrange(0,width,8):
        d = image_data[y:y+8,x:x+8]
        # further processing on the 8x8-block

I can't help to wonder if there is a way to vectorize this operation or another approach using numpy/scipy that I can use instead? An iterator of some kind?

A MWE¹:

#!/usr/bin/env python

import sys
import numpy as np
from scipy.fftpack import dct, idct
import scipy.misc
import matplotlib.pyplot as plt

def dctdemo(coeffs=1):
    unzig = np.array([
         0,  1,  8, 16,  9,  2,  3, 10,
        17, 24, 32, 25, 18, 11,  4, 5,
        12, 19, 26, 33, 40, 48, 41, 34,
        27, 20, 13,  6,  7, 14, 21, 28,
        35, 42, 49, 56, 57, 50, 43, 36,
        29, 22, 15, 23, 30, 37, 44, 51,
        58, 59, 52, 45, 38, 31, 39, 46,
        53, 60, 61, 54, 47, 55, 62, 63])

    lena = scipy.misc.lena()
    width, height = lena.shape

    # reconstructed
    rec = np.zeros(lena.shape, dtype=np.int64)

    # Can this part be vectorized?
    for y in xrange(0,height,8):
        for x in xrange(0,width,8):
            d = lena[y:y+8,x:x+8].astype(np.float)
            D = dct(dct(d.T, norm='ortho').T, norm='ortho').reshape(64)
            Q = np.zeros(64, dtype=np.float)
            Q[unzig[:coeffs]] = D[unzig[:coeffs]]
            Q = Q.reshape([8,8])
            q = np.round(idct(idct(Q.T, norm='ortho').T, norm='ortho'))
            rec[y:y+8,x:x+8] = q.astype(np.int64)

    plt.imshow(rec, cmap='gray')
    plt.show()

if __name__ == '__main__':
    try:
        c = int(sys.argv[1])
    except ValueError:
        sys.exit()
    else:
        if 1 <= int(sys.argv[1]) <= 64:
            dctdemo(int(sys.argv[1]))

Footnotes:

Actual application: https://github.com/figgis/dctdemo

Currently trying to fiddle it so that the function I suggested works with your code... — YXD
– YXD, Commented May 7, 2014 at 20:58

YXD · Accepted Answer · 2014-05-07 21:02:54Z

5

There's a function view_as_windows for this in Scikit Image

http://scikit-image.org/docs/dev/api/skimage.util.html#view-as-windows

Unfortunately I will have to finish this answer another time, but you can grab the windows in a form that you can pass to dct with:

from skimage.util import view_as_windows
# your code...
d = view_as_windows(lena.astype(np.float), (8, 8)).reshape(-1, 8, 8)
dct(d, axis=0)

edited May 7, 2014 at 21:02

answered May 7, 2014 at 20:43

YXD

32.6k15 gold badges79 silver badges117 bronze badges

Sign up to request clarification or add additional context in comments.

Comments

eickenberg · Accepted Answer · 2014-05-07 22:16:01Z

3

There is a function called extract_patches in the scikit-learn feature extraction routines. You need to specify a patch_size and an extraction_step. The result will be a view on your image as patches, which may overlap. The resulting array is 4D, the first 2 index the patch, and the last two index the pixels of the patch. Try this

from sklearn.feature_extraction.image import extract_patches
patches = extract_patches(image_data, patch_size=(8, 8), extraction_step=(4, 4))

This gives (8, 8) size patches that overlap by half.

Note that up until now this uses no extra memory, because it is implemented using stride tricks. You can force a copy by reshaping

patches = patches.reshape(-1, 8, 8)

which will basically yield a list of patches.

answered May 7, 2014 at 22:16

eickenberg

14.4k1 gold badge52 silver badges52 bronze badges

1 Comment

Jred Over a year ago

Is there a way to retain the partial patches at the edge of an image?

Collectives™ on Stack Overflow

numpy/scipy, loop over subarrays

2 Answers 2

Comments

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related