Java tutorial
/* * Copyright (c) 2012 Diamond Light Source Ltd. * * All rights reserved. This program and the accompanying materials * are made available under the terms of the Eclipse Public License v1.0 * which accompanies this distribution, and is available at * http://www.eclipse.org/legal/epl-v10.html */ package uk.ac.diamond.scisoft.analysis.diffraction; import java.util.ArrayList; import java.util.Arrays; import java.util.HashSet; import java.util.List; import java.util.TreeSet; import javax.measure.unit.NonSI; import javax.vecmath.Vector3d; import org.apache.commons.math3.optim.SimpleBounds; import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer; import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics; import org.eclipse.dawnsci.analysis.api.diffraction.DetectorProperties; import org.eclipse.dawnsci.analysis.api.diffraction.DiffractionCrystalEnvironment; import org.eclipse.dawnsci.analysis.api.fitting.IConicSectionFitFunction; import org.eclipse.dawnsci.analysis.api.fitting.IConicSectionFitter; import org.eclipse.dawnsci.analysis.api.monitor.IMonitor; import org.eclipse.dawnsci.analysis.api.roi.IROI; import org.eclipse.dawnsci.analysis.dataset.impl.BooleanDataset; import org.eclipse.dawnsci.analysis.dataset.impl.BooleanIterator; import org.eclipse.dawnsci.analysis.dataset.impl.Comparisons; import org.eclipse.dawnsci.analysis.dataset.impl.Dataset; import org.eclipse.dawnsci.analysis.dataset.impl.DatasetFactory; import org.eclipse.dawnsci.analysis.dataset.impl.DatasetUtils; import org.eclipse.dawnsci.analysis.dataset.impl.DoubleDataset; import org.eclipse.dawnsci.analysis.dataset.impl.Stats; import org.eclipse.dawnsci.analysis.dataset.roi.CircularROI; import org.eclipse.dawnsci.analysis.dataset.roi.EllipticalFitROI; import org.eclipse.dawnsci.analysis.dataset.roi.EllipticalROI; import org.eclipse.dawnsci.analysis.dataset.roi.PointROI; import org.eclipse.dawnsci.analysis.dataset.roi.PolylineROI; import org.eclipse.dawnsci.analysis.dataset.roi.RectangularROI; import org.slf4j.Logger; import org.slf4j.LoggerFactory; import org.uncommons.maths.combinatorics.CombinationGenerator; import uk.ac.diamond.scisoft.analysis.crystallography.HKL; import uk.ac.diamond.scisoft.analysis.fitting.Fitter; import uk.ac.diamond.scisoft.analysis.fitting.Generic1DFitter; import uk.ac.diamond.scisoft.analysis.fitting.functions.IdentifiedPeak; import uk.ac.diamond.scisoft.analysis.fitting.functions.Polynomial; import uk.ac.diamond.scisoft.analysis.roi.ROIProfile; /** * Utilities to fit powder rings */ public class PowderRingsUtils { private static Logger logger = LoggerFactory.getLogger(PowderRingsUtils.class); private static final double FULL_CIRCLE = 2.0 * Math.PI; private static final double ARC_LENGTH = 8; private static final double RADIAL_DELTA = 10; private static final int MAX_POINTS = 200; private static final int PEAK_SMOOTHING = 15; private static final double MAX_FWHM_FACTOR = 2; private static final double RING_SEPARATION = 4; public static PolylineROI findPOIsNearCircle(IMonitor mon, Dataset image, BooleanDataset mask, CircularROI circle) { return findPOIsNearCircle(mon, image, mask, circle, ARC_LENGTH, RADIAL_DELTA, MAX_POINTS); } public static PolylineROI findPOIsNearCircle(IMonitor mon, Dataset image, BooleanDataset mask, CircularROI circle, double arcLength, double radialDelta, int maxPoints) { return findPOIsNearEllipse(mon, image, mask, new EllipticalROI(circle), arcLength, radialDelta, maxPoints); } public static PolylineROI findPOIsNearEllipse(IMonitor mon, Dataset image, BooleanDataset mask, EllipticalROI ellipse) { return findPOIsNearEllipse(mon, image, mask, ellipse, ARC_LENGTH, RADIAL_DELTA, MAX_POINTS); } /** * Find a set of points of interests near given ellipse from an image. * <p> * The ellipse is divided into sub-areas and these POIs are considered to * be the locations of maximum pixel intensities found within those sub-areas. * @param mon * @param image * @param mask (can be null) * @param ellipse * @param arcLength step size along arc in pixels * @param radialDelta +/- value to define area to search * @param maxPoints maximum number of points to return * @return polyline ROI */ public static PolylineROI findPOIsNearEllipse(IMonitor mon, Dataset image, BooleanDataset mask, EllipticalROI ellipse, double arcLength, double radialDelta, int maxPoints) { if (image.getRank() != 2) { logger.error("Dataset must have two dimensions"); throw new IllegalArgumentException("Dataset must have two dimensions"); } if (mask != null && !image.isCompatibleWith(mask)) { logger.error("Mask must match image shape"); throw new IllegalArgumentException("Mask must match image shape"); } final int[] shape = image.getShape(); final int h = shape[0]; final int w = shape[1]; if (ellipse.containsPoint(-1, -1) && ellipse.containsPoint(-1, h + 1) && ellipse.containsPoint(w + 1, h + 1) && ellipse.containsPoint(w + 1, -1)) { throw new IllegalArgumentException("Ellipse does not intersect image!"); } final double aj = ellipse.getSemiAxis(0); final double an = ellipse.getSemiAxis(1); if (an < arcLength) { logger.error("Ellipse/circle is too small"); throw new IllegalArgumentException("Ellipse/circle is too small"); } final double xc = ellipse.getPointX(); final double yc = ellipse.getPointY(); final double ang = ellipse.getAngle(); final double ca = Math.cos(ang); final double sa = Math.sin(ang); final double pdelta = arcLength / aj; // change in angle double rdelta = radialDelta; // semi-width of annulus of interest if (rdelta < 1) { logger.warn("Radial delta was set too low: setting to 1"); rdelta = 1; } final double rsj = aj - rdelta; final double rej = aj + rdelta; final double rsn = an - rdelta; final double ren = an + rdelta; final int imax = (int) Math.ceil(FULL_CIRCLE / pdelta); logger.debug("Major semi-axis = [{}, {}]; {}", new Object[] { rsj, rej, imax }); final int[] start = new int[2]; final int[] stop = new int[2]; final int[] step = new int[] { 1, 1 }; HashSet<PointROI> pointSet = new HashSet<PointROI>(); for (int i = 0; i < imax; i++) { double p = i * pdelta; double cp = Math.cos(p); double sp = Math.sin(p); Dataset sub; final int[] beg = new int[] { (int) (yc + rsj * sa * cp + rsn * ca * sp), (int) (xc + rsj * ca * cp - rsn * sa * sp) }; final int[] end = new int[] { (int) (yc + rej * sa * cp + ren * ca * sp), (int) (xc + rej * ca * cp - ren * sa * sp) }; start[0] = Math.max(0, Math.min(beg[0], end[0])); stop[0] = Math.min(h, Math.max(beg[0], end[0])); if (start[0] == stop[0]) { if (stop[0] == h) { start[0]--; } else { stop[0]++; } } else if (start[0] > stop[0] || start[0] >= h) { continue; } else { stop[0] = Math.min(Math.max(stop[0], start[0] + (int) radialDelta), h); } start[1] = Math.max(0, Math.min(beg[1], end[1])); stop[1] = Math.min(w, Math.max(beg[1], end[1])); if (start[1] == stop[1]) { if (stop[1] == w) { start[1]--; } else { stop[1]++; } } else if (start[1] > stop[1] || start[1] >= w) { continue; } else { stop[1] = Math.min(Math.max(stop[1], start[1] + (int) radialDelta), w); } sub = image.getSlice(start, stop, step); // TODO ensure slice has peaky data double iqr = (Double) Stats.iqr(sub); double low = (Double) sub.mean() + 0.5 * iqr; if (sub.max().doubleValue() < low) { logger.info("Discard sub at {} ({}): {}; {}; {} [{}] => {} cf {}", new Object[] { Arrays.toString(start), sub.getShape(), sub.min(), sub.mean(), sub.max(), low, 0.5 * iqr, sub.stdDeviation() }); continue; } int[] pos = sub.maxPos(); pos[0] += start[0]; pos[1] += start[1]; if (mask != null) { if (mask.get(pos)) { Dataset sorted = DatasetUtils.sort(sub.flatten(), null); int l = sorted.getSize() - 1; do { double x = sorted.getElementDoubleAbs(l); pos = sub.getNDPosition(DatasetUtils.findIndexEqualTo(sub, x)); pos[0] += start[0]; pos[1] += start[1]; } while (!mask.get(pos) && --l >= 0); if (l < 0) { logger.warn("Could not find unmasked value for slice!"); } else { //add 0.5 to make pointROI at centre of pixel // pointSet.add(new PointROI(pos[1], pos[0])); pointSet.add(new PointROI(pos[1] + 0.5, pos[0] + 0.5)); } } } else { // System.err.printf("Slice: %s, %s has max at %s\n", Arrays.toString(start), Arrays.toString(stop), Arrays.toString(pos)); //add 0.5 to make pointROI at centre of pixel pointSet.add(new PointROI(pos[1] + 0.5, pos[0] + 0.5)); // pointSet.add(new PointROI(pos[1], pos[0])); } if (mon != null) mon.worked(1); } // analyse pixel values int n = pointSet.size(); double[] values = new double[n]; int i = 0; for (PointROI p : pointSet) { int[] pos = p.getIntPoint(); values[i++] = image.getDouble(pos[1], pos[0]); } DoubleDataset pixels = new DoubleDataset(values); // System.err.println(pixels.toString(true)); // threshold with population stats from maxima logger.debug("Stats: {} {} {} {}", new Object[] { pixels.min(), pixels.mean(), pixels.max(), Arrays.toString(Stats.quantile(pixels, 0.25, 0.5, 0.75)) }); double threshold; if (n > maxPoints) { threshold = Stats.quantile(pixels, 1 - maxPoints / (double) n); logger.debug("Threshold: {} setting for highest {}", threshold, maxPoints); } else { double iqr = (Double) Stats.iqr(pixels); threshold = (Double) pixels.mean() - 2 * iqr; double min = pixels.min().doubleValue(); while (threshold < min) { threshold += iqr; } logger.debug("Threshold: {} setting by mean - 2IQR", threshold); } PolylineROI polyline = new PolylineROI(); if (threshold > (Double) pixels.min()) { for (PointROI p : pointSet) { int[] pos = p.getIntPoint(); double v = image.getDouble(pos[1], pos[0]); if (v >= threshold) { // System.err.printf("Adding %f %s\n", v, Arrays.toString(pos)); polyline.insertPoint(p); // } else { // System.err.println("Rejecting " + p + " = " + v); } } } else { for (PointROI p : pointSet) { polyline.insertPoint(p); } } if (mon != null) mon.worked(polyline.getNumberOfPoints()); logger.debug("Used {} of {} pixels", polyline.getNumberOfPoints(), pointSet.size()); return polyline; } public static EllipticalFitROI fitAndTrimOutliers(IMonitor mon, PolylineROI points, boolean circleOnly) { return fitAndTrimOutliers(mon, points, RADIAL_DELTA, circleOnly); } /** * Fit an ellipse to given points, trim points if they fall outside given distance * and re-fit * @param points * @param trimDelta trim distance * @param circleOnly if true, then fit a circle * @return fitted ellipse */ public static EllipticalFitROI fitAndTrimOutliers(IMonitor mon, PolylineROI points, double trimDelta, boolean circleOnly) { try { EllipticalFitROI efroi = new EllipticalFitROI(points, circleOnly); PolylineROI cpts = points; int n = cpts.getNumberOfPoints(); IConicSectionFitter f = efroi.getFitter(); IConicSectionFitFunction fn = f.getFitFunction(null, null); Dataset d = DatasetUtils.convertToDataset(fn.calcDistanceSquared(f.getParameters())); // find outliers double h = trimDelta * trimDelta; double ds = d.max().doubleValue(); logger.debug("Range: [0, {}] cf [{}, {}, {}]", new Object[] { h, d.min(), d.mean(), d.max() }); if (ds < h) return efroi; BooleanDataset b = Comparisons.lessThanOrEqualTo(d, h); BooleanIterator it = d.getBooleanIterator(b, true); PolylineROI npts = new PolylineROI(); while (it.hasNext()) { npts.insertPoint(cpts.getPoint(it.index)); } int m = npts.getNumberOfPoints(); if (m < n) { logger.debug("Found some outliers: {}/{}", n - m, n); efroi.setPoints(npts); if (mon != null) mon.worked(m); } return efroi; } catch (Exception e) { logger.error("Problem with trimming: {}", e); throw new IllegalArgumentException("Problem!: ", e); } } /** * Find other ellipses from given ellipse and image. * <p> * This is done by looking at the box profile along spokes from the * given centre (from a minimum distance of the minor semi-axis) and finding peaks. * Then the distance out to those peaks is used to search for more POIs and so more ellipses * @param image * @param mask (can be null) * @param roi initial ellipse * @return list of ellipses */ public static List<EllipticalROI> findOtherEllipses(IMonitor mon, Dataset image, BooleanDataset mask, EllipticalROI roi) { return findOtherEllipses(mon, image, mask, roi, roi.getSemiAxis(1), RADIAL_DELTA, ARC_LENGTH, 0.3 * RADIAL_DELTA, MAX_POINTS); } /** * Find other ellipses from given ellipse and image. * <p> * This is done by looking at the box profile along spokes from the * given centre and finding peaks. Then the distance out to those peaks is used * to search for more POIs and so more ellipses * @param image * @param mask (can be null) * @param roi initial ellipse * @param radialMin * @param radialDelta * @param arcLength * @param trimDelta * @param maxPoints * @return list of ellipses */ public static List<EllipticalROI> findOtherEllipses(IMonitor mon, Dataset image, BooleanDataset mask, EllipticalROI roi, double radialMin, double radialDelta, double arcLength, double trimDelta, int maxPoints) { if (image.getRank() != 2) { logger.error("Dataset must have two dimensions"); throw new IllegalArgumentException("Dataset must have two dimensions"); } if (mask != null && !image.isCompatibleWith(mask)) { logger.error("Mask must match image shape"); throw new IllegalArgumentException("Mask must match image shape"); } // explore all corners final int[] shape = image.getShape(); final int h = shape[0]; final int w = shape[1]; double[] ec = roi.getPoint(); TreeSet<Double> majors = new TreeSet<Double>(); // TODO farm this out across several threads findMajorAxes(mon, majors, image, mask, roi, radialMin, radialDelta, ec, 0 - ec[0], 0 - ec[1]); // TL findMajorAxes(mon, majors, image, mask, roi, radialMin, radialDelta, ec, w - ec[0], 0 - ec[1]); // TR findMajorAxes(mon, majors, image, mask, roi, radialMin, radialDelta, ec, w - ec[0], h - ec[1]); // BR findMajorAxes(mon, majors, image, mask, roi, radialMin, radialDelta, ec, 0 - ec[0], h - ec[1]); // BL findMajorAxes(mon, majors, image, mask, roi, radialMin, radialDelta, ec, 0, h - ec[1]); // T findMajorAxes(mon, majors, image, mask, roi, radialMin, radialDelta, ec, w - ec[0], 0); // R findMajorAxes(mon, majors, image, mask, roi, radialMin, radialDelta, ec, 0, 0 - ec[1]); // B findMajorAxes(mon, majors, image, mask, roi, radialMin, radialDelta, ec, 0 - ec[0], 0); // L // and finally find POIs List<EllipticalROI> ells = new ArrayList<EllipticalROI>(); double major = roi.getSemiAxis(0); double aspect = roi.getSemiAxis(0) / roi.getSemiAxis(1); double last = Double.NEGATIVE_INFINITY; for (double a : majors) { System.err.println("Current " + a + ", last " + last); if (a < last) { System.err.println("Dropped as less than last"); continue; } if (Math.abs(a - last) < RING_SEPARATION) { // omit close rings last = a; System.err.println("Dropped as too close"); continue; } if (Math.abs(a - major) < RING_SEPARATION) { last = major; System.err.println("Add original"); ells.add(roi); } else { EllipticalROI er = new EllipticalROI(a, a / aspect, roi.getAngle(), ec[0], ec[1]); try { PolylineROI polyline = findPOIsNearEllipse(mon, image, mask, er, arcLength, 0.8 * radialDelta, maxPoints); if (polyline.getNumberOfPoints() > 2) { er = fitAndTrimOutliers(mon, polyline, trimDelta, roi.isCircular()); double emaj = er.getSemiAxis(0); if (Math.abs(emaj - last) < RING_SEPARATION) { // omit close rings last = a; System.err.println("Dropped as fit is too close"); continue; } double[] c = er.getPointRef(); if (Math.hypot(c[0] - ec[0], c[1] - ec[1]) > 8 * radialDelta) { last = a; // omit fits with far-off centres System.err.println("Dropped as centre is far-off"); continue; } if (Math.abs(emaj - major) < RING_SEPARATION) { System.err.println("Add fit that is close to original"); } last = Math.max(a, emaj); ells.add(er); } else { logger.warn("Could not find enough points at {}", er); } } catch (IllegalArgumentException e) { logger.debug("Problem with {}", er, e); last = a; } if (mon != null) mon.worked(1); } } return ells; } /** * Find major axes by looking along thick line given by relative coordinates to centre for * maximum intensity values * @param mon * @param axes * @param image * @param mask * @param roi * @param offset minimum position of peaks * @param width of line * @param centre * @param dx * @param dy */ private static void findMajorAxes(IMonitor mon, TreeSet<Double> axes, Dataset image, Dataset mask, EllipticalROI roi, double offset, double width, double[] centre, double dx, double dy) { RectangularROI rroi = new RectangularROI(); rroi.setPoint(centre); rroi.setAngle(Math.atan2(dy, dx)); rroi.setLengths(Math.hypot(dx, dy), width); rroi.translate(0, -0.5); rroi.setClippingCompensation(true); Dataset profile = ROIProfile.maxInBox(image, mask, rroi)[0]; List<IdentifiedPeak> peaks = Generic1DFitter .findPeaks(DatasetFactory.createRange(profile.getSize(), Dataset.INT), profile, PEAK_SMOOTHING); if (mon != null) mon.worked(profile.getSize()); System.err.printf("\n"); DescriptiveStatistics stats = new DescriptiveStatistics(); int[] pb = new int[1]; int[] pe = new int[1]; for (IdentifiedPeak p : peaks) { if (p.getPos() < offset) { continue; } pb[0] = (int) p.getMinXVal(); pe[0] = (int) p.getMaxXVal(); p.setArea((Double) profile.getSlice(pb, pe, null).sum()); stats.addValue(p.getArea()); System.err.printf("P %f A %f W %f H %f\n", p.getPos(), p.getArea(), p.getFWHM(), p.getHeight()); } double area = stats.getMean() + 0.4 * (stats.getPercentile(75) - stats.getPercentile(25)); logger.debug("Area: {}", stats); logger.debug("Minimum threshold: {}", area); double majorFactor = roi.getSemiAxis(0) / roi.getDistance(rroi.getAngle()); double maxFWHM = MAX_FWHM_FACTOR * width; for (IdentifiedPeak p : peaks) { double l = p.getPos(); if (l < offset) { continue; } // System.err.println(p); // filter on area and FWHM if (p.getFWHM() > maxFWHM) { continue; } if (p.getArea() < area) { break; } axes.add(l * majorFactor); } if (mon != null) mon.worked(peaks.size()); } /** * Fit ellipses to a single detector with/without fixing wavelength. * <p> * Combinations of spacings are used to fit the ellipses (if there are more ellipses than spacings then * the outer ellipses are used) * @param mon * @param detector * @param env * @param ellipses * @param spacings * a list of possible spacings * @param fixedWavelength * @return q-space */ public static QSpace fitEllipsesToQSpace(IMonitor mon, DetectorProperties detector, DiffractionCrystalEnvironment env, List<EllipticalROI> ellipses, List<HKL> spacings, boolean fixedWavelength) { int n = ellipses.size(); double dmax = spacings.get(0).getD().doubleValue(NonSI.ANGSTROM); { double rmin = detector.getDetectorDistance() * Math.tan(2.0 * Math.asin(0.5 * env.getWavelength() / dmax)) / detector.getVPxSize(); int l = 0; for (EllipticalROI e : ellipses) { if (e.getSemiAxis(0) > rmin) break; l++; } if (l >= n) { throw new IllegalArgumentException("Maybe all rings are too small!"); } if (l > 0) { logger.debug("Discarding first {} rings", l); ellipses = ellipses.subList(l, n); n = ellipses.size(); } } if (n >= spacings.size()) { // always allow a choice to be made int l = n - spacings.size(); logger.warn("The number of d-spacings ({}) should be greater than or equal to {}: using outer rings", l, n - 1); ellipses = ellipses.subList(l + 1, n); n = ellipses.size(); } logger.debug("Using {} rings:", n); boolean allCircles = true; for (int i = 0; i < n; i++) { EllipticalROI e = ellipses.get(i); logger.debug(" {}", e); if (allCircles && !e.isCircular()) { allCircles = false; } } DetectorFitFunction f; if (allCircles) { logger.debug("All rings are circular"); f = createQFitFunction4(ellipses, detector, env.getWavelength(), fixedWavelength); } else { f = createQFitFunction7(ellipses, detector, env.getWavelength(), fixedWavelength); } logger.debug("Init: {}", f.getInitial()); // set up a combination generator for all List<Double> s = new ArrayList<Double>(); for (int i = 0, imax = spacings.size(); i < imax; i++) { HKL d = spacings.get(i); s.add(d.getD().doubleValue(NonSI.ANGSTROM)); } CombinationGenerator<Double> gen = new CombinationGenerator<Double>(s, n); if (mon != null) { mon.worked(1); } logger.debug("There are {} combinations", gen.getTotalCombinations()); double min = Double.POSITIVE_INFINITY; MultivariateOptimizer opt = FittingUtils.createOptimizer(f.getN()); // opt.setMaxEvaluations(2000); List<Double> fSpacings = null; int i = 0; for (List<Double> list : gen) { // find combination that minimizes residuals f.setSpacings(list); double res = FittingUtils.optimize(f, opt, min); if (res < min) { min = res; fSpacings = list; } // System.err.printf("."); if (mon != null) { mon.worked(10); if (mon.isCancelled()) return null; } if (i++ == 100) { // System.err.printf("\n"); i = 0; } } // System.err.printf("\n"); if (fSpacings == null || f.getParameters() == null) { logger.warn("Problem with fitting - as could not find a single fit!"); return null; } logger.debug("Parameters: w {}, D {}, e {} (min {})", new Object[] { f.getWavelength(), f.getDistance(), f.getNormalAngles(), min }); logger.debug("Spacings used: {}", fSpacings); f.setSpacings(fSpacings); logger.debug("Residual value: {}", f.value(f.getParameters())); QSpace q = new QSpace(f.getDetectorProperties().get(0).clone(), new DiffractionCrystalEnvironment(f.getWavelength())); q.setResidual(f.value(f.getParameters())); return q; } /** * Fit ellipses to a single detector with/without fixing wavelength. * <p> * The list of ROIs should be of equal size to the list of d-spacings and can have null * entries if there are no corresponding figures on the detector * @param mon * @param detector * @param env * @param rois * @param spacings * a list of spacings * @param fixedWavelength * @return q-space */ public static QSpace fitAllEllipsesToQSpace(IMonitor mon, DetectorProperties detector, DiffractionCrystalEnvironment env, List<? extends IROI> rois, List<HKL> spacings, boolean fixedWavelength) { int n = rois.size(); if (n != spacings.size()) { // always allow a choice to be made throw new IllegalArgumentException("Number of ellipses should be equal to spacings"); } // build up list of ellipses and spacings boolean allCircles = true; List<EllipticalROI> ellipses = new ArrayList<EllipticalROI>(); List<Double> s = new ArrayList<Double>(); for (int i = 0; i < n; i++) { IROI r = rois.get(i); if (r instanceof CircularROI) { r = new EllipticalROI((CircularROI) r); } else if (!(r instanceof EllipticalROI)) continue; EllipticalROI e = (EllipticalROI) r; ellipses.add(e); HKL d = spacings.get(i); s.add(d.getD().doubleValue(NonSI.ANGSTROM)); logger.debug(" {}", e); if (allCircles && !e.isCircular()) { allCircles = false; } } n = ellipses.size(); DetectorFitFunction f; if (allCircles) { logger.debug("All rings are circular"); f = createQFitFunction4(ellipses, detector, env.getWavelength(), fixedWavelength); } else { f = createQFitFunction7(ellipses, detector, env.getWavelength(), fixedWavelength); } logger.debug("Init: {}", f.getInitial()); if (mon != null) { mon.worked(1); if (mon.isCancelled()) return null; } MultivariateOptimizer opt = FittingUtils.createOptimizer(f.getN()); f.setSpacings(s); double res = FittingUtils.optimize(f, opt, Double.POSITIVE_INFINITY); if (mon != null) { mon.worked(10); if (mon.isCancelled()) return null; } if (f.getParameters() == null) { logger.warn("Problem with fitting - as could not find a single fit!"); } logger.debug("Parameters: w {}, D {}, e {} (min {})", new Object[] { f.getWavelength(), f.getDistance(), f.getNormalAngles(), res }); logger.debug("Spacings used: {}", s); logger.debug("Residual value: {}", f.value(f.getParameters())); QSpace q = new QSpace(f.getDetectorProperties().get(0).clone(), new DiffractionCrystalEnvironment(f.getWavelength())); q.setResidual(f.value(f.getParameters())); return q; } private static double calcBaseRollAngle(List<EllipticalROI> ellipses) { int n = ellipses.size(); double base = 0; for (int i = 0; i < n; i++) { double a = ellipses.get(i).getAngle(); base += a > Math.PI ? a - FULL_CIRCLE : a; // ensure angle is in range +/- pi } return base / n; } /** * Fit all ellipses for a list of detectors to a single wavelength that is initially fixed * <p> * Each list of ROIs should be of equal size to the list of d-spacings and can have null * entries if there are no corresponding figures on the detector * @param mon * @param lDetectors * @param env * @param lROIs list containing lists of ROIs for each detector * @param spacings * a list of spacings * @param postFitChange if true then linearly fit wavelength after detector fits * @return a list of q-spaces */ public static List<QSpace> fitAllEllipsesToAllQSpacesAtFixedWavelength(IMonitor mon, List<DetectorProperties> lDetectors, DiffractionCrystalEnvironment env, List<List<? extends IROI>> lROIs, List<HKL> spacings, boolean postFitChange) { int n = lDetectors.size(); if (n != lROIs.size()) { throw new IllegalArgumentException("Number of detectors must match number of lists of ROIs"); } List<QSpace> qs = new ArrayList<QSpace>(); for (int i = 0; i < n; i++) { DetectorProperties dp = lDetectors.get(i); List<? extends IROI> rois = lROIs.get(i); QSpace q = null; try { q = fitAllEllipsesToQSpace(mon, dp, env, rois, spacings, true); } catch (IllegalArgumentException e) { logger.warn("Problem in calibrating image: {}", i, e); } qs.add(q); } if (!postFitChange) return qs; List<Double> odist = new ArrayList<Double>(); List<Double> ndist = new ArrayList<Double>(); for (int i = 0; i < n; i++) { DetectorProperties dp = lDetectors.get(i); QSpace q = qs.get(i); odist.add(dp.getDetectorDistance()); ndist.add(q.getDetectorProperties().getDetectorDistance()); } if (odist.size() < 3) { logger.warn("Need to use three or more images"); return qs; } Polynomial p; try { p = Fitter.polyFit(new Dataset[] { DatasetFactory.createFromList(odist) }, DatasetFactory.createFromList(ndist), 1e-15, 1); } catch (Exception e) { logger.error("Problem with fit", e); return qs; } logger.debug("Straight line fit: 0 {}, 1 {}", p.getParameter(0), p.getParameter(1)); double l = env.getWavelength() * p.getParameterValue(0); for (int i = 0; i < n; i++) { qs.get(i).setDiffractionCrystalEnvironment(new DiffractionCrystalEnvironment(l)); } return qs; } /** * Fit all ellipses for a list of detectors to a single wavelength. * <p> * Each list of ROIs should be of equal size to the list of d-spacings and can have null * entries if there are no corresponding figures on the detector * @param mon * @param lDetectors * @param env * @param lROIs list containing lists of ROIs for each detector * @param spacings * a list of spacings * @return a list of q-spaces */ public static List<QSpace> fitAllEllipsesToAllQSpaces(IMonitor mon, List<DetectorProperties> lDetectors, DiffractionCrystalEnvironment env, List<List<? extends IROI>> lROIs, List<HKL> spacings) { int m = lDetectors.size(); if (lROIs.size() != m) { throw new IllegalArgumentException( "Number of lists of ROIs should be equal to number of detectors/images"); } // build up list of lists of ellipses and flattened list of spacings List<List<EllipticalROI>> lEllipses = new ArrayList<List<EllipticalROI>>(); List<Double> s = new ArrayList<Double>(); for (int j = 0; j < m; j++) { List<? extends IROI> rois = lROIs.get(j); int n = rois.size(); if (n != spacings.size()) { // always allow a choice to be made throw new IllegalArgumentException("Number of ellipses should be equal to spacings"); } List<EllipticalROI> ellipses = new ArrayList<EllipticalROI>(); for (int i = 0; i < n; i++) { IROI r = rois.get(i); if (r instanceof CircularROI) { r = new EllipticalROI((CircularROI) r); } else if (!(r instanceof EllipticalROI)) continue; EllipticalROI e = (EllipticalROI) r; ellipses.add(e); HKL d = spacings.get(i); s.add(d.getD().doubleValue(NonSI.ANGSTROM)); logger.debug(" {}", e); } lEllipses.add(ellipses); } DetectorFitFunction f = createQFitFunctionForAllImages(lEllipses, lDetectors, env.getWavelength()); logger.debug("Init: {}", f.getInitial()); if (mon != null) { mon.worked(1); if (mon.isCancelled()) return null; } MultivariateOptimizer opt = FittingUtils.createOptimizer(f.getN()); f.setSpacings(s); double res = FittingUtils.optimize(f, opt, Double.POSITIVE_INFINITY); if (mon != null) { mon.worked(10); if (mon.isCancelled()) return null; } if (f.getParameters() == null) { logger.warn("Problem with fitting - as could not find a single fit!"); } logger.debug("Parameters: w {}, D {}, e {} (min {})", new Object[] { f.getWavelength(), f.getDistances(), f.getNormalAngles(), res }); logger.debug("Spacings used: {}", s); logger.debug("Residual value: {}", f.value(f.getParameters())); List<QSpace> qs = new ArrayList<QSpace>(); DiffractionCrystalEnvironment nEnv = new DiffractionCrystalEnvironment(f.getWavelength()); for (DetectorProperties d : f.getDetectorProperties()) { QSpace q = new QSpace(d.clone(), nEnv.clone()); q.setResidual(f.value(f.getParameters())); qs.add(q); } return qs; } /** * Create function which uses 6/7 parameters: wavelength (Angstrom), detector origin (mm), orientation angles (degrees) */ static DetectorFitFunction createQFitFunction7(List<EllipticalROI> ellipses, DetectorProperties dp, double wavelength, boolean fixedWavelength) { int n = ellipses.size(); double[][] known = new double[n][FitFunctionBase.nC]; double[] weight = new double[n]; double base = -calcBaseRollAngle(ellipses); logger.debug("Mean roll angle: {}", Math.toDegrees(base)); for (int i = 0; i < n; i++) { EllipticalROI e = ellipses.get(i); weight[i] = e.getSemiAxis(0); int j = 0; double a = base - e.getAngle(); for (double off : FitFunctionBase.angles) { double[] pt = e.getPoint(a + off); known[i][j++] = pt[0]; known[i][j++] = pt[1]; } } DetectorFitFunction f; Vector3d o = dp.getOrigin(); double[] a = dp.getNormalAnglesInDegrees(); if (fixedWavelength) { f = new QSpaceFitFixedWFunction7(known, weight, dp.getVPxSize(), wavelength); f.setInitial(o.getX(), o.getY(), o.getZ(), a[0], a[1], a[2]); } else { f = new QSpaceFitFunction7(known, weight, dp.getVPxSize()); f.setInitial(wavelength, o.getX(), o.getY(), o.getZ(), a[0], a[1], a[2]); } f.setBaseRollAngle(base); return f; } /** * Create function which uses 3/4 parameters: wavelength (Angstrom), detector origin (mm) */ static DetectorFitFunction createQFitFunction4(List<EllipticalROI> ellipses, DetectorProperties dp, double wavelength, boolean fixedWavelength) { int n = ellipses.size(); double[][] known = new double[n][FitFunctionBase.nC]; double[] weight = new double[n]; for (int i = 0; i < n; i++) { EllipticalROI e = ellipses.get(i); weight[i] = e.getSemiAxis(0); int j = 0; for (double off : FitFunctionBase.angles) { double[] pt = e.getPoint(off); known[i][j++] = pt[0]; known[i][j++] = pt[1]; } } DetectorFitFunction f; Vector3d o = dp.getOrigin(); if (fixedWavelength) { f = new QSpaceFitFixedWFunction4(known, weight, dp.getVPxSize(), wavelength); f.setInitial(o.getX(), o.getY(), o.getZ()); } else { f = new QSpaceFitFunction4(known, weight, dp.getVPxSize()); f.setInitial(wavelength, o.getX(), o.getY(), o.getZ()); } f.setBaseRollAngle(ellipses.get(0).getAngle()); return f; } /** * Create function which uses 6N+1 parameters: wavelength (Angstrom), and per image, detector origin (mm), orientation angles (degrees) */ static DetectorFitFunction createQFitFunctionForAllImages(List<List<EllipticalROI>> lEllipses, List<DetectorProperties> lDP, double wavelength) { int m = lEllipses.size(); if (lDP.size() != m) { throw new IllegalArgumentException( "Number of lists of ellipses should be equal to number of detectors"); } double[][][] allKnowns = new double[m][][]; double[][] allWeights = new double[m][]; double[] bases = new double[m]; double[] init = new double[6 * m + 1]; int l = 0; init[l++] = wavelength; for (int k = 0; k < m; k++) { List<EllipticalROI> ellipses = lEllipses.get(k); DetectorProperties dp = lDP.get(k); double dist = dp.getBeamCentreDistance(); int n = ellipses.size(); double[][] known = new double[n][FitFunctionBase.nC]; allKnowns[k] = known; double[] weight = new double[n]; allWeights[k] = weight; double base = -calcBaseRollAngle(ellipses); bases[k] = base; logger.debug("Mean roll angle: {}", Math.toDegrees(base)); for (int i = 0; i < n; i++) { EllipticalROI e = ellipses.get(i); weight[i] = e.getSemiAxis(0) / dist; int j = 0; double a = base - e.getAngle(); for (double off : FitFunctionBase.angles) { double[] pt = e.getPoint(a + off); known[i][j++] = pt[0]; known[i][j++] = pt[1]; } } Vector3d o = dp.getOrigin(); init[l++] = o.getX(); init[l++] = o.getY(); init[l++] = o.getZ(); double[] a = dp.getNormalAnglesInDegrees(); init[l++] = a[0]; init[l++] = a[1]; init[l++] = a[2]; } DetectorFitFunction f = new QSpacesFitFunction(allKnowns, allWeights, lDP.get(0).getVPxSize()); f.setInitial(init); f.setBaseRollAngles(bases); return f; } interface DetectorFitFunction extends FittingUtils.FitFunction { public double[] getNormalAngles(); /** * @return wavelength (in Angstrom) */ public double getWavelength(); /** * @return distance (in mm) */ public double getDistance(); public double[] getDistances(); /** * @param spacings (in Angstrom) */ public void setSpacings(List<Double> spacings); public void setSigma(double[] sigma); /** * @param baseAngle (in radians) */ public void setBaseRollAngle(double baseAngle); /** * @param baseAngles (in radians) */ public void setBaseRollAngles(double[] baseAngles); /** * * @return list of detector properties fitted */ public List<DetectorProperties> getDetectorProperties(); } static abstract class FitFunctionBase implements DetectorFitFunction { private double[] initial; protected double[] spacing; // in Angstroms protected int nR; // number of rings to fit // protected int nV; // number of calculated internally values protected double[] angle; protected double[] base; protected double[] parameters; protected SimpleBounds bounds; protected double[] sigma; protected int n; // number of parameters protected double maxWlen; // maximum wavelength allowed, in Angstroms protected static final double WAVE_MIN = 5e-2; // 0.05A protected static final double WAVE_MAX = 10; // 10.0A protected static final double DIST_MIN = 10; // (in mm) protected static final double DIST_MAX = 2e5; // (in mm) protected static final double SIGMA_WAVE = 2e-3; // 0.002A protected static final double SIGMA_POSN = 3; // (in mm) protected static final double SIGMA_ANG = 6; // (in degrees) // points are selected from given offset angles to fit to on each ring // protected static final double[] angles = {0, RIGHT_ANGLE, Math.PI, -RIGHT_ANGLE}; // offset angles protected static final double[] angles = { 0, Math.PI / 3, 2 * Math.PI / 3, Math.PI, -2 * Math.PI / 3, -Math.PI / 3 }; // offset angles // protected static final double[] angles = {0, Math.PI/6, Math.PI/3, Math.PI/2, 2*Math.PI/3, 5*Math.PI/6, Math.PI, -5*Math.PI/6, -2*Math.PI/3, -Math.PI/2, -Math.PI/3, -Math.PI/6}; // offset angles protected static final int nC = 2 * angles.length; // number of coordinate values per ring @Override public int getN() { return n; } @Override public void setParameters(double[] arg) { parameters = arg; } @Override public double[] getParameters() { return parameters; } @Override public double[] getSigma() { return sigma; } @Override public void setSigma(double[] sigma) { this.sigma = sigma; } @Override public SimpleBounds getBounds() { return bounds; } @Override public void setInitial(double... init) { initial = init; } @Override public double[] getInitial() { return initial; } @Override public void setBaseRollAngle(double baseAngle) { base[0] = baseAngle; } @Override public double[] getNormalAngles() { return null; } @Override public void setSpacings(List<Double> spacings) { double min = Double.POSITIVE_INFINITY; for (int i = 0; i < nR; i++) { spacing[i] = spacings.get(i); if (spacing[i] < min) { min = spacing[i]; } } // wavelength must be smaller than twice the smallest spacing /* * l/(2d_i) = sin(2t) < 1 * l < 2 d_i * => l < 2 d_min = l_max for any solution of Bragg's law */ maxWlen = 2 * min; } @Override public void setBaseRollAngles(double[] baseAngles) { base = baseAngles; } } /** * LS function uses 7 parameters: wavelength (Angstrom), detector origin (mm), orientation angles (degrees) */ static class QSpaceFitFunction7 extends FitFunctionBase { protected DetectorProperties dp; private double[][] target; private double[] weight; public QSpaceFitFunction7(double[][] known, double[] weight, double pix) { n = 7; nR = known.length; target = known; spacing = new double[nR]; bounds = new SimpleBounds( new double[] { WAVE_MIN, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, DIST_MIN, -90, -90, -180 }, new double[] { WAVE_MAX, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, DIST_MAX, 90, 90, 180 }); sigma = new double[] { SIGMA_WAVE, SIGMA_POSN, SIGMA_POSN, SIGMA_POSN, SIGMA_ANG, SIGMA_ANG, SIGMA_ANG }; dp = new DetectorProperties(); dp.setBeamVector(new Vector3d(0, 0, 1)); dp.setHPxSize(pix); dp.setVPxSize(pix); dp.setOrigin(new Vector3d(0, 0, 200)); base = new double[1]; if (weight.length != nR) { throw new IllegalArgumentException("Weight array should have length to match number of rings"); } this.weight = weight; double w = 0; // normalise weights for (int i = 0; i < nR; i++) { w += weight[i] * nC; } for (int i = 0; i < nR; i++) { weight[i] /= w; } } @Override public double getWavelength() { return parameters[0]; } protected void setDetector(int off, double... arg) { dp.setNormalAnglesInDegrees(arg[off + 3], arg[off + 4], arg[off + 5]); dp.setOrigin(new Vector3d(arg[off], arg[off + 1], arg[off + 2])); } @Override public double getDistance() { setDetector(1, parameters); return dp.getDetectorDistance(); } @Override public double[] getDistances() { return new double[] { getDistance() }; } @Override public double[] getNormalAngles() { setDetector(1, parameters); return dp.getNormalAnglesInDegrees(); } @Override public List<DetectorProperties> getDetectorProperties() { setDetector(1, parameters); List<DetectorProperties> l = new ArrayList<DetectorProperties>(); l.add(dp); return l; } double calcValue(double wlen, double... arg) { if (wlen > maxWlen) { return Double.POSITIVE_INFINITY; } setDetector(0, arg); double s = 0; boolean any = false; for (int i = 0; i < nR; i++) { IROI r; try { r = DSpacing.conicFromDSpacing(dp, wlen, spacing[i]); } catch (Exception e) { return Double.POSITIVE_INFINITY; } if (r instanceof EllipticalROI) { any = true; EllipticalROI ell = (EllipticalROI) r; double a = base[0] - ell.getAngle(); double[] pt = target[i]; int j = 0; double t = 0; for (double off : angles) { double[] pv = ell.getPoint(a + off); double u = pv[0] - pt[j++]; t += u * u; u = pv[1] - pt[j++]; t += u * u; } s += t * weight[i]; } } // System.err.println(s + "; w=" + wlen + "; d=" + dp); // XXX return any ? s : Double.POSITIVE_INFINITY; } @Override public double value(double[] arg) { return calcValue(arg[0], arg[1], arg[2], arg[3], arg[4], arg[5], arg[6]); } } /** * LS function uses 6 parameters: detector origin (mm), orientation angles (degrees) */ static class QSpaceFitFixedWFunction7 extends QSpaceFitFunction7 { protected double lambda; public QSpaceFitFixedWFunction7(double[][] known, double[] weight, double pix, double wavelength) { super(known, weight, pix); n = 6; int bl = bounds.getLower().length; bounds = new SimpleBounds(Arrays.copyOfRange(bounds.getLower(), 1, bl), Arrays.copyOfRange(bounds.getUpper(), 1, bl)); sigma = Arrays.copyOfRange(sigma, 1, sigma.length); lambda = wavelength; } @Override public double getWavelength() { return lambda; } @Override public double getDistance() { setDetector(0, parameters); return dp.getDetectorDistance(); } @Override public double[] getNormalAngles() { setDetector(0, parameters); return dp.getNormalAnglesInDegrees(); } @Override public List<DetectorProperties> getDetectorProperties() { setDetector(0, parameters); List<DetectorProperties> l = new ArrayList<DetectorProperties>(); l.add(dp); return l; } @Override public double value(double[] arg) { return calcValue(lambda, arg); } } /** * LS function uses 4 parameters: wavelength (Angstrom), detector origin (mm) */ static class QSpaceFitFunction4 extends QSpaceFitFunction7 { public QSpaceFitFunction4(double[][] known, double[] weight, double pix) { super(known, weight, pix); n = 4; bounds = new SimpleBounds( new double[] { WAVE_MIN, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, DIST_MIN }, new double[] { WAVE_MAX, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, DIST_MAX }); sigma = new double[] { SIGMA_WAVE, SIGMA_POSN, SIGMA_POSN, SIGMA_POSN }; } @Override protected void setDetector(int off, double... arg) { dp.setNormalAnglesInDegrees(0, 0, Math.toDegrees(base[0])); // need to correct for roll dp.setOrigin(new Vector3d(arg[off], arg[off + 1], arg[off + 2])); } @Override public double[] getNormalAngles() { return new double[3]; } @Override public double value(double[] arg) { return calcValue(arg[0], arg[1], arg[2], arg[3]); } } /** * LS function uses 3 parameters: detector origin (mm) */ static class QSpaceFitFixedWFunction4 extends QSpaceFitFunction4 { protected double lambda; public QSpaceFitFixedWFunction4(double[][] known, double[] weight, double pix, double wavelength) { super(known, weight, pix); n = 3; int bl = bounds.getLower().length; bounds = new SimpleBounds(Arrays.copyOfRange(bounds.getLower(), 1, bl), Arrays.copyOfRange(bounds.getUpper(), 1, bl)); sigma = Arrays.copyOfRange(sigma, 1, sigma.length); lambda = wavelength; } @Override public double getWavelength() { return lambda; } @Override public double getDistance() { setDetector(0, parameters); return dp.getDetectorDistance(); } @Override public List<DetectorProperties> getDetectorProperties() { setDetector(0, parameters); List<DetectorProperties> l = new ArrayList<DetectorProperties>(); l.add(dp); return l; } @Override public double value(double[] arg) { return calcValue(lambda, arg); } } /** * LS function uses 6I+1 parameters: wavelength (Angstrom), per image, detector origin (mm), orientation angles (degrees) */ static class QSpacesFitFunction extends FitFunctionBase { protected List<DetectorProperties> dps; private double[][][] target; private double[][] weight; public QSpacesFitFunction(double[][][] known, double[][] weight, double pix) { target = known; int m = known.length; n = 6 * m + 1; nR = 0; double[] lb = new double[n]; double[] ub = new double[n]; sigma = new double[n]; dps = new ArrayList<DetectorProperties>(); int j = 0; lb[j] = WAVE_MIN; ub[j] = WAVE_MAX; sigma[j++] = SIGMA_WAVE; DetectorProperties dp; for (int k = 0; k < m; k++) { nR += known[k].length; lb[j] = Double.NEGATIVE_INFINITY; ub[j] = Double.POSITIVE_INFINITY; sigma[j++] = SIGMA_POSN; lb[j] = Double.NEGATIVE_INFINITY; ub[j] = Double.POSITIVE_INFINITY; sigma[j++] = SIGMA_POSN; lb[j] = DIST_MIN; ub[j] = DIST_MAX; sigma[j++] = SIGMA_POSN; lb[j] = -90; ub[j] = 90; sigma[j++] = SIGMA_ANG; lb[j] = -90; ub[j] = 90; sigma[j++] = SIGMA_ANG; lb[j] = -180; ub[j] = 180; sigma[j++] = SIGMA_ANG; dp = new DetectorProperties(); dp.setBeamVector(new Vector3d(0, 0, 1)); dp.setHPxSize(pix); dp.setVPxSize(pix); dp.setOrigin(new Vector3d(0, 0, 200)); dps.add(dp); } if (weight.length != m) { throw new IllegalArgumentException("Weight array should have length to match number of images"); } this.weight = weight; double w = 0; // normalise weights for (int k = 0; k < m; k++) { double[] wgt = weight[k]; int n = wgt.length; if (n != known[k].length) { throw new IllegalArgumentException("Weight array should have length to match number of rings"); } for (int i = 0; i < n; i++) { w += wgt[i] * nC; } } for (int k = 0; k < m; k++) { double[] wgt = weight[k]; int n = wgt.length; for (int i = 0; i < n; i++) { wgt[i] /= w; } } spacing = new double[nR]; bounds = new SimpleBounds(lb, ub); } @Override public double getWavelength() { return parameters[0]; } protected void setDetector(int off, double... arg) { int i = off; for (DetectorProperties dp : dps) { int j = i; i += 3; dp.setNormalAnglesInDegrees(arg[i++], arg[i++], arg[i++]); dp.setOrigin(new Vector3d(arg[j++], arg[j++], arg[j++])); } } @Override public double getDistance() { setDetector(1, parameters); return dps.get(0).getDetectorDistance(); } @Override public double[] getDistances() { setDetector(1, parameters); int m = dps.size(); double[] ds = new double[m]; for (int k = 0; k < m; k++) { ds[k] = dps.get(k).getDetectorDistance(); } return ds; } @Override public double[] getNormalAngles() { setDetector(1, parameters); int m = dps.size(); double[] as = new double[3 * m]; int i = 0; for (int k = 0; k < m; k++) { double[] a = dps.get(k).getNormalAnglesInDegrees(); as[i++] = a[0]; as[i++] = a[1]; as[i++] = a[2]; } return as; } @Override public List<DetectorProperties> getDetectorProperties() { setDetector(1, parameters); return dps; } double calcValue(double wlen, double... arg) { if (wlen > maxWlen) { return Double.POSITIVE_INFINITY; } setDetector(0, arg); double s = 0; boolean any = false; int m = target.length; int i = 0; for (int k = 0; k < m; k++) { double[][] tgt = target[k]; double[] wgt = weight[k]; int nr = tgt.length; DetectorProperties dp = dps.get(k); for (int l = 0; l < nr; l++) { IROI r; try { r = DSpacing.conicFromDSpacing(dp, wlen, spacing[i + l]); } catch (Exception e) { return Double.POSITIVE_INFINITY; } if (r instanceof EllipticalROI) { any = true; EllipticalROI ell = (EllipticalROI) r; double a = base[k] - ell.getAngle(); double[] pt = tgt[l]; int j = 0; double t = 0; for (double off : angles) { double[] pv = ell.getPoint(a + off); double u = pv[0] - pt[j++]; t += u * u; u = pv[1] - pt[j++]; t += u * u; } s += t * wgt[l]; } } // System.err.println(s + "; w=" + wlen + "; d=" + dp); // XXX i += nr; } return any ? s : Double.POSITIVE_INFINITY; } @Override public double value(double[] arg) { double[] a = Arrays.copyOfRange(arg, 1, arg.length); return calcValue(arg[0], a); } } }