View Javadoc

1   /**
2    * Copyright (c) 2011, The University of Southampton and the individual contributors.
3    * All rights reserved.
4    *
5    * Redistribution and use in source and binary forms, with or without modification,
6    * are permitted provided that the following conditions are met:
7    *
8    *   * 	Redistributions of source code must retain the above copyright notice,
9    * 	this list of conditions and the following disclaimer.
10   *
11   *   *	Redistributions in binary form must reproduce the above copyright notice,
12   * 	this list of conditions and the following disclaimer in the documentation
13   * 	and/or other materials provided with the distribution.
14   *
15   *   *	Neither the name of the University of Southampton nor the names of its
16   * 	contributors may be used to endorse or promote products derived from this
17   * 	software without specific prior written permission.
18   *
19   * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
20   * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
21   * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
22   * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
23   * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
24   * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
25   * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
26   * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27   * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
28   * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29   */
30  package org.openimaj.image.pixel.statistics;
31  
32  import java.util.Arrays;
33  
34  import org.apache.commons.math.DimensionMismatchException;
35  import org.apache.commons.math.stat.descriptive.MultivariateSummaryStatistics;
36  import org.openimaj.image.MBFImage;
37  import org.openimaj.image.pixel.sampling.FLineSampler;
38  import org.openimaj.math.geometry.line.Line2d;
39  import org.openimaj.math.geometry.point.Point2d;
40  import org.openimaj.math.geometry.point.Point2dImpl;
41  import org.openimaj.math.util.FloatArrayStatsUtils;
42  import org.openimaj.util.array.ArrayUtils;
43  
44  import Jama.Matrix;
45  
46  /**
47   * An {@link MBFStatisticalPixelProfileModel} is a statistical model of pixels
48   * from an {@link MBFImage} sampled along a line.
49   * 
50   * The model allows for various sampling strategies (see {@link FLineSampler})
51   * and uses the mean and covariance as its internal state.
52   * 
53   * The model is updateable, but does not hold on to previously seen samples to
54   * reduce memory usage.
55   * 
56   * Internally, the model is one-dimensional, and is created by stacking the
57   * samples from each image band.
58   * 
59   * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
60   */
61  @SuppressWarnings("deprecation")
62  public class MBFStatisticalPixelProfileModel implements PixelProfileModel<MBFImage> {
63  	private MultivariateSummaryStatistics statistics;
64  	private int nsamples;
65  	private FLineSampler sampler;
66  
67  	private double[] mean;
68  	private Matrix invCovar;
69  
70  	/**
71  	 * Construct a new {@link MBFStatisticalPixelProfileModel} with the given
72  	 * number of samples per line, and the given sampling strategy.
73  	 * 
74  	 * @param nsamples
75  	 *            number of samples
76  	 * @param sampler
77  	 *            line sampling strategy
78  	 */
79  	public MBFStatisticalPixelProfileModel(int nsamples, FLineSampler sampler) {
80  		this.nsamples = nsamples;
81  		this.statistics = new MultivariateSummaryStatistics(nsamples, true);
82  		this.sampler = sampler;
83  	}
84  
85  	private float[] normaliseSamples(float[] samples) {
86  		final float sum = FloatArrayStatsUtils.sum(samples);
87  
88  		for (int i = 0; i < samples.length; i++) {
89  			samples[i] /= sum;
90  		}
91  
92  		return samples;
93  	}
94  
95  	@Override
96  	public void updateModel(MBFImage image, Line2d line) {
97  		final float[] samples = extractNormalisedStacked(image, line);
98  
99  		try {
100 			statistics.addValue(ArrayUtils.convertToDouble(samples));
101 		} catch (final DimensionMismatchException e) {
102 			throw new RuntimeException(e);
103 		}
104 
105 		invCovar = null;
106 		mean = null;
107 	}
108 
109 	private float[][] extract(MBFImage image, Line2d line, int numSamples) {
110 		final float[][] samples = new float[image.numBands()][];
111 
112 		for (int i = 0; i < image.numBands(); i++) {
113 			samples[i] = sampler.extractSamples(line, image.getBand(i), numSamples);
114 		}
115 
116 		return samples;
117 	}
118 
119 	/**
120 	 * Extract normalised stacked samples b1b1b1bb2b2b2b3b3b3...
121 	 * 
122 	 * @param image
123 	 * @param line
124 	 * @return
125 	 */
126 	private float[] extractNormalisedStacked(MBFImage image, Line2d line) {
127 		final float[] samples = new float[nsamples * image.numBands()];
128 
129 		for (int i = 0; i < image.numBands(); i++) {
130 			final float[] s = sampler.extractSamples(line, image.getBand(i), nsamples);
131 			System.arraycopy(s, 0, samples, i * nsamples, nsamples);
132 		}
133 
134 		return normaliseSamples(samples);
135 	}
136 
137 	/**
138 	 * @return the mean of the model
139 	 */
140 	public double[] getMean() {
141 		if (mean == null) {
142 			mean = statistics.getMean();
143 			invCovar = new Matrix(statistics.getCovariance().getData()).inverse();
144 		}
145 
146 		return mean;
147 	}
148 
149 	/**
150 	 * @return the covariance of the model
151 	 */
152 	public Matrix getCovariance() {
153 		if (mean == null) {
154 			mean = statistics.getMean();
155 			invCovar = new Matrix(statistics.getCovariance().getData()).inverse();
156 		}
157 		return new Matrix(statistics.getCovariance().getData());
158 	}
159 
160 	/**
161 	 * @return the inverse of the covariance matrix
162 	 */
163 	public Matrix getInverseCovariance() {
164 		if (mean == null) {
165 			mean = statistics.getMean();
166 			invCovar = new Matrix(statistics.getCovariance().getData()).inverse();
167 		}
168 		return invCovar;
169 	}
170 
171 	/**
172 	 * Compute the Mahalanobis distance of the given vector to the internal
173 	 * model. The vector must have the same size as the number of samples given
174 	 * during construction.
175 	 * 
176 	 * @param vector
177 	 *            the vector
178 	 * @return the computed Mahalanobis distance
179 	 */
180 	public float computeMahalanobis(float[] vector) {
181 		if (mean == null) {
182 			mean = statistics.getMean();
183 			try {
184 				invCovar = new Matrix(statistics.getCovariance().getData()).inverse();
185 			} catch (final RuntimeException e) {
186 				invCovar = Matrix.identity(nsamples, nsamples);
187 			}
188 		}
189 
190 		final double[] meanCentered = new double[mean.length];
191 		for (int i = 0; i < mean.length; i++) {
192 			meanCentered[i] = vector[i] - mean[i];
193 		}
194 
195 		final Matrix mct = new Matrix(new double[][] { meanCentered });
196 		final Matrix mc = mct.transpose();
197 
198 		final Matrix dist = mct.times(invCovar).times(mc);
199 
200 		return (float) dist.get(0, 0);
201 	}
202 
203 	/**
204 	 * Compute the Mahalanobis distance of a vector of samples extracted along a
205 	 * line in the given image to the internal model.
206 	 * 
207 	 * @param image
208 	 *            the image to sample
209 	 * @param line
210 	 *            the line to sample along
211 	 * @return the computed Mahalanobis distance
212 	 */
213 	public float computeMahalanobis(MBFImage image, Line2d line) {
214 		final float[] samples = extractNormalisedStacked(image, line);
215 		return computeMahalanobis(samples);
216 	}
217 
218 	/**
219 	 * Extract numSamples samples from the line in the image and then compare
220 	 * this model at each overlapping position starting from the first sample at
221 	 * the beginning of the line.
222 	 * 
223 	 * numSamples must be bigger than the number of samples used to construct
224 	 * the model. In addition, callers are responsible for ensuring the sampling
225 	 * rate between the new samples and the model is equal.
226 	 * 
227 	 * @param image
228 	 *            the image to sample
229 	 * @param line
230 	 *            the line to sample along
231 	 * @param numSamples
232 	 *            the number of samples to make
233 	 * @return an array of the computed Mahalanobis distances at each offset
234 	 */
235 	public float[] computeMahalanobisWindowed(MBFImage image, Line2d line, int numSamples) {
236 		final float[][] samples = extract(image, line, numSamples);
237 		return computeMahalanobisWindowed(samples);
238 	}
239 
240 	@Override
241 	public Point2d computeNewBest(MBFImage image, Line2d line, int numSamples) {
242 		final float[] resp = computeMahalanobisWindowed(image, line, numSamples);
243 
244 		final int minIdx = ArrayUtils.minIndex(resp);
245 		final int offset = (numSamples - nsamples) / 2;
246 
247 		if (resp[offset] == resp[minIdx]) // prefer the centre over another
248 											// value if same response
249 			return line.calculateCentroid();
250 
251 		// the sample line might be different, so we need to measure relative to
252 		// it...
253 		line = this.sampler.getSampleLine(line, image.getBand(0), numSamples);
254 
255 		final float x = line.begin.getX();
256 		final float y = line.begin.getY();
257 		final float dxStep = (line.end.getX() - x) / (numSamples - 1);
258 		final float dyStep = (line.end.getY() - y) / (numSamples - 1);
259 
260 		return new Point2dImpl(x + (minIdx + offset) * dxStep, y + (minIdx + offset) * dyStep);
261 	}
262 
263 	@Override
264 	public float computeMovementDistance(MBFImage image, Line2d line, int numSamples, Point2d pt) {
265 		final Line2d sampleLine = sampler.getSampleLine(line, image.getBand(0), numSamples);
266 
267 		return (float) (2 * Line2d.distance(sampleLine.calculateCentroid(), pt) / sampleLine.calculateLength());
268 	}
269 
270 	/**
271 	 * Compare this model at each overlapping position of the given vector
272 	 * starting from the first sample and return the distance for each overlap.
273 	 * 
274 	 * The length of the vector must be bigger than the number of samples used
275 	 * to construct the model. In addition, callers are responsible for ensuring
276 	 * the sampling rate between the new samples and the model is equal.
277 	 * 
278 	 * @param vector
279 	 *            array of samples, one vector per band
280 	 * @return an array of the computed Mahalanobis distances at each offset
281 	 */
282 	public float[] computeMahalanobisWindowed(float[][] vector) {
283 		final int maxShift = vector.length - nsamples + 1;
284 
285 		final float[] responses = new float[maxShift];
286 		float[] samples = new float[nsamples * vector.length];
287 		for (int i = 0; i < maxShift; i++) {
288 			for (int j = 0; j < vector.length; j++) {
289 				System.arraycopy(vector[j], i, samples, nsamples * j, nsamples);
290 			}
291 			samples = normaliseSamples(samples);
292 			responses[i] = computeMahalanobis(samples);
293 		}
294 
295 		return responses;
296 	}
297 
298 	@Override
299 	public String toString() {
300 		return "\nPixelProfileModel[\n" +
301 				"\tcount = " + statistics.getN() + "\n" +
302 				"\tmean = " + Arrays.toString(statistics.getMean()) + "\n" +
303 				"\tcovar = " + statistics.getCovariance() + "\n" +
304 				"]";
305 	}
306 
307 	/**
308 	 * @return the number of samples used along each profile line
309 	 */
310 	public int getNumberSamples() {
311 		return nsamples;
312 	}
313 
314 	/**
315 	 * @return the sampler used by the model to extract samples along profiles
316 	 */
317 	public FLineSampler getSampler() {
318 		return sampler;
319 	}
320 
321 	@Override
322 	public float computeCost(MBFImage image, Line2d line) {
323 		return computeMahalanobis(image, line);
324 	}
325 }