https://doi.org/10.5201/ipol.2016.172
Revision edce3fbb7ce8b72fd316f77f13ff3f1d75091a52 authored by Software Heritage on 19 February 2016, 00:00:00 UTC, committed by Software Heritage on 07 November 2016, 00:00:00 UTC
0 parent
Tip revision: edce3fbb7ce8b72fd316f77f13ff3f1d75091a52 authored by Software Heritage on 19 February 2016, 00:00:00 UTC
ipol: Deposit 1282 in collection ipol
ipol: Deposit 1282 in collection ipol
Tip revision: edce3fb
smoothness.h
// This program is free software: you can use, modify and/or redistribute it
// under the terms of the simplified BSD License. You should have received a
// copy of this license along this program. If not, see
// <http://www.opensource.org/licenses/bsd-license.html>.
//
// Copyright (C) 2012, Javier Sánchez Pérez <jsanchez@dis.ulpgc.es>
// Copyright (C) 2014, Agustín Salgado de la Nuez <asalgado@dis.ulpgc.es>
// Copyright (C) 2014, 2015 Nelson Monzón López <nmonzon@ctim.es>
// All rights reserved.
#ifndef SMOOTHNESS_H
#define SMOOTHNESS_H
#include <algorithm>
#include <vector>
#define EPSILON 0.001
#define XI 0.05
#define TAU 0.94
#define BETA 0.001
/**
*
* Compute the coefficients of the robust functional (smoothness term)
*
**/
void psi_smooth(
const float *ux, //gradient of x component of the optical flow
const float *uy, //gradient of x component of the optical flow
const float *vx, //gradient of y component of the optical flow
const float *vy, //gradient of y component of the optical flow
const float *expo, //exponential smoothing factor
const int size_flow,
float *psi //output coefficients
)
{
//compute 1/(sqrt(ux²+uy²+vx²+vy²+e²) in each pixel
#pragma omp parallel for
for(int i = 0; i < size_flow; i++)
{
const float du = expo[i] * ux[i] * ux[i] + expo[i] * uy[i] * uy[i];
const float dv = expo[i] * vx[i] * vx[i] + expo[i] * vy[i] * vy[i];
const float normFlow = du + dv;
psi[i] = expo[i] / sqrt(normFlow + EPSILON * EPSILON);
}
}
void max_gradients(
const float *Ix, // Computed Image 1 derivative in x
const float *Iy, // Computed Image 1 derivative in y
const int size, // Total image size (height * weight * nchannels)
const int nz, // nº channels
float *maximum_gradients_per_pixel // vector with the maximum gradients per pixel
){
int index_flow = 0;
for(int index_image = 0; index_image < size; index_image+=nz){
maximum_gradients_per_pixel[index_flow] = sqrt(Ix[index_image]*Ix[index_image] + Iy[index_image]*Iy[index_image]);
for(int index_multichannel = 1; index_multichannel < nz; index_multichannel++){
const int real_index = index_image + index_multichannel;
const float gradient_in_this_pixel = sqrt(Ix[real_index]*Ix[real_index] + Iy[real_index]*Iy[real_index]);
if (maximum_gradients_per_pixel[index_flow] < gradient_in_this_pixel)
maximum_gradients_per_pixel[index_flow] = gradient_in_this_pixel;
}
index_flow++;
}
return;
}
/**
* Calculate the lambda optimum using the maximum gradient from all the multi-channel image
* It also return the maximum gradient per pixel
**/
float lambda_optimum_using_maximum_gradient_per_pixel(
const float *Ix, // Computed Image 1 derivative in x
const float *Iy, // Computed Image 1 derivative in y
const int size, // Total image size (height * weight * nchannels)
const int size_flow, // Total flow size (height * weight)
const int nz, // nº channels
const float alpha, // smoothness weight
float *lambda_per_pixel, // local lambda per pixel
float *maximum_gradients_per_pixel // vector with the maximum gradients per pixel
)
{
int index_flow = 0;
for(int index_image = 0; index_image < size; index_image+=nz){
maximum_gradients_per_pixel[index_flow] = sqrt(Ix[index_image]*Ix[index_image] + Iy[index_image]*Iy[index_image]);
for(int index_multichannel = 1; index_multichannel < nz; index_multichannel++){
const int real_index = index_image + index_multichannel;
const float gradient_in_this_pixel = sqrt(Ix[real_index]*Ix[real_index] + Iy[real_index]*Iy[real_index]);
if (maximum_gradients_per_pixel[index_flow] < gradient_in_this_pixel)
maximum_gradients_per_pixel[index_flow] = gradient_in_this_pixel;
}
lambda_per_pixel[index_flow] = (-log(XI) + log(alpha)) / maximum_gradients_per_pixel[index_flow];
index_flow++;
}
std::vector<float> gradients_ordered(size_flow);
#pragma omp parallel for
for (int index_flow = 0; index_flow < size_flow; index_flow++)
gradients_ordered[index_flow] = maximum_gradients_per_pixel[index_flow];
std::sort(gradients_ordered.begin(), gradients_ordered.end());
const float c = -log(XI) + log(alpha);
int pos_ref = TAU * size_flow;
float lambda_pi = c / gradients_ordered[pos_ref-1];
while ((pos_ref < size_flow) && (c/2 > gradients_ordered[pos_ref-1]))
pos_ref++;
if (pos_ref == size_flow)
lambda_pi = 0;
else
lambda_pi = (c / gradients_ordered[pos_ref-1]);
return lambda_pi;
}
/**
** Calculate the exponential values.
**
**/
void exponential_calculation(
const float *Ix, // Computed Image 1 derivative in x
const float *Iy, // Computed Image 1 derivative in y
const int size_flow, // size of the flow field
const int size, // size of the multi-channel image
const int nz, // nº of image channels
const float alpha, // smoothness weight
const float lambda, // Coeffient for decreasing function
const int method_type, // (1 = DF, 2 = DF_BETA, 3 = DF_AUTO)
float *expo // e^(lambda * DI)
)
{
float *maximum_gradients_per_pixel = new float[size_flow];
switch(method_type){
case 1: case 2:{
float beta = 0;
if (method_type == 2) // For Exponential Beta Approximation
beta = BETA;
max_gradients(Ix, Iy, size, nz, maximum_gradients_per_pixel);
#pragma omp parallel for
for(int index_flow = 0; index_flow < size_flow; index_flow++){
expo[index_flow] = exp(- lambda * maximum_gradients_per_pixel[index_flow]) + beta;
}
break;
}
case 3:{
float *lambda_per_pixel = new float[size_flow];
float lambda_omega = lambda_optimum_using_maximum_gradient_per_pixel(Ix, Iy, size, size_flow, nz, alpha, lambda_per_pixel, maximum_gradients_per_pixel);
#pragma omp parallel for
for(int index_flow = 0; index_flow < size_flow; index_flow++){
float lambda_pi = lambda_omega;
if(lambda_omega > lambda_per_pixel[index_flow])
lambda_pi = lambda_per_pixel[index_flow];
expo[index_flow] = exp(- lambda_pi * maximum_gradients_per_pixel[index_flow]);
}
delete [] lambda_per_pixel;
break;
}
} // End switch
delete [] maximum_gradients_per_pixel;
}
#endif
Computing file changes ...