Eulerian Video Magnification

eulerian_vid_mag.py

#!/usr/bin/env python

import cv
import os
import sys
from skimage import io
from skimage.transform import pyramid_expand
from skimage import img_as_ubyte
import numpy as np

from yiq import yiq2rgb, rgb2yiq

ALPHA=25
MAX_FRAMES=301

# we need to hack the skimage version of pyramid laplacian to
# * bail out when it can't maintain its upscale factor constantly
# * return Ln = gn

def pyramid_laplacian(image, max_layer=-1, downscale=2, sigma=None, order=1,
                      mode='reflect', cval=0):
    """Yield images of the laplacian pyramid formed by the input image.

    Each layer contains the difference between the downsampled and the
    downsampled, smoothed image::

        layer = resize(prev_layer) - smooth(resize(prev_layer))

    Note that the first image of the pyramid will be the difference between the
    original, unscaled image and its smoothed version. The total number of
    images is `max_layer + 1`. In case all layers are computed, the last image
    is either a one-pixel image or the image where the reduction does not
    change its shape.

    Parameters
    ----------
    image : array
        Input image.
    max_layer : int
        Number of layers for the pyramid. 0th layer is the original image.
        Default is -1 which builds all possible layers.
    downscale : float, optional
        Downscale factor.
    sigma : float, optional
        Sigma for gaussian filter. Default is `2 * downscale / 6.0` which
        corresponds to a filter mask twice the size of the scale factor that
        covers more than 99% of the gaussian distribution.
    order : int, optional
        Order of splines used in interpolation of downsampling. See
        `scipy.ndimage.map_coordinates` for detail.
    mode :  {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional
        The mode parameter determines how the array borders are handled, where
        cval is the value when mode is equal to 'constant'.
    cval : float, optional
        Value to fill past edges of input if mode is 'constant'.

    Returns
    -------
    pyramid : generator
        Generator yielding pyramid layers as float images.

    References
    ----------
    .. [1] http://web.mit.edu/persci/people/adelson/pub_pdfs/pyramid83.pdf
    .. [2] http://sepwww.stanford.edu/~morgan/texturematch/paper_html/node3.html

    """

    from skimage.transform.pyramids import _smooth,_check_factor, img_as_float, resize
    import math

    _check_factor(downscale)

    # cast to float for consistent data type in pyramid
    image = img_as_float(image)

    if sigma is None:
        # automatically determine sigma which covers > 99% of distribution
        sigma = 2 * downscale / 6.0

    layer = 0
    rows = image.shape[0]
    cols = image.shape[1]

    smoothed_image = _smooth(image, sigma, mode, cval)
    yield image - smoothed_image

    # build downsampled images until max_layer is reached or downscale process
    # does not change image size
    while layer != max_layer:
        layer += 1

        out_rows = math.ceil(rows / float(downscale))
        out_cols = math.ceil(cols / float(downscale))

        if (out_rows/float(rows) != out_cols/float(cols) and
            out_rows/float(rows) != 1):
            yield smoothed_image
            break

        resized_image = resize(smoothed_image, (out_rows, out_cols),
                               order=order, mode=mode, cval=cval)
        smoothed_image = _smooth(resized_image, sigma, mode, cval)

        prev_rows = rows
        prev_cols = cols
        rows = resized_image.shape[0]
        cols = resized_image.shape[1]

        # no change to previous pyramid layer
        if prev_rows == rows and prev_cols == cols:
            break

        yield resized_image - smoothed_image

def generatePyramidSet(frameSet, frameCount):
    zeroFrame = frameSet[0]

    laplacian_gen = pyramid_laplacian(zeroFrame)
    pyramidSet = []

    levelCount = 0
    for l in laplacian_gen:
        levelCount = levelCount + 1
        pyramidSet.append(np.empty(l.shape + (frameCount,),
                                   dtype='float64'))

    print "generating Laplacian pyramids"
    for frameNumber in xrange(frameCount):
        print frameNumber, '/', frameCount

        frame = frameSet[frameNumber]
        laplacian_gen = pyramid_laplacian(frame)

        level = 0
        for l in laplacian_gen:
            pyramidSet[level][:,:,:,frameNumber] = l
            level = level + 1

    return pyramidSet

def filterAndAmplify(frameRate, pyramidSet):

    # courtesy http://stackoverflow.com/a/12233959/122022
    from scipy.signal import butter, lfilter
    from scipy.signal import freqz

    def butter_bandpass(lowcut, highcut, fs, order=3):
        nyq = 0.5 * fs
        low = lowcut / nyq
        high = highcut / nyq
        b, a = butter(order, [low, high], btype='band')
        return b, a

    order = 3
    fs = frameRate

    # consider 30 to 180 bpm
    # 0.5 to 3Hz
    lowcut = 0.5
    highcut = 3

    # consider 50 to 60 bpm
    lowcut=0.83
    highcut=1

    # for each level of the pyramid, and for each coordinate, filter the signal
    levelCount = len(pyramidSet)

    for level in xrange(4,levelCount-1):
        pSet = pyramidSet[level]
        b, a = butter_bandpass(lowcut, highcut, fs, order=3)
        for x in xrange(pSet.shape[0]):
            print level, levelCount, x, pSet.shape[0]
            for y in xrange(pSet.shape[1]):
                for z in xrange(pSet.shape[2]):
                    p_filtered = lfilter(b, a, pSet[x,y,z,:])
                    pSet[x,y,z,:] = pSet[x,y,z,:] + ALPHA*p_filtered

def saveVideo(frameRate, frameCount, frameSet, pyramidSet):

    print "reconstructing video"

    zeroFrame = frameSet[0]

    levelCount = len(pyramidSet)

    for frameNumber in xrange(frameCount):
        print frameNumber
        # g[-2] = g[-2] + g[-1]
        pyramidSet[-2][:,:,:,frameNumber] = pyramidSet[-2][:,:,:,frameNumber] + \
                                             pyramidSet[-1][:,:,:,frameNumber]

        for i in xrange(levelCount-3, -1, -1):
            upscale = pyramidSet[i].shape[0]/float(pyramidSet[i+1].shape[0])
            pyramidSet[i][:,:,:,frameNumber] = pyramidSet[i][:,:,:,frameNumber] + \
                pyramid_expand(pyramidSet[i+1][:,:,:,frameNumber], upscale=upscale)


        frame = pyramidSet[0][:,:,:,frameNumber]

        clippedFrame = frame.clip(0,1)

        io.imsave('out_frames/amplify_%4.4d.png' % frameNumber, clippedFrame)

file = 'capture.avi'

v = io.Video(file)
frameSet    = v.get_collection()
frameCount  = min(len(frameSet), MAX_FRAMES)
videoTimeMS = v.duration()
frameRate   = float(len(frameSet))/videoTimeMS*1000


if __name__ == "__main__":
    pyramidSet = generatePyramidSet(frameSet, frameCount)

    print "applying ideal bandpass filters"
    filterAndAmplify(frameRate, pyramidSet)

    saveVideo(frameRate, frameCount, frameSet, pyramidSet)