001: /*
002: * Geotools2 - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2002-2005, Geotools Project Managment Committee (PMC)
005: *
006: * This library is free software; you can redistribute it and/or
007: * modify it under the terms of the GNU Lesser General Public
008: * License as published by the Free Software Foundation;
009: * version 2.1 of the License.
010: *
011: * This library is distributed in the hope that it will be useful,
012: * but WITHOUT ANY WARRANTY; without even the implied warranty of
013: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
014: * Lesser General Public License for more details.
015: */
016: package org.geotools.referencing.operation.builder;
017:
018: import org.geotools.referencing.operation.matrix.GeneralMatrix;
019: import org.geotools.referencing.operation.transform.ProjectiveTransform;
020: import org.opengis.referencing.FactoryException;
021: import org.opengis.referencing.cs.CartesianCS;
022: import org.opengis.referencing.operation.MathTransform;
023: import org.opengis.geometry.DirectPosition;
024: import org.opengis.geometry.MismatchedDimensionException;
025: import org.opengis.geometry.MismatchedReferenceSystemException;
026: import java.util.List;
027: import javax.vecmath.MismatchedSizeException;
028:
029: /**
030: * Builds {@linkplain MathTransform
031: * MathTransform} setup as Projective transformation from a list of
032: * {@linkplain org.geotools.referencing.operation.builder.MappedPosition
033: * MappedPosition}. The calculation uses least square method. The Projective
034: * transform equation: (2D). The calculation uses least square method.
035: * Projective transform equation:<pre> [ x'] [ m00 m01 m02 ] [ x ]
036: * [ y'] = [ m10 m11 m12 ] [ y ]
037: * [ 1 ] [ m20 m21 1 ] [ 1 ] x' = m * x
038: * </pre>In the case that we have more identical points we can write it
039: * like this (in Matrix):<pre>
040: * [ x'<sub>1</sub> ] [ x<sub>1</sub> y<sub>1</sub> 1 0 0 0 -x'x -x'y] [ m00 ]
041: * [ x'<sub>2</sub> ] [ x<sub>2</sub> y<sub>2</sub> 1 0 0 0 -x'x -x'y] [ m01 ]
042: * [ . ] [ . ] [ m02 ]
043: * [ . ] [ . ] * [ m10 ]
044: * [ x'<sub>n</sub> ] = [ x<sub>n</sub> y<sub>n</sub> 1 0 0 0 -x'x -x'y] [ m11 ]
045: * [ y'<sub>1</sub> ] [ 0 0 0 x<sub>1</sub> y<sub>1</sub> 1 -y'x -y'y] [ m12 ]
046: * [ y'<sub>2</sub> ] [ 0 0 0 x<sub>2</sub> y<sub>2</sub> 1 -y'x -y'y] [ m20 ]
047: * [ . ] [ . ] [ m21 ]
048: * [ . ] [ . ]
049: * [ y'<sub>n</sub> ] [ 0 0 0 x<sub>n</sub> y<sub>n</sub> 1 -y'x -y'y]
050: * x' = A*m </pre>Using the least square method we get this result:
051: * <pre><blockquote>
052: * m = (A<sup>T</sup>PA)<sup>-1</sup> A<sup>T</sup>Px' </blockquote> </pre>
053: *
054: * @author Jan Jezek
055: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/builder/ProjectiveTransformBuilder.java $
056: * @version $Id: ProjectiveTransformBuilder.java 28982 2008-01-28 16:27:33Z acuster $
057: * @since 2.4
058: */
059: public class ProjectiveTransformBuilder extends MathTransformBuilder {
060: /** Matrix of derivations */
061: protected GeneralMatrix A;
062:
063: /** Matrix of wights */
064: protected GeneralMatrix P = null;
065:
066: /** Matrix of target values */
067: protected GeneralMatrix X;
068:
069: protected ProjectiveTransformBuilder() {
070: }
071:
072: /**
073: * Creates ProjectiveTransformBuilder for the set of properties.
074: *
075: *
076: * @param vectors list of {@linkplain MappedPosition
077: * MappedPosition}
078: * @throws MismatchedSizeException
079: * if the number of properties is not set properly.
080: * @throws MismatchedDimensionException
081: * if the dimension of properties is not set properly.
082: * @throws MismatchedReferenceSystemException
083: * -if there is mismatch in coordinate system in
084: * {@linkplain MappedPosition MappedPosition}
085: */
086: public ProjectiveTransformBuilder(List vectors)
087: throws MismatchedSizeException,
088: MismatchedDimensionException,
089: MismatchedReferenceSystemException {
090: super .setMappedPositions(vectors);
091: }
092:
093: /**
094: * Returns the minimum number of points required by this builder,
095: * which is 4 by default. Subclasses like {@linkplain AffineTransformBuilder
096: * affine transform builders} will reduce this minimum.
097: *
098: * @return minimum number of points required by this builder, which is 4 by
099: * default.
100: */
101: public int getMinimumPointCount() {
102: return 4;
103: }
104:
105: /**
106: * Returns the required coordinate system type, which is
107: * {@linkplain CartesianCS cartesian CS}.
108: *
109: * @return required coordinate system type
110: */
111: public Class /* <? extends CartesianCS> */getCoordinateSystemType() {
112: return CartesianCS.class;
113: }
114:
115: /**
116: * Fills P matrix for m = (A<sup>T</sup>PA)<sup>-1</sup>
117: * A<sup>T</sup>Px' equation
118: *
119: * @throws MissingInfoException if accuracy is not defined.
120: */
121: protected void fillPMatrix() throws MissingInfoException {
122: this .P = new GeneralMatrix(getMappedPositions().size() * 2,
123: getMappedPositions().size() * 2);
124:
125: for (int i = 0; i < getMappedPositions().size(); i = i + 2) {
126: if (Double.compare((((MappedPosition) getMappedPositions()
127: .get(i)).getAccuracy()), Double.NaN) == 0) {
128: throw new MissingInfoException(
129: "Accuracy has to be defined for all points");
130: }
131:
132: // weight for x
133: P.setElement(i, i,
134: 1 / ((MappedPosition) getMappedPositions().get(i))
135: .getAccuracy());
136: // weight for y
137: P.setElement(i + 1, i + 1,
138: 1 / ((MappedPosition) getMappedPositions().get(i))
139: .getAccuracy());
140: }
141: }
142:
143: /**
144: * Fills A matrix for m = (A<sup>T</sup>PA)<sup>-1</sup>
145: * A<sup>T</sup>Px' equation
146: */
147: protected void fillAMatrix() {
148: final DirectPosition[] sourcePoints = getSourcePoints();
149: final DirectPosition[] targetPoints = getTargetPoints();
150: A = new GeneralMatrix(2 * sourcePoints.length, 8);
151:
152: int numRow = 2 * sourcePoints.length;
153:
154: // fill first half of matrix
155: for (int j = 0; j < ((2 * sourcePoints.length) / 2); j++) {
156: double xs = sourcePoints[j].getCoordinates()[0];
157: double ys = sourcePoints[j].getCoordinates()[1];
158: double xd = targetPoints[j].getCoordinates()[0];
159:
160: A.setRow(j, new double[] { xs, ys, 1, 0, 0, 0, -xd * xs,
161: -xd * ys });
162: }
163:
164: // fill second half
165: for (int j = numRow / 2; j < numRow; j++) {
166: double xs = sourcePoints[j - (numRow / 2)].getCoordinates()[0];
167: double ys = sourcePoints[j - (numRow / 2)].getCoordinates()[1];
168: double yd = targetPoints[j - (numRow / 2)].getCoordinates()[1];
169:
170: A.setRow(j, new double[] { 0, 0, 0, xs, ys, 1, -yd * xs,
171: -yd * ys });
172: }
173: }
174:
175: /**
176: * Fills x' matrix for m = (A<sup>T</sup>PA)<sup>-1</sup>
177: * A<sup>T</sup>Px' equation
178: */
179: protected void fillXMatrix() {
180: X = new GeneralMatrix(2 * getTargetPoints().length, 1);
181:
182: int numRow = X.getNumRow();
183:
184: // Creates X matrix
185: for (int j = 0; j < (numRow / 2); j++) {
186: X
187: .setElement(j, 0, getTargetPoints()[j]
188: .getCoordinates()[0]);
189: }
190:
191: for (int j = numRow / 2; j < numRow; j++) {
192: X.setElement(j, 0, getTargetPoints()[j - (numRow / 2)]
193: .getCoordinates()[1]);
194: }
195: }
196:
197: /**
198: * Switch whether to include weights into the calculation. Weights are derived from each point accuracy.
199: * Weight p = 1 / accuracy<sup>2<sup>.
200: * @param include if true then the weights will be included onto the calculation. False is default.
201: * @throws FactoryException if all or some of the {@linkplain #setMappedPositions(List) points} does not have accuracy setup properly.
202: */
203: public void includeWeights(boolean include)
204: throws MissingInfoException {
205: this .P = new GeneralMatrix(getMappedPositions().size() * 2,
206: getMappedPositions().size() * 2);
207:
208: if (include) {
209: fillPMatrix();
210: } else {
211: for (int j = 0; j < getMappedPositions().size(); j++) {
212: P.setElement(j, j, 1);
213: }
214: }
215: }
216:
217: /**
218: * Calculates the parameters using the least square method. The
219: * equation:
220: * <pre><blockquote>
221: * m = (A<sup>T</sup>A)<sup>-1</sup> A<sup>T</sup>x'
222: * </blockquote> </pre>
223: *
224: * @return m matrix.
225: */
226: protected double[] calculateLSM() {
227: fillAMatrix();
228: // fillPMatrix();
229: fillXMatrix();
230:
231: if (P == null) {
232: try {
233: includeWeights(false);
234: } catch (FactoryException e) {
235: // should never reach here - weights are not included
236: }
237: }
238:
239: GeneralMatrix AT = (GeneralMatrix) A.clone();
240: AT.transpose();
241:
242: GeneralMatrix ATP = new GeneralMatrix(AT.getNumRow(), P
243: .getNumCol());
244: GeneralMatrix ATPA = new GeneralMatrix(AT.getNumRow(), A
245: .getNumCol());
246: GeneralMatrix ATPX = new GeneralMatrix(AT.getNumRow(), 1);
247: GeneralMatrix x = new GeneralMatrix(A.getNumCol(), 1);
248: ATP.mul(AT, P); // ATP
249: ATPA.mul(ATP, A); // ATPA
250: ATPX.mul(ATP, X); // ATPX
251: ATPA.invert();
252: x.mul(ATPA, ATPX);
253: ATPA.invert();
254:
255: x.transpose();
256:
257: return x.getElements()[0];
258: }
259:
260: /**
261: * Returns the matrix of parameters for Projective transformation.
262: * This method should by override for the special cases like affine or
263: * similar transformation. The M matrix looks like this:<pre>
264: *
265: * [ m00 m01 m02 ]
266: * [ m10 m11 m12 ]
267: * [ m20 m21 1 ]
268: * </pre>
269: *
270: * @return Matrix M
271: */
272: protected GeneralMatrix getProjectiveMatrix() {
273: GeneralMatrix M = new GeneralMatrix(3, 3);
274:
275: // double[] param = generateMMatrix();
276: double[] param = calculateLSM();
277: double[] m0 = { param[0], param[1], param[2] };
278: double[] m1 = { param[3], param[4], param[5] };
279: double[] m2 = { param[6], param[7], 1 };
280: M.setRow(0, m0);
281: M.setRow(1, m1);
282: M.setRow(2, m2);
283:
284: return M;
285: }
286:
287: protected MathTransform computeMathTransform() {
288: return ProjectiveTransform.create(getProjectiveMatrix());
289: }
290: }
|