adjust_phase.jl

function adjust_phase(phase_delta::ImagePyramid, corrected_phase_delta::ImagePyramid)
    adjusted_phase = ImagePyramid(phase_delta)

    update_subband!(adjusted_phase, 0, adjust_phase(subband(phase_delta, 0), subband(corrected_phase_delta, 0)))

    for l = 1:pyramid1.num_levels
        for o = 1:pyramid1.num_orientations
            update_subband!(adjusted_phase, l, adjust_phase(subband(phase_delta, l, orientation=o), subband(corrected_phase_delta, l, orientation=o)), orientation=o)
        end
    end

    return adjusted_phase
end

function adjust_phase(phase_delta::Array, corrected_phase_delta)
    adjusted_phase_delta = copy(phase_delta)

    for n = 1:length(adjusted_phase_delta)
        while adjusted_phase_delta[n] > (corrected_phase_delta[n] + π)
            adjusted_phase_delta[n] -= 2*π
        end

        while adjusted_phase_delta[n] < (corrected_phase_delta[n] - π)
            adjusted_phase_delta[n] += 2*π
        end
    end

    return adjusted_phase_delta
end

blend_and_interpolate.jl

function blend_and_interpolate(pyramid1::ImagePyramid, pyramid2::ImagePyramid, phase_delta::ImagePyramid, alpha)
    blended_pyramid = ImagePyramid(pyramid1)

    for l = 1:pyramid1.num_levels
        for o = 1:pyramid1.num_orientations
            new_band = ((1 - alpha) * abs(subband(pyramid1, l, orientation=o)) + alpha * abs(subband(pyramid2, l, orientation=o))) .* exp(complex(0,1) * (angle(subband(pyramid1, l, orientation=o)) + alpha * subband(phase_delta, l, orientation=o)))
            update_subband!(blended_pyramid, l, new_band, orientation=o)
        end
    end

    # blend low frequency residuals
    update_subband!(blended_pyramid, pyramid1.num_levels+1, subband(pyramid1, pyramid1.num_levels+1) * (1 - alpha) + alpha * subband(pyramid2, pyramid2.num_levels+1))

    # choose one high frequency residual
    if (alpha > 0.5)
        update_subband!(blended_pyramid, 0, subband(pyramid2, 0))
    else
        update_subband!(blended_pyramid, 0, subband(pyramid1, 0))
    end

    return blended_pyramid

end

generate_output.jl

phase_delta = phase_difference(pyramid1, pyramid2)
corrected_phase_delta = shift_correction(phase_delta, shift_limit=1)
adjusted_phase_delta = adjust_phase(phase_delta, corrected_phase_delta)

alpha = 0.5

println("Generating image with alpha $(alpha)")
newpyr = blend_and_interpolate(pyramid1, pyramid2, adjusted_phase_delta, alpha)
newim = toimage(newpyr) # conver the pyramid back to an image

newLabIm = zeros(im1)
newLabIm[:,:,1] = newim
# linearly interpolate the chrominance channels
newLabIm[:,:,2] = (1 - alpha) * im1[:,:,2] + alpha * im2[:,:,2]
newLabIm[:,:,3] = (1 - alpha) * im1[:,:,3] + alpha * im2[:,:,3]

# convert back to Julia Image class with the correct colorspace
newLabIm = convert(Image, newLabIm)
newLabIm.properties["colorspace"] = "Lab"
newLabIm = convert(Image{RGB}, newLabIm)
save("interpolated_frame_$(alpha).png", newLabIm')

loading.jl

println("Loading images")

im1 = load("frame_0.jpg")
im1 = convert(Image{Lab}, float32(im1))

im2 = load("frame_1.jpg")
im2 = convert(Image{Lab}, float32(im2))

im1 = convert(Array, separate(im1))
im2 = convert(Array, separate(im2))

L1 = im1[:,:,1]
L2 = im2[:,:,1]

println("Converting images to complex steerable pyramids")

pyramid1 = ImagePyramid(L1, ComplexSteerablePyramid(), scale=0.5^0.25)
pyramid2 = ImagePyramid(L2, ComplexSteerablePyramid(), scale=0.5^0.25)

phase_difference.jl

function phase_difference(pyramid1::ImagePyramid, pyramid2::ImagePyramid)
    phase_bands = Dict{Int, Union{Array,Dict{Int, Array}}}()

    phase_bands[0] = angle(subband(pyramid2, 0)) - angle(subband(pyramid1, 0))

    for l = 1:pyramid1.num_levels
        this_band = Dict{Int, Array}()

        for o = 1:pyramid1.num_orientations
            this_band[o] = angle(subband(pyramid2, l, orientation=o)) - angle(subband(pyramid1, l, orientation=o))
        end

        phase_bands[l] = copy(this_band)
    end

    return ImagePyramid(phase_bands, pyramid1.scale, PhasePyramid(), pyramid1.num_levels, pyramid1.num_orientations)
end

setup.jl

using Images
using ColorTypes
using Pyramids
using Interpolations

type PhasePyramid <: PyramidType
end

shift_correction.jl

function shift_correction(phi_delta::ImagePyramid; shift_limit=0.5)
    # create a new pyramid of the same shape and type as our input pyramid
    corrected_phase_delta = ImagePyramid(phi_delta)

    # iterate through each level and orientation of the complex steerable pyramid
    for level = corrected_phase_delta.num_levels:-1:1
        for orientation = 1:corrected_phase_delta.num_orientations
        
            # grab the phase delta for this level and orientation
            phi_l = subband(corrected_phase_delta, level, orientation=orientation)
            dims_l = size(phi_l)
            
            # set a limit on the maximum phase correction that is allowed (Eq. 11)
            phi_limit = (π / (phi_delta.scale^(phi_delta.num_levels - level))) * shift_limit

            # if this isn't the broadest spatial scale (which we must use as the "ground truth" phase), then correct phase
            if level < phi_delta.num_levels
                # get the next larger spatial scale
                phi_lp1 = subband(corrected_phase_delta, level+1, orientation=orientation)
                dims_lp1 = size(phi_lp1)
                
                # interpolate it up to the size of the spatial scale we're corrected (bilinear interpolation)
                phi_lp1_itp = interpolate(phi_lp1, BSpline(Linear()), OnGrid())
               
                # iterate over each pixel in the pyramid level
                for r = 1:dims_l[1]
                    for c = 1:dims_l[2]

                        r_lp1 = ((r-1)/(dims_l[1]-1)) * (dims_lp1[1]-1) + 1
                        c_lp1 = ((c-1)/(dims_l[2]-1)) * (dims_lp1[2]-1) + 1

                        # unwrap (subtract or add 2π until the phase of the level to be corrected matches the larger spatial scale level as closely as possible)
                        while phi_l[r,c] > (phi_lp1_itp[r_lp1,c_lp1]/phi_delta.scale + π)
                            phi_l[r,c] -= 2*π
                        end

                        while phi_l[r,c] < (phi_lp1_itp[r_lp1,c_lp1]/phi_delta.scale - π)
                            phi_l[r,c] += 2*π
                        end

                        # compute a confidence score (ambiguity of π -> least confident, ambiguity of 0 -> most confident)
                        confidence = atan2(sin(phi_l[r,c] - phi_lp1_itp[r_lp1,c_lp1]/phi_delta.scale), cos(phi_l[r,c] - phi_lp1_itp[r_lp1,c_lp1]/phi_delta.scale))

                        # if confidence is beyond a threshold, discard the current level's phase information
                        if (abs(confidence) > π/2)
                            phi_l[r,c] = phi_lp1_itp[r_lp1,c_lp1]/phi_delta.scale
                        end
                       
                        # if the magnitude of the phase is too large, discard the current level's phase information 
                        if abs(phi_l[r,c]) > phi_limit
                            phi_l[r,c] = phi_lp1_itp[r_lp1,c_lp1]/phi_delta.scale
                        end

                    end
                end
            else
                phi_l[abs(phi_l) .> phi_limit] = 0
            end

            # update the subband of the new pyramid
            update_subband!(corrected_phase_delta, level, phi_l, orientation=orientation)
        end
    end

    return corrected_phase_delta
end