List of usage examples for org.apache.commons.math3.special BesselJ value
public double value(double x) throws MathIllegalArgumentException, ConvergenceException
From source file:ch.epfl.leb.sass.models.psfs.internal.GibsonLanniPSF.java
/** * Computes a digital representation of the PSF. * // ww w .j a v a2s . c om * @param z The stage displacement. * @return An image stack of the PSF. **/ private void computeDigitalPSF(double z) { // Has a PSF has already been computed for this emitter's z-plane? Long zDiscrete = getNearestZPlane(this.eZ); if (interpolators.containsKey(zDiscrete)) { // PSF already computed for this z-plane, so don't do anything. return; } // Otherwise, compute the PSF and store the result in the hash map double x0 = (this.sizeX - 1) / 2.0D; double y0 = (this.sizeY - 1) / 2.0D; double xp = x0; double yp = y0; //ImageStack stack = new ImageStack(this.sizeX, this.sizeY); int maxRadius = (int) Math .round(Math.sqrt((this.sizeX - x0) * (this.sizeX - x0) + (this.sizeY - y0) * (this.sizeY - y0))) + 1; double[] r = new double[maxRadius * this.oversampling]; double[] h = new double[r.length]; double a = 0.0D; double b = Math.min(1.0D, this.ns / this.NA); double k0 = 2 * Math.PI / this.wavelength; double factor1 = this.MINWAVELENGTH / this.wavelength; double factor = factor1 * this.NA / 1.4; double deltaRho = (b - a) / (this.numSamples - 1); // basis construction double rho = 0.0D; double am = 0.0; double[][] Basis = new double[this.numSamples][this.numBasis]; BesselJ bj0 = new BesselJ(0); BesselJ bj1 = new BesselJ(1); for (int m = 0; m < this.numBasis; m++) { am = (3 * m + 1) * factor; // am = (3 * m + 1); for (int rhoi = 0; rhoi < this.numSamples; rhoi++) { rho = rhoi * deltaRho; Basis[rhoi][m] = bj0.value(am * rho); } } // compute the function to be approximated double ti = 0.0D; double OPD = 0; double W = 0; double[][] Coef = new double[1][this.numBasis * 2]; double[][] Ffun = new double[this.numSamples][2]; // Oil thickness. ti = (this.ti0 + z); double sqNA = this.NA * this.NA; double rhoNA2; for (int rhoi = 0; rhoi < this.numSamples; rhoi++) { rho = rhoi * deltaRho; rhoNA2 = rho * rho * sqNA; // OPD in the sample OPD = this.eZ * Math.sqrt(this.ns * this.ns - rhoNA2); // OPD in the immersion medium OPD += ti * Math.sqrt(this.ni * this.ni - rhoNA2) - this.ti0 * Math.sqrt(this.ni0 * this.ni0 - rhoNA2); // OPD in the coverslip OPD += this.tg * Math.sqrt(this.ng * this.ng - rhoNA2) - this.tg0 * Math.sqrt(this.ng0 * this.ng0 - rhoNA2); W = k0 * OPD; Ffun[rhoi][0] = Math.cos(W); Ffun[rhoi][1] = Math.sin(W); } // solve the linear system // begin....... (Using Common Math) RealMatrix coefficients = new Array2DRowRealMatrix(Basis, false); RealMatrix rhsFun = new Array2DRowRealMatrix(Ffun, false); DecompositionSolver solver = null; if (this.solverName.equals("svd")) { // slower but more accurate solver = new SingularValueDecomposition(coefficients).getSolver(); } else { // faster, less accurate solver = new QRDecomposition(coefficients).getSolver(); } RealMatrix solution = solver.solve(rhsFun); Coef = solution.getData(); // end....... double[][] RM = new double[this.numBasis][r.length]; double beta = 0.0D; double rm = 0.0D; for (int n = 0; n < r.length; n++) { r[n] = (n * 1.0 / this.oversampling); beta = k0 * this.NA * r[n] * this.resPSF; for (int m = 0; m < this.numBasis; m++) { am = (3 * m + 1) * factor; rm = am * bj1.value(am * b) * bj0.value(beta * b) * b; rm = rm - beta * b * bj0.value(am * b) * bj1.value(beta * b); RM[m][n] = rm / (am * am - beta * beta); } } // obtain one component double maxValue = 0.0D; for (int n = 0; n < r.length; n++) { double realh = 0.0D; double imgh = 0.0D; for (int m = 0; m < this.numBasis; m++) { realh = realh + RM[m][n] * Coef[m][0]; imgh = imgh + RM[m][n] * Coef[m][1]; } h[n] = realh * realh + imgh * imgh; } // Interpolate the PSF onto a 2D grid of physical coordinates double[] pixel = new double[this.sizeX * this.sizeY]; for (int x = 0; x < this.sizeX; x++) { for (int y = 0; y < this.sizeY; y++) { // Distance in pixels from the center of the array double rPixel = Math.sqrt((x - xp) * (x - xp) + (y - yp) * (y - yp)); // Find the index of the PSF h array that matches this distance int index = (int) Math.floor(rPixel * this.oversampling); // Perform a linear interpolation from h onto the current point. // (1 / this.oversampling) is the distance between samples in // the h array in units of pixels in the final output array. double value = h[index] + (h[(index + 1)] - h[index]) * (rPixel - r[index]) * this.oversampling; pixel[(x + this.sizeX * y)] = value; } } // Compute the constant that normalizes the PSF to its area double normConst = 0; for (int x = 0; x < this.sizeX; x++) { for (int y = 0; y < this.sizeY; y++) { normConst += pixel[x + this.sizeX * y] * this.resPSF * this.resPSF; } } // Compute the (discrete) cumulative distribution function. // First, compute the sums in the y-direction. double[] CDF = new double[this.sizeX * this.sizeY]; double sum = 0; for (int x = 0; x < this.sizeX; x++) { sum = 0; for (int y = 0; y < this.sizeY; y++) { // Normalize to the integrated area pixel[x + this.sizeX * y] /= normConst; // Increment the sum in the y-direction and the CDF sum += pixel[x + this.sizeX * y] * this.resPSF; CDF[x + this.sizeX * y] = sum; } } // Next, compute the sums in the x-direction. for (int y = 0; y < this.sizeY; y++) { sum = 0; for (int x = 0; x < this.sizeX; x++) { // Increment the sum in the x-direction sum += CDF[x + y * this.sizeX] * this.resPSF; CDF[x + y * this.sizeX] = sum; } } // Construct the x-coordinates double[] mgridX = new double[this.sizeX]; for (int x = 0; x < this.sizeX; x++) { mgridX[x] = (x - 0.5 * (this.sizeX - 1)) * this.resPSF; } // Construct the y-coordinates double[] mgridY = new double[this.sizeY]; for (int y = 0; y < this.sizeY; y++) { mgridY[y] = (y - 0.5 * (this.sizeY - 1)) * this.resPSF; } //stack.addSlice(new FloatProcessor(this.sizeX, this.sizeY, pixel)); //stack.addSlice(new FloatProcessor(this.sizeX, this.sizeY, CDF)); // Reshape CDF for interpolation double[][] rCDF = new double[this.sizeY][this.sizeX]; for (int x = 0; x < this.sizeX; x++) { for (int y = 0; y < this.sizeY; y++) { rCDF[y][x] = CDF[x + y * this.sizeX]; } } // Compute the interpolating spline for this PSF. this.interpCDF = new PiecewiseBicubicSplineInterpolatingFunction(mgridX, mgridY, rCDF); interpolators.put(zDiscrete, new PiecewiseBicubicSplineInterpolatingFunction(mgridX, mgridY, rCDF)); }