package edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.registration.math;
import ij.IJ;
/**
* BSpline library.
* 2009 Ignacio Arganda-Carreras and Arrate Munoz-Barrutia
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation (http://www.gnu.org/licenses/gpl.txt )
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*
*/
/**
* This class implements basic operations to work with cubic B-splines.
*
* @author Ignacio Arganda-Carreras - ignacio.arganda@gmail.com
*/
public class BSpline
{
// --------------------------------------------------------------
/**
* <p>
* Ref: A. Muñoz-Barrutia, T. Blu, M. Unser, "Least-Squares Image Resizing
* Using Finite Differences," IEEE Transactions on Image Processing, vol.
* 10, n. 9, September 2001.
* Use:
* out(x) = input(x/scale); scale > 1 expansion; a < 1 reduction.
* Use n=n2
* if (n1==-1) the standard method (interpolation and resampling) is
* applied.
* If (n=n1=n2), the orthogonal projection is computed. The error is the
* minimum in the least-squares sense.
* If ((n1 > -1) && (n1<n2)), the oblique projection is computed. The error
* is slightly greater than the case above.
* </p>
*
* @param input input signal
* @param n degree of the interpolation spline. It can vary from 0 to 3.
* @param n1 degree of the analysis spline. It can vary from -1 to 3.
* @param n2 degree of the synthesis spline. It can vary from 0 to 3.
* @param scale zoom factor (scale > 1 expansion; scale < 1 reduction)
* @param coefOrSamples if working with coefficients coefOrSamples is true; if working with samples coefOrSamples is false
* @return out signal
*/
public static double[] resize1D(
double[] input,
int n,
int n1,
int n2,
double scale,
boolean coefOrSamples)
{
int nxIn = input.length;
double bb = Math.round( ((nxIn-1) * scale) * 1 / scale);
while ((nxIn-1-bb) != 0)
{
nxIn = nxIn + 1;
bb = Math.round(((nxIn-1)*scale)*1/scale);
}
// Output signal
int nxOut = (int) Math.round((nxIn-1)*scale)+1;
double[] out = new double[nxOut];
// From samples to coefficients
double[] coef = null;
if (coefOrSamples == false)
{
if (n < 2)
coef = input;
else
{
double[] pole = tableOfPoles(n);
double tolerance = 1e-11;
coef = convertToInterpCoef(input, pole, tolerance);
}
}
// (n1+1)-running sums
double[] integ = new double[input.length];
double med = 0;
if (n1 == -1)
{
integ = input;
med = 0;
}
if (n1 == 0)
{
med = integSA(input, integ);
}
if (n1 == 1)
{
med = integSA(input, integ);
integ = integAS(integ);
}
if (n1 == 2)
{
med = integSA(input, integ);
integ = integAS(integ);
integSA(integ, integ);
}
if (n1 == 3)
{
med = integSA(input, integ);
integ = integAS(integ);
integSA(integ, integ);
integ = integAS(integ);
}
// geometric transformation and resampling
double [] resize = resampling(integ, nxIn, nxOut, scale, n1, n+n1+1);
//(n1+1)-centered finite differences
double [] fd = null;
if (n1 == -1)
fd=resize;
else if (n1 == 0)
fd = finDiffAS(resize);
else if (n1 == 1)
{
fd=finDiffSA(resize);
fd=finDiffAS(fd);
}
else if (n1==2)
{
fd=finDiffAS(resize);
fd=finDiffSA(fd);
fd=finDiffAS(fd);
}
else if (n1==3)
{
fd=finDiffSA(resize);
fd=finDiffAS(fd);
fd=finDiffSA(fd);
fd=finDiffAS(fd);
}
//Recover the output size and add the mean
final int n11 = n1+1;
final int val1 = (n11*0.5 == (int) Math.floor(n11*0.5)) ? (int) (n11*0.5) : (int) Math.floor(n11*0.5)+1;
//final int val2 = (int) Math.floor(n11*0.5);
double[] fdShort = new double[nxOut];
for(int i = 0 ; i < nxOut ; i ++)
fdShort[i] = fd[val1 + i] + med;
//fdShort = [fdShort(1+val2:nxOut) fdShort(val2:-1:1)];
// post-filtering - q=a^(-1)
double[] coefFull = null;
if ( (n1+n2+1) < 2 )
coefFull = fdShort;
else
{
double[] pole = tableOfPoles(n1+n2+1);
double tolerance = 1e-11;
coefFull=convertToInterpCoef(fdShort, pole, tolerance);
}
// from coeficients to samples
if (coefOrSamples==true)
{
if (n2 < 2)
out = coefFull;
else
{
double [] samples = tableOfSamples(n2);
out = convertToSamples(coefFull, samples);
}
}
return out;
} // end method resize1D
// --------------------------------------------------------------
/**
* Convert the spline coefficients to samples
*
* @param coef input coefficients values
* @param samples output samples
* @return
*/
public static double[] convertToSamples(double[] coef, double[] samples)
{
final int NbCoef = coef.length;
final int kn = 2*NbCoef-2;
double[] out = new double[NbCoef];
for (int k = 0; k < NbCoef; k++)
{
double yaux = coef[k] * samples[0];
for (int i = 1; i < samples.length-1; i ++)
{
int k1 = k-i;
int k2 = k+i;
if (k1 < 0)
{
k1 = -k1;
if (k1 >= NbCoef)
k1 = kn-k1;
}
if (k2 >= NbCoef)
{
k2 = 2*(NbCoef-1)-k2;
if (k2 >= NbCoef)
k2 = kn-k2;
}
yaux = yaux + coef[k1] + coef[k2] * samples[i];
}
out[k] = yaux;
}
return out;
} // end method convertToSamples
// --------------------------------------------------------------
/**
* Finite Difference. Symmetric input boundary conditions. Anti-symmetric
* output boundary conditions.
*
* @param c input signal
* @return output signal
*/
private static double[] finDiffSA(double[] c)
{
int N = c.length;
double[] y = new double[N];
for(int i = 0; i < N-1; i++)
y[i] = c[i] - c[i+1];
y[N-1] = c[N-1] - c[N-2];
for(int i = 0; i < N; i++)
y[i] *= -1;
return y;
} // end method finDiffSA
// --------------------------------------------------------------
/**
* Finite Difference. Anti-symmetric input boundary conditions. Symmetric
* output boundary conditions.
*
* @param c input signal
* @return output signal
*/
public static double[] finDiffAS(double[] c)
{
int N = c.length;
double[] y = new double[N];
for (int i = 0; i < N-1; i ++)
y[i+1] = c[i+1] - c [i];
y[0]=2*c[0];
return y;
} // end method finDiffAS
// --------------------------------------------------------------
/**
* Geometric transform and resampling
*
* @param input input signal
* @param nInput length of the input signal
* @param nOut length of the resize signal
* @param scale expansion/reduction
* @param n1 degree of the analysis spline
* @param nt (n+n1+1)
* @return re-sampled signal
*/
public static double[] resampling(double[] input, int nInput, int nOut,
double scale, int n1, int nt)
{
final boolean even = ( (n1+1)*0.5 == Math.floor((n1+1)*0.5) ) ? true : false;
int m = factorial(nt);
double[] temp = new double[nt + 1];
final int val1 = (int) (even ? (n1+1)*0.5 : (int) Math.floor((n1+1)*0.5)+1);
double newx0 = even ? (nt+1)*0.5-val1/scale : 0.5*(1/scale-1) + (nt+1)*0.5-val1/scale;
final double aa = Math.pow(scale, n1+1);
final double val2 = Math.floor((n1+1)*0.5);
double[] out = new double[ nOut+n1+1 ];
for(int k = -val1; k < (nOut + val2); k++)
{
final int val = (int) Math.ceil(-(nt+1) + newx0);
for (int s = 0; s <= nt; s++)
temp[s] = Math.pow(-newx0+(nt+1)+val+s, nt)/m;
for (int s = 0; s <= nt; s++)
for(int l = nt+1; l > 0; l--)
temp[l] -= temp[l-1];
if(even)
{
for(int i = 0; i <= nt; i++)
{
if (i + val < 0)
out[k+val1] += input[-i-val]*temp[i];
else
{
if (i+val > nInput-1)
out[k+val1] += input[2*(nInput-1)-i-val]*temp[i];
else
out[k+val1] += input[i+val] * temp[i];
}
}
}
else
{
for(int i = 0; i <= nt; i++)
{
if (i+val < 0)
out[k+val1] -= input[-i-val]*temp[i];
else
if (i+val > nInput-1)
out[k+val1] -= input[2*(nInput-1)-i-val]*temp[i];
else
out[k+val1] += input[i+val]*temp[i];
}
}
newx0 += 1/scale;
}
for(int i = 0; i < out.length; i++)
out[i] *= aa;
return out;
}// end method resampling
//-----------------------------------------------------------------------------
/**
* Factorial
* @param n
* @return n factorial
*/
public static int factorial(int n)
{
int fact = 0;
switch(n)
{
case 0:
fact = 1;
break;
case 1:
fact = 1;
break;
case 2:
fact = 2;
break;
case 3:
fact = 6;
break;
case 4:
fact = 24;
break;
case 5:
fact = 120;
break;
case 6:
fact = 720;
break;
case 7:
fact = 5040;
break;
default:
fact = 1;
for (int i=0 ; i < n; i++)
fact *= i;
}
return fact;
} // end factorial
// --------------------------------------------------------------
/**
* Begin computation of the running sums. The input boundary conditions are
* anti-symmetric. The output boundary conditions are symmetric.
*
* @param c input coefficients
* @return sums
*/
public static double[] integAS(double[] c)
{
int N = c.length;
double [] y = new double[N];
y[0] = c[0];
y[1] = 0;
for(int i = 2; i < N; i++)
y[i] = y[i-1] - c[i-1];
for(int i = 0; i < y.length; i++)
y[i] *= -1;
return y;
}
// --------------------------------------------------------------
/**
* Begin computation of the running sums. The input boundary conditions are
* symmetric. The output boundary conditions are anti-symmetric.
*
* @param c coefficients
* @param z sums
*
* @return mean value on the period (2*length(c)-2)
*/
public static double integSA(double[] c, double[] z)
{
double m = mean(c);
double[] coeffs = c.clone();
for(int i = 0; i < coeffs.length; i++)
coeffs[i] -= m;
z[0] = coeffs[0] * 0.5;
for(int i = 1; i < coeffs.length; i++)
z[i] = coeffs[i] + z[i-1];
return m;
} // end method integSA
// --------------------------------------------------------------
/**
* Calculate mean value on the period (2*length(A)-2)
* @param A
* @return mean value on the period (2*length(A)-2)
*/
public static double mean(double[] A)
{
double sum = 0;
for(int i = 0; i < A.length; i++)
sum += A[i];
double sumValue = 2 * sum - A[0] - A[A.length-1];
return sumValue/(2*A.length-2);
}// end method mean
// --------------------------------------------------------------
/**
* Convert to interpolation coefficients
* @param input
* @param poles
* @param tolerance
* @return
*/
public static double[] convertToInterpCoef(double[] input, double[] poles, double tolerance)
{
double lambda =1.0;
int dataLength = input.length;
double [] coef = new double[dataLength];
// special case required by mirror boundaries
if (dataLength == 1)
return coef;
// compute the overall gain
for(int k = 0 ; k < poles.length; k++)
lambda = lambda * (1.0 - poles[k]) * (1.0 - 1.0 / poles[k]);
// apply the gain
for(int n = 0; n < dataLength; n ++)
coef[n] = input[n] * lambda;
// loop over all poles
for (int k = 0 ; k < poles.length; k++)
{
// causal initialization
coef[0] = initialCausalCoefficient(coef, poles[k], tolerance);
// causal recursion
for (int n = 1; n < dataLength; n++)
coef[n] = coef[n] + poles[k] * coef[n - 1];
// anti-causal initialization
coef[dataLength-1] = initialAntiCausalCoefficient(coef, poles[k]);
// anti-causal recursion
for (int n = dataLength - 2 ; n > 0; n--)
coef[n] = poles[k] * coef[n + 1] - coef[n];
}
return coef;
} // end method convertToInterpCoef
// --------------------------------------------------------------
/**
* Get initial anti-causal coefficients
*
* @param c coefficients
* @param z pole
* @return value of the initial anti-causal coefficient
*/
public static double initialAntiCausalCoefficient(double[] c, double z)
{
return ((z / (z * z - 1.0)) * (z * c[c.length- 2] + c[c.length-1]));
} //end method initialAntiCausalCoefficient
// --------------------------------------------------------------
/**
* Initialization of causal coefficients for mirror boundary conditions
*
* @param c coefficients
* @param z pole
* @param tolerance accuracy
* @return value of the initial causal coefficient for the given accuracy
*/
public static double initialCausalCoefficient(double[] c, double z, double tolerance)
{
double sum = 0;
final int dataLength = c.length;
int horizon = dataLength;
if (tolerance > 0.0)
horizon = (int) (Math.ceil(Math.log(tolerance))/Math.log(Math.abs(z)));
if (horizon < dataLength)
{
// accelerated loop
double zn = z;
sum = c[0];
for (int n = 1; n < horizon; n++)
{
sum = sum + zn * c[n];
zn = zn * z;
}
}
else
{
// full loop
double zn = z;
double iz = 1.0 / z;
double z2n = Math.pow(z, (dataLength - 1));
sum = c[0] + z2n * c[dataLength-1];
z2n = z2n * z2n * iz;
for (int n = 1; n < dataLength - 1; n++)
{
sum = sum + (zn + z2n) * c[n];
zn = zn * z;
z2n = z2n * iz;
}
sum = (sum / (1.0 - zn * zn));
}
return sum;
} // end method initialCausalCoefficient
// --------------------------------------------------------------
/**
* Recover the poles from a lookup table
*
* @param splineDegree degree of the spline
* @return Pole values
*/
public static double[] tableOfPoles(int splineDegree)
{
double[] pole = null;
switch (splineDegree)
{
case 2:
pole = new double[1];
pole[0] = Math.sqrt(8.0) - 3.0;
break;
case 3:
pole = new double[1];
pole[0] = Math.sqrt(3.0) - 2.0;
break;
case 4:
pole = new double[2];
pole[0] = Math.sqrt(664.0 - Math.sqrt(438976.0)) + Math.sqrt(304.0) - 19.0;
pole[1] = Math.sqrt(664.0 + Math.sqrt(438976.0)) - Math.sqrt(304.0) - 19.0;
break;
case 5:
pole = new double[2];
pole[0] = Math.sqrt(135.0 / 2.0 - Math.sqrt(17745.0 / 4.0)) + Math.sqrt(105.0 / 4.0) - 13.0 / 2.0;
pole[1] = Math.sqrt(135.0 / 2.0 + Math.sqrt(17745.0 / 4.0)) - Math.sqrt(105.0 / 4.0) - 13.0 / 2.0;
break;
case 6:
pole = new double[3];
pole[0] = -0.48829458930304475513011803888378906211227916123938;
pole[1] = -0.081679271076237512597937765737059080653379610398148;
pole[2] = -0.0014141518083258177510872439765585925278641690553467;
break;
case 7:
pole = new double[3];
pole[0] = -0.53528043079643816554240378168164607183392315234269;
pole[1] = -0.12255461519232669051527226435935734360548654942730;
pole[2] = -0.0091486948096082769285930216516478534156925639545994;
break;
default:
IJ.error("Invalid spline degree");
}
return pole;
}// end method tableOfPoles
// --------------------------------------------------------------
/**
* Recover the B-spline samples values from a lookup table
*
* @param splineDegree degree of the spline
* @return Value of the samples
*/
public static double[] tableOfSamples(int splineDegree)
{
double[] samples = null;
switch (splineDegree)
{
case 0:
case 1:
samples = null;
break;
case 2:
samples = new double[2];
samples[0] = 6/8;
samples[1] = 1/8;
break;
case 3:
samples = new double[2];
samples[0] = 4/6;
samples[1] = 1/6;
break;
case 4:
samples = new double[3];
samples[0] = 230/384;
samples[1] = 76/384;
samples[2] = 1/384;
break;
case 5:
samples = new double[3];
samples[0] = 66/120;
samples[1] = 26/120;
samples[2] = 1/120;
break;
case 6:
samples = new double[4];
samples[0] = 23548/46080;
samples[1] = 10543/46080;
samples[2] = 722/46080;
samples[3] = 1/46080;
break;
case 7:
samples = new double[4];
samples[0] = 2416/5040;
samples[1] = 1191/5040;
samples[2] = 120/5040;
samples[3] = 1/5040;
break;
default:
IJ.error("Invalid spline degree");
}
return samples;
}// end method tableOfSamples
} // end class BSpline