edu.mines.jtk.dsp.SteerablePyramid Maven / Gradle / Ivy
Show all versions of edu-mines-jtk Show documentation
/****************************************************************************
Copyright 2008, 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.*;
/**
* Steerable pyramid for 2D and 3D images. Includes creation of the steerable
* pyramid, estimation of local orientation and dimensionality attributes,
* application of steering weights and scaling to enhance locally linear
* features, and construction of the filtered image.
*
* This implementation of the steerable pyramid transform performs a
* multi-scale, multi-orientation decomposition of an input image through
* application of radial and directional filters in wavenumber domain.
* The basis steerable filter amplitudes are proportional to cos^2(theta).
* Three basis orientations are used for 2D, and six orientations are used
* for 3D images. Radial filters are used to partition the data into
* 1-octave bands, with a cosine taper. Images are subsampled for each
* pyramid level which greatly reduces processing effort for lower
* wavenumbers.
*
* Directionally-filtered basis images are used to estimate local orientation
* and dimensionality. Preprocessing, which includes averaging in space
* and scale domains, is applied for these estimates. Steering weights can be
* calculated and applied. Scaling and thresholding can also be applied,
* based on a local dimensionality attribute. For 3D images, processing can
* be applied to enhance either locally-linear or locally-planar features.
*
* The number of pyramid levels to use is calculated from the size of the
* input image, assuming a minimum basis image dimension of 9 samples in x,y
* or z. The input image is padded before it is transformed to wavenumber
* domain, where the filters are applied. The main reason for this padding is
* to avoid losing first or last samples when subsampling. We like the number
* of samples to be such that, for number of samples n in x, y, or z,
* (n-1)/2+1 will always yield an integer.
*
* The format of the steerable pyramid is a 4-dimensional array for 2D, and a
* 5-dimensional array for 3D. The first dimension is level number and second
* dimension is basis filter orientation. Below these are either 2D or 3D
* arrays.
* To illustrate for 2D:
*
* [0][0][0][0] to [0][0][n2][n1] = level 0, theta 0
*
* [0][1][0][0] to [0][1][n2][n1] = level 0, theta PI/3
*
* [0][2][0][0] to [0][2][n2][n1] = level 0, theta 2*PI/3
*
* [1][0][0][0] to [1][0][(n2-1)/2+1][(n1-1)/2+1] = level 1, theta 0
*
* [1][1][0][0] to [1][1][(n2-1)/2+1][(n1-1)/2+1] = level 1, theta PI/3
*
* [1][2][0][0] to [1][2][(n2-1)/2+1][(n1-1)/2+1] = level 1, theta 2*PI/3
*
* ...
*
* [NLEVEL-1][2][0][0] to
* [NLEVEL-1][2][(n2-1)/(2^(NLEVEL-1))+1][(n1-1)/(2^(NLEVEL-1))+1]
* = level N-1, theta 2*PI/3
*
* [NLEVEL][0][0][0] to [NLEVEL][0][(n2-1)/(2^NLEVEL)+1][(n1-1)/(2^NLEVEL)+1]
* = residual low-wavenumber image
*
* The 3D steerable pyramid array is the same except that it is arrays of
* arrays of 3D, rather than 2D arrays.
*
* @author John Mathewson, Colorado School of Mines
* @version 2008.12.01
*/
public class SteerablePyramid {
/**
* Construct a steerable pyramid with default cutoff wavenumbers
* used in the radial low-pass filters. Default values are:
* ka=0.60 and kb=1.00, where ka and kb are wavenumbers at start and
* end of taper (Amp(ka)=1.0, Amp(kb)=0.0).
*/
public SteerablePyramid() {
ka = 0.60;
kb = 1.00;
}
/**
* Construct a steerable pyramid with specified cutoff wavenumbers
* used in the radial low-pass filters.
* @param ka wavenumber at start of taper. Amp(ka)=1.
* @param kb wavenumber at end of taper. Amp(ka)=0.
*/
public SteerablePyramid(double ka,double kb) {
this.ka = ka;
this.kb = kb;
}
/**
* Creates a steerable pyramid representation of an input 2D image.
* @param x input 2D image.
* @return array containing steerable pyramid representation of the input
* 2D image.
*/
public float[][][][] makePyramid(float[][] x) {
nx2 = x.length;
nx1 = x[0].length;
// Compute number of levels in pyramid from size of input image. Also
// determine dimensions n1 and n2 for the finest-sampled pyramid level
// that will allow us to subsample each pyramid level without losing the
// last sample. In our pyramid images we will carry this number of
// samples and copy the original number of samples only for final output.
nlev = 1;
int nlev2 = 1;
n1 = 9;
n2 = 9;
while (nx1>n1) {
n1 = (n1-1)*2+1;
nlev += 1;
}
while (nx2>n2) {
n2 = (n2-1)*2+1;
nlev2 += 1;
}
if (nlev>nlev2) {
nlev = nlev2;
}
/**
* Create output 4-dimensional array, consisting of:
* Basis images: nlev*ndir 2D sub-arrays; for each pyramid level images
* are subsampled, so they are 25% size of images in preceding level.
* Low-wavenumber image: Single image, residual low-wavenumber energy.
*/
float[][][][] spyr = new float[nlev+1][1][1][1];
for (int lev=0; lev0) {
cf = ftForward(lev,spyr[lev][0]);
}
makePyramidLevel(lev,cf,spyr);
}
return spyr;
}
/**
* Creates a steerable pyramid representation of an input 3D image.
* @param x input 3D image.
* @return array containing steerable pyramid representation of the input
* 3D image.
*/
public float[][][][][] makePyramid(float[][][] x) {
nx3 = x.length;
nx2 = x[0].length;
nx1 = x[0][0].length;
// Compute number of levels in pyramid from size of input image. Also
// determine dimensions n1,n2,n3 for the finest-sampled pyramid level
// that will allow us to subsample each pyramid level without losing the
// last sample. In our pyramid images we will carry this number of
// samples and copy the original number of samples only for final output.
nlev = 1;
int nlev2 = 1;
int nlev3 = 1;
n1 = 9;
n2 = 9;
n3 = 9;
while (nx1>n1) {
n1 = (n1-1)*2+1;
nlev += 1;
}
while (nx2>n2) {
n2 = (n2-1)*2+1;
nlev2 += 1;
}
while (nx3>n3) {
n3 = (n3-1)*2+1;
nlev3 += 1;
}
if (nlev>nlev2) {
nlev = nlev2;
}
if (nlev>nlev3) {
nlev = nlev3;
}
/**
* Create output 5-dimensional array, consisting of:
* Basis images: nlev*ndir 3D sub-arrays; for each pyramid level images
* are subsampled, so they are 12.5% size of images in preceding level.
* Low-wavenumber image: Single image, residual low-wavenumber energy.
*/
float[][][][][] spyr = new float[nlev+1][1][1][1][1];
for (int lev=0; lev0) {
cf = ftForward(lev,spyr[lev][0]);
}
makePyramidLevel(lev,cf,spyr);
}
return spyr;
}
/**
* Sums all basis images from an input 2D steerable pyramid to create a
* filtered output image. Optionally, residual low-wavenumber image can
* be zeroed.
* @param keeplow if true:keep low-wavenumber energy, if false: zero it.
* @param spyr input 2D steerable pyramid.
* @return array containing output filtered 2D image.
*/
public float[][] sumPyramid(boolean keeplow,float[][][][] spyr) {
int lev;
// Optionally zero the low-wavenumber image.
if (!keeplow) {
zero(spyr[nlev][0]);
}
// Sum basis images to create filtered image. Sinc interpolation is
// performed on subsampled images prior to summing adjacent levels.
int lfactor = (int)pow(2.0,(double)(nlev-1));
int nl2 = (n2-1)/lfactor+1;
int nl1 = (n1-1)/lfactor+1;
for (int i=0; i
* The format of the output attributes is a 4-dimensional array.
* The first dimension is level number and second dimension is type
* of attribute. Below these are 2D arrays:
*
* [0][0][0][0] to [0][0][n2][n1] = level 0, theta attribute (radians)
*
* [0][1][0][0] to [0][1][n2][n1] = level 0, linearity attribute
*
* [1][0][0][0] to [1][0][(n2-1)/2+1][(n1-1)/2+1] = level 1,
* theta attribute (radians)
*
* [1][1][0][0] to [1][1][(n2-1)/2+1][(n1-1)/2+1] = level 1,
* linearity attribute
*
* ...
*
* [NLEVEL-1][1][0][0] to
* [NLEVEL-1][1][(n2-1)/(2^(NLEVEL-1))+1][(n1-1)/(2^(NLEVEL-1))+1]
* = level N-1, linearity attribute
* @param sigma half-width of 2D Gaussian smoothing filter.
* @param spyr input 2D steerable pyramid.
* @return array containing local orientation estimate and linearity
* attribute for all sample locations in every pyramid level. Array
* consists of numlevels*2 2D sub-arrays. For each level the first
* sub-array contains orientation theta in radians, and the second
* contains linearity attribute ranging from 0 to 1.
*/
public float[][][][] estimateAttributes(double sigma,float[][][][] spyr) {
double sigmaa = 2.0*sigma;
double sigmac = 0.5*sigma;
double a0,a1,a2;
double e0,e1;
double[][] feout;
int i1a,i2a,i1c,i2c;
// Allocate output 4-dimensional array.
float[][][][] attr = new float[nlev][2][1][1];
for (int lev=0; lev=0) {
float[][][] pqja = pqjShiftSmooth(sigmaa,leva,spyr);
for (int i2b=0; i2b
* In 3D we have a choice of filtering to enhance locally-planar or
* locally-linear image features. There is a parameter in this
* method to select one of these choices. If enhancement of planar
* features is selected the output orientation attributes define the
* normal to locally-planar features, and the dimensionality attribute
* is a measure of planarity. If enhancement of locally-linear
* features is selected the output orientation attributes define the
* orientation of a locally-linear feature, and the dimensionality
* attribute is a local measure of linearity.
*
* The format of the output attributes is a 5-dimensional array.
* The first dimension is level number and second dimension is type
* of attribute. Below these are 3D arrays:
*
* [0][0][0][0][0] to [0][0][n3][n2][n1] = level 0, direction cosine a
*
* [0][0][0][0][0] to [0][0][n3][n2][n1] = level 0, direction cosine b
*
* [0][0][0][0][0] to [0][0][n3][n2][n1] = level 0, direction cosine c
*
* [0][1][0][0][0] to [0][1][n3][n2][n1] = level 0, dimensionality attribute
*
* These are repeated for all levels, subsampled for every successive level.
* @param forlinear true: apply to enhance locally linear, false: apply for
* planar.
* @param sigma half-width of 3D Gaussian smoothing filter.
* @param spyr input 3D steerable pyramid.
* @return array containing local orientation estimate and dimensionality
* attribute for all sample locations in every pyramid level. Array
* consists of numlevels*4 3D sub-arrays. For each level the first three
* sub-arrays contain direction cosines, and the fourth contains
* dimensionality attribute ranging from 0 to 1.
*/
public float[][][][][] estimateAttributes(boolean forlinear,double sigma,
float[][][][][] spyr) {
double sigmaa = 2.0*sigma;
double sigmac = 0.5*sigma;
double[] f = zerodouble(NDIR3);
double[][] abcf = zerodouble(4,3);
int i1a,i2a,i3a,i1c,i2c,i3c;
// Parameters to select estimation for locally planar or linear features
int abcindx = 2;
int e0indx = 2;
int e1indx = 1;
statelinear = forlinear;
if (forlinear) {
abcindx = 0;
e0indx = 1;
e1indx = 0;
}
// Allocate output 5-dimensional array.
float[][][][][] attr = new float[nlev][4][1][1][1];
for (int lev=0; lev=0) {
float[][][][] pqja = pqjShiftSmooth(sigmaa,leva,spyr);
for (int i3b=0; i3b1&&linpowr<99) {
scal = pow(attr[lev][3][i3][i2][i1],linpowr);
}
else if(linpowr==99) {
scal = 1.0f/(1.0f+exp(k*(thresh-attr[lev][3][i3][i2][i1])));
}
spyr[lev][dir][i3][i2][i1] *= scal*wi;
}
}
}
}
}
}
/**
* Applies steering weights and scaling or thresholding based on linearity
* attribute to the basis images in the input 2D steerable pyramid.
* Scaling options include:
* No scaling (linpowr=0).
* Linearity raised to a power (1 ≤ linpowr ≤ 98).
* Sigmoidal thresholding (linpowr=99) Typical values for thresholding
* parameters are k=50.0, thresh=0.5.
* @param linpowr linearity power and scaling type switch.
* @param k sigmoidal thresholding steepness.
* @param thresh threshold.
* @param attr input array containing local orientation and linearity.
* @param spyr input/output 2D steerable pyramid.
*/
public void steerScale(int linpowr,float k, float thresh,
float[][][][] attr, float[][][][] spyr) {
float w0,w1,w2;
float scal = 0.0f;
double theta;
int nl2,nl1;
for (int lev=0; lev1&&linpowr<99) {
scal = pow(attr[lev][1][i2][i1],linpowr);
}
else if(linpowr==99) {
scal = 1.0f/(1.0f+exp(k*(thresh-attr[lev][1][i2][i1])));
}
spyr[lev][0][i2][i1] *= scal*w0;
spyr[lev][1][i2][i1] *= scal*w1;
spyr[lev][2][i2][i1] *= scal*w2;
}
}
}
}
///////////////////////////////////////////////////////////////////////////
// private
private static final double THETA0 = 0.0*PI/3.0;
private static final double THETA1 = 1.0*PI/3.0;
private static final double THETA2 = 2.0*PI/3.0;
private static final double ONE_THIRD = 1.0/3.0;
private static final float TWO_THIRDS = 2.0f/3.0f;
private static final double SQRT3 = sqrt(3.0);
private static final int NDIR2=3;
private static final int NDIR3=6;
// Smallest significant change to a, b, or c (for convergence test).
private static final double ABC_SMALL = 1.0e-2;
// Smallest significant determinant (denominator in Newton's method).
private static final double DET_SMALL = 1.0e-40;
// Other constants.
private static final double COS_PIO3 = cos(PI/3.0);
private static final double SIN_PIO3 = sin(PI/3.0);
private int nlev,nx1,nx2,nx3,n1,n2,n3;
private boolean statelinear;
double ka,kb;
/**
* Make a single 2D pyramid level consisting of three directionally-filtered
* basis images and a subsampled lower-wavenumber image. (The low-wavenumber
* image is the input for creation of the next level). The output is
* written into the appropriate part of the 2D steerable pyramid.
* @param lev level number.
* @param cf input image in wavenumber domain (complex array).
* @param spyr input/output 2D steerable pyramid.
*/
private void makePyramidLevel(int lev,float[][] cf,float[][][][] spyr) {
int lfactor = (int)pow(2.0,(double)lev);
int nl2 = (n2-1)/lfactor+1;
int nl1 = (n1-1)/lfactor+1;
float[][] clo1 = copy(cf);
applyRadial(ka/2.0,kb/2.0,clo1);
sub(cf,clo1,cf);
ftInverse(lev,0,clo1,spyr);
int ml2 = (nl2-1)/2+1;
int ml1 = (nl1-1)/2+1;
copy(ml1,ml2,0,0,2,2,spyr[lev][0],0,0,1,1,spyr[lev+1][0]);
for (int dir=0; dirk2) {
cf[i2][ir] = 0.0f;
cf[i2][ii] = 0.0f;
}
else if (wd>k1&&wd= k2) {
cf[i3][i2][ir] = 0.0f;
cf[i3][i2][ii] = 0.0f;
}
else if (wd>k1&&wd
* Given three values output from the three steering filters, this method
* computes two angles and corresponding steered output values such that
* the derivative of steered output with respect to angle is zero.
* Returned angles are in the range [0,PI].
* @param f0 the value for steering filter with angle 0*PI/3.
* @param f1 the value for steering filter with angle 1*PI/3.
* @param f2 the value for steering filter with angle 2*PI/3.
* @return array {{theta1,value1},{theta2,value2}}. The angle theta1
* corresponds to value1 with maximum absolute value; the angle theta2
* corresponds to value2 with minimum absolute value. Both theta1 and
* theta2 are in the range [0,PI] and the absolute difference in these
* angles is PI/2. (This method by Dave Hale 2007.12.19.)
*/
private static double[][] findExtrema(double f0, double f1, double f2) {
double theta1 = 0.5*(PI+atan2(SQRT3*(f1-f2),(2.0*f0-f1-f2)));
double value1 = eval0(f0,f1,f2,theta1);
double theta2 = modPi(theta1+0.5*PI);
double value2 = eval0(f0,f1,f2,theta2);
if (abs(value1)>=abs(value2)) {
return new double[][]{{theta1,value1},{theta2,value2}};
} else {
return new double[][]{{theta2,value2},{theta1,value1}};
}
}
// Returns the specified angle modulo PI.
// The returned angle is in the range [0,PI].
private static double modPi(double theta) {
return theta-floor(theta/PI)*PI;
}
// Steered output (zeroth derivative).
private static double eval0(double f0, double f1, double f2, double theta) {
double k0 = ONE_THIRD*(1.0+2.0*cos(2.0*(theta-THETA0)));
double k1 = ONE_THIRD*(1.0+2.0*cos(2.0*(theta-THETA1)));
double k2 = ONE_THIRD*(1.0+2.0*cos(2.0*(theta-THETA2)));
return k0*f0+k1*f1+k2*f2;
}
/**
* Finds three critical points of the 3D steering function.
* The function is
*
* f(a,b,c) = ((a+b)*(a+b)-c*c)*f[0] + ((a-b)*(a-b)-c*c)*f[1] +
* ((a+c)*(a+c)-b*b)*f[2] + ((a-c)*(a-c)-b*b)*f[3] +
* ((b+c)*(b+c)-a*a)*f[4] + ((b-c)*(b-c)-a*a)*f[5]
*
* where (a,b,c) are components of a unit-vector.
* @param f array[6] of function values.
* @param abcf array of critical points; if null, a new array is returned.
* @return array{{a0,b0,c0,f0},{a1,b1,c1,f1},{a2,b2,c2,f2}} of critical
* points. If a non-null array is specified, that array will be returned
* with possibly modified values. (This method by Dave Hale 2008.04.24.,
* revised 2008.09.23.)
*/
private static double[][] findCriticalPoints(double[] f, double[][] abcf) {
// Coefficients for evaluating the function f' = (f-fsum)/2.
// The function f' can be evaluated more efficiently and has
// the same critical points (a,b,c) as the specified function f.
double fsum = f[0]+f[1]+f[2]+f[3]+f[4]+f[5];
double fab = f[0]-f[1];
double fac = f[2]-f[3];
double fbc = f[4]-f[5];
double faa = f[4]+f[5];
double fbb = f[2]+f[3];
double fcc = f[0]+f[1];
// Initialize (a,b,c).
double a,b,c;
{
double haa,hbb,hcc,hab,hac,hbc;
// Determinant for a = 1.0, b = c = 0.0.
hbb = 2.0*(fbb-faa);
hcc = 2.0*(fcc-faa);
hbc = -fbc;
double da = hbb*hcc-hbc*hbc;
double dda = da*da;
// Determinant for b = 1.0, a = c = 0.0.
haa = 2.0*(faa-fbb);
hcc = 2.0*(fcc-fbb);
hac = -fac;
double db = haa*hcc-hac*hac;
double ddb = db*db;
// Determinant for c = 1.0, a = b = 0.0.
haa = 2.0*(faa-fcc);
hbb = 2.0*(fbb-fcc);
hab = -fab;
double dc = haa*hbb-hab*hab;
double ddc = dc*dc;
// Choose the point with largest curvature (determinant-squared).
if (dda>=ddb && dda>=ddc) {
a = 1.0; b = 0.0; c = 0.0;
} else if (ddb>=ddc) {
a = 0.0; b = 1.0; c = 0.0;
} else {
a = 0.0; b = 0.0; c = 1.0;
}
}
// Search for a critical point using Newton's method.
for (int niter=0; niter<50; ++niter) {
// We have three cases, depending on which component of the
// vector (a,b,c) has largest magnitude. For example, if that
// component is c, we express c as c(a,b) so that f is a
// function of f(a,b). We then compute a Newton update (da,db),
// constrained so that after the update a*a+b*b <= 1.0. This
// last condition is necessary so that a*a+b*b+c*c = 1.0 after
// updating (a,b,c). We eliminate the component with largest
// magnitude (in this case, c), because that choice permits the
// largest updates (da,db).
double da,db,dc;
double aa = a*a;
double bb = b*b;
double cc = c*c;
// If the component c has largest magnitude, ...
if (aa<=cc && bb<=cc) {
double aoc = a/c;
double boc = b/c;
double aocs = aoc*aoc;
double bocs = boc*boc;
double faamcc2 = (faa-fcc)*2.0;
double fbbmcc2 = (fbb-fcc)*2.0;
double aocsp1 = aocs+1.0;
double bocsp1 = bocs+1.0;
double ga = a*(faamcc2+boc*fbc)-c*(1.0-aocs)*fac-b*fab;
double gb = b*(fbbmcc2+aoc*fac)-c*(1.0-bocs)*fbc-a*fab;
double haa = faamcc2+boc*aocsp1*fbc+aoc*(3.0+aocs)*fac;
double hbb = fbbmcc2+aoc*bocsp1*fac+boc*(3.0+bocs)*fbc;
double hab = boc*aocsp1*fac+aoc*bocsp1*fbc-fab;
double det = haa*hbb-hab*hab;
if (det<=DET_SMALL && -det<=DET_SMALL) det = DET_SMALL;
double odet = 1.0/det;
da = odet*(hbb*ga-hab*gb);
db = odet*(haa*gb-hab*ga);
dc = 0.0;
for (double at=a-da,bt=b-db; at*at+bt*bt>=1.0; at=a-da,bt=b-db) {
da *= 0.5;
db *= 0.5;
}
a -= da;
b -= db;
c = (c>=0.0)?sqrt(1.0-a*a-b*b):-sqrt(1.0-a*a-b*b);
}
// Else if the component b has largest magnitude, ...
else if (aa<=bb) {
double aob = a/b;
double cob = c/b;
double aobs = aob*aob;
double cobs = cob*cob;
double faambb2 = (faa-fbb)*2.0;
double fccmbb2 = (fcc-fbb)*2.0;
double aobsp1 = aobs+1.0;
double cobsp1 = cobs+1.0;
double ga = a*(faambb2+cob*fbc)-b*(1.0-aobs)*fab-c*fac;
double gc = c*(fccmbb2+aob*fab)-b*(1.0-cobs)*fbc-a*fac;
double haa = faambb2+cob*aobsp1*fbc+aob*(3.0+aobs)*fab;
double hcc = fccmbb2+aob*cobsp1*fab+cob*(3.0+cobs)*fbc;
double hac = cob*aobsp1*fab+aob*cobsp1*fbc-fac;
double det = haa*hcc-hac*hac;
if (det<=DET_SMALL && -det<=DET_SMALL) det = DET_SMALL;
double odet = 1.0/det;
da = odet*(hcc*ga-hac*gc);
db = 0.0;
dc = odet*(haa*gc-hac*ga);
for (double at=a-da,ct=c-dc; at*at+ct*ct>=1.0; at=a-da,ct=c-dc) {
da *= 0.5;
dc *= 0.5;
}
a -= da;
c -= dc;
b = (b>=0.0)?sqrt(1.0-a*a-c*c):-sqrt(1.0-a*a-c*c);
}
// Else if the component a has largest magnitude, ...
else {
double boa = b/a;
double coa = c/a;
double boas = boa*boa;
double coas = coa*coa;
double fbbmaa2 = (fbb-faa)*2.0;
double fccmaa2 = (fcc-faa)*2.0;
double boasp1 = boas+1.0;
double coasp1 = coas+1.0;
double gb = b*(fbbmaa2+coa*fac)-a*(1.0-boas)*fab-c*fbc;
double gc = c*(fccmaa2+boa*fab)-a*(1.0-coas)*fac-b*fbc;
double hbb = fbbmaa2+coa*boasp1*fac+boa*(3.0+boas)*fab;
double hcc = fccmaa2+boa*coasp1*fab+coa*(3.0+coas)*fac;
double hbc = coa*boasp1*fab+boa*coasp1*fac-fbc;
double det = hbb*hcc-hbc*hbc;
if (det<=DET_SMALL && -det<=DET_SMALL) det = DET_SMALL;
double odet = 1.0/det;
da = 0.0;
db = odet*(hcc*gb-hbc*gc);
dc = odet*(hbb*gc-hbc*gb);
for (double bt=b-db,ct=c-dc; bt*bt+ct*ct>=1.0; bt=b-db,ct=c-dc) {
db *= 0.5;
dc *= 0.5;
}
b -= db;
c -= dc;
a = (a>=0.0)?sqrt(1.0-b*b-c*c):-sqrt(1.0-b*b-c*c);
}
// Test for convergence.
if (daf1) {
double at = a0; a0 = a1; a1 = at;
double bt = b0; b0 = b1; b1 = bt;
double ct = c0; c0 = c1; c1 = ct;
double ft = f0; f0 = f1; f1 = ft;
}
if (f1>f2) {
double at = a1; a1 = a2; a2 = at;
double bt = b1; b1 = b2; b2 = bt;
double ct = c1; c1 = c2; c2 = ct;
double ft = f1; f1 = f2; f2 = ft;
}
if (f0>f1) {
double at = a0; a0 = a1; a1 = at;
double bt = b0; b0 = b1; b1 = bt;
double ct = c0; c0 = c1; c1 = ct;
double ft = f0; f0 = f1; f1 = ft;
}
// Results.
if (abcf==null) abcf = new double[3][4];
abcf[0][0] = a0; abcf[0][1] = b0; abcf[0][2] = c0; abcf[0][3] = f0;
abcf[1][0] = a1; abcf[1][1] = b1; abcf[1][2] = c1; abcf[1][3] = f1;
abcf[2][0] = a2; abcf[2][1] = b2; abcf[2][2] = c2; abcf[2][3] = f2;
return abcf;
}
/**
* Applies forward 2D Fourier transform.
* @param level level number.
* @param x input 2D image.
* @return a 2D complex array, which is the forward 2D Fourier transform of
* the input image.
*/
private float[][] ftForward(int level, float[][] x) {
FftReal fft1;
FftComplex fft2;
int ny2 = x.length;
int ny1 = x[0].length;
int mpad = round(20.0f/(1.0f+(float)level));
int lfactor = (int)pow(2.0,(double)level);
int nl2 = (n2-1)/lfactor+1;
int nl1 = (n1-1)/lfactor+1;
int nf1 = FftReal.nfftSmall(nl1+mpad*2);
int nf1c = nf1/2+1;
int nf2 = FftComplex.nfftSmall(nl2+mpad*2);
float[][] xr = zerofloat(nf1,nf2);
copy(ny1,ny2,0,0,x,mpad,mpad,xr);
float[][] cx = czerofloat(nf1c,nf2);
fft1 = new FftReal(nf1);
fft2 = new FftComplex(nf2);
fft1.realToComplex1(1,nf2,xr,cx);
flipSign(2,cx);
fft2.complexToComplex2(1,nf1c,cx,cx);
return cx;
}
/**
* Applies forward 3D Fourier transform.
* @param level level number.
* @param x input 3D image.
* @return a 3D complex array, which is the forward 3D Fourier transform of
* the input image.
*/
private float[][][] ftForward(int level,float[][][] x) {
FftReal fft1;
FftComplex fft2;
FftComplex fft3;
int ny3 = x.length;
int ny2 = x[0].length;
int ny1 = x[0][0].length;
int mpad = round(20.0f/(1.0f+(float)level));
int lfactor = (int)pow(2.0,(double)level);
int nl3 = (n3-1)/lfactor+1;
int nl2 = (n2-1)/lfactor+1;
int nl1 = (n1-1)/lfactor+1;
int nf3 = FftComplex.nfftSmall(nl3+mpad*2);
int nf2 = FftComplex.nfftSmall(nl2+mpad*2);
int nf1 = FftReal.nfftSmall(nl1+mpad*2);
int nf1c = nf1/2+1;
float[][][] xr = zerofloat(nf1,nf2,nf3);
copy(ny1,ny2,ny3,0,0,0,x,mpad,mpad,mpad,xr);
float[][][] cx = czerofloat(nf1c,nf2,nf3);
fft1 = new FftReal(nf1);
fft2 = new FftComplex(nf2);
fft3 = new FftComplex(nf3);
fft1.realToComplex1(1,nf2,nf3,xr,cx);
flipSign(2, cx);
fft2.complexToComplex2(1,nf1c,nf3,cx,cx);
flipSign(3, cx);
fft3.complexToComplex3(1,nf1c,nf2,cx,cx);
return cx;
}
/**
* Applies inverse 2D Fourier transform to an input wavenumber-domain image.
* The output space-domain image is written to the appropriate part of the
* steerable pyramid array.
* @param lev level number.
* @param dir basis filter orientation index for the input image.
* @param cf input image in wavenumber domain (complex array).
* @param spyr input/output 2D steerable pyramid.
*/
private void ftInverse(int lev,int dir,
float[][] cf,float spyr[][][][]) {
FftReal fft1;
FftComplex fft2;
int nf2 = cf.length;
int nf1c = cf[0].length/2;
int nf1 = (nf1c-1)*2;
int mpad = round(20.0f/(1.0f+(float)lev));
int lfactor = (int)pow(2.0,(double)lev);
int nl2 = (n2-1)/lfactor+1;
int nl1 = (n1-1)/lfactor+1;
fft1 = new FftReal(nf1);
fft2 = new FftComplex(nf2);
fft2.complexToComplex2(-1,nf1c,cf,cf);
flipSign(2,cf);
fft2.scale(nf1c,nf2,cf);
fft1.complexToReal1(-1,nf2,cf,cf);
fft1.scale(nf1,nf2,cf);
copy(nl1,nl2,mpad,mpad,1,1,cf,0,0,1,1,spyr[lev][dir]);
}
/**
* Applies inverse 3D Fourier transform to an input wavenumber-domain image.
* The output space-domain image is written to the appropriate part of the
* steerable pyramid array.
* @param lev level number.
* @param dir basis filter orientation index for the input image.
* @param cf input image in wavenumber domain (complex array).
* @param spyr input/output 3D steerable pyramid.
*/
private void ftInverse(int lev,int dir,
float[][][] cf,float spyr[][][][][]) {
FftReal fft1;
FftComplex fft2;
FftComplex fft3;
int nf3 = cf.length;
int nf2 = cf[0].length;
int nf1c = cf[0][0].length/2;
int nf1 = (nf1c-1)*2;
int mpad = round(20.0f/(1.0f+(float)lev));
int lfactor = (int)pow(2.0,(double)lev);
int nl3 = (n3-1)/lfactor+1;
int nl2 = (n2-1)/lfactor+1;
int nl1 = (n1-1)/lfactor+1;
fft1 = new FftReal(nf1);
fft2 = new FftComplex(nf2);
fft3 = new FftComplex(nf3);
fft3.complexToComplex3(-1,nf1c,nf2,cf,cf);
flipSign(3, cf);
fft3.scale(nf1c,nf2,nf3,cf);
fft2.complexToComplex2(-1,nf1c,nf3,cf,cf);
flipSign(2, cf);
fft2.scale(nf1c,nf2,nf3,cf);
fft1.complexToReal1(-1,nf2,nf3,cf,cf);
fft1.scale(nf1,nf2,nf3,cf);
copy(nl1,nl2,nl3,mpad,mpad,mpad,1,1,1,cf,
0,0,0,1,1,1,spyr[lev][dir]);
}
/**
* Multiplies every other sample in input space-domain image by -1.
* Applied before forward Fourier transform where it has the effect of
* shifting the image by -PI in the wavenumber domain. This is done to
* center the image after Fourier transform. The sign-reversal is backed
* out after inverse Fourier transform.
* @param a direction of Fourier transform (1=x1 direction,2=x2).
* @param x input/output 2D array.
*/
private void flipSign(int a, float[][] x) {
int ir,ii;
int nx2 = x.length;
int nx1 = x[0].length;
if (a==1) {
for (int i2=0; i2