List of usage examples for org.apache.commons.math3.linear RealVector setEntry
public abstract void setEntry(int index, double value) throws OutOfRangeException;
From source file:edu.stanford.cfuller.colocalization3d.Colocalization3DMain.java
private String formatPositionData(java.util.List<ImageObject> imageObjects, Correction c) { StringBuilder sb = new StringBuilder(); StringBuilder sbSingle = new StringBuilder(); String currentImageName = ""; File lastFile = null;// ww w . ja v a 2 s .co m for (int o = 0; o < imageObjects.size(); o++) { ImageObject currObject = imageObjects.get(o); String objectImageName = currObject.getImageID(); if (!currentImageName.equals(objectImageName)) { if (!(sbSingle.length() == 0)) { sb.append(sbSingle.toString()); sbSingle = new StringBuilder(); } sb.append(objectImageName + "\n"); currentImageName = objectImageName; } RealVector correction = null; try { correction = c.correctPosition( currObject.getPositionForChannel(c.getReferenceChannelIndex()).getEntry(0), currObject.getPositionForChannel(c.getReferenceChannelIndex()).getEntry(1)); } catch (UnableToCorrectException e) { continue; } sbSingle.append(currObject.getLabel()); sbSingle.append(" "); sbSingle.append(currObject.getPositionForChannel(c.getReferenceChannelIndex()).getEntry(0)); sbSingle.append(" "); sbSingle.append(currObject.getPositionForChannel(c.getReferenceChannelIndex()).getEntry(1)); sbSingle.append(" "); sbSingle.append(currObject.getPositionForChannel(c.getReferenceChannelIndex()).getEntry(2)); sbSingle.append(" "); sbSingle.append(currObject.getPositionForChannel(c.getCorrectionChannelIndex()).getEntry(0)); sbSingle.append(" "); sbSingle.append(currObject.getPositionForChannel(c.getCorrectionChannelIndex()).getEntry(1)); sbSingle.append(" "); sbSingle.append(currObject.getPositionForChannel(c.getCorrectionChannelIndex()).getEntry(2)); sbSingle.append(" "); RealVector corrPos = null; RealVector modCorrection = correction; if (this.parameters.hasKeyAndTrue("flip_channels_at_end")) { modCorrection = correction.mapMultiply(-1.0); } if (this.parameters.hasKeyAndTrue("inverted_z_axis")) { modCorrection = correction.mapMultiply(1.0); modCorrection.setEntry(2, -1.0 * modCorrection.getEntry(2)); } corrPos = currObject.getPositionForChannel(c.getCorrectionChannelIndex()).subtract(modCorrection); sbSingle.append(corrPos.getEntry(0)); sbSingle.append(" "); sbSingle.append(corrPos.getEntry(1)); sbSingle.append(" "); sbSingle.append(corrPos.getEntry(2)); sbSingle.append(" "); sbSingle.append( currObject.getFitParametersByChannel().get(c.getReferenceChannelIndex()).getAmplitude()); sbSingle.append(" "); sbSingle.append( currObject.getFitParametersByChannel().get(c.getCorrectionChannelIndex()).getAmplitude()); sbSingle.append("\n"); } sb.append(sbSingle.toString()); return sb.toString(); }
From source file:agent.perceptors.RestrictedVisionPerceptor.java
@Override protected synchronized void setField(int guess, int start, int end, char[] msg) { String[] str = new String(msg, start, (end - start - 1)).split(" "); RealVector v = null; switch (guess) { case 1:/*from w w w.jav a2s.c o m*/ return; case 2: return; case 3: v = ball; break; case 4: v = flag1L; break; case 5: v = flag1R; break; case 6: v = flag2L; break; case 7: v = flag2R; break; case 8: v = goalPost1L; break; case 9: v = goalPost1R; break; case 10: v = goalPost2L; break; case 11: v = goalPost2R; break; } if (v == null) { throw new UnsupportedOperationException("Invalid guess!"); } double r = Double.parseDouble(str[1]); double a = Double.parseDouble(str[2]); double b = Double.parseDouble(str[3]); v.setEntry(0, r * sin(a) * cos(b)); v.setEntry(1, r * sin(a) * sin(b)); v.setEntry(2, r * cos(a)); // for (int i = 0; i < 3; i++) { // v.setEntry(i, Double.parseDouble(str[i+1])); // } }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.BisquareLinearFit.java
/** * Calculates the leverages of data points for least squares fitting (assuming equal variances). * * @param indVarValues The values of the independent variable used for the fitting. * @return a RealVector containing a leverage value for each independent variable value. *//*from w ww. j a v a 2 s.co m*/ protected RealVector calculateLeverages(RealVector indVarValues) { RealMatrix indVarMatrix = null; if (this.noIntercept) { indVarMatrix = new Array2DRowRealMatrix(indVarValues.getDimension(), 1); } else { indVarMatrix = new Array2DRowRealMatrix(indVarValues.getDimension(), 2); } indVarMatrix.setColumnVector(0, indVarValues); if (!this.noIntercept) { indVarMatrix.setColumnVector(1, indVarValues.mapMultiply(0).mapAdd(1)); } RealVector leverages = new ArrayRealVector(indVarValues.getDimension()); QRDecomposition xQR = new QRDecomposition(indVarMatrix); RealMatrix xR = xQR.getR(); int smallerDim = xR.getRowDimension() < xR.getColumnDimension() ? xR.getRowDimension() : xR.getColumnDimension(); RealMatrix xRSq = xR.getSubMatrix(0, smallerDim - 1, 0, smallerDim - 1); QRDecomposition xRQR = new QRDecomposition(xRSq); RealMatrix xRInv = xRQR.getSolver().getInverse(); RealMatrix xxRInv = indVarMatrix.multiply(xRInv); for (int i = 0; i < indVarValues.getDimension(); i++) { double sum = 0; for (int j = 0; j < xxRInv.getColumnDimension(); j++) { sum += Math.pow(xxRInv.getEntry(i, j), 2); } leverages.setEntry(i, sum); } return leverages; }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.NelderMeadMinimizer.java
/** * Runs the minimization of the specified function starting from an initial simplex. * @param f The ObjectiveFunction to be minimized. * @param initialSimplex The initial simplex to use to start optimization, as might be returned from {@link #generateInitialSimplex} * @return The parameters at the function minimum in the same order as specified for each point on the simplex. *//* w w w . j a v a 2 s .c o m*/ public RealVector optimize(ObjectiveFunction f, RealMatrix initialSimplex) { RealMatrix currentSimplex = initialSimplex.copy(); double currTolVal = 1.0e6; RealVector values = new ArrayRealVector(initialSimplex.getRowDimension(), 0.0); RealVector centerOfMass = new ArrayRealVector(initialSimplex.getColumnDimension(), 0.0); boolean shouldEvaluate = false; long iterCounter = 0; while (Math.abs(currTolVal) > this.tol) { int maxIndex = 0; int minIndex = 0; double maxValue = -1.0 * Double.MAX_VALUE; double minValue = Double.MAX_VALUE; double secondMaxValue = -1.0 * Double.MAX_VALUE; centerOfMass.mapMultiplyToSelf(0.0); if (shouldEvaluate) { for (int i = 0; i < currentSimplex.getRowDimension(); i++) { RealVector currRow = currentSimplex.getRowVector(i); values.setEntry(i, f.evaluate(currRow)); } } for (int i = 0; i < currentSimplex.getRowDimension(); i++) { double currValue = values.getEntry(i); if (currValue < minValue) { minValue = currValue; minIndex = i; } if (currValue > maxValue) { secondMaxValue = maxValue; maxValue = currValue; maxIndex = i; } else if (currValue > secondMaxValue) { secondMaxValue = currValue; } } for (int i = 0; i < currentSimplex.getRowDimension(); i++) { if (i == maxIndex) continue; centerOfMass = centerOfMass.add(currentSimplex.getRowVector(i)); } centerOfMass.mapDivideToSelf(currentSimplex.getRowDimension() - 1); RealVector oldPoint = currentSimplex.getRowVector(maxIndex); RealVector newPoint = centerOfMass.subtract(oldPoint).mapMultiplyToSelf(a).add(centerOfMass); // newpoint = COM + a*(COM-oldpoint) double newValue = f.evaluate(newPoint); if (newValue < secondMaxValue) { // success if (newValue < minValue) { // best found so far //expansion RealVector expPoint = centerOfMass.subtract(oldPoint).mapMultiplyToSelf(g).add(centerOfMass); double expValue = f.evaluate(expPoint); if (expValue < newValue) { currentSimplex.setRowVector(maxIndex, expPoint); currTolVal = 2.0 * (expValue - maxValue) / (1.0e-20 + expValue + maxValue); values.setEntry(maxIndex, expValue); shouldEvaluate = false; continue; } } //reflection currentSimplex.setRowVector(maxIndex, newPoint); currTolVal = 2.0 * (newValue - maxValue) / (1.0e-20 + newValue + maxValue); values.setEntry(maxIndex, newValue); shouldEvaluate = false; continue; } //contraction RealVector conPoint = centerOfMass.subtract(oldPoint).mapMultiplyToSelf(r).add(centerOfMass); double conValue = f.evaluate(conPoint); if (conValue < maxValue) { currentSimplex.setRowVector(maxIndex, conPoint); currTolVal = 2.0 * (conValue - maxValue) / (1.0e-20 + conValue + maxValue); values.setEntry(maxIndex, conValue); shouldEvaluate = false; continue; } //reduction for (int i = 0; i < currentSimplex.getRowDimension(); i++) { if (i == minIndex) continue; currentSimplex.setRowVector(i, currentSimplex.getRowVector(i).subtract(currentSimplex.getRowVector(minIndex)) .mapMultiplyToSelf(s).add(currentSimplex.getRowVector(minIndex))); } double redValue = f.evaluate(currentSimplex.getRowVector(maxIndex)); currTolVal = 2.0 * (redValue - maxValue) / (1.0e-20 + redValue + maxValue); shouldEvaluate = true; if (iterCounter++ > 100000) { System.out.println("stalled? tol: " + currTolVal + " minValue: " + minValue); } } double minValue = Double.MAX_VALUE; RealVector minVector = null; for (int i = 0; i < currentSimplex.getRowDimension(); i++) { values.setEntry(i, f.evaluate(currentSimplex.getRowVector(i))); if (values.getEntry(i) < minValue) { minValue = values.getEntry(i); minVector = currentSimplex.getRowVector(i); } } return minVector; }
From source file:edu.stanford.cfuller.imageanalysistools.filter.ConvexHullByLabelFilter.java
/** * Applies the convex hull filter to the supplied mask. * @param im The Image to process-- a mask whose regions will be replaced by their filled convex hulls. *//* w w w . j a va2s . com*/ @Override public void apply(WritableImage im) { RelabelFilter RLF = new RelabelFilter(); RLF.apply(im); Histogram h = new Histogram(im); java.util.Hashtable<Integer, java.util.Vector<Integer>> xLists = new java.util.Hashtable<Integer, java.util.Vector<Integer>>(); java.util.Hashtable<Integer, java.util.Vector<Integer>> yLists = new java.util.Hashtable<Integer, java.util.Vector<Integer>>(); java.util.Vector<Integer> minValues = new java.util.Vector<Integer>(h.getMaxValue() + 1); java.util.Vector<Integer> minIndices = new java.util.Vector<Integer>(h.getMaxValue() + 1); for (int i = 0; i < h.getMaxValue() + 1; i++) { minValues.add(im.getDimensionSizes().get(ImageCoordinate.X)); minIndices.add(0); } for (ImageCoordinate i : im) { int value = (int) im.getValue(i); if (value == 0) continue; if (!xLists.containsKey(value)) { xLists.put(value, new java.util.Vector<Integer>()); yLists.put(value, new java.util.Vector<Integer>()); } xLists.get(value).add(i.get(ImageCoordinate.X)); yLists.get(value).add(i.get(ImageCoordinate.Y)); if (i.get(ImageCoordinate.X) < minValues.get(value)) { minValues.set(value, i.get(ImageCoordinate.X)); minIndices.set(value, xLists.get(value).size() - 1); } } java.util.Vector<Integer> hullPointsX = new java.util.Vector<Integer>(); java.util.Vector<Integer> hullPointsY = new java.util.Vector<Integer>(); ImageCoordinate ic = ImageCoordinate.createCoordXYZCT(0, 0, 0, 0, 0); for (int k = 1; k < h.getMaxValue() + 1; k++) { hullPointsX.clear(); hullPointsY.clear(); java.util.Vector<Integer> xList = xLists.get(k); java.util.Vector<Integer> yList = yLists.get(k); int minIndex = (int) minIndices.get(k); //start at the leftmost point int currentIndex = minIndex; int currentX = xList.get(currentIndex); int currentY = yList.get(currentIndex); hullPointsX.add(currentX); hullPointsY.add(currentY); org.apache.commons.math3.linear.RealVector angles = new org.apache.commons.math3.linear.ArrayRealVector( xList.size()); Vector3D currentVector = new Vector3D(0, -1, 0); java.util.HashSet<Integer> visited = new java.util.HashSet<Integer>(); do { visited.add(currentIndex); int maxIndex = 0; double maxAngle = -2 * Math.PI; double dist = Double.MAX_VALUE; for (int i = 0; i < xList.size(); i++) { if (i == currentIndex) continue; Vector3D next = new Vector3D(xList.get(i) - xList.get(currentIndex), yList.get(i) - yList.get(currentIndex), 0); double angle = Vector3D.angle(currentVector, next); angles.setEntry(i, angle); if (angle > maxAngle) { maxAngle = angle; maxIndex = i; dist = next.getNorm(); } else if (angle == maxAngle) { double tempDist = next.getNorm(); if (tempDist < dist) { dist = tempDist; maxAngle = angle; maxIndex = i; } } } currentX = xList.get(maxIndex); currentY = yList.get(maxIndex); currentVector = new Vector3D(xList.get(currentIndex) - currentX, yList.get(currentIndex) - currentY, 0); hullPointsX.add(currentX); hullPointsY.add(currentY); currentIndex = maxIndex; } while (!visited.contains(currentIndex)); //hull vertices have now been determined .. need to fill in the lines //between them so I can apply a fill filter //approach: x1, y1 to x0, y0: //start at min x, min y, go to max x, max y // if x_i, y_i = x0, y0 + slope to within 0.5 * sqrt(2), then add to hull double eps = Math.sqrt(2); for (int i = 0; i < hullPointsX.size() - 1; i++) { int x0 = hullPointsX.get(i); int y0 = hullPointsY.get(i); int x1 = hullPointsX.get(i + 1); int y1 = hullPointsY.get(i + 1); int xmin = (x0 < x1) ? x0 : x1; int ymin = (y0 < y1) ? y0 : y1; int xmax = (x0 > x1) ? x0 : x1; int ymax = (y0 > y1) ? y0 : y1; x1 -= x0; y1 -= y0; double denom = (x1 * x1 + y1 * y1); for (int x = xmin; x <= xmax; x++) { for (int y = ymin; y <= ymax; y++) { int rel_x = x - x0; int rel_y = y - y0; double projLength = (x1 * rel_x + y1 * rel_y) / denom; double projPoint_x = x1 * projLength; double projPoint_y = y1 * projLength; if (Math.hypot(rel_x - projPoint_x, rel_y - projPoint_y) < eps) { ic.set(ImageCoordinate.X, x); ic.set(ImageCoordinate.Y, y); im.setValue(ic, k); } } } } } ic.recycle(); FillFilter ff = new FillFilter(); ff.apply(im); }
From source file:edu.stanford.cfuller.colocalization3d.correction.PositionCorrector.java
/** * Creates a correction from a set of objects whose positions should be the same in each channel. * * @param imageObjects A Vector containing all the ImageObjects to be used for the correction * or in the order it appears in a multiwavelength image file. * @return A Correction object that can be used to correct the positions of other objects based upon the standards provided. *///from w ww . j a v a 2 s.c om public Correction getCorrection(java.util.List<ImageObject> imageObjects) { int referenceChannel = this.parameters.getIntValueForKey(REF_CH_PARAM); int channelToCorrect = this.parameters.getIntValueForKey(CORR_CH_PARAM); if (!this.parameters.hasKeyAndTrue(DET_CORR_PARAM)) { try { return Correction.readFromDisk(FileUtils.getCorrectionFilename(this.parameters)); } catch (java.io.IOException e) { java.util.logging.Logger .getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME) .severe("Exception encountered while reading correction from disk: "); e.printStackTrace(); } catch (ClassNotFoundException e) { java.util.logging.Logger .getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME) .severe("Exception encountered while reading correction from disk: "); e.printStackTrace(); } return null; } int numberOfPointsToFit = this.parameters.getIntValueForKey(NUM_POINT_PARAM); RealMatrix correctionX = new Array2DRowRealMatrix(imageObjects.size(), numberOfCorrectionParameters); RealMatrix correctionY = new Array2DRowRealMatrix(imageObjects.size(), numberOfCorrectionParameters); RealMatrix correctionZ = new Array2DRowRealMatrix(imageObjects.size(), numberOfCorrectionParameters); RealVector distanceCutoffs = new ArrayRealVector(imageObjects.size(), 0.0); RealVector ones = new ArrayRealVector(numberOfPointsToFit, 1.0); RealVector distancesToObjects = new ArrayRealVector(imageObjects.size(), 0.0); RealMatrix allCorrectionParametersMatrix = new Array2DRowRealMatrix(numberOfPointsToFit, numberOfCorrectionParameters); for (int i = 0; i < imageObjects.size(); i++) { RealVector ithPos = imageObjects.get(i).getPositionForChannel(referenceChannel); for (int j = 0; j < imageObjects.size(); j++) { double d = imageObjects.get(j).getPositionForChannel(referenceChannel).subtract(ithPos).getNorm(); distancesToObjects.setEntry(j, d); } //the sorting becomes a bottleneck once the number of points gets large //reverse comparator so we can use the priority queue and get the max element at the head Comparator<Double> cdReverse = new Comparator<Double>() { public int compare(Double o1, Double o2) { if (o1.equals(o2)) return 0; if (o1 > o2) return -1; return 1; } }; PriorityQueue<Double> pq = new PriorityQueue<Double>(numberOfPointsToFit + 2, cdReverse); double maxElement = Double.MAX_VALUE; for (int p = 0; p < numberOfPointsToFit + 1; p++) { pq.add(distancesToObjects.getEntry(p)); } maxElement = pq.peek(); for (int p = numberOfPointsToFit + 1; p < distancesToObjects.getDimension(); p++) { double value = distancesToObjects.getEntry(p); if (value < maxElement) { pq.poll(); pq.add(value); maxElement = pq.peek(); } } double firstExclude = pq.poll(); double lastDist = pq.poll(); double distanceCutoff = (lastDist + firstExclude) / 2.0; distanceCutoffs.setEntry(i, distanceCutoff); RealVector xPositionsToFit = new ArrayRealVector(numberOfPointsToFit, 0.0); RealVector yPositionsToFit = new ArrayRealVector(numberOfPointsToFit, 0.0); RealVector zPositionsToFit = new ArrayRealVector(numberOfPointsToFit, 0.0); RealMatrix differencesToFit = new Array2DRowRealMatrix(numberOfPointsToFit, imageObjects.get(0).getPositionForChannel(referenceChannel).getDimension()); int toFitCounter = 0; for (int j = 0; j < imageObjects.size(); j++) { if (distancesToObjects.getEntry(j) < distanceCutoff) { xPositionsToFit.setEntry(toFitCounter, imageObjects.get(j).getPositionForChannel(referenceChannel).getEntry(0)); yPositionsToFit.setEntry(toFitCounter, imageObjects.get(j).getPositionForChannel(referenceChannel).getEntry(1)); zPositionsToFit.setEntry(toFitCounter, imageObjects.get(j).getPositionForChannel(referenceChannel).getEntry(2)); differencesToFit.setRowVector(toFitCounter, imageObjects.get(j) .getVectorDifferenceBetweenChannels(referenceChannel, channelToCorrect)); toFitCounter++; } } RealVector x = xPositionsToFit.mapSubtractToSelf(ithPos.getEntry(0)); RealVector y = yPositionsToFit.mapSubtractToSelf(ithPos.getEntry(1)); allCorrectionParametersMatrix.setColumnVector(0, ones); allCorrectionParametersMatrix.setColumnVector(1, x); allCorrectionParametersMatrix.setColumnVector(2, y); allCorrectionParametersMatrix.setColumnVector(3, x.map(new Power(2))); allCorrectionParametersMatrix.setColumnVector(4, y.map(new Power(2))); allCorrectionParametersMatrix.setColumnVector(5, x.ebeMultiply(y)); DecompositionSolver solver = (new QRDecomposition(allCorrectionParametersMatrix)).getSolver(); RealVector cX = solver.solve(differencesToFit.getColumnVector(0)); RealVector cY = solver.solve(differencesToFit.getColumnVector(1)); RealVector cZ = solver.solve(differencesToFit.getColumnVector(2)); correctionX.setRowVector(i, cX); correctionY.setRowVector(i, cY); correctionZ.setRowVector(i, cZ); } Correction c = new Correction(correctionX, correctionY, correctionZ, distanceCutoffs, imageObjects, referenceChannel, channelToCorrect); return c; }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.DifferentialEvolutionMinimizer.java
/** * Performs a minimization of a function starting with a given population. * /*from w w w . ja va 2s . co m*/ * @param population The population of parameters to start from, one population entry per row, one parameter per column. * @param f The function to be minimized. * @param parameterLowerBounds The lower bounds of each parameter. This must have the same size as the column dimension of the population. Generated parameter values less than these values will be discarded. * @param parameterUpperBounds The upper bounds of each paraemter. This must have the same size as the column dimension of the population. Generated parameter values greater than these values will be discarded. * @param populationSize The size of the population of parameters sets. This must be equal to the row dimension of the population. * @param scaleFactor Factor controlling the scale of crossed over entries during crossover events. * @param maxIterations The maximum number of iterations to allow before returning a result. * @param crossoverFrequency The frequency of crossover from 0 to 1. At any given parameter, this is the probability of initiating a crossover as well as the probability of ending one after it has started. * @param tol Relative function value tolerance controlling termination; if the maximal and minimal population values differ by less than this factor times the maximal value, optimization will terminate. * @return The parameter values at the minimal function value found. */ public RealVector minimizeWithPopulation(RealMatrix population, ObjectiveFunction f, RealVector parameterLowerBounds, RealVector parameterUpperBounds, int populationSize, double scaleFactor, int maxIterations, double crossoverFrequency, double tol) { int numberOfParameters = parameterLowerBounds.getDimension(); double currMinValue = Double.MAX_VALUE; double currMaxValue = -1.0 * Double.MAX_VALUE; int iterationCounter = maxIterations; double mutationProb = 0.01; // int totalIterations =0; RealVector values = new ArrayRealVector(populationSize); boolean[] valuesChanged = new boolean[populationSize]; java.util.Arrays.fill(valuesChanged, true); computeValues(f, population, values, valuesChanged); RealVector newVec = new ArrayRealVector(numberOfParameters); RealMatrix newPopulation = new Array2DRowRealMatrix(populationSize, numberOfParameters); while (iterationCounter > 0) { for (int i = 0; i < populationSize; i++) { int i1 = RandomGenerator.getGenerator().randInt(populationSize); int i2 = RandomGenerator.getGenerator().randInt(populationSize); int i3 = RandomGenerator.getGenerator().randInt(populationSize); newVec.mapMultiplyToSelf(0.0); boolean inBounds = true; boolean isCrossingOver = false; for (int j = 0; j < numberOfParameters; j++) { if ((RandomGenerator.rand() < crossoverFrequency) ^ isCrossingOver) { if (!isCrossingOver) { isCrossingOver = true; } newVec.setEntry(j, scaleFactor * (population.getEntry(i2, j) - population.getEntry(i1, j)) + population.getEntry(i3, j)); } else { if (isCrossingOver) { isCrossingOver = false; } newVec.setEntry(j, population.getEntry(i, j)); } //random 10% range +/- mutation if ((RandomGenerator.rand() < mutationProb)) { double magnitude = 0.2 * (parameterUpperBounds.getEntry(j) - parameterLowerBounds.getEntry(j)); newVec.setEntry(j, newVec.getEntry(j) + (RandomGenerator.rand() - 0.5) * magnitude); } if (newVec.getEntry(j) < parameterLowerBounds.getEntry(j) || newVec.getEntry(j) > parameterUpperBounds.getEntry(j)) { inBounds = false; } } double functionValue = Double.MAX_VALUE; if (inBounds) functionValue = f.evaluate(newVec); //if (inBounds) System.out.printf("in bounds candidate value: %1.2f old value: %1.2f \n", functionValue, values.getEntry(i)); if (functionValue < values.getEntry(i)) { newPopulation.setRowVector(i, newVec); valuesChanged[i] = true; values.setEntry(i, functionValue); } else { newPopulation.setRowVector(i, population.getRowVector(i)); valuesChanged[i] = false; } } population = newPopulation; double tempMinValue = Double.MAX_VALUE; double tempMaxValue = -1.0 * Double.MAX_VALUE; // double averageValue = 0; for (int i = 0; i < values.getDimension(); i++) { double value = values.getEntry(i); // averageValue += value; if (value < tempMinValue) { tempMinValue = value; } if (value > tempMaxValue) { tempMaxValue = value; } } // averageValue /= values.getDimension(); currMinValue = tempMinValue; currMaxValue = tempMaxValue; //LoggingUtilities.getLogger().info("Iteration counter: " + Integer.toString(totalIterations) + " best score: " + currMinValue + " worst score: " + currMaxValue + " average score: " + averageValue); if (Math.abs(currMaxValue - currMinValue) < Math.abs(tol * currMaxValue) + Math.abs(tol * currMinValue)) { iterationCounter--; } else { iterationCounter = 1; } // totalIterations++; } for (int i = 0; i < populationSize; i++) { valuesChanged[i] = true; } computeValues(f, population, values, valuesChanged); double tempMinValue = Double.MAX_VALUE; int tempMinIndex = 0; for (int i = 0; i < populationSize; i++) { if (values.getEntry(i) < tempMinValue) { tempMinValue = values.getEntry(i); tempMinIndex = i; } } RealVector output = new ArrayRealVector(numberOfParameters); for (int i = 0; i < numberOfParameters; i++) { output.setEntry(i, population.getEntry(tempMinIndex, i)); } return output; }
From source file:edu.stanford.cfuller.imageanalysistools.filter.LocalMaximumSeparabilityThresholdingFilter.java
/** * Applies the filter to an Image, optionally turning on adaptive determination of the increment size used to find the threshold, and specifying a size for the threshold determination increment. * Turning on adaptive determination of the increment will generally make the threshold slightly less optimal, but can sometimes speed up the filtering, especially * for images with a large dynamic range. * <p>/* ww w.j av a 2 s . co m*/ * The increment size specified (in greylevels) will be used to determine the threshold only if adaptive determination is off; otherwise this parameter will be ignored. * @param im The Image to be thresholded; this will be overwritten by the thresholded Image. * @param adaptiveincrement true to turn on adaptive determination of the threshold increment; false to turn it off and use the default value * @param increment the increment size (in greylevels) to use for determining the threshold; must be positive. */ public void apply_ext(WritableImage im, boolean adaptiveincrement, int increment) { Histogram h = new Histogram(im); int thresholdValue = 0; final int numSteps = 1000; double best_eta = 0; int best_index = Integer.MAX_VALUE; int nonzerocounts = h.getTotalCounts() - h.getCounts(0); double meannonzero = h.getMeanNonzero(); ArrayRealVector omega_v = new ArrayRealVector(h.getMaxValue()); ArrayRealVector mu_v = new ArrayRealVector(h.getMaxValue()); ArrayRealVector eta_v = new ArrayRealVector(h.getMaxValue()); int c = 0; if (adaptiveincrement) { increment = (int) ((h.getMaxValue() - h.getMinValue() + 1) * 1.0 / numSteps); if (increment < 1) increment = 1; } for (int k = h.getMinValue(); k < h.getMaxValue() + 1; k += increment) { if (k == 0) continue; omega_v.setEntry(c, h.getCumulativeCounts(k) * 1.0 / nonzerocounts); if (c == 0) { mu_v.setEntry(c, k * omega_v.getEntry(c)); } else { mu_v.setEntry(c, mu_v.getEntry(c - 1) + k * h.getCounts(k) * 1.0 / nonzerocounts); for (int i = k - increment + 1; i < k; i++) { mu_v.setEntry(c, mu_v.getEntry(c) + h.getCounts(i) * i * 1.0 / nonzerocounts); } } double omega = omega_v.getEntry(c); double mu = mu_v.getEntry(c); if (omega > 1e-8 && 1 - omega > 1e-8) { double eta = omega * (1 - omega) * Math.pow((meannonzero - mu) / (1 - omega) - mu / omega, 2); eta_v.setEntry(c, eta); if (eta >= best_eta) { best_eta = eta; best_index = k; } } else { eta_v.setEntry(c, 0); } c++; } int orig_method_best_index = best_index; c = 1; ArrayList<Integer> maxima = new ArrayList<Integer>(); Map<Integer, Integer> k_by_c = new HashMap<Integer, Integer>(); Map<Integer, Integer> c_by_k = new HashMap<Integer, Integer>(); for (int k = h.getMinValue() + 1; k < h.getMaxValue(); k += increment) { //detect if this is a local maximum k_by_c.put(c, k); c_by_k.put(k, c); int lastEntryNotEqual = c - 1; int nextEntryNotEqual = c + 1; while (lastEntryNotEqual > 0 && eta_v.getEntry(lastEntryNotEqual) == eta_v.getEntry(c)) { --lastEntryNotEqual; } while (nextEntryNotEqual < (eta_v.getDimension() - 1) && eta_v.getEntry(nextEntryNotEqual) == eta_v.getEntry(c)) { ++nextEntryNotEqual; } if (eta_v.getEntry(c) > eta_v.getEntry(lastEntryNotEqual) && eta_v.getEntry(c) > eta_v.getEntry(nextEntryNotEqual)) { maxima.add(k); } c++; } //now that we have maxima, try doing a gaussian fit to find the positions. If there's only one, we need to guess at a second RealVector parameters = new ArrayRealVector(6, 0.0); int position0 = 0; int position1 = h.getMaxValue(); if (maxima.size() > 1) { double best_max = 0; double second_best_max = 0; int best_pos = 0; int second_best_pos = 0; for (Integer k : maxima) { int ck = c_by_k.get(k); double eta_k = eta_v.getEntry(ck); if (eta_k > best_max) { second_best_max = best_max; second_best_pos = best_pos; best_max = eta_k; best_pos = ck; } else if (eta_k > second_best_max) { second_best_max = eta_k; second_best_pos = ck; } } position0 = best_pos; position1 = second_best_pos; } else { position0 = c_by_k.get(maxima.get(0)); position1 = (eta_v.getDimension() - position0) / 2 + position0; } //make sure that position 1 is larger than position 0 if (position1 < position0) { int temp = position0; position0 = position1; position1 = temp; } double s = (position1 - position0) / 4.0; parameters.setEntry(0, eta_v.getEntry(position0));//*Math.sqrt(2*Math.PI)*s); parameters.setEntry(1, position0); parameters.setEntry(2, s); parameters.setEntry(3, eta_v.getEntry(position1));//*Math.sqrt(2*Math.PI)*s); parameters.setEntry(4, position1); parameters.setEntry(5, s); DoubleGaussianObjectiveFunction dgof = new DoubleGaussianObjectiveFunction(); dgof.setEta(eta_v); NelderMeadMinimizer nmm = new NelderMeadMinimizer(); RealVector result = nmm.optimize(dgof, parameters); best_index = (int) result.getEntry(4); if (k_by_c.containsKey(best_index)) { best_index = k_by_c.get(best_index); } else { //fall back to the normal global maximum if the fitting seems to have found an invalid value. best_index = orig_method_best_index; } thresholdValue = best_index; if (thresholdValue == Integer.MAX_VALUE) { thresholdValue = 0; } for (ImageCoordinate coord : im) { if (im.getValue(coord) < thresholdValue) im.setValue(coord, 0); } }
From source file:edu.stanford.cfuller.colocalization3d.correction.PositionCorrector.java
/** * Creates an in situ cellular aberration correction. This will take a specified group of datasets from the parameters file and perform a robust linar fit between * the aberrations in two pairs of channels (also specified in the parameters file). * * @return a RealVector with three elements: the slope of the linear fit in the x, y, and z directions. *///www. j a v a 2s. com public RealVector determineInSituAberrationCorrection() { int referenceChannel = this.parameters.getIntValueForKey(REF_CH_PARAM); int inSituAberrCorrChannel = this.parameters.getIntValueForKey(IN_SITU_ABERR_SECOND_CH_PARAM); int measurementChannel = this.parameters.getIntValueForKey(CORR_CH_PARAM); List<ImageObject> iobjsForInSituAberrCorr = null; try { iobjsForInSituAberrCorr = FileUtils.readInSituAberrCorrPositionData(this.parameters); } catch (java.io.IOException e) { java.util.logging.Logger .getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME) .severe("Exception encountered while reading data for in situ aberration correction from disk: " + e.getMessage()); } catch (ClassNotFoundException e) { java.util.logging.Logger .getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME) .severe("Exception encountered while reading data for in situ aberration correction from disk: " + e.getMessage()); } RealVector xForCorr = new org.apache.commons.math3.linear.ArrayRealVector(iobjsForInSituAberrCorr.size(), 0.0); RealVector xForExpt = new org.apache.commons.math3.linear.ArrayRealVector(iobjsForInSituAberrCorr.size(), 0.0); RealVector yForCorr = new org.apache.commons.math3.linear.ArrayRealVector(iobjsForInSituAberrCorr.size(), 0.0); RealVector yForExpt = new org.apache.commons.math3.linear.ArrayRealVector(iobjsForInSituAberrCorr.size(), 0.0); RealVector zForCorr = new org.apache.commons.math3.linear.ArrayRealVector(iobjsForInSituAberrCorr.size(), 0.0); RealVector zForExpt = new org.apache.commons.math3.linear.ArrayRealVector(iobjsForInSituAberrCorr.size(), 0.0); for (int i = 0; i < iobjsForInSituAberrCorr.size(); i++) { ImageObject currObj = iobjsForInSituAberrCorr.get(i); RealVector corrDiff = currObj .getCorrectedVectorDifferenceBetweenChannels(referenceChannel, inSituAberrCorrChannel) .ebeMultiply(this.pixelToDistanceConversions); RealVector exptDiff = currObj .getCorrectedVectorDifferenceBetweenChannels(referenceChannel, measurementChannel) .ebeMultiply(this.pixelToDistanceConversions); xForCorr.setEntry(i, corrDiff.getEntry(0)); yForCorr.setEntry(i, corrDiff.getEntry(1)); zForCorr.setEntry(i, corrDiff.getEntry(2)); xForExpt.setEntry(i, exptDiff.getEntry(0)); yForExpt.setEntry(i, exptDiff.getEntry(1)); zForExpt.setEntry(i, exptDiff.getEntry(2)); } BisquareLinearFit bslf = new BisquareLinearFit(); bslf.disableIntercept(); RealVector xFit = bslf.fit(xForCorr, xForExpt); RealVector yFit = bslf.fit(yForCorr, yForExpt); RealVector zFit = bslf.fit(zForCorr, zForExpt); RealVector allSlopes = new org.apache.commons.math3.linear.ArrayRealVector(3, 0.0); allSlopes.setEntry(0, xFit.getEntry(0)); allSlopes.setEntry(1, yFit.getEntry(0)); allSlopes.setEntry(2, zFit.getEntry(0)); return allSlopes; }
From source file:edu.stanford.cfuller.colocalization3d.correction.PositionCorrector.java
/** * Determines the target registration error for a correction by successively leaving out each ImageObject in a set used to make a correction, * calculating a correction from the remaining objects, and assessing the error in correcting the object left out. * /*w w w . j av a 2 s . c o m*/ * @param imageObjects A Vector containing all the ImageObjects to be used for the correction * or in the order it appears in a multiwavelength image file. * @return The average value of the error over all objects. */ public double determineTRE(java.util.List<ImageObject> imageObjects) { int referenceChannel = this.parameters.getIntValueForKey(REF_CH_PARAM); int channelToCorrect = this.parameters.getIntValueForKey(CORR_CH_PARAM); RealVector treVector = new ArrayRealVector(imageObjects.size(), 0.0); RealVector treXYVector = new ArrayRealVector(imageObjects.size(), 0.0); java.util.Deque<TREThread> startedThreads = new java.util.LinkedList<TREThread>(); int maxThreads = 1; if (this.parameters.hasKey(THREAD_COUNT_PARAM)) { maxThreads = this.parameters.getIntValueForKey(THREAD_COUNT_PARAM); } final int threadWaitTime_ms = 1000; for (int removeIndex = 0; removeIndex < imageObjects.size(); removeIndex++) { if (removeIndex % 10 == 0) { java.util.logging.Logger .getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME) .finer("calulating TRE: point " + (removeIndex + 1) + " of " + imageObjects.size()); } TREThread nextFit = new TREThread(imageObjects, referenceChannel, channelToCorrect, removeIndex, this); if (startedThreads.size() < maxThreads) { startedThreads.add(nextFit); nextFit.start(); } else { TREThread next = startedThreads.poll(); try { next.join(threadWaitTime_ms); } catch (InterruptedException e) { e.printStackTrace(); } while (next.isAlive()) { startedThreads.add(next); next = startedThreads.poll(); try { next.join(threadWaitTime_ms); } catch (InterruptedException e) { e.printStackTrace(); } } treVector.setEntry(next.getRemoveIndex(), next.getTre()); treXYVector.setEntry(next.getRemoveIndex(), next.getTreXY()); startedThreads.add(nextFit); nextFit.start(); } } java.util.List<Integer> unsuccessful_TRE = new java.util.ArrayList<Integer>(); while (startedThreads.size() > 0) { TREThread next = startedThreads.poll(); try { next.join(); if (next.getSuccess()) { treVector.setEntry(next.getRemoveIndex(), next.getTre()); } else { unsuccessful_TRE.add(next.getRemoveIndex()); } } catch (InterruptedException e) { e.printStackTrace(); } } RealVector treVector_mod = new ArrayRealVector(treVector.getDimension() - unsuccessful_TRE.size()); RealVector treXYVector_mod = new ArrayRealVector(treVector_mod.getDimension()); int c = 0; //unsuccessful TRE calculation results when there is incomplete coverage in the correction dataset for (int i = 0; i < treVector.getDimension(); ++i) { if (!unsuccessful_TRE.contains(i)) { treVector_mod.setEntry(c, treVector.getEntry(i)); treXYVector_mod.setEntry(c, treXYVector.getEntry(i)); ++c; } } treVector = treVector_mod; treXYVector = treXYVector_mod; double tre = treVector.getL1Norm() / treVector.getDimension(); double xy_tre = (treXYVector.getL1Norm() / treXYVector.getDimension()); java.util.logging.Logger.getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME) .info("TRE: " + tre); java.util.logging.Logger.getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME) .info("x-y TRE: " + xy_tre); return tre; }