Fast Image Pyramid Creation and Reconstruction in Python

image_pyramid.py

"""
image_pyramid.py
Fast Image Pyramid Creation and Reconstruction in Python

Copyright (c) 2014 Jack Doerner

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""

import numpy
import scipy.signal
from image_pyramid_filterTaps import filterTaps, filterTapsDirect

def buildImagePyramid( im ):
	# Returns a multi-scale pyramid of im. pyr is the pyramid concatenated as a
	# column vector while pind is the size of each level. im is expected to be
	# a grayscale two dimenionsal image in either single floating
	# point precision.
	#
	# Based on matlab code povided provided by Neal Wadhwa
	# in the supplimentary material of:
	#
	# Riesz Pyramids for Fast Phase-Based Video Magnification
	# Neal Wadhwa, Michael Rubinstein, Fredo Durand and William T. Freeman
	# Computational Photography (ICCP), 2014 IEEE International Conference on
	
	# Get the filter taps
	(bL, bH, t) = filterTaps()
	bL = bL.reshape((1, 1, bL.size))
	bL = 2 * bL # To make up for the energy lost during downsampling
	bH = bH.reshape((1, 1, bH.size))

	pyr = []
	pind = []

	while (numpy.amin(im.shape) >= 10): # Stop building the pyramid when the image is too small
		Y = numpy.zeros((im.shape[0], im.shape[1], bL.size)) 
		Y[:,:,0] = im
		# We apply the McClellan transform repeated to the image
		for k in range(1,bL.size):
			previousFiltered = Y[:,:,k-1]
			Y[:,:,k] = scipy.signal.convolve2d(previousFiltered, t, mode='same', boundary='symm')
		
		# Use Y to compute lo and highpass filtered image
		lopassedIm = numpy.sum(Y * bL,axis=2)
		hipassedIm = numpy.sum(Y * bH,axis=2)
		
		# Add highpassed image to the pyramid
		pyr.append(hipassedIm)
		pind.append(im.shape)

		# Downsample lowpassed image
		lopassedIm = lopassedIm[0:lopassedIm.shape[0]:2,0:lopassedIm.shape[1]:2]
		
		# Recurse on the lowpassed image
		im = lopassedIm

	# Add a residual level for the remaining low frequencies
	pyr.append(im)
	pind.append(im.shape)

	return (pyr, pind)

def reconstructImagePyramid( pyr, pind ):
	# Collapases a multi-scale pyramid of and returns the reconstructed image.
	# pyr is a column vector, in which each level of the pyramid is
	# concatenated, pind is the size of each level. 
	#
	# Based on matlab code povided provided by Neal Wadhwa
	#
	# Get the filter taps
	# Because we won't be simultaneously lowpassing/highpassing anything and 
	# most of the computational savings comes from the simultaneous application 
	# of the filters, we use the direct form of the filters rather the 
	# McClellan transform form

	(directL, directH) = filterTapsDirect()
	directL = 2*directL # To make up for the energy lost during downsampling

	nLevels = len(pind)
	lo = pyr[nLevels -1]
	for k in range(nLevels-1,0,-1):
		upsz = pind[k-1]
		# Upsample the lowest level
		lowest = numpy.zeros(upsz)
		lowest[::2,::2 ] = lo
		# Lowpass it with reflective boundary conditions
		lowest = scipy.signal.convolve2d(lowest, directL, mode='same', boundary='symm')
		# Get the next level    
		nextLevel = pyr[k-1]
		# Highpass the level and add it to lowest level to form a new lowest level
		lowest = lowest + scipy.signal.convolve2d(nextLevel, directH, mode='same', boundary='symm')
		lo = lowest
	return lo

image_pyramid_filterTaps.py

import numpy
import scipy.signal

hL =  numpy.array([-0.0209, -0.0219, 0.0900, 0.2723, 0.3611, 0.2723, 0.09, -0.0219, -0.0209])
hH =  numpy.array([0.0099, 0.0492, 0.1230, 0.2020, -0.7633, 0.2020, 0.1230, 0.0492, 0.0099])
# McClellan Transform
t = numpy.array([[1.0/8, 1.0/4, 1.0/8], [1.0/4, -1.0/2, 1.0/4], [1.0/8, 1.0/4, 1.0/8]])

def filterTaps():
	# Returns the lowpass and highpass filters specified in the supplementary
	# materials of "Riesz Pyramid for Fast Phase-Based Video Magnification"
	#
	# Based on matlab code povided provided by Neal Wadhwa
	# in the supplimentary material of:
	#
	# Riesz Pyramids for Fast Phase-Based Video Magnification
	# Neal Wadhwa, Michael Rubinstein, Fredo Durand and William T. Freeman
	# Computational Photography (ICCP), 2014 IEEE International Conference on
	#
	# hL and hH are the one dimenionsal filters designed by our optimization
	# bL and bH are the corresponding Chebysheve polynomials
	# t is the 3x3 McClellan transform matrix
	# directL and directH are the direct forms of the 2d filters

	# These are computed using Chebyshev polynomials, see filterToChebyCoeff
	# for more details
	bL = filterToChebyCoeff(hL)
	bH = filterToChebyCoeff(hH)

	return (bL, bH, t)

def filterTapsDirect():
	(bL, bH, t) = filterTaps()
	directL = filterTo2D(bL, t)
	directH = filterTo2D(bH, t)
	return (directL,directH)

 
# Returns the Chebyshev polynomial coefficients corresponding to a 
# symmetric 1D filter
def filterToChebyCoeff(taps):
	# taps should be an odd symmetric filter
	M = taps.size#
	N = int((M+1)/2) # Number of unique entries
	
	# Compute frequency response
	# g(1) + g(2)*cos(\omega) + g(3) \cos(2\omega) + ...
	g = numpy.empty((N,))
	g[0] = taps[N-1]
	g[1:N] = taps[N:]*2
	
	# Only need five polynomials for our filters
	ChebyshevPolynomial = []
	ChebyshevPolynomial.append(numpy.asarray([0, 0, 0, 0, 1]))
	ChebyshevPolynomial.append(numpy.asarray([0, 0, 0, 1, 0]))
	ChebyshevPolynomial.append(numpy.asarray([0, 0, 2, 0, -1]))
	ChebyshevPolynomial.append(numpy.asarray([0, 4, 0, -3, 0]))
	ChebyshevPolynomial.append(numpy.asarray([8, 0, -8, 0, 1]))
	
	# Now, convert frequency response to polynomials form
	# b(1) + b(2)\cos(\omega) + b(3) \cos(\omega)^2 + ...
	b = numpy.zeros((N,))
	for k in range(N):
		p = ChebyshevPolynomial[k]
		b = b + g[k]*p

	return b[::-1]

def filterTo2D(chebyshevPolyCoefficients, mcClellanTransform):
	ctr = chebyshevPolyCoefficients.size
	N = 2*ctr-1
	
	# Initial an impulse and then filter it
	X = numpy.zeros((N,N))
	X[ctr-1, ctr-1] = 1

	Y = numpy.zeros((X.shape[0], X.shape[1], chebyshevPolyCoefficients.size)) 
	Y[:,:,0] = X
	# We apply the McClellan transform repeated to the image
	for k in range(1,chebyshevPolyCoefficients.size):
		Y[:,:,k] = scipy.signal.convolve2d(Y[:,:,k-1],mcClellanTransform, mode='same', boundary='symm')

	# Take a linear combination of these to get the full 2D response
	chebyshevPolyCoefficients = chebyshevPolyCoefficients.reshape((1, 1, chebyshevPolyCoefficients.size))    
	return numpy.sum(Y * chebyshevPolyCoefficients,axis=2)

Thank you