Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fast Image Pyramid Creation and Reconstruction in Python
"""
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
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.