an implemention of automated rig removal

optical_flow.py

"""
The rest of the example code that I ganked.

// 3. Draw arrows on second image to show movement
for (int x=0; x<frame2->width; x=x+10) {
  for (int y=0; y<frame2->height; y=y+10) {
    int vel_x_here = (int)cvGetReal2D( vel_x, y, x);
    int vel_y_here = (int)cvGetReal2D( vel_y, y, x);
    cvLine( frame2, cvPoint(x, y), cvPoint(x+vel_x_here,
y+vel_y_here), cvScalarAll(255));
  }
}
		
// 4. Display image with arrows and wait for key to be pressed
cvShowImage("result", frame2 );
cvWaitKey(0);
"""
from cv import CalcOpticalFlowLK, CreateImage, GetReal2D, IPL_DEPTH_32F, LoadImage, \
               CV_LOAD_IMAGE_GRAYSCALE
import csv

# load images
frame1 = LoadImage('/Users/kobold/Desktop/forest_start/forest 001.png', CV_LOAD_IMAGE_GRAYSCALE)
frame2 = LoadImage('/Users/kobold/Desktop/forest_start/forest 002.png', CV_LOAD_IMAGE_GRAYSCALE)

# calculate optical flow
vel_x = CreateImage((frame1.width, frame1.height), IPL_DEPTH_32F, 1)
vel_y = CreateImage((frame1.width, frame1.height), IPL_DEPTH_32F, 1)
CalcOpticalFlowLK(frame2, frame1, (5, 5), vel_x, vel_y)

# dump the pixel velocities
x_writer = csv.writer(open('vel_x.csv', 'wb'))
y_writer = csv.writer(open('vel_y.csv', 'wb'))
for y in xrange(frame1.height):
    x_writer.writerow([GetReal2D(vel_x, y, x) for x in xrange(frame1.width)])
    y_writer.writerow([GetReal2D(vel_y, y, x) for x in xrange(frame1.width)])

probabilities.py

"""
how do we update the occlusion data?
where does the weight field come from,

how is the estimate worked into the main estimation--there's overlap in terms

do we work from back to front

where the fuck does lambda_o in equation 6 come from?
"""

SIGMA_E = 0.01 # allowance of acceleration
SIGMA_V = 0.01 # allowance of acceleration
LAMBDA = 2.0
ALPHA = 0.0331 # ~ 3 * SIGMA_V according to section 3.4


#
# Helpers
#

def neighborhood(x_r):
    """The function S_n. Returns the eight points surrounding x_r."""
    x, y = x_r
    return [(x - 1, y - 1), (x, y - 1), (x + 1, y - 1),
            (x - 1, y),                 (x + 1, y),
            (x - 1, y + 1), (x, y + 1), (x + 1, y + 1)]

def lambda_(s, x_r):
    """Discourages 'smoothness' over too large a range."""
    return LAMBDA / magnitude(s - x_r)

#
# Probability model components
#

def p_l(x_r, w_n, w_n_1, d_h, I_n, I_n_1):
    """
    x_r: x_r
        a rig site
    w_n: w_n
        continuous matrix. 1 indicates data available, 0 indicates data missing.
        in the rig area w(x_r) = 0, i.e. this is the "not-rig" matrix
    w_n_1: w_n-1
        continuous matrix. 1 indicates data available, 0 indicates data missing.
        in the rig area w(x_r) = 0, i.e. this is the "not-rig" matrix
    d_h: d^h_n,n-1
        vector matrix estimating the motion of the hidden area
    I_n: I_n
        frame at n
    I_n_1: I_n-1
        frame at n-1
    """
    # the motion compensated site x_r + d^h_n,n-1(x_r)
    x_r_prime = x_r + d_h(x_r)

    temp = (1. / (2. * (SIGMA_E ** 2.))) *
           w_n(x_r) *
           w_n_1(x_r_prime) *
           ((I_n(x_r) - I_n_1(x_r_prime)) ** 2)
           
    return exp(-temp)


def p_t(x_r, occlusion, data_available, d_estimate, d_prev):
    """
    x_r:
        a rig site
    occlusion: o_n,n-1
        binary matrix. 1 indicates data that that point in the frame n does not
        exist in the frame n-1. 0 indicates no discontinuity
    data_available: w_n-1
        continuous matrix. 1 indicates data available, 0 indicates data missing.
        in the rig area w(x_r) = 0, i.e. this is the "not-rig" matrix
    d_estimate: d^h_n,n-1
        vector matrix estimating the motion of the hidden area
    d_prev: d_n-1,n-2
        vector matrix with the motion mapping from frame n-1 to frame n-2
    """
    # the motion compensated site x_r + d^h_n,n-1(x_r)
    x_r_prime = x_r + d_estimate(x_r)
    
    temp = (1. / (SIGMA_V ** 2.)) *
           (1. - occlusion(x_r)) *
           data_available(x_r_prime) *
           (magnitude(d_estimate(x_r) - d_prev(x_r_prime)) ** 2.)
    
    return exp(-temp)


def p_s(x_r, d_estimate, d):
    """
    x_r: rig site
    d_estimate: vector matrix estimating motion 
        of the hidden area, as in other functions
    d: matrix of motion vectors (maybe - section 4.1)

    """
    temp = sum(lambda_(s, x_r) * (magnitude(d_estimate(x_r) - d(s)) ** 2.)
               for s in neighborhood(x_r))
    
    return exp(-temp)
    

def p_so(x_r, occlusion):
    """
    x_r:
        a rig site
    occlusion: o_n,n-1
        binary matrix. 1 indicates data that that point in the frame n does not
        exist in the frame n-1. 0 indicates no discontinuity
    """
    occ = sum(lambda_(s, x_r) * magnitude(occlusion(x_r) - occlusion(s))
              for s in neighborhood(x_r))
    
    # TODO: is this sum term *actually* the sum of all the occlusion values?
    penalty = ALPHA * sum(occlusion)
    
    return exp(-acc) * exp(-penalty)