Vladislav Samsonov
9 years ago
parent
54e746bebd
commit
7a18e78865
6 changed files with 149 additions and 251 deletions
@ -1,143 +0,0 @@ |
||||
/*M///////////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
|
||||
//
|
||||
// By downloading, copying, installing or using the software you agree to this license.
|
||||
// If you do not agree to this license, do not download, install,
|
||||
// copy or use the software.
|
||||
//
|
||||
//
|
||||
// License Agreement
|
||||
// For Open Source Computer Vision Library
|
||||
//
|
||||
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
|
||||
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
|
||||
// Third party copyrights are property of their respective owners.
|
||||
//
|
||||
// Redistribution and use in source and binary forms, with or without modification,
|
||||
// are permitted provided that the following conditions are met:
|
||||
//
|
||||
// * Redistribution's of source code must retain the above copyright notice,
|
||||
// this list of conditions and the following disclaimer.
|
||||
//
|
||||
// * Redistribution's in binary form must reproduce the above copyright notice,
|
||||
// this list of conditions and the following disclaimer in the documentation
|
||||
// and/or other materials provided with the distribution.
|
||||
//
|
||||
// * The name of the copyright holders may not be used to endorse or promote products
|
||||
// derived from this software without specific prior written permission.
|
||||
//
|
||||
// This software is provided by the copyright holders and contributors "as is" and
|
||||
// any express or implied warranties, including, but not limited to, the implied
|
||||
// warranties of merchantability and fitness for a particular purpose are disclaimed.
|
||||
// In no event shall the Intel Corporation or contributors be liable for any direct,
|
||||
// indirect, incidental, special, exemplary, or consequential damages
|
||||
// (including, but not limited to, procurement of substitute goods or services;
|
||||
// loss of use, data, or profits; or business interruption) however caused
|
||||
// and on any theory of liability, whether in contract, strict liability,
|
||||
// or tort (including negligence or otherwise) arising in any way out of
|
||||
// the use of this software, even if advised of the possibility of such damage.
|
||||
//
|
||||
//M*/
|
||||
|
||||
#include "precomp.hpp" |
||||
|
||||
namespace cv |
||||
{ |
||||
namespace optflow |
||||
{ |
||||
|
||||
class OpticalFlowBlockMatching : public DenseOpticalFlow |
||||
{ |
||||
protected: |
||||
int windowSize; |
||||
int blockSize; |
||||
|
||||
inline float submatrixAbsDiff( int x0, int y0, const Mat &I0, int x1, int y1, const Mat &I1 ) const; |
||||
|
||||
public: |
||||
OpticalFlowBlockMatching() : windowSize( 3 ), blockSize( 8 ){}; |
||||
|
||||
void calc( InputArray I0, InputArray I1, InputOutputArray flow ); |
||||
void collectGarbage(); |
||||
}; |
||||
|
||||
inline float OpticalFlowBlockMatching::submatrixAbsDiff( int x0, int y0, const Mat &I0, int x1, int y1, |
||||
const Mat &I1 ) const |
||||
{ |
||||
float error = 0; |
||||
const Size size = I0.size(); |
||||
for ( int i = -windowSize; i <= windowSize; ++i ) |
||||
{ |
||||
if ( i + y0 < 0 || i + y0 >= size.height || i + y1 < 0 || i + y1 >= size.height ) |
||||
{ |
||||
error += 1; |
||||
continue; |
||||
} |
||||
const Vec3f *I0X = I0.ptr<Vec3f>( i + y0 ); |
||||
const Vec3f *I1X = I1.ptr<Vec3f>( i + y1 ); |
||||
for ( int j = -windowSize; j <= windowSize; ++j ) |
||||
{ |
||||
if ( j + x0 < 0 || j + x0 >= size.width || j + x1 < 0 || j + x1 >= size.width ) |
||||
{ |
||||
error += 1; |
||||
continue; |
||||
} |
||||
const Vec3f diff = I0X[j + x0] - I1X[j + x1]; |
||||
error += abs( diff[0] ); |
||||
error += abs( diff[1] ); |
||||
error += abs( diff[2] ); |
||||
} |
||||
} |
||||
return error; |
||||
} |
||||
|
||||
void OpticalFlowBlockMatching::calc( InputArray I0, InputArray I1, InputOutputArray flow_out ) |
||||
{ |
||||
CV_Assert( I0.channels() == 3 ); |
||||
CV_Assert( I1.channels() == 3 ); |
||||
Size size = I0.size(); |
||||
CV_Assert( size == I1.size() ); |
||||
|
||||
flow_out.create( size, CV_32FC2 ); |
||||
Mat flow = flow_out.getMat(); |
||||
Mat from = I0.getMat(); |
||||
Mat to = I1.getMat(); |
||||
|
||||
from.convertTo( from, CV_32FC3, 1.0 / 255.0 ); |
||||
to.convertTo( to, CV_32FC3, 1.0 / 255.0 ); |
||||
|
||||
const float distNormalize = blockSize * sqrt( 2 ); |
||||
|
||||
for ( int y0 = 0; y0 < size.height; ++y0 ) |
||||
{ |
||||
Vec2f *flowX = flow.ptr<Vec2f>( y0 ); |
||||
const int yEnd = std::min( size.height - 1, y0 + blockSize ); |
||||
|
||||
for ( int x0 = 0; x0 < size.width; ++x0 ) |
||||
{ |
||||
float minDiff = 1e10; |
||||
Vec2f du( 0, 0 ); |
||||
const int xEnd = std::min( size.width - 1, x0 + blockSize ); |
||||
for ( int y1 = std::max( 0, y0 - blockSize ); y1 <= yEnd; ++y1 ) |
||||
for ( int x1 = std::max( 0, x0 - blockSize ); x1 <= xEnd; ++x1 ) |
||||
{ |
||||
const float distance = sqrt( ( x0 - x1 ) * ( x0 - x1 ) + ( y0 - y1 ) * ( y0 - y1 ) ) / distNormalize; |
||||
const float kernel = 1.0 + 0.5 * distance; |
||||
const float diff = kernel * submatrixAbsDiff( x0, y0, from, x1, y1, to ); |
||||
if ( diff < minDiff ) |
||||
{ |
||||
minDiff = diff; |
||||
du = Vec2f( ( x1 - x0 ), ( y1 - y0 ) ); |
||||
} |
||||
} |
||||
flowX[x0] = du; |
||||
} |
||||
} |
||||
} |
||||
|
||||
void OpticalFlowBlockMatching::collectGarbage() {} |
||||
|
||||
Ptr<DenseOpticalFlow> createOptFlow_BlockMatching() { return makePtr<OpticalFlowBlockMatching>(); } |
||||
} |
||||
} |
||||
@ -1,123 +1,166 @@ |
||||
import os, sys |
||||
#!/usr/bin/env python |
||||
|
||||
import os |
||||
import sys |
||||
import numpy as np |
||||
import cv2 |
||||
import struct |
||||
import argparse |
||||
from math import sqrt |
||||
|
||||
basis_size = (14, 18) |
||||
lambd = 0.2 |
||||
argparser = argparse.ArgumentParser( |
||||
description='''Use this script to generate prior for using with PCAFlow. |
||||
Basis size here must match corresponding parameter in the PCAFlow. |
||||
Gamma should be selected experimentally.''') |
||||
|
||||
argparser.add_argument('-f', |
||||
'--files', |
||||
nargs='+', |
||||
help='List of optical flow .flo files for learning', |
||||
required=True) |
||||
argparser.add_argument('-o', |
||||
'--output', |
||||
help='Output file for prior', |
||||
required=True) |
||||
argparser.add_argument('--width', |
||||
type=int, |
||||
help='Size of the basis first dimension', |
||||
required=True, |
||||
default=18) |
||||
argparser.add_argument('--height', |
||||
type=int, |
||||
help='Size of the basis second dimension', |
||||
required=True, |
||||
default=14) |
||||
argparser.add_argument( |
||||
'-g', |
||||
'--gamma', |
||||
type=float, |
||||
help='Amount of regularization. The greater this parameter, the bigger will be an impact of the regularization.', |
||||
required=True) |
||||
args = argparser.parse_args() |
||||
|
||||
basis_size = (args.height, args.width) |
||||
gamma = args.gamma |
||||
|
||||
|
||||
def find_flo(pp): |
||||
f = [] |
||||
for p in pp: |
||||
if os.path.isfile(p): |
||||
f.append(p) |
||||
else: |
||||
for root, subdirs, files in os.walk(p): |
||||
f += map(lambda x: os.path.join(root, x), filter(lambda x: x.split('.')[-1] == 'flo', files)) |
||||
return list(set(f)) |
||||
f = [] |
||||
for p in pp: |
||||
if os.path.isfile(p): |
||||
f.append(p) |
||||
else: |
||||
for root, subdirs, files in os.walk(p): |
||||
f += map(lambda x: os.path.join(root, x), |
||||
filter(lambda x: x.split('.')[-1] == 'flo', files)) |
||||
return list(set(f)) |
||||
|
||||
|
||||
def load_flo(flo): |
||||
with open(flo) as f: |
||||
magic = np.fromfile(f, np.float32, count=1) |
||||
if 202021.25 != magic: |
||||
print('Magic number incorrect. Invalid .flo file') |
||||
else: |
||||
w = np.fromfile(f, np.int32, count=1) |
||||
h = np.fromfile(f, np.int32, count=1) |
||||
print('Reading %dx%d flo file %s' % (w, h, flo)) |
||||
data = np.fromfile(f, np.float32, count=2*w*h) |
||||
# Reshape data into 3D array (columns, rows, bands) |
||||
flow = np.resize(data, (h, w, 2)) |
||||
return flow[:,:,0], flow[:,:,1] |
||||
magic = np.fromfile(f, np.float32, count=1)[0] |
||||
if 202021.25 != magic: |
||||
print('Magic number incorrect. Invalid .flo file') |
||||
else: |
||||
w = np.fromfile(f, np.int32, count=1)[0] |
||||
h = np.fromfile(f, np.int32, count=1)[0] |
||||
print('Reading %dx%d flo file %s' % (w, h, flo)) |
||||
data = np.fromfile(f, np.float32, count=2 * w * h) |
||||
# Reshape data into 3D array (columns, rows, bands) |
||||
flow = np.reshape(data, (h, w, 2)) |
||||
return flow[:, :, 0], flow[:, :, 1] |
||||
|
||||
|
||||
def get_w(m): |
||||
s = m.shape |
||||
w = cv2.dct(m) |
||||
w *= 2.0 / sqrt(s[0] * s[1]) |
||||
#w[0,0] *= 0.5 |
||||
w[:,0] *= sqrt(0.5) |
||||
w[0,:] *= sqrt(0.5) |
||||
w = w[0:basis_size[0],0:basis_size[1]].transpose().flatten() |
||||
return w |
||||
s = m.shape |
||||
w = cv2.dct(m) |
||||
w *= 2.0 / sqrt(s[0] * s[1]) |
||||
#w[0,0] *= 0.5 |
||||
w[:, 0] *= sqrt(0.5) |
||||
w[0, :] *= sqrt(0.5) |
||||
w = w[0:basis_size[0], 0:basis_size[1]].transpose().flatten() |
||||
return w |
||||
|
||||
|
||||
w1 = [] |
||||
w2 = [] |
||||
|
||||
for flo in find_flo(sys.argv[1:]): |
||||
x,y = load_flo(flo) |
||||
w1.append(get_w(x)) |
||||
w2.append(get_w(y)) |
||||
for flo in find_flo(args.files): |
||||
x, y = load_flo(flo) |
||||
w1.append(get_w(x)) |
||||
w2.append(get_w(y)) |
||||
|
||||
w1mean = sum(w1) / len(w1) |
||||
w2mean = sum(w2) / len(w2) |
||||
|
||||
for i in xrange(len(w1)): |
||||
w1[i] -= w1mean |
||||
w1[i] -= w1mean |
||||
for i in xrange(len(w2)): |
||||
w2[i] -= w2mean |
||||
w2[i] -= w2mean |
||||
|
||||
Q1 = sum([w1[i].reshape(-1,1).dot(w1[i].reshape(1,-1)) for i in xrange(len(w1))]) / len(w1) |
||||
Q2 = sum([w2[i].reshape(-1,1).dot(w2[i].reshape(1,-1)) for i in xrange(len(w2))]) / len(w2) |
||||
Q1 = sum([w1[i].reshape(-1, 1).dot(w1[i].reshape(1, -1)) |
||||
for i in xrange(len(w1))]) / len(w1) |
||||
Q2 = sum([w2[i].reshape(-1, 1).dot(w2[i].reshape(1, -1)) |
||||
for i in xrange(len(w2))]) / len(w2) |
||||
Q1 = np.matrix(Q1) |
||||
Q2 = np.matrix(Q2) |
||||
|
||||
if len(w1) > 1: |
||||
while True: |
||||
try: |
||||
L1 = np.linalg.cholesky(Q1) |
||||
break |
||||
except np.linalg.linalg.LinAlgError: |
||||
mev = min(np.linalg.eig(Q1)[0]).real |
||||
assert(mev < 0) |
||||
print('Q1', mev) |
||||
if -mev < 1e-6: |
||||
mev = -1e-6 |
||||
Q1 += (-mev*1.000001) * np.identity(Q1.shape[0]) |
||||
|
||||
while True: |
||||
try: |
||||
L2 = np.linalg.cholesky(Q2) |
||||
break |
||||
except np.linalg.linalg.LinAlgError: |
||||
mev = min(np.linalg.eig(Q2)[0]).real |
||||
assert(mev < 0) |
||||
print('Q2', mev) |
||||
if -mev < 1e-6: |
||||
mev = -1e-6 |
||||
Q2 += (-mev*1.000001) * np.identity(Q2.shape[0]) |
||||
while True: |
||||
try: |
||||
L1 = np.linalg.cholesky(Q1) |
||||
break |
||||
except np.linalg.linalg.LinAlgError: |
||||
mev = min(np.linalg.eig(Q1)[0]).real |
||||
assert (mev < 0) |
||||
print('Q1', mev) |
||||
if -mev < 1e-6: |
||||
mev = -1e-6 |
||||
Q1 += (-mev * 1.000001) * np.identity(Q1.shape[0]) |
||||
|
||||
while True: |
||||
try: |
||||
L2 = np.linalg.cholesky(Q2) |
||||
break |
||||
except np.linalg.linalg.LinAlgError: |
||||
mev = min(np.linalg.eig(Q2)[0]).real |
||||
assert (mev < 0) |
||||
print('Q2', mev) |
||||
if -mev < 1e-6: |
||||
mev = -1e-6 |
||||
Q2 += (-mev * 1.000001) * np.identity(Q2.shape[0]) |
||||
else: |
||||
L1 = np.identity(Q1.shape[0]) |
||||
L2 = np.identity(Q2.shape[0]) |
||||
L1 = np.identity(Q1.shape[0]) |
||||
L2 = np.identity(Q2.shape[0]) |
||||
|
||||
L1 = np.linalg.inv(L1) * lambd |
||||
L2 = np.linalg.inv(L2) * lambd |
||||
L1 = np.linalg.inv(L1) * gamma |
||||
L2 = np.linalg.inv(L2) * gamma |
||||
|
||||
assert(L1.shape == L2.shape) |
||||
assert(L1.shape[0] == L1.shape[1]) |
||||
assert (L1.shape == L2.shape) |
||||
assert (L1.shape[0] == L1.shape[1]) |
||||
|
||||
f = open('cov.dat', 'wb') |
||||
f = open(args.output, 'wb') |
||||
|
||||
f.write(struct.pack('I', L1.shape[0])) |
||||
f.write(struct.pack('I', L1.shape[1])) |
||||
|
||||
for i in xrange(L1.shape[0]): |
||||
for j in xrange(L1.shape[1]): |
||||
f.write(struct.pack('f', L1[i,j])) |
||||
for j in xrange(L1.shape[1]): |
||||
f.write(struct.pack('f', L1[i, j])) |
||||
|
||||
for i in xrange(L2.shape[0]): |
||||
for j in xrange(L2.shape[1]): |
||||
f.write(struct.pack('f', L2[i,j])) |
||||
for j in xrange(L2.shape[1]): |
||||
f.write(struct.pack('f', L2[i, j])) |
||||
|
||||
b1 = L1.dot(w1mean.reshape(-1,1)) |
||||
b2 = L2.dot(w2mean.reshape(-1,1)) |
||||
b1 = L1.dot(w1mean.reshape(-1, 1)) |
||||
b2 = L2.dot(w2mean.reshape(-1, 1)) |
||||
|
||||
assert(L1.shape[0] == b1.shape[0]) |
||||
assert (L1.shape[0] == b1.shape[0]) |
||||
|
||||
for i in xrange(b1.shape[0]): |
||||
f.write(struct.pack('f', b1[i,0])) |
||||
f.write(struct.pack('f', b1[i, 0])) |
||||
|
||||
for i in xrange(b2.shape[0]): |
||||
f.write(struct.pack('f', b2[i,0])) |
||||
f.write(struct.pack('f', b2[i, 0])) |
||||
|
||||
f.close() |
||||
f.close() |
||||
|
||||
Loading…
Reference in new issue