package gdsc.smlm.model;
import gdsc.core.utils.DoubleEquality;
import gdsc.core.utils.Maths;
import java.util.Arrays;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.integration.SimpsonIntegrator;
import org.apache.commons.math3.analysis.integration.UnivariateIntegrator;
import org.apache.commons.math3.analysis.interpolation.SplineInterpolator;
import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math3.random.RandomDataGenerator;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator;
import org.apache.commons.math3.util.FastMath;
/**
* Contains methods for generating models of a Point Spread Function using a Airy pattern.
* <p>
* Out-of-focus regions are computed using a width spreading of the Airy pattern. A true diffraction model for
* out-of-focus regions is not implemented.
*/
public class AiryPSFModel extends PSFModel
{
private double zeroW0, zeroW1;
private double w0;
private double w1;
private double zDepth = 0;
private int ring = 2;
private boolean singlePixelApproximation = false;
private int minSamplesPerDimension = 2;
private int maxSamplesPerDimension = 50;
// Used for the random sampling of the Airy function
private static int SAMPLE_RINGS = 4;
private static double[] r;
private static double[] sum;
private static PolynomialSplineFunction spline;
/**
* The zeros of J1(x) corresponding to the rings of the Airy pattern
*/
public static double[] RINGS = { 0, 3.8317, 7.0156, 10.1735, 13.3237, 16.4706 };
/**
* The Airy power corresponding to the rings of the Airy pattern
*/
public static double[] POWER;
static
{
POWER = new double[RINGS.length];
for (int i = 1; i < POWER.length; i++)
POWER[i] = AiryPattern.power(RINGS[i]);
}
/**
* @param w0
* The Airy width for dimension 0
* @param w1
* The Airy width for dimension 1
*/
public AiryPSFModel(double w0, double w1)
{
super();
this.zeroW0 = w0;
this.zeroW1 = w1;
}
/**
* @param randomGenerator
* @param w0
* The Airy width for dimension 0
* @param w1
* The Airy width for dimension 1
*/
public AiryPSFModel(RandomGenerator randomGenerator, double w0, double w1)
{
super(randomGenerator);
this.zeroW0 = w0;
this.zeroW1 = w1;
}
/**
* @param randomGenerator
* @param w0
* The Airy width for dimension 0
* @param w1
* The Airy width for dimension 1
* @param zDepth
* the Z-depth where the 3D PSF is 1.5x the width (1.5 x FWHM)
*/
public AiryPSFModel(RandomGenerator randomGenerator, double w0, double w1, double zDepth)
{
super(randomGenerator);
this.zeroW0 = w0;
this.zeroW1 = w1;
setzDepth(zDepth);
}
/**
* @param randomDataGenerator
* @param w0
* The Airy width for dimension 0
* @param w1
* The Airy width for dimension 1
*/
public AiryPSFModel(RandomDataGenerator randomDataGenerator, double w0, double w1)
{
super(randomDataGenerator);
this.zeroW0 = w0;
this.zeroW1 = w1;
}
/**
* @param randomDataGenerator
* @param w0
* The Airy width for dimension 0
* @param w1
* The Airy width for dimension 1
* @param zDepth
* the Z-depth where the 3D PSF is twice the width (2 x FWHM)
*/
public AiryPSFModel(RandomDataGenerator randomDataGenerator, double w0, double w1, double zDepth)
{
super(randomDataGenerator);
this.zeroW0 = w0;
this.zeroW1 = w1;
setzDepth(zDepth);
}
/*
* (non-Javadoc)
*
* @see gdsc.smlm.model.PSFModel#create3D(float[], int, int, double, double, double, double, boolean)
*/
public double create3D(float[] data, final int width, final int height, final double sum, double x0, double x1,
double x2, boolean poissonNoise)
{
if (sum == 0)
return 0;
final double scale = createWidthScale(x2);
try
{
final double d = airy2D(data, width, height, sum, x0, x1, scale * zeroW0, scale * zeroW1, poissonNoise);
return d;
}
catch (IllegalArgumentException e)
{
//System.out.println(e.getMessage());
return 0;
}
}
/*
* (non-Javadoc)
*
* @see gdsc.smlm.model.PSFModel#create3D(double[], int, int, double, double, double, double, boolean)
*/
public double create3D(double[] data, final int width, final int height, final double sum, double x0, double x1,
double x2, boolean poissonNoise)
{
if (sum == 0)
return 0;
final double scale = createWidthScale(x2);
try
{
return airy2D(data, width, height, sum, x0, x1, scale * zeroW0, scale * zeroW1, poissonNoise);
}
catch (IllegalArgumentException e)
{
//System.out.println(e.getMessage());
return 0;
}
}
/**
* Generate a scale so that at the configured zDepth the scale is 1.5.
*
* @param z
* @return The scale
*/
private double createWidthScale(double z)
{
if (zDepth == 0) // Not 3D data
return 1;
// PSF fitting on data from the GDSC microscope show that the PSF width spread can be modelled
// by a simple quadratic up to 1.5 times the width:
// width = 1 + z^2 / 2
// = 1.5 @ z=1
z /= zDepth; // Scale so z=1 at the configured z-depth
return 1.0 + z * z * 0.5;
}
/**
* Construct a Airy pattern on the provided data. Only evaluates the function up to the configured dark ring.
*
* @param data
* The data (can be null)
* @param width
* The data width
* @param height
* The data height
* @param sum
* The integral
* @param x0
* The centre in dimension 0
* @param x1
* The centre in dimension 1
* @param w0
* The Airy width for dimension 0
* @param w1
* The Airy width for dimension 1
* @param poissonNoise
* Add Poisson noise
* @return The total sum added to the image (useful when poissonNoise is added)
*/
public double airy2D(float[] data, final int width, final int height, final double sum, double x0, double x1,
double w0, double w1, boolean poissonNoise)
{
if (sum == 0)
return 0;
// Parameter check
if (width < 1)
throw new IllegalArgumentException("Width cannot be less than 1");
if (height < 1)
throw new IllegalArgumentException("Height cannot be less than 1");
if (data == null)
data = new float[width * height];
else if (data.length < width * height)
throw new IllegalArgumentException("Data length cannot be smaller than width * height");
w0 = Math.abs(w0);
w1 = Math.abs(w1);
// The second zero (dark ring of an Airy pattern is at 7.0156 of the width
final int x0min = clip((int) (x0 - RINGS[ring] * w0), width);
final int x1min = clip((int) (x1 - RINGS[ring] * w1), height);
final int x0max = clip((int) Math.ceil(x0 + RINGS[ring] * w0), width);
final int x1max = clip((int) Math.ceil(x1 + RINGS[ring] * w1), height);
final int x0range = x0max - x0min;
final int x1range = x1max - x1min;
// min should always be less than max
if (x0range < 1)
throw new IllegalArgumentException("Dimension 0 range not within data bounds");
if (x1range < 1)
throw new IllegalArgumentException("Dimension 1 range not within data bounds");
// Shift centre to origin and compute gaussian
double[] gauss = airy2D(x0range, x1range, sum, x0 - x0min, x1 - x1min, w0, w1);
return insert(data, x0min, x1min, x0max, x1max, width, gauss, poissonNoise);
}
/**
* Construct a Airy pattern on the provided data. Only evaluates the function up to the configured dark ring.
*
* @param data
* The data (can be null)
* @param width
* The data width
* @param height
* The data height
* @param sum
* The integral
* @param x0
* The centre in dimension 0
* @param x1
* The centre in dimension 1
* @param w0
* The Airy width for dimension 0
* @param w1
* The Airy width for dimension 1
* @param poissonNoise
* Add Poisson noise
* @return The total sum added to the image (useful when poissonNoise is added)
*/
public double airy2D(double[] data, final int width, final int height, final double sum, double x0, double x1,
double w0, double w1, boolean poissonNoise)
{
if (sum == 0)
return 0;
// Parameter check
if (width < 1)
throw new IllegalArgumentException("Width cannot be less than 1");
if (height < 1)
throw new IllegalArgumentException("Height cannot be less than 1");
if (data == null)
data = new double[width * height];
else if (data.length < width * height)
throw new IllegalArgumentException("Data length cannot be smaller than width * height");
w0 = Math.abs(w0);
w1 = Math.abs(w1);
// The second zero (dark ring of an Airy pattern is at 7.0156 of the width
final int x0min = clip((int) (x0 - RINGS[ring] * w0), width);
final int x1min = clip((int) (x1 - RINGS[ring] * w1), height);
final int x0max = clip((int) Math.ceil(x0 + RINGS[ring] * w0), width);
final int x1max = clip((int) Math.ceil(x1 + RINGS[ring] * w1), height);
final int x0range = x0max - x0min;
final int x1range = x1max - x1min;
// min should always be less than max
if (x0range < 1)
throw new IllegalArgumentException("Dimension 0 range not within data bounds");
if (x1range < 1)
throw new IllegalArgumentException("Dimension 1 range not within data bounds");
// Shift centre to origin and compute gaussian
double[] gauss = airy2D(x0range, x1range, sum, x0 - x0min, x1 - x1min, w0, w1);
return insert(data, x0min, x1min, x0max, x1max, width, gauss, poissonNoise);
}
/**
* Construct a Airy pattern on the provided data. Only evaluates the function up to the configured dark ring.
*
* @param x0range
* The maximum range in dimension 0 (width)
* @param x1range
* The maximum range in dimension 1 (height)
* @param sum
* The integral
* @param x0
* The centre in dimension 0
* @param x1
* The centre in dimension 1
* @param w0
* The Airy width for dimension 0
* @param w1
* The Airy width for dimension 1
* @return The data (packed in yx order, length = x0range * x1range)
*/
public double[] airy2D(int x0range, int x1range, double sum, double x0, double x1, double w0, double w1)
{
w0 = Math.abs(w0);
w1 = Math.abs(w1);
this.w0 = w0;
this.w1 = w1;
// Limit to nth dark ring
final double limit = RINGS[ring] * RINGS[ring];
double[] data = new double[x0range * x1range];
// Store if the Airy pattern has been clipped
boolean clipped = (x0 - RINGS[ring] * w0 < 0) || (x1 - RINGS[ring] * w1 < 0) ||
(x0 + RINGS[ring] * w0 > x0range) || (x1 + RINGS[ring] * w0 > x1range);
// Offset by pixel centres by 0.5
x0 -= 0.5;
x1 -= 0.5;
// Pre-compute the Airy intensity used for interpolation.
// Find the maximum distance from the centre to the edge of the image (normalised using the widths)
final double max = Maths.max(x0 / w0, x1 / w1, (x0range - x0) / w0, (x1range - x1) / w1);
// Find the maximum distance needed to evaluate the Airy pattern
final double maxD = FastMath.min(RINGS[ring], Math.sqrt(2 * max * max));
// Limit the total samples used for interpolation but always sample at least every pixel:
final double samplesPerPixel = FastMath.max(200 / maxD, 1);
final int maxR = (int) Math.ceil(maxD * samplesPerPixel);
double[] radius = new double[maxR + 1];
double[] intensity = new double[maxR + 1];
for (int r = 0; r <= maxR; r++)
{
// TODO - To simulate out of focus planes the intensity function can be pre-computed using
// a different equation, e.g. Born-Wolf model.
// Note that the pixel width for evaluation (e.g. the dark rings) would need to be calculated
// and a different normalisation factor would have to be calculated for clipped data. This may be
// achieved by pre-calculation of widths and normalisation factors for different z-depths.
intensity[r] = AiryPattern.intensity(r / samplesPerPixel);
radius[r] = r / samplesPerPixel;
}
double integral = 0;
// Pre-calculate x offset
double[] d0 = new double[x0range];
double[] d02 = new double[x0range];
for (int x = 0; x < x0range; x++)
{
d0[x] = (x - x0) / w0;
d02[x] = d0[x] * d0[x];
}
if (singlePixelApproximation)
{
// Single point approximation
for (int y = 0, i = 0; y < x1range; y++)
{
double d1 = (y - x1) / w1;
d1 *= d1;
for (int x = 0; x < x0range; x++, i++)
{
final double distance2 = d02[x] + d1;
if (distance2 < limit)
{
final double a = intensity(d02[x], d1, limit, samplesPerPixel, intensity, radius);
data[i] = a;
integral += a;
}
}
}
}
else
{
// Integration using Simpson's composite interval
// Set the number of subintervals adaptively, i.e. for small widths use more samples per pixel.
double nPixels = Math.PI * maxD * maxD;
double number = Math.sqrt(1000 / nPixels); // Approx 1000 (or so) samples across the image
final int nSubintervals = FastMath.max(minSamplesPerDimension,
FastMath.min(maxSamplesPerDimension, (int) Math.ceil(number * 0.5) * 2));
final double range0 = 0.5 / w0;
final double range1 = 0.5 / w1;
// Allow any point of the square pixel to be within the limit
final double pixelLimit = limit + Math.sqrt(0.5);
final double rescale = w0 * w1;
for (int y = 0, i = 0; y < x1range; y++)
{
final double d1 = (y - x1) / w1;
final double d12 = d1 * d1;
for (int x = 0; x < x0range; x++, i++)
{
final double distance2 = d02[x] + d12;
if (distance2 < pixelLimit)
{
final double a = integral(d0[x] - range0, d0[x] + range0, d1 - range1, d1 + range1, limit,
samplesPerPixel, intensity, radius, nSubintervals) * rescale;
data[i] = a;
integral += a;
}
}
}
}
//System.out.printf("Integral = %g (nPixels=%g) w0=%g, power = %g (norm = %g) = %g (clipped=%b)\n", integral,
// Math.PI * RINGS[ring] * w0 * RINGS[ring] * w1, w0, POWER[ring], integral / POWER[ring], (4 * Math.PI *
// w0 * w1), clipped);
// We must normalise the integral we calculated to the correct power of the Airy pattern,
// i.e. make the function we calculated a probability density that sums to 1.
if (clipped)
{
// Analysis has shown on unclipped data that the integral up to the nth ring is:
// integral ~ POWER[ring] * (Math.PI * 4 * w0 * w1)
// i.e. the full power of the Airy pattern is (Math.PI * 4 * w0 * w1)
sum *= 1.0 / (4 * Math.PI * w0 * w1);
}
else
{
// The integral we calculated corresponds to the power at the nth ring
sum *= POWER[ring] / integral;
}
for (int i = 0; i < data.length; i++)
data[i] *= sum;
return data;
}
/**
* Calculate the intensity of the Airy pattern at the given distances by interpolation using the lookup table
*
* @param d0
* squared distance in dimension 0
* @param d1
* squared distance in dimension 1
* @param limit
* The squared distance limit of the Airy pattern
* @param samplesPerPixel
* The number of samples per pixel of the pattern
* @param intensity
* The Airy intensity at the provided radii
* @param radius
* The radii
* @return The intensity
*/
private double intensity(final double d0, final double d1, final double limit, final double samplesPerPixel,
final double[] intensity, final double[] radius)
{
final double distance2 = d0 + d1;
if (distance2 < limit)
{
final double r = Math.sqrt(distance2);
//return AiryPattern.intensity(r);
// Interpolate the intensity at this pixel
final int index = (int) (r * samplesPerPixel);
return intensity[index] + (intensity[index + 1] - intensity[index]) * (r - radius[index]) * samplesPerPixel;
}
return 0;
}
/**
* Calculate the intensity of the Airy pattern between the specified ranges using the composite Simpson's rule
*
* @param ax
* Lower limit of x
* @param bx
* Upper limit of x
* @param ay
* Lower limit of y
* @param by
* Upper limit of y
* @param limit
* The squared distance limit of the Airy pattern
* @param samplesPerPixel
* The number of samples per pixel of the pattern
* @param intensity
* The Airy intensity at the provided radii
* @param radius
* The radii
* @param N
* The number of subintervals
* @return
*/
private double integral(final double ax, final double bx, final double ay, final double by, final double limit,
final double samplesPerPixel, final double[] intensity, final double[] radius, final int N)
{
final double h = (bx - ax) / N;
// TODO - The upper and lower bounds can be pre-computed since they are used for each pixel boundary
double s = integral(ax * ax, ay, by, limit, samplesPerPixel, intensity, radius, N) +
integral(bx * bx, ay, by, limit, samplesPerPixel, intensity, radius, N);
for (int n = 1; n < N; n += 2)
{
final double x = ax + n * h;
s += 4 * integral(x * x, ay, by, limit, samplesPerPixel, intensity, radius, N);
}
for (int n = 2; n < N; n += 2)
{
final double x = ax + n * h;
s += 2 * integral(x * x, ay, by, limit, samplesPerPixel, intensity, radius, N);
}
return s * h / 3;
}
/**
* Calculate the intensity of the Airy pattern between the specified ranges using the composite Simpson's rule
*
* @param x2
* The squared x distance
* @param ay
* Lower limit of y
* @param by
* Upper limit of y
* @param limit
* The squared distance limit of the Airy pattern
* @param samplesPerPixel
* The number of samples per pixel of the pattern
* @param intensity
* The Airy intensity at the provided radii
* @param radius
* The radii
* @param N
* The number of subintervals
* @return
*/
private double integral(final double x2, final double ay, final double by, final double limit,
final double samplesPerPixel, final double[] intensity, final double[] radius, final int N)
{
final double h = (by - ay) / N;
// TODO - The upper and lower bounds can be pre-computed since they are used for each pixel boundary
double s = intensity(x2, ay * ay, limit, samplesPerPixel, intensity, radius) +
intensity(x2, by * by, limit, samplesPerPixel, intensity, radius);
for (int n = 1; n < N; n += 2)
{
final double y = ay + n * h;
s += 4 * intensity(x2, y * y, limit, samplesPerPixel, intensity, radius);
}
for (int n = 2; n < N; n += 2)
{
final double y = ay + n * h;
s += 2 * intensity(x2, y * y, limit, samplesPerPixel, intensity, radius);
}
return s * h / 3;
}
private int clip(int x, int max)
{
if (x < 0)
x = 0;
if (x > max)
x = max;
return x;
}
/**
* @return the Z-depth where the 3D PSF is 1.5x the width (1.5 x FWHM)
*/
public double getzDepth()
{
return zDepth;
}
/**
* @param zDepth
* the Z-depth where the 3D PSF is 1.5x the width (1.5 x FWHM)
*/
public void setzDepth(double zDepth)
{
this.zDepth = Math.abs(zDepth);
}
/**
* @return The width in dimension 0 for the last drawn Airy pattern
*/
public double getW0()
{
return w0;
}
/**
* @return The width in dimension 1 for the last drawn Airy pattern
*/
public double getW1()
{
return w1;
}
/*
* (non-Javadoc)
*
* @see gdsc.smlm.model.PSFModel#getFwhm()
*/
public double getFwhm()
{
// Use the Gaussian approximation to the Airy pattern to get the FWHM
// Gaussian SD = 1.323 * Airy width
// Gaussian FWHM = SD * 2.35 = SD * (2 * Math.sqrt(2.0 * Math.log(2.0)))
return AiryPattern.FACTOR * (w0 + w1) * Math.sqrt(2.0 * Math.log(2.0));
}
/**
* @return the ring limit for the calculated Airy pattern
*/
public int getRing()
{
return ring;
}
/**
* Set the limit of the Airy pattern, defined by the dark rings where the pattern is zero. Allowed values are 1-5.
*
* @param ring
* the ring limit for the calculated Airy pattern
*/
public void setRing(int ring)
{
if (ring < RINGS.length && ring > 1)
this.ring = ring;
}
/**
* @return True if the Airy pattern is evaluated once per pixel, otherwise use Simpson's integration
*/
public boolean isSinglePixelApproximation()
{
return singlePixelApproximation;
}
/**
* @param singlePixelApproximation
* True if the Airy pattern is evaluated once per pixel, otherwise use Simpson's integration
*/
public void setSinglePixelApproximation(boolean singlePixelApproximation)
{
this.singlePixelApproximation = singlePixelApproximation;
}
/**
* @return The minimum number of samples per dimension for Simpson's integration over each pixel
*/
public int getMinSamplesPerDimension()
{
return minSamplesPerDimension;
}
/**
* Set the minimum number of samples per dimension for Simpson's integration over each pixel. Must be above 0 and is
* set to the next even number.
*
* @param n
* The minimum number of samples per dimension for Simpson's integration over each pixel
*/
public void setMinSamplesPerDimension(int n)
{
if (n >= 2)
this.minSamplesPerDimension = (n % 2 == 0) ? n : n + 1;
}
/**
* @return The maximum number of samples per dimension for Simpson's integration over each pixel
*/
public int getMaxSamplesPerDimension()
{
return maxSamplesPerDimension;
}
/**
* Set the maximum number of samples per dimension for Simpson's integration over each pixel. Must be above 0 and is
* set to the next even number.
*
* @param n
* The maximum number of samples per dimension for Simpson's integration over each pixel
*/
public void setMaxSamplesPerDimension(int n)
{
if (n >= 2)
this.maxSamplesPerDimension = (n % 2 == 0) ? n : n + 1;
}
@Override
public int sample3D(float[] data, int width, int height, int n, double x0, double x1, double x2)
{
if (n <= 0)
return insertSample(data, width, height, null, null);
final double scale = createWidthScale(x2);
double[][] sample = sample(n, x0, x1, scale * zeroW0, scale * zeroW1);
return insertSample(data, width, height, sample[0], sample[1]);
}
@Override
public int sample3D(double[] data, int width, int height, int n, double x0, double x1, double x2)
{
if (n <= 0)
return insertSample(data, width, height, null, null);
final double scale = createWidthScale(x2);
double[][] sample = sample(n, x0, x1, scale * zeroW0, scale * zeroW1);
return insertSample(data, width, height, sample[0], sample[1]);
}
/**
* Sample from an Airy distribution
*
* @param n
* The number of samples
* @param x0
* The centre in dimension 0
* @param x1
* The centre in dimension 1
* @param w0
* The Airy width for dimension 0
* @param w1
* The Airy width for dimension 1
* @return The sample x and y values
*/
public double[][] sample(final int n, final double x0, final double x1, final double w0, final double w1)
{
this.w0 = w0;
this.w1 = w1;
if (spline == null)
createAiryDistribution();
double[] x = new double[n];
double[] y = new double[n];
final RandomGenerator random = rand.getRandomGenerator();
UnitSphereRandomVectorGenerator vg = new UnitSphereRandomVectorGenerator(2, random);
int c = 0;
for (int i = 0; i < n; i++)
{
final double p = random.nextDouble();
if (p > POWER[SAMPLE_RINGS])
{
// TODO - We could add a simple interpolation here using a spline from AiryPattern.power()
continue;
}
final double r = spline.value(p);
// Convert to xy using a random vector generator
final double[] v = vg.nextVector();
x[c] = v[0] * r * w0 + x0;
y[c] = v[1] * r * w1 + x1;
c++;
}
if (c < n)
{
x = Arrays.copyOf(x, c);
y = Arrays.copyOf(y, c);
}
return new double[][] { x, y };
}
private static synchronized void createAiryDistribution()
{
if (spline != null)
return;
final double relativeAccuracy = 1e-4;
final double absoluteAccuracy = 1e-8;
final int minimalIterationCount = 3;
final int maximalIterationCount = 32;
UnivariateIntegrator integrator = new SimpsonIntegrator(relativeAccuracy, absoluteAccuracy,
minimalIterationCount, maximalIterationCount);
UnivariateFunction f = new UnivariateFunction()
{
public double value(double x)
{
// The pattern profile is in one dimension.
// Multiply by the perimeter of a circle to convert to 2D volume then normalise by 4 pi
//return AiryPattern.intensity(x) * 2 * Math.PI * x / (4 * Math.PI);
return AiryPattern.intensity(x) * 0.5 * x;
}
};
// Integrate up to a set number of dark rings
int samples = 1000;
final double step = RINGS[SAMPLE_RINGS] / samples;
double to = 0, from = 0;
r = new double[samples + 1];
sum = new double[samples + 1];
for (int i = 1; i < sum.length; i++)
{
from = to;
r[i] = to = step * i;
sum[i] = integrator.integrate(2000, f, from, to) + sum[i - 1];
}
if (DoubleEquality.relativeError(sum[samples], POWER[SAMPLE_RINGS]) > 1e-3)
throw new RuntimeException("Failed to create the Airy distribution");
SplineInterpolator si = new SplineInterpolator();
spline = si.interpolate(sum, r);
}
}