boofcv.alg.flow.BroxWarpingSpacial Maven / Gradle / Ivy
/*
* Copyright (c) 2011-2016, Peter Abeles. All Rights Reserved.
*
* This file is part of BoofCV (http://boofcv.org).
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package boofcv.alg.flow;
import boofcv.abst.filter.derivative.ImageGradient;
import boofcv.abst.filter.derivative.ImageHessian;
import boofcv.alg.interpolate.InterpolatePixelS;
import boofcv.alg.misc.ImageMiscOps;
import boofcv.alg.misc.PixelMath;
import boofcv.factory.filter.derivative.FactoryDerivative;
import boofcv.struct.image.GrayF32;
import boofcv.struct.image.ImageGray;
import boofcv.struct.pyramid.ImagePyramid;
import java.util.Arrays;
/**
*
* Dense optical flow which adheres to a brightness constancy assumption, a gradient constancy
* assumption, and a discontinuity-preserving spatio-temporal smoothness constraint. Based on the
* work of Brox [2] with implementation details taken from [1].
*
*
*
* -
* Javier Sánchez Pérez, Nelson Monzón López, and Agustín Salgado de la Nuez, Robust Optical Flow Estimation,
* Image Processing On Line, 3 (2013), pp. 252–270. http://dx.doi.org/10.5201/ipol.2013.21
* -
* Thomas Brox, Andr´es Bruhn, Nils Papenberg, and Joachim Weickert. High accuracy optical
* flow estimation based on a theory for warping. In T. Pajdla and J. Matas, editors, European
* Conference on Computer Vision (ECCV), volume 3024 of Lecture Notes in Computer Science,
* pages 25–36, Prague, Czech Republic, May 2004. Springer.
*
*
*
* @author Peter Abeles
*/
public class BroxWarpingSpacial extends DenseFlowPyramidBase {
// regularization term
private static final double EPSILON = 0.001;
// brightness error weighting factor
protected float alpha;
// gradient error weighting factor
protected float gamma;
// relaxation parameter for SOR 0 < w < 2. Recommended default is 1.9
private float SOR_RELAXATION;
// number of iterations for inner and outer loops
private int numOuter;
private int numInner;
// maximum number of iterations for SOR
private int maxIterationsSor;
// convergence tolerance for SOR
private float convergeTolerance;
// derivative of first image
private GrayF32 deriv1X = new GrayF32(1,1);
private GrayF32 deriv1Y = new GrayF32(1,1);
// derivatives of second image
private GrayF32 deriv2X = new GrayF32(1,1);
private GrayF32 deriv2Y = new GrayF32(1,1);
private GrayF32 deriv2XX = new GrayF32(1,1);
private GrayF32 deriv2YY = new GrayF32(1,1);
private GrayF32 deriv2XY = new GrayF32(1,1);
private ImageGradient gradient = FactoryDerivative.three(GrayF32.class, GrayF32.class);
private ImageHessian hessian = FactoryDerivative.hessianThree(GrayF32.class);
// flow estimation at the start of the iteration
protected GrayF32 flowU = new GrayF32(1,1); // flow along x-axis
protected GrayF32 flowV = new GrayF32(1,1); // flow along y-axis
// storage for the warped flow
protected GrayF32 warpImage2 = new GrayF32(1,1);
protected GrayF32 warpDeriv2X = new GrayF32(1,1);
protected GrayF32 warpDeriv2Y = new GrayF32(1,1);
protected GrayF32 warpDeriv2XX = new GrayF32(1,1);
protected GrayF32 warpDeriv2YY = new GrayF32(1,1);
protected GrayF32 warpDeriv2XY = new GrayF32(1,1);
// derivative of flow
protected GrayF32 derivFlowUX = new GrayF32(1,1);
protected GrayF32 derivFlowUY = new GrayF32(1,1);
protected GrayF32 derivFlowVX = new GrayF32(1,1);
protected GrayF32 derivFlowVY = new GrayF32(1,1);
protected GrayF32 psiSmooth = new GrayF32(1,1);
protected GrayF32 psiData = new GrayF32(1,1);
protected GrayF32 psiGradient = new GrayF32(1,1);
// divergence for u,v,d
protected GrayF32 divU = new GrayF32(1,1);
protected GrayF32 divV = new GrayF32(1,1);
protected GrayF32 divD = new GrayF32(1,1);
// motion increments for optical flow
protected GrayF32 du = new GrayF32(1,1);
protected GrayF32 dv = new GrayF32(1,1);
/**
* Configures flow estimation
*
* @param config Configuration parameters
* @param interp Interpolation for image flow between image layers and warping. Overrides selection in config.
*/
public BroxWarpingSpacial(ConfigBroxWarping config, InterpolatePixelS interp)
{
super(config.pyrScale,config.pyrSigma,config.pyrMaxLayers,interp);
this.alpha = config.alpha;
this.gamma = config.gamma;
this.SOR_RELAXATION = config.SOR_RELAXATION;
this.numOuter = config.numOuter;
this.numInner = config.numInner;
this.maxIterationsSor = config.maxIterationsSor;
this.convergeTolerance = config.convergeToleranceSor;
}
/**
* Computes dense optical flow from the provided image pyramid. Image gradient for each layer should be
* computed directly from the layer images.
*
* @param image1 Pyramid of first image
* @param image2 Pyramid of second image
*/
public void process(ImagePyramid image1 , ImagePyramid image2 )
{
// Process the pyramid from low resolution to high resolution
boolean first = true;
for( int i = image1.getNumLayers()-1; i >= 0; i-- ) {
GrayF32 layer1 = image1.getLayer(i);
GrayF32 layer2 = image2.getLayer(i);
resizeForLayer(layer1.width,layer2.height);
// compute image derivatives
gradient.process(layer1,deriv1X,deriv1Y);
gradient.process(layer2,deriv2X,deriv2Y);
hessian.process(deriv2X,deriv2Y,deriv2XX,deriv2YY,deriv2XY);
if( !first ) {
// interpolate initial flow from previous layer
interpolateFlowScale(layer1.width, layer1.height);
} else {
// for the very first layer there is no information on flow so set everything to 0
first = false;
flowU.reshape(layer1.width, layer1.height);
flowV.reshape(layer1.width, layer1.height);
ImageMiscOps.fill(flowU,0);
ImageMiscOps.fill(flowV,0);
}
// compute flow for this layer
processLayer(layer1,layer2,deriv1X,deriv1Y,deriv2X,deriv2Y,deriv2XX,deriv2YY,deriv2XY);
}
}
/**
* Resize images for the current layer being processed
*/
protected void resizeForLayer( int width , int height ) {
deriv1X.reshape(width,height);
deriv1Y.reshape(width,height);
deriv2X.reshape(width,height);
deriv2Y.reshape(width,height);
deriv2XX.reshape(width,height);
deriv2YY.reshape(width,height);
deriv2XY.reshape(width,height);
warpImage2.reshape(width,height);
warpDeriv2X.reshape(width,height);
warpDeriv2Y.reshape(width,height);
warpDeriv2XX.reshape(width,height);
warpDeriv2YY.reshape(width,height);
warpDeriv2XY.reshape(width,height);
derivFlowUX.reshape(width,height);
derivFlowUY.reshape(width,height);
derivFlowVX.reshape(width,height);
derivFlowVY.reshape(width,height);
psiData.reshape(width,height);
psiGradient.reshape(width,height);
psiSmooth.reshape(width,height);
divU.reshape(width,height);
divV.reshape(width,height);
divD.reshape(width,height);
du.reshape(width,height);
dv.reshape(width,height);
}
/**
* Provides an initial estimate for the flow by interpolating values from the previous layer.
*/
protected void interpolateFlowScale(int widthNew, int heightNew) {
// warping isn't done until later so use those images as temporary storage
GrayF32 enlargedU = warpDeriv2X;
GrayF32 enlargedV = warpDeriv2Y;
// use the previous low resolution flow estimate to initialize the new image
interpolateFlowScale(flowU, enlargedU);
interpolateFlowScale(flowV, enlargedV);
flowU.reshape(widthNew, heightNew);
flowV.reshape(widthNew, heightNew);
// save the initial flow values
flowU.setTo(enlargedU);
flowV.setTo(enlargedV);
}
/**
* Computes the flow for a layer using Taylor series expansion and Successive Over-Relaxation linear solver.
* Flow estimates from previous layers are feed into this by setting initFlow and flow to their values.
*/
protected void processLayer(GrayF32 image1 , GrayF32 image2 ,
GrayF32 deriv1X , GrayF32 deriv1Y,
GrayF32 deriv2X , GrayF32 deriv2Y,
GrayF32 deriv2XX , GrayF32 deriv2YY, GrayF32 deriv2XY) {
int N = image1.width*image1.height;
int stride = image1.stride;
// outer Taylor expansion iterations
for( int indexOuter = 0; indexOuter < numOuter; indexOuter++ ) {
// warp the image and the first + second derivatives
warpImageTaylor(image2, flowU, flowV, warpImage2);
warpImageTaylor(deriv2X, flowU, flowV, warpDeriv2X);
warpImageTaylor(deriv2Y, flowU, flowV, warpDeriv2Y);
warpImageTaylor(deriv2XX, flowU, flowV, warpDeriv2XX);
warpImageTaylor(deriv2YY, flowU, flowV, warpDeriv2YY);
warpImageTaylor(deriv2XY, flowU, flowV, warpDeriv2XY);
gradient.process(flowU,derivFlowUX,derivFlowUY);
gradient.process(flowV,derivFlowVX,derivFlowVY);
computePsiSmooth(derivFlowUX,derivFlowUY,derivFlowVX,derivFlowVY,psiSmooth);
computeDivUVD(flowU, flowV,psiSmooth,divU,divV,divD);
// initialize the motion increments to zero
Arrays.fill(du.data,0,N,0);
Arrays.fill(dv.data,0,N,0);
for( int indexInner = 0; indexInner < numInner; indexInner++ ) {
computePsiDataPsiGradient(image1, image2,
deriv1X, deriv1Y,
deriv2X, deriv2Y, deriv2XX, deriv2YY, deriv2XY,
du, dv, psiData, psiGradient);
float error;
int iter = 0;
do {
// inner SOR iteration.
error = 0;
// inner portion
for (int y = 1; y < image1.height - 1; y++) {
int i = y * image1.width + 1;
for (int x = 1; x < image1.width - 1; x++, i++) {
error += iterationSor(image1, deriv1X, deriv1Y, i, i + 1, i - 1, i + stride, i - stride);
}
}
// border regions require special treatment
int y0 = 0;
int y1 = image1.height-1;
for (int x = 0; x < image1.width; x++ ) {
error += iterationSor(image1, deriv1X, deriv1Y,
s(x, y0), s(x + 1, y0), s(x - 1, y0), s(x, y0 - 1), s(x, y0 + 1));
error += iterationSor(image1, deriv1X, deriv1Y,
s(x, y1), s(x + 1, y1), s(x - 1, y1), s(x, y1 - 1), s(x, y1 + 1));
}
int x0 = 0;
int x1 = image1.width-1;
for (int y = 1; y < image1.height - 1; y++) {
error += iterationSor(image1, deriv1X, deriv1Y,
s(x0, y), s(x0 - 1, y), s(x0 + 1, y), s(x0, y - 1), s(x0, y + 1));
error += iterationSor(image1, deriv1X, deriv1Y,
s(x1, y), s(x1 - 1, y), s(x1 + 1, y), s(x1, y - 1), s(x1, y + 1));
}
} while (error > convergeTolerance * image1.width * image1.height && ++iter < maxIterationsSor);
}
// update the flow with the motion increments
PixelMath.add(flowU,du, flowU);
PixelMath.add(flowV,dv, flowV);
}
}
/**
* Inner SOR iteration step
*
* @param i Index of target pixel at (x,y)
* @param ipx (x+1,y)
* @param imx (x-1,y)
* @param ipy (x,y+1)
* @param imy (x,y-1)
*/
private float iterationSor(GrayF32 image1, GrayF32 deriv1X, GrayF32 deriv1Y,
int i, int ipx, int imx, int ipy, int imy) {
float w = SOR_RELAXATION;
// these variables could be precomputed once. See equation 11
float psid = psiData.data[i];
float psig = gamma*psiGradient.data[i];
float di = warpImage2.data[i] - image1.data[i];
float dx2 = warpDeriv2X.data[i];
float dy2 = warpDeriv2Y.data[i];
float dxx2 = warpDeriv2XX.data[i];
float dyy2 = warpDeriv2YY.data[i];
float dxy2 = warpDeriv2XY.data[i];
float Au = -psid*di*dx2 + alpha*divU.data[i]
- psig*((dx2 - deriv1X.data[i])*warpDeriv2XX.data[i] + (dy2 - deriv1Y.data[i])*warpDeriv2XY.data[i]);
float Av = -psid*di*dy2 + alpha*divV.data[i]
- psig*((dx2 - deriv1X.data[i])*warpDeriv2XY.data[i] + (dy2 - deriv1Y.data[i])*warpDeriv2YY.data[i]);
float Du = psid*dx2*dx2 + psig*(dxx2*dxx2 + dxy2*dxy2) + alpha*divD.data[i];
float Dv = psid*dy2*dy2 + psig*(dyy2*dyy2 + dxy2*dxy2) + alpha*divD.data[i];
float D = psid*dx2*dy2 + psig*(dxx2 + dyy2)*dxy2;
// update the change in flow
float psi_index = psiSmooth.data[i];
float coef0 = 0.5f*(psiSmooth.data[ipx] + psi_index);
float coef1 = 0.5f*(psiSmooth.data[imx] + psi_index);
float coef2 = 0.5f*(psiSmooth.data[ipy] + psi_index);
float coef3 = 0.5f*(psiSmooth.data[imy] + psi_index);
float div_du = coef0 * du.data[ipx] + coef1 * du.data[imx] + coef2 * du.data[ipy] + coef3 * du.data[imy];
float div_dv = coef0 * dv.data[ipx] + coef1 * dv.data[imx] + coef2 * dv.data[ipy] + coef3 * dv.data[imy] ;
final float dui = du.data[i];
final float dvi = dv.data[i];
du.data[i] = (1f-w)*dui + w*(Au - D*dvi + alpha*div_du)/Du;
dv.data[i] = (1f-w)*dvi + w*(Av - D*du.data[i] + alpha*div_dv)/Dv;
return (du.data[i] - dui) * (du.data[i] - dui) + (dv.data[i] - dvi) * (dv.data[i] - dvi);
}
/**
* Equation 5. Psi_s
*/
private void computePsiSmooth(GrayF32 ux , GrayF32 uy , GrayF32 vx , GrayF32 vy ,
GrayF32 psiSmooth ) {
int N = derivFlowUX.width * derivFlowUX.height;
for( int i = 0; i < N; i++ ) {
float vux = ux.data[i];
float vuy = uy.data[i];
float vvx = vx.data[i];
float vvy = vy.data[i];
float mu = vux*vux + vuy*vuy;
float mv = vvx*vvx + vvy*vvy;
psiSmooth.data[i] = (float)(1.0/(2.0*Math.sqrt(mu + mv + EPSILON*EPSILON)));
}
}
/**
* Compute Psi-data using equation 6 and approximation in equation 5
*/
protected void computePsiDataPsiGradient(GrayF32 image1, GrayF32 image2,
GrayF32 deriv1x, GrayF32 deriv1y,
GrayF32 deriv2x, GrayF32 deriv2y,
GrayF32 deriv2xx, GrayF32 deriv2yy, GrayF32 deriv2xy,
GrayF32 du, GrayF32 dv,
GrayF32 psiData, GrayF32 psiGradient ) {
int N = image1.width * image1.height;
for( int i = 0; i < N; i++ ) {
float du_ = du.data[i];
float dv_ = dv.data[i];
// compute Psi-data
float taylor2 = image2.data[i] + deriv2x.data[i]*du_ + deriv2y.data[i]*dv_;
float v = taylor2 - image1.data[i];
psiData.data[i] = (float)(1.0/(2.0*Math.sqrt(v*v + EPSILON*EPSILON)));
// compute Psi-gradient
float dIx = deriv2x.data[i] + deriv2xx.data[i]*du_ + deriv2xy.data[i]*dv_ - deriv1x.data[i];
float dIy = deriv2y.data[i] + deriv2xy.data[i]*du_ + deriv2yy.data[i]*dv_ - deriv1y.data[i];
float dI2 = dIx*dIx + dIy*dIy;
psiGradient.data[i] = (float)(1.0/(2.0*Math.sqrt(dI2 + EPSILON*EPSILON)));
}
}
/**
* Computes the divergence for u,v, and d. Equation 8 and Equation 10.
*/
private void computeDivUVD(GrayF32 u , GrayF32 v , GrayF32 psi ,
GrayF32 divU , GrayF32 divV , GrayF32 divD ) {
final int stride = psi.stride;
// compute the inside pixel
for (int y = 1; y < psi.height-1; y++) {
// index of the current pixel
int index = y*stride + 1;
for (int x = 1; x < psi.width-1; x++ , index++) {
float psi_index = psi.data[index];
float coef0 = 0.5f*(psi.data[index+1] + psi_index);
float coef1 = 0.5f*(psi.data[index-1] + psi_index);
float coef2 = 0.5f*(psi.data[index+stride] + psi_index);
float coef3 = 0.5f*(psi.data[index-stride] + psi_index);
float u_index = u.data[index];
divU.data[index] = coef0*(u.data[index+1] - u_index) + coef1*(u.data[index-1] - u_index) +
coef2*(u.data[index+stride] - u_index) + coef3*(u.data[index-stride] - u_index);
float v_index = v.data[index];
divV.data[index] = coef0*(v.data[index+1] - v_index) + coef1*(v.data[index-1] - v_index) +
coef2*(v.data[index+stride] - v_index) + coef3*(v.data[index-stride] - v_index);
divD.data[index] = coef0 + coef1 + coef2 + coef3;
}
}
// handle the image borders
for( int x = 0; x < psi.width; x++ ) {
computeDivUVD_safe(x,0,u,v,psi,divU,divV,divD);
computeDivUVD_safe(x,psi.height-1,u,v,psi,divU,divV,divD);
}
for( int y = 1; y < psi.height-1; y++ ) {
computeDivUVD_safe(0,y,u,v,psi,divU,divV,divD);
computeDivUVD_safe(psi.width-1,y,u,v,psi,divU,divV,divD);
}
}
protected void computeDivUVD_safe(int x , int y ,
GrayF32 u , GrayF32 v , GrayF32 psi ,
GrayF32 divU , GrayF32 divV , GrayF32 divD ) {
int index = u.getIndex(x, y);
int index_px = s(x + 1, y);
int index_mx = s(x - 1, y);
int index_py = s(x, y + 1);
int index_my = s(x, y - 1);
float psi_index = psi.data[index];
float coef0 = 0.5f*(psi.data[index_px] + psi_index);
float coef1 = 0.5f*(psi.data[index_mx] + psi_index);
float coef2 = 0.5f*(psi.data[index_py] + psi_index);
float coef3 = 0.5f*(psi.data[index_my] + psi_index);
float u_index = u.data[index];
divU.data[index] = coef0*(u.data[index_px] - u_index) + coef1*(u.data[index_mx] - u_index) +
coef2*(u.data[index_py] - u_index) + coef3*(u.data[index_my] - u_index);
float v_index = v.data[index];
divV.data[index] = coef0*(v.data[index_px] - v_index) + coef1*(v.data[index_mx] - v_index) +
coef2*(v.data[index_py] - v_index) + coef3*(v.data[index_my] - v_index);
divD.data[index] = coef0 + coef1 + coef2 + coef3;
}
protected int s( int x , int y ) {
if( x < 0 ) x = 0;
else if( x >= warpImage2.width ) x = warpImage2.width-1;
if( y < 0 ) y = 0;
else if( y >= warpImage2.height ) y = warpImage2.height-1;
return warpImage2.getIndex(x,y);
}
public GrayF32 getFlowX() {
return flowU;
}
public GrayF32 getFlowY() {
return flowV;
}
}