001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2003-2006, GeoTools Project Managment Committee (PMC)
005: * (C) 2003, Institut de Recherche pour le Développement
006: * (C) 1998, Pêches et Océans Canada
007: *
008: * This library is free software; you can redistribute it and/or
009: * modify it under the terms of the GNU Lesser General Public
010: * License as published by the Free Software Foundation;
011: * version 2.1 of the License.
012: *
013: * This library is distributed in the hope that it will be useful,
014: * but WITHOUT ANY WARRANTY; without even the implied warranty of
015: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
016: * Lesser General Public License for more details.
017: */
018: package org.geotools.math;
019:
020: // J2SE and vecmath dependencies
021: import java.io.Serializable;
022: import javax.vecmath.MismatchedSizeException;
023: import javax.vecmath.Point3d;
024:
025: // OpenGIS dependencies
026: import org.opengis.util.Cloneable;
027:
028: /**
029: * Equation of a plane in a three-dimensional space (<var>x</var>,<var>y</var>,<var>z</var>).
030: * The plane equation is expressed by {@link #c}, {@link #cx} and {@link #cy} coefficients as
031: * below:
032: *
033: * <blockquote>
034: * <var>z</var>(<var>x</var>,<var>y</var>) = <var>c</var> +
035: * <var>cx</var>*<var>x</var> + <var>cy</var>*<var>y</var>
036: * </blockquote>
037: *
038: * Those coefficients can be set directly, or computed by a linear regression of this plane
039: * through a set of three-dimensional points.
040: *
041: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/metadata/src/main/java/org/geotools/math/Plane.java $
042: * @version $Id: Plane.java 22443 2006-10-27 20:47:22Z desruisseaux $
043: * @author Martin Desruisseaux
044: * @author Howard Freeland
045: *
046: * @since 2.0
047: */
048: public class Plane implements Cloneable, Serializable {
049: /**
050: * Serial number for compatibility with different versions.
051: */
052: private static final long serialVersionUID = 2956201711131316723L;
053:
054: /**
055: * The <var>c</var> coefficient for this plane. This coefficient appears in the place equation
056: * <var><strong>c</strong></var>+<var>cx</var>*<var>x</var>+<var>cy</var>*<var>y</var>.
057: */
058: public double c;
059:
060: /**
061: * The <var>cx</var> coefficient for this plane. This coefficient appears in the place equation
062: * <var>c</var>+<var><strong>cx</strong></var>*<var>x</var>+<var>cy</var>*<var>y</var>.
063: */
064: public double cx;
065:
066: /**
067: * The <var>cy</var> coefficient for this plane. This coefficient appears in the place equation
068: * <var>c</var>+<var>cx</var>*<var>x</var>+<var><strong>cy</strong></var>*<var>y</var>.
069: */
070: public double cy;
071:
072: /**
073: * Construct a new plane. All coefficients are set to 0.
074: */
075: public Plane() {
076: }
077:
078: /**
079: * Compute the <var>z</var> value for the specified (<var>x</var>,<var>y</var>) point.
080: * The <var>z</var> value is computed using the following equation:
081: *
082: * <blockquote><code>
083: * z(x,y) = c + cx*x + cy*y
084: * </code></blockquote>
085: *
086: * @param x The <var>x</var> value.
087: * @param y The <var>y</var> value.
088: * @return The <var>z</var> value.
089: */
090: public final double z(final double x, final double y) {
091: return c + cx * x + cy * y;
092: }
093:
094: /**
095: * Compute the <var>y</var> value for the specified (<var>x</var>,<var>z</var>) point.
096: * The <var>y</var> value is computed using the following equation:
097: *
098: * <blockquote><code>
099: * y(x,z) = (z - (c+cx*x)) / cy
100: * </code></blockquote>
101: *
102: * @param x The <var>x</var> value.
103: * @param z The <var>y</var> value.
104: * @return The <var>y</var> value.
105: */
106: public final double y(final double x, final double z) {
107: return (z - (c + cx * x)) / cy;
108: }
109:
110: /**
111: * Compute the <var>x</var> value for the specified (<var>y</var>,<var>z</var>) point.
112: * The <var>x</var> value is computed using the following equation:
113: *
114: * <blockquote><code>
115: * x(y,z) = (z - (c+cy*y)) / cx
116: * </code></blockquote>
117: *
118: * @param y The <var>x</var> value.
119: * @param z The <var>y</var> value.
120: * @return The <var>x</var> value.
121: */
122: public final double x(final double y, final double z) {
123: return (z - (c + cy * y)) / cx;
124: }
125:
126: /**
127: * Computes the plane's coefficients from the specified points. Three points
128: * are enough for determining exactly the plan, providing that the points are
129: * not colinear.
130: *
131: * @throws ArithmeticException If the three points are colinear.
132: */
133: public void setPlane(final Point3d P1, final Point3d P2,
134: final Point3d P3) throws ArithmeticException {
135: final double m00 = P2.x * P3.y - P3.x * P2.y;
136: final double m01 = P3.x * P1.y - P1.x * P3.y;
137: final double m02 = P1.x * P2.y - P2.x * P1.y;
138: final double det = m00 + m01 + m02;
139: if (det == 0) {
140: throw new ArithmeticException("Points are colinear");
141: }
142: c = ((m00) * P1.z + (m01) * P2.z + (m02) * P3.z) / det;
143: cx = ((P2.y - P3.y) * P1.z + (P3.y - P1.y) * P2.z + (P1.y - P2.y)
144: * P3.z)
145: / det;
146: cy = ((P3.x - P2.x) * P1.z + (P1.x - P3.x) * P2.z + (P2.x - P1.x)
147: * P3.z)
148: / det;
149: }
150:
151: /**
152: * Compute the plane's coefficients from a set of points. This method use
153: * a linear regression in the least-square sense. Result is undertermined
154: * if all points are colinear.
155: *
156: * @param x vector of <var>x</var> coordinates
157: * @param y vector of <var>y</var> coordinates
158: * @param z vector of <var>z</var> values
159: *
160: * @throws MismatchedSizeException if <var>x</var>, <var>y</var> and <var>z</var>
161: * don't have the same length.
162: */
163: public void setPlane(final double[] x, final double[] y,
164: final double[] z) throws MismatchedSizeException {
165: final int N = x.length;
166: if (N != y.length || N != z.length) {
167: throw new MismatchedSizeException();
168: }
169: double sum_x = 0;
170: double sum_y = 0;
171: double sum_z = 0;
172: double sum_xx = 0;
173: double sum_yy = 0;
174: double sum_xy = 0;
175: double sum_zx = 0;
176: double sum_zy = 0;
177: for (int i = 0; i < N; i++) {
178: final double xi = x[i];
179: final double yi = y[i];
180: final double zi = z[i];
181: sum_x += xi;
182: sum_y += yi;
183: sum_z += zi;
184: sum_xx += xi * xi;
185: sum_yy += yi * yi;
186: sum_xy += xi * yi;
187: sum_zx += zi * xi;
188: sum_zy += zi * yi;
189: }
190: /*
191: * ( sum_zx - sum_z*sum_x ) = cx*(sum_xx - sum_x*sum_x) + cy*(sum_xy - sum_x*sum_y)
192: * ( sum_zy - sum_z*sum_y ) = cx*(sum_xy - sum_x*sum_y) + cy*(sum_yy - sum_y*sum_y)
193: */
194: final double ZX = sum_zx - sum_z * sum_x / N;
195: final double ZY = sum_zy - sum_z * sum_y / N;
196: final double XX = sum_xx - sum_x * sum_x / N;
197: final double XY = sum_xy - sum_x * sum_y / N;
198: final double YY = sum_yy - sum_y * sum_y / N;
199: final double den = (XY * XY - XX * YY);
200:
201: cy = (ZX * XY - ZY * XX) / den;
202: cx = (ZY * XY - ZX * YY) / den;
203: c = (sum_z - (cx * sum_x + cy * sum_y)) / N;
204: }
205:
206: /**
207: * Returns a string representation of this plane.
208: * The string will contains the plane's equation, as below:
209: *
210: * <blockquote>
211: * <var>z</var>(<var>x</var>,<var>y</var>) = {@link #c} +
212: * {@link #cx}*<var>x</var> + {@link #cy}*<var>y</var>
213: * </blockquote>
214: */
215: public String toString() {
216: final StringBuffer buffer = new StringBuffer("z(x,y)= ");
217: if (c == 0 && cx == 0 && cy == 0) {
218: buffer.append(0);
219: } else {
220: if (c != 0) {
221: buffer.append(c);
222: buffer.append(" + ");
223: }
224: if (cx != 0) {
225: buffer.append(cx);
226: buffer.append("*x");
227: if (cy != 0) {
228: buffer.append(" + ");
229: }
230: }
231: if (cy != 0) {
232: buffer.append(cy);
233: buffer.append("*y");
234: }
235: }
236: return buffer.toString();
237: }
238:
239: /**
240: * Compare this plane with the specified object for equality.
241: */
242: public boolean equals(final Object object) {
243: if (object != null && getClass().equals(object.getClass())) {
244: final Plane that = (Plane) object;
245: return Double.doubleToLongBits(this .c) == Double
246: .doubleToLongBits(that.c)
247: && Double.doubleToLongBits(this .cy) == Double
248: .doubleToLongBits(that.cx)
249: && Double.doubleToLongBits(this .cx) == Double
250: .doubleToLongBits(that.cy);
251: } else {
252: return false;
253: }
254: }
255:
256: /**
257: * Returns a hash code value for this plane.
258: */
259: public int hashCode() {
260: final long code = Double.doubleToLongBits(c)
261: + 37
262: * (Double.doubleToLongBits(cy) + 37 * (Double
263: .doubleToLongBits(cx)));
264: return (int) code ^ (int) (code >>> 32);
265: }
266:
267: /**
268: * Returns a clone of this plane.
269: */
270: public Object clone() {
271: try {
272: return super .clone();
273: } catch (CloneNotSupportedException exception) {
274: throw new AssertionError(exception);
275: }
276: }
277: }
|