001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021package org.biojava.nbio.structure.geometry;
022
023import javax.vecmath.Matrix4d;
024import javax.vecmath.Point3d;
025import javax.vecmath.Vector3d;
026
027import org.biojava.nbio.structure.jama.Matrix;
028
029/**
030 * Utility operations on Point3d.
031 * 
032 * @author Aleix Lafita
033 * @since 5.0.0
034 *
035 */
036public class CalcPoint {
037
038        /** Prevent instantiation */
039        private CalcPoint() {
040        }
041
042        /**
043         * Center a cloud of points. This means subtracting the {@lin
044         * #centroid(Point3d[])} of the cloud to each point.
045         * 
046         * @param x
047         *            array of points. Point objects will be modified
048         */
049        public static void center(Point3d[] x) {
050                Point3d center = centroid(x);
051                center.negate();
052                translate(new Vector3d(center), x);
053        }
054
055        /**
056         * Calculate the centroid of the point cloud.
057         * 
058         * @param x
059         *            array of points. Point objects will not be modified
060         * @return centroid as Point3d
061         */
062        public static Point3d centroid(Point3d[] x) {
063                Point3d center = new Point3d();
064                for (Point3d p : x) {
065                        center.add(p);
066                }
067                center.scale(1.0 / x.length);
068                return center;
069        }
070
071        /**
072         * Transform all points with a 4x4 transformation matrix.
073         * 
074         * @param rotTrans
075         *            4x4 transformation matrix
076         * @param x
077         *            array of points. Point objects will be modified
078         */
079        public static void transform(Matrix4d rotTrans, Point3d[] x) {
080                for (Point3d p : x) {
081                        rotTrans.transform(p);
082                }
083        }
084
085        /**
086         * Translate all points with a translation vector.
087         * 
088         * @param trans
089         *            the translation vector to apply
090         * @param x
091         *            array of points. Point objects will be modified
092         */
093        public static void translate(Vector3d trans, Point3d[] x) {
094                for (Point3d p : x) {
095                        p.add(trans);
096                }
097        }
098
099        /**
100         * Clone an array of points.
101         * 
102         * @param x
103         *            original array of points. Point objects will not be modified
104         * @return new array of points, identical clone of x
105         */
106        public static Point3d[] clonePoint3dArray(Point3d[] x) {
107                Point3d[] clone = new Point3d[x.length];
108                for (int i = 0; i < x.length; i++) {
109                        clone[i] = new Point3d(x[i]);
110                }
111                return clone;
112        }
113
114        /*
115         * Peter can you document this method? TODO
116         * 
117         * @param moved
118         * @param fixed
119         * @return
120         */
121        public static Matrix formMatrix(Point3d[] a, Point3d[] b) {
122                double xx = 0.0, xy = 0.0, xz = 0.0;
123                double yx = 0.0, yy = 0.0, yz = 0.0;
124                double zx = 0.0, zy = 0.0, zz = 0.0;
125
126                for (int i = 0; i < a.length; i++) {
127                        xx += a[i].x * b[i].x;
128                        xy += a[i].x * b[i].y;
129                        xz += a[i].x * b[i].z;
130                        yx += a[i].y * b[i].x;
131                        yy += a[i].y * b[i].y;
132                        yz += a[i].y * b[i].z;
133                        zx += a[i].z * b[i].x;
134                        zy += a[i].z * b[i].y;
135                        zz += a[i].z * b[i].z;
136                }
137
138                double[][] f = new double[4][4];
139                f[0][0] = xx + yy + zz;
140                f[0][1] = zy - yz;
141                f[1][0] = f[0][1];
142                f[1][1] = xx - yy - zz;
143                f[0][2] = xz - zx;
144                f[2][0] = f[0][2];
145                f[1][2] = xy + yx;
146                f[2][1] = f[1][2];
147                f[2][2] = yy - zz - xx;
148                f[0][3] = yx - xy;
149                f[3][0] = f[0][3];
150                f[1][3] = zx + xz;
151                f[3][1] = f[1][3];
152                f[2][3] = yz + zy;
153                f[3][2] = f[2][3];
154                f[3][3] = zz - xx - yy;
155
156                return new Matrix(f);
157        }
158
159        /**
160         * Returns the TM-Score for two superimposed sets of coordinates Yang Zhang
161         * and Jeffrey Skolnick, PROTEINS: Structure, Function, and Bioinformatics
162         * 57:702–710 (2004)
163         * 
164         * @param x
165         *            coordinate set 1
166         * @param y
167         *            coordinate set 2
168         * @param lengthNative
169         *            total length of native sequence
170         * @return
171         */
172        public static double TMScore(Point3d[] x, Point3d[] y, int lengthNative) {
173                
174                if (x.length != y.length) {
175                        throw new IllegalArgumentException(
176                                        "Point arrays are not of the same length.");
177                }
178                
179                double d0 = 1.24 * Math.cbrt(x.length - 15.0) - 1.8;
180                double d0Sq = d0 * d0;
181
182                double sum = 0;
183                for (int i = 0; i < x.length; i++) {
184                        sum += 1.0 / (1.0 + x[i].distanceSquared(y[i]) / d0Sq);
185                }
186
187                return sum / lengthNative;
188        }
189
190        /*
191         * Needs documentation!
192         * 
193         * @param x
194         * 
195         * @param y
196         * 
197         * @return
198         */
199        public static double GTSlikeScore(Point3d[] x, Point3d[] y) {
200                
201                if (x.length != y.length) {
202                        throw new IllegalArgumentException(
203                                        "Point arrays are not of the same length.");
204                }
205                
206                int contacts = 0;
207
208                for (Point3d px : x) {
209                        double minDist = Double.MAX_VALUE;
210
211                        for (Point3d py : y) {
212                                minDist = Math.min(minDist, px.distanceSquared(py));
213                        }
214
215                        if (minDist > 64)
216                                continue;
217                        contacts++;
218
219                        if (minDist > 16)
220                                continue;
221                        contacts++;
222
223                        if (minDist > 4)
224                                continue;
225                        contacts++;
226
227                        if (minDist > 1)
228                                continue;
229                        contacts++;
230                }
231
232                return contacts * 25.0 / x.length;
233        }
234
235        /**
236         * Calculate the RMSD of two point arrays, already superposed.
237         * 
238         * @param x
239         *            array of points superposed to y
240         * @param y
241         *            array of points superposed to x
242         * @return RMSD
243         */
244        public static double rmsd(Point3d[] x, Point3d[] y) {
245                
246                if (x.length != y.length) {
247                        throw new IllegalArgumentException(
248                                        "Point arrays are not of the same length.");
249                }
250                
251                double sum = 0.0;
252                for (int i = 0; i < x.length; i++) {
253                        sum += x[i].distanceSquared(y[i]);
254                }
255                return Math.sqrt(sum / x.length);
256        }
257
258        /*
259         * Needs documentation!
260         * 
261         * @param x
262         * 
263         * @param y
264         * 
265         * @return
266         */
267        public static double rmsdMin(Point3d[] x, Point3d[] y) {
268                
269                if (x.length != y.length) {
270                        throw new IllegalArgumentException(
271                                        "Point arrays are not of the same length.");
272                }
273                
274                double sum = 0.0;
275                for (int i = 0; i < x.length; i++) {
276                        double minDist = Double.MAX_VALUE;
277                        for (int j = 0; j < y.length; j++) {
278                                minDist = Math.min(minDist, x[i].distanceSquared(y[j]));
279                        }
280                        sum += minDist;
281                }
282                return Math.sqrt(sum / x.length);
283        }
284
285        /*
286         * Needs documentation!
287         * 
288         * @param x
289         * 
290         * @param y
291         * 
292         * @param maxDistance
293         * 
294         * @return
295         */
296        public static int contacts(Point3d[] x, Point3d[] y, double maxDistance) {
297                int contacts = 0;
298                for (int i = 0; i < x.length; i++) {
299                        double minDist = Double.MAX_VALUE;
300                        for (int j = 0; j < y.length; j++) {
301                                minDist = Math.min(minDist, x[i].distanceSquared(y[j]));
302                        }
303                        if (minDist < maxDistance * maxDistance) {
304                                contacts++;
305                        }
306                }
307                return contacts;
308        }
309
310}