edu.mines.jtk.dsp.LocalShiftFinder Maven / Gradle / Ivy
Show all versions of edu-mines-jtk Show documentation
/****************************************************************************
Copyright 2006, Colorado School of Mines and others.
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 edu.mines.jtk.dsp;
import static edu.mines.jtk.util.ArrayMath.*;
import edu.mines.jtk.util.Check;
/**
* Estimates displacement vector fields for two images. For example, given
* two 2-D images f(x1,x2) and g(x1,x2), a shift finder estimates two vector
* components of displacement u1(x1,x2) and u2(x1,x2) such that
* f(x1,x2) ~ g(x1+u1(x1,x2),x2+u2(x1,x2)).
*
* Like the images f and g, the components of displacement are sampled
* functions of coordinates x1 and x2. That is, displacements may vary
* from sample to sample. The components u1 and u2 represent displacements
* in the x1 and x2 coordinate directions, respectively.
*
* This shift finder estimates each component of displacement using local
* cross-correlations. For each image sample, the estimated shift is that
* which yields the maximum correlation coefficient. This coefficient is
* found by quadratic interpolation of correlation functions sampled at
* integer lags.
*
* The peak (maximum) correlation coefficient may be used to measure
* quality of an estimated shift. However, because a correlation function
* may have more than one peak (local maxima), a better measure of quality
* may be the difference between the coefficients for the correlation peak
* and next highest peak. Both the peak coefficient and this difference may
* be computed with the shifts.
*
* Methods are provided to find and compensate for each component of shift
* sequentially. As each component is found, that component can be removed
* from the image g before estimating another component. For example, again
* for 2-D images f(x1,x2) and g(x1,x2), we might first estimate u1(x1,x2).
* If we then compute an image h(x1,x2) = g(x1+u1(x1,x2),x2), we can use
* f(x1,x2) and h(x1,x2) to estimate u2(x1,x2). By repeating this process
* sequentially, we obtain estimates for both u1(x1,x2) and u2(x1,x2) such
* that f(x1,x2) ~ g(x1+u1(x1,x2),x2+u2(x1,x2)).
*
* Methods are also provided to whiten 2-D and 3-D images before estimating
* displacement vectors. This (spectral) whitening improves estimates of
* displacements parallel to image features that may be otherwise poorly
* resolved. Whitening is performed with local prediction error filters
* computed from local auto-correlations.
*
* @author Dave Hale, Colorado School of Mines
* @version 2006.11.18
*/
public class LocalShiftFinder {
/**
* Construct a shift estimator with specified parameters.
* When applied to multi-dimensional arrays, the estimator has the
* same correlation window half-width for all dimensions.
* @param sigma the correlation window half-width; must not be less than 1.
*/
public LocalShiftFinder(double sigma) {
this(sigma,sigma,sigma);
}
/**
* Construct a shift estimator with specified parameters.
* When applied to multi-dimensional arrays, the estimator has half-width
* sigma1 for the 1st dimension and half-width sigma2 for 2nd and higher
* dimensions.
* @param sigma1 correlaton window half-width for 0st dimension;
* must not be less than 1.
* @param sigma2 correlation window half-width for 2nd and higher
* dimensions; must not be less than 1.
*/
public LocalShiftFinder(double sigma1, double sigma2) {
this(sigma1,sigma2,sigma2);
}
/**
* Construct a shift estimator with specified parameters.
* When applied to multi-dimensional arrays, the estimator has half-width
* sigma1 for the 1st dimension, half-width sigma2 for the 2nd dimension,
* and half-width sigma3 for 3rd and higher dimensions.
* @param sigma1 correlation window half-width for 1st dimension;
* must not be less than 1.
* @param sigma2 correlation window half-width for 2nd dimension;
* must not be less than 1.
* @param sigma3 correlation window half-width for 3rd and higher
* dimensions; must not be less than 1.
*/
public LocalShiftFinder(double sigma1, double sigma2, double sigma3) {
Check.argument(sigma1>=1.0,"sigma1>=1.0");
Check.argument(sigma2>=1.0,"sigma2>=1.0");
Check.argument(sigma3>=1.0,"sigma3>=1.0");
_lcfSimple = new LocalCorrelationFilter(
LocalCorrelationFilter.Type.SIMPLE,
LocalCorrelationFilter.Window.GAUSSIAN,
sigma1,sigma2,sigma3);
_si = new SincInterpolator();
_si.setExtrapolation(SincInterpolator.Extrapolation.CONSTANT);
}
/**
* Enables or disables interpolation of displacements when shifting.
* The default is to interpolate displacements. This is the most
* accurate method when sequentially applying non-constant shifts.
* @param enable true, to enable interpolation; false, to disable.
*/
public void setInterpolateDisplacements(boolean enable) {
_interpolateDisplacements = enable;
}
/**
* Finds shifts in the 1st (and only) dimension.
* @param min1 the minimum shift.
* @param max1 the maximum shift.
* @param f the input array f.
* @param g the input array g.
* @param u output array of shifts.
*/
public void find1(
int min1, int max1, float[] f, float[] g, float[] u)
{
findShifts(min1,max1,f,g,u,null,null);
}
/**
* Finds shifts in the 1st (and only) dimension.
* Also computes peak correlation coefficients and differences between
* the peak and next-highest-peak coeffcients.
* @param min1 the minimum shift.
* @param max1 the maximum shift.
* @param f the input array f.
* @param g the input array g.
* @param u output array of shifts.
* @param c output array of peak correlation coefficients.
* @param d output array of differences, peak minus next-highest-peak.
*/
public void find1(
int min1, int max1, float[] f, float[] g,
float[] u, float[] c, float[] d)
{
findShifts(min1,max1,f,g,u,c,d);
}
/**
* Finds shifts in the 1st dimension.
* @param min1 the minimum shift.
* @param max1 the maximum shift.
* @param f the input array f.
* @param g the input array g.
* @param u output array of shifts.
*/
public void find1(
int min1, int max1, float[][] f, float[][] g, float[][] u)
{
findShifts(1,min1,max1,f,g,u);
}
/**
* Finds shifts in the 2nd dimension.
* @param min2 the minimum shift.
* @param max2 the maximum shift.
* @param f the input array f.
* @param g the input array g.
* @param u output array of shifts.
*/
public void find2(
int min2, int max2, float[][] f, float[][] g, float[][] u)
{
findShifts(2,min2,max2,f,g,u);
}
/**
* Finds shifts in the 1st dimension.
* @param min1 the minimum shift.
* @param max1 the maximum shift.
* @param f the input array f.
* @param g the input array g.
* @param u output array of shifts.
*/
public void find1(
int min1, int max1, float[][][] f, float[][][] g, float[][][] u)
{
findShifts(1,min1,max1,f,g,u);
}
/**
* Finds shifts in the 2nd dimension.
* @param min2 the minimum shift.
* @param max2 the maximum shift.
* @param f the input array f.
* @param g the input array g.
* @param u output array of shifts.
*/
public void find2(
int min2, int max2, float[][][] f, float[][][] g, float[][][] u)
{
findShifts(2,min2,max2,f,g,u);
}
/**
* Finds shifts in the 3rd dimension.
* @param min3 the minimum shift.
* @param max3 the maximum shift.
* @param f the input array f.
* @param g the input array g.
* @param u output array of shifts.
*/
public void find3(
int min3, int max3, float[][][] f, float[][][] g, float[][][] u)
{
findShifts(3,min3,max3,f,g,u);
}
/**
* Applies specified shift in the 1st (and only) dimension.
* @param du input array of changes to displacements in 1st dimension.
* @param u1 input/output array of displacements in 1st dimension.
* @param h input/output array of image samples.
*/
public void shift1(float[] du, float[] u1, float[] h) {
int n1 = h.length;
float[] xu1 = new float[n1];
float[] u1a = u1;
float[] u1b = new float[n1];
float[] ha = h;
float[] hb = new float[n1];
for (int i1=0; i10.0) {
double l22 = sqrt(d22);
double v1 = b1/l11;
double v2 = (b2-l21*v1)/l22;
x2 = v2/l22;
x1 = (v1-l21*x2)/l11;
}
float a1 = (float)x1;
float a2 = (float)x2;
s[i2][i1] = f[i2][i1]
- a1*f[i2][i1-1]
- a2*f[i2-1][i1];
}
}
if (sigma>=1.0) {
RecursiveGaussianFilter rgf = new RecursiveGaussianFilter(sigma);
rgf.apply0X(s,t);
rgf.applyX0(t,g);
} else {
copy(s,g);
}
}
/**
* Applies local prediction-error (spectal whitening) filters.
* The input and output arrays f and g can be the same array.
* Smooths the output with a Gaussian filter with half-width sigma = 1.0.
* @param f the input array.
* @param g the output array.
*/
public void whiten(float[][][] f, float[][][] g) {
whiten(1.0,f,g);
}
/**
* Applies local prediction-error (spectal whitening) filters.
* The input and output arrays f and g can be the same array.
* @param sigma half-width of Gaussian smoothing applied after whitening;
* less than one for no smoothing.
* @param f the input array.
* @param g the output array.
*/
public void whiten(double sigma, float[][][] f, float[][][] g) {
int n1 = f[0][0].length;
int n2 = f[0].length;
int n3 = f.length;
float[][][] r000 = new float[n3][n2][n1];
float[][][] rpm0 = new float[n3][n2][n1];
float[][][] rp0m = new float[n3][n2][n1];
float[][][] r0pm = new float[n3][n2][n1];
float[][][] rm00 = new float[n3][n2][n1];
float[][][] r0m0 = new float[n3][n2][n1];
float[][][] r00m = new float[n3][n2][n1];
float[][][] s = rm00;
float[][][] t = r0m0;
_lcfSimple.setInputs(f,f);
_lcfSimple.correlate( 0, 0, 0,r000);
_lcfSimple.correlate( 1,-1, 0,rpm0);
_lcfSimple.correlate( 1, 0,-1,rp0m);
_lcfSimple.correlate( 0, 1,-1,r0pm);
_lcfSimple.correlate(-1, 0, 0,rm00);
_lcfSimple.correlate( 0,-1, 0,r0m0);
_lcfSimple.correlate( 0, 0,-1,r00m);
for (int i3=0; i30.0) {
double l22 = sqrt(d22);
double l32 = (a32-l31*l21)/l22;
double d33 = a33-l31*l31-l32*l32;
if (d33>0.0) {
double l33 = sqrt(d33);
double v1 = b1/l11;
double v2 = (b2-l21*v1)/l22;
double v3 = (b3-l31*v1-l32*v2)/l33;
x3 = v3/l33;
x2 = (v2-l32*x3)/l22;
x1 = (v1-l21*x2-l31*x3)/l11;
}
}
float a1 = (float)x1;
float a2 = (float)x2;
float a3 = (float)x3;
s[i3][i2][i1] = f[i3][i2][i1]
- a1*f[i3][i2][i1-1]
- a2*f[i3][i2-1][i1]
- a3*f[i3-1][i2][i1];
}
}
}
if (sigma>=1.0) {
RecursiveGaussianFilter rgf = new RecursiveGaussianFilter(sigma);
rgf.apply0XX(s,t);
rgf.applyX0X(t,s);
rgf.applyXX0(s,g);
} else {
copy(s,g);
}
}
///////////////////////////////////////////////////////////////////////////
// private
private LocalCorrelationFilter _lcfSimple;
private SincInterpolator _si;
private boolean _interpolateDisplacements = true;
private void findShifts(
int min, int max, float[] f, float[] g, float[] u, float[] c, float[] d)
{
int n1 = f.length;
// Initially shifts, correlations, and differences are zero.
zero(u);
if (c!=null)
zero(c);
else
c = zerofloat(n1);
if (d!=null)
zero(d);
// Arrays to contain cross-correlations for three consecutive lags.
float[][] c3 = new float[3][n1];
// Correlate for min lag.
LocalCorrelationFilter lcf = _lcfSimple;
lcf.setInputs(f,g);
int lag1 = min;
lcf.correlate(lag1,c3[1]);
lcf.normalize(lag1,c3[1]);
// For all lags in range [min,max], ...
for (int lag=min; lag<=max; ++lag) {
// Arrays ca, cb, and cc will contain three cross-correlations. For
// first and last lags, buffers a and c are the same. In other words,
// assume that correlation values are symmetric about the min and max
// lags scanned. This assumption enables local maxima to occur at the
// specified min and max lags, but forces displacements to lie within
// the range [min,max].
int i = lag-min;
float[] ca = (lag>min)?c3[(i )%3]:c3[(i+2)%3];
float[] cb = c3[(i+1)%3];
float[] cc = (lag=ai && bi>=ci) {
double c0 = bi;
double c1 = 0.5*(ci-ai);
double c2 = 0.5*(ci+ai)-bi;
double up = (c2<0.0)?-0.5*c1/c2:0.0;
double cp = c0+up*(c1+up*c2);
if (cp>c[i1]) {
if (d!=null) d[i1] = (float)cp-c[i1];
c[i1] = (float)cp;
u[i1] = (float)(lag+up);
}
}
}
}
}
private void findShifts(
int min, int max, float[] f, float[] g, float[] u)
{
int n1 = f.length;
// Default shifts are zero.
zero(u);
// Arrays to contain cross-correlations for three consecutive lags.
float[][] c = new float[3][n1];
// Array for current correlation maximum values.
float[] cmax = new float[n1];
// Correlate for min lag.
LocalCorrelationFilter lcf = _lcfSimple;
lcf.setInputs(f,g);
int lag1 = min;
lcf.correlate(lag1,c[1]);
lcf.normalize(lag1,c[1]);
// For all lags in range [min,max], ...
for (int lag=min; lag<=max; ++lag) {
// Arrays ca, cb, and cc will contain three cross-correlations. For
// first and last lags, buffers a and c are the same. In other words,
// assume that correlation values are symmetric about the min and max
// lags scanned. This assumption enables local maxima to occur at the
// specified min and max lags, but forces displacements to lie within
// the range [min,max].
int i = lag-min;
float[] ca = (lag>min)?c[(i )%3]:c[(i+2)%3];
float[] cb = c[(i+1)%3];
float[] cc = (lag=ai && bi>=ci) {
double c0 = bi;
double c1 = 0.5*(ci-ai);
double c2 = 0.5*(ci+ai)-bi;
double up = (c2<0.0)?-0.5*c1/c2:0.0;
double cp = c0+up*(c1+up*c2);
if (cp>cmax[i1]) {
cmax[i1] = (float)cp;
u[i1] = (float)(lag+up);
}
}
}
}
}
private void findShifts(
int dim, int min, int max, float[][] f, float[][] g, float[][] u)
{
int n1 = f[0].length;
int n2 = f.length;
// Default shifts are zero.
zero(u);
// Arrays to contain cross-correlations for three consecutive lags.
float[][][] c = new float[3][n2][n1];
// Array for current correlation maximum values.
float[][] cmax = new float[n2][n1];
// Correlate for min lag.
LocalCorrelationFilter lcf = _lcfSimple;
lcf.setInputs(f,g);
int lag1 = (dim==1)?min:0;
int lag2 = (dim==2)?min:0;
lcf.correlate(lag1,lag2,c[1]);
lcf.normalize(lag1,lag2,c[1]);
// For all lags in range [min,max], ...
for (int lag=min; lag<=max; ++lag) {
// Arrays ca, cb, and cc will contain three cross-correlations. For
// first and last lags, buffers a and c are the same. In other words,
// assume that correlation values are symmetric about the min and max
// lags scanned. This assumption enables local maxima to occur at the
// specified min and max lags, but forces displacements to lie within
// the range [min,max].
int i = lag-min;
float[][] ca = (lag>min)?c[(i )%3]:c[(i+2)%3];
float[][] cb = c[(i+1)%3];
float[][] cc = (lag=ai && bi>=ci) {
double c0 = bi;
double c1 = 0.5*(ci-ai);
double c2 = 0.5*(ci+ai)-bi;
double up = (c2<0.0)?-0.5*c1/c2:0.0;
double cp = c0+up*(c1+up*c2);
if (cp>cmax[i2][i1]) {
cmax[i2][i1] = (float)cp;
u[i2][i1] = (float)(lag+up);
}
}
}
}
}
}
private void findShifts(
int dim, int min, int max, float[][][] f, float[][][] g, float[][][] u)
{
int n1 = f[0][0].length;
int n2 = f[0].length;
int n3 = f.length;
// Default shifts are zero.
zero(u);
// Arrays to contain cross-correlations for three consecutive lags.
float[][][][] c = new float[3][n3][n2][n1];
// Array for current correlation maximum values.
float[][][] cmax = new float[n3][n2][n1];
// Correlate for min lag.
LocalCorrelationFilter lcf = _lcfSimple;
lcf.setInputs(f,g);
int lag1 = (dim==1)?min:0;
int lag2 = (dim==2)?min:0;
int lag3 = (dim==3)?min:0;
lcf.correlate(lag1,lag2,lag3,c[1]);
lcf.normalize(lag1,lag2,lag3,c[1]);
// For all lags in range [min,max], ...
for (int lag=min; lag<=max; ++lag) {
// Arrays ca, cb, and cc will contain three cross-correlations. For
// first and last lags, buffers a and c are the same. In other words,
// assume that correlation values are symmetric about the min and max
// lags scanned. This assumption enables local maxima to occur at the
// specified min and max lags, but forces displacements to lie within
// the range [min,max].
int i = lag-min;
float[][][] ca = (lag>min)?c[(i )%3]:c[(i+2)%3];
float[][][] cb = c[(i+1)%3];
float[][][] cc = (lag=ai && bi>=ci) {
double c0 = bi;
double c1 = 0.5*(ci-ai);
double c2 = 0.5*(ci+ai)-bi;
double up = (c2<0.0)?-0.5*c1/c2:0.0;
double cp = c0+up*(c1+up*c2);
if (cp>cmax[i3][i2][i1]) {
cmax[i3][i2][i1] = (float)cp;
u[i3][i2][i1] = (float)(lag+up);
}
}
}
}
}
}
}
}