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.contact; 022 023 024import java.util.ArrayList; 025import java.util.Arrays; 026import java.util.Collection; 027import java.util.List; 028 029import javax.vecmath.Point3d; 030 031import org.biojava.nbio.structure.Atom; 032import org.biojava.nbio.structure.Calc; 033 034 035/** 036 * A grid to be used for calculating atom contacts through a spatial hashing algorithm. 037 * <p> 038 * The grid is composed of cells of size of the cutoff so that the distances that need to be calculated 039 * are reduced to those within each cell and to the neighbouring cells. 040 * <p> 041 * Usage, for generic 3D points: 042 * <pre> 043 * Point3d[] points = ...; 044 * Grid grid = new Grid(8.0); 045 * grid.addCoords(points); 046 * List<Contact> contacts = getIndicesContacts(); 047 * </pre> 048 * Usage, for atoms: 049 * <pre> 050 * Atom[] atoms = ...; 051 * Grid grid = new Grid(8.0); 052 * grid.addCoords(atoms); 053 * AtomContactSet contacts = getAtomContacts(); 054 * </pre> 055 * 056 * @author Jose Duarte 057 * 058 */ 059public class Grid { 060 061 /** 062 * The scale: we use units of hundredths of Angstroms (thus cutoffs can be specified with a maximum precision of 0.01A) 063 */ 064 private static final int SCALE=100; 065 066 private GridCell[][][] cells; 067 068 private double cutoff; 069 private int cellSize; 070 071 private Point3d[] iAtoms; 072 private Point3d[] jAtoms; 073 074 private Atom[] iAtomObjects; 075 private Atom[] jAtomObjects; 076 077 // the bounds in int grid coordinates 078 private int[] bounds; 079 080 // the i and j bounding boxes in original double coordinates 081 private BoundingBox ibounds; 082 private BoundingBox jbounds; 083 084 private boolean noOverlap; // if the 2 sets of atoms are found not to overlap then this is set to true 085 086 /** 087 * Creates a <code>Grid</code>, the cutoff is in the same units as the coordinates 088 * (Angstroms if they are atom coordinates) and can 089 * be specified to a precision of 0.01. 090 * @param cutoff 091 */ 092 public Grid(double cutoff) { 093 this.cutoff = cutoff; 094 this.cellSize = (int) Math.floor(cutoff*SCALE); 095 this.noOverlap = false; 096 } 097 098 private int getFloor(double number) { 099 return (cellSize*((int)Math.floor(number*SCALE/cellSize))); 100 } 101 102 private int xintgrid2xgridindex(int xgridDim) { 103 return (xgridDim-bounds[0])/cellSize; 104 } 105 106 private int yintgrid2ygridindex(int ygridDim) { 107 return (ygridDim-bounds[1])/cellSize; 108 } 109 110 private int zintgrid2zgridindex(int zgridDim) { 111 return (zgridDim-bounds[2])/cellSize; 112 } 113 114 /** 115 * Adds the i and j atoms and fills the grid. Their bounds will be computed. 116 * Subsequent call to {@link #getIndicesContacts()} or {@link #getAtomContacts()} will produce the interatomic contacts. 117 * @param iAtoms 118 * @param jAtoms 119 */ 120 public void addAtoms(Atom[] iAtoms, Atom[] jAtoms) { 121 addAtoms(iAtoms, null, jAtoms, null); 122 } 123 124 /** 125 * Adds the i and j atoms and fills the grid, passing their bounds (array of size 6 with x,y,z minima and x,y,z maxima) 126 * This way the bounds don't need to be recomputed. 127 * Subsequent call to {@link #getIndicesContacts()} or {@link #getAtomContacts()} will produce the interatomic contacts. 128 * @param iAtoms 129 * @param icoordbounds 130 * @param jAtoms 131 * @param jcoordbounds 132 */ 133 public void addAtoms(Atom[] iAtoms, BoundingBox icoordbounds, Atom[] jAtoms, BoundingBox jcoordbounds) { 134 this.iAtoms = Calc.atomsToPoints(iAtoms); 135 this.iAtomObjects = iAtoms; 136 137 if (icoordbounds!=null) { 138 this.ibounds = icoordbounds; 139 } else { 140 this.ibounds = new BoundingBox(this.iAtoms); 141 } 142 143 this.jAtoms = Calc.atomsToPoints(jAtoms); 144 this.jAtomObjects = jAtoms; 145 146 if (jAtoms==iAtoms) { 147 this.jbounds=ibounds; 148 } else { 149 if (jcoordbounds!=null) { 150 this.jbounds = jcoordbounds; 151 } else { 152 this.jbounds = new BoundingBox(this.jAtoms); 153 154 } 155 } 156 157 fillGrid(); 158 } 159 160 /** 161 * Adds a set of atoms, subsequent call to {@link #getIndicesContacts()} or {@link #getAtomContacts()} will produce the interatomic contacts. 162 * The bounding box of the atoms will be computed based on input array 163 * @param atoms 164 */ 165 public void addAtoms(Atom[] atoms) { 166 addAtoms(atoms, (BoundingBox) null); 167 } 168 169 /** 170 * Adds a set of atoms, subsequent call to {@link #getIndicesContacts()} or {@link #getAtomContacts()} will produce the interatomic contacts. 171 * The bounds calculated elsewhere can be passed, or if null they are computed. 172 * @param atoms 173 * @param bounds 174 */ 175 public void addAtoms(Atom[] atoms, BoundingBox bounds) { 176 this.iAtoms = Calc.atomsToPoints(atoms); 177 this.iAtomObjects = atoms; 178 179 if (bounds!=null) { 180 this.ibounds = bounds; 181 } else { 182 this.ibounds = new BoundingBox(iAtoms); 183 } 184 185 this.jAtoms = null; 186 this.jAtomObjects = null; 187 this.jbounds = null; 188 189 fillGrid(); 190 } 191 192 /** 193 * Adds the i and j coordinates and fills the grid. Their bounds will be computed. 194 * Subsequent call to {@link #getIndicesContacts()} will produce the 195 * contacts, i.e. the set of points within distance cutoff. 196 * 197 * Subsequent calls to method {@link #getAtomContacts()} will produce a NullPointerException 198 * since this only adds coordinates and no atom information. 199 * @param iAtoms 200 * @param jAtoms 201 */ 202 public void addCoords(Point3d[] iAtoms, Point3d[] jAtoms) { 203 addCoords(iAtoms, null, jAtoms, null); 204 } 205 206 /** 207 * Adds the i and j coordinates and fills the grid, passing their bounds (array of size 6 with x,y,z minima and x,y,z maxima) 208 * This way the bounds don't need to be recomputed. 209 * Subsequent call to {@link #getIndicesContacts()} will produce the 210 * contacts, i.e. the set of points within distance cutoff. 211 * 212 * Subsequent calls to method {@link #getAtomContacts()} will produce a NullPointerException 213 * since this only adds coordinates and no atom information. 214 * @param iAtoms 215 * @param icoordbounds 216 * @param jAtoms 217 * @param jcoordbounds 218 */ 219 public void addCoords(Point3d[] iAtoms, BoundingBox icoordbounds, Point3d[] jAtoms, BoundingBox jcoordbounds) { 220 this.iAtoms = iAtoms; 221 this.iAtomObjects = null; 222 223 if (icoordbounds!=null) { 224 this.ibounds = icoordbounds; 225 } else { 226 this.ibounds = new BoundingBox(this.iAtoms); 227 } 228 229 this.jAtoms = jAtoms; 230 this.jAtomObjects = null; 231 232 if (jAtoms==iAtoms) { 233 this.jbounds=ibounds; 234 } else { 235 if (jcoordbounds!=null) { 236 this.jbounds = jcoordbounds; 237 } else { 238 this.jbounds = new BoundingBox(this.jAtoms); 239 240 } 241 } 242 243 fillGrid(); 244 } 245 246 /** 247 * Adds a set of coordinates, subsequent call to {@link #getIndicesContacts()} will produce the 248 * contacts, i.e. the set of points within distance cutoff. 249 * The bounding box of the atoms will be computed based on input array. 250 * Subsequent calls to method {@link #getAtomContacts()} will produce a NullPointerException 251 * since this only adds coordinates and no atom information. 252 * @param atoms 253 */ 254 public void addCoords(Point3d[] atoms) { 255 addCoords(atoms, (BoundingBox) null); 256 } 257 258 /** 259 * Adds a set of coordinates, subsequent call to {@link #getIndicesContacts()} will produce the 260 * contacts, i.e. the set of points within distance cutoff. 261 * The bounds calculated elsewhere can be passed, or if null they are computed. 262 * Subsequent calls to method {@link #getAtomContacts()} will produce a NullPointerException 263 * since this only adds coordinates and no atom information. 264 * @param atoms 265 * @param bounds 266 */ 267 public void addCoords(Point3d[] atoms, BoundingBox bounds) { 268 this.iAtoms = atoms; 269 this.iAtomObjects = null; 270 271 if (bounds!=null) { 272 this.ibounds = bounds; 273 } else { 274 this.ibounds = new BoundingBox(iAtoms); 275 } 276 277 this.jAtoms = null; 278 this.jAtomObjects = null; 279 this.jbounds = null; 280 281 fillGrid(); 282 } 283 284 /** 285 * Creates the grid based on the boundaries defined by all atoms given (iAtoms and jAtoms) 286 * and places the atoms in their corresponding grid cells. 287 * Checks also if the i and j grid overlap, i.e. the enclosing bounds of 288 * the 2 grids (i and j) are no more than one cell size apart. If they don't 289 * overlap then they are too far apart so there's nothing to calculate, we set 290 * the noOverlap flag and then {@link #getIndicesContacts()} will do no calculation at all. 291 */ 292 private void fillGrid() { 293 294 if (jbounds!=null && !ibounds.overlaps(jbounds, cutoff)) { 295 //System.out.print("-"); 296 noOverlap = true; 297 return; 298 } 299 300 findFullGridIntBounds(); 301 302 cells = new GridCell[1+(bounds[3]-bounds[0])/cellSize] 303 [1+(bounds[4]-bounds[1])/cellSize] 304 [1+(bounds[5]-bounds[2])/cellSize]; 305 306 int i = 0; 307 for (Point3d atom:iAtoms) { 308 309 int xind = xintgrid2xgridindex(getFloor(atom.x)); 310 int yind = yintgrid2ygridindex(getFloor(atom.y)); 311 int zind = zintgrid2zgridindex(getFloor(atom.z)); 312 if (cells[xind][yind][zind]==null) { 313 cells[xind][yind][zind] = new GridCell(this); 314 } 315 cells[xind][yind][zind].addIindex(i); 316 i++; 317 } 318 319 if (jAtoms==null) return; 320 321 int j = 0; 322 for (Point3d atom:jAtoms) { 323 324 int xind = xintgrid2xgridindex(getFloor(atom.x)); 325 int yind = yintgrid2ygridindex(getFloor(atom.y)); 326 int zind = zintgrid2zgridindex(getFloor(atom.z)); 327 if (cells[xind][yind][zind]==null) { 328 cells[xind][yind][zind] = new GridCell(this); 329 } 330 cells[xind][yind][zind].addJindex(j); 331 j++; 332 } 333 334 } 335 336 /** 337 * Calculates an int array of size 6 into member variable bounds: 338 * - elements 0,1,2: minimum x,y,z of the iAtoms and jAtoms 339 * - elements 3,4,5: maximum x,y,z of the iAtoms and jAtoms 340 */ 341 private void findFullGridIntBounds() { 342 int[] iIntBounds = getIntBounds(ibounds); 343 344 bounds = new int[6]; 345 if (jbounds==null) { 346 bounds = iIntBounds; 347 } else { 348 int[] jIntBounds = getIntBounds(jbounds); 349 bounds[0] = Math.min(iIntBounds[0],jIntBounds[0]); 350 bounds[1] = Math.min(iIntBounds[1],jIntBounds[1]); 351 bounds[2] = Math.min(iIntBounds[2],jIntBounds[2]); 352 bounds[3] = Math.max(iIntBounds[3],jIntBounds[3]); 353 bounds[4] = Math.max(iIntBounds[4],jIntBounds[4]); 354 bounds[5] = Math.max(iIntBounds[5],jIntBounds[5]); 355 } 356 } 357 358 /** 359 * Returns an int array of size 6 : 360 * - elements 0,1,2: minimum x,y,z (in grid int coordinates) of the given atoms 361 * - elements 3,4,5: maximum x,y,z (in grid int coordinates) of the given atoms 362 * @return 363 */ 364 private int[] getIntBounds(BoundingBox coordbounds) { 365 int[] bs = new int[6]; 366 bs[0] = getFloor(coordbounds.xmin); 367 bs[1] = getFloor(coordbounds.ymin); 368 bs[2] = getFloor(coordbounds.zmin); 369 bs[3] = getFloor(coordbounds.xmax); 370 bs[4] = getFloor(coordbounds.ymax); 371 bs[5] = getFloor(coordbounds.zmax); 372 return bs; 373 } 374 375 /** 376 * Returns all contacts, i.e. all atoms that are within the cutoff distance. 377 * If both iAtoms and jAtoms are defined then contacts are between iAtoms and jAtoms, 378 * if jAtoms is null, then contacts are within the iAtoms. 379 * @return 380 */ 381 public AtomContactSet getAtomContacts() { 382 383 AtomContactSet contacts = new AtomContactSet(cutoff); 384 385 List<Contact> list = getIndicesContacts(); 386 387 if (jAtomObjects == null) { 388 for (Contact cont : list) { 389 contacts.add(new AtomContact(new Pair<Atom>(iAtomObjects[cont.getI()],iAtomObjects[cont.getJ()]),cont.getDistance())); 390 } 391 392 } else { 393 for (Contact cont : list) { 394 contacts.add(new AtomContact(new Pair<Atom>(iAtomObjects[cont.getI()],jAtomObjects[cont.getJ()]),cont.getDistance())); 395 } 396 } 397 398 return contacts; 399 } 400 401 /** 402 * Returns all contacts, i.e. all atoms that are within the cutoff distance. 403 * If both iAtoms and jAtoms are defined then contacts are between iAtoms and jAtoms, 404 * if jAtoms is null, then contacts are within the iAtoms. 405 * @return 406 * @deprecated use {@link #getAtomContacts()} instead 407 */ 408 @Deprecated 409 public AtomContactSet getContacts() { 410 return getAtomContacts(); 411 } 412 413 /** 414 * Returns all contacts, i.e. all atoms that are within the cutoff distance, as simple Contact objects containing the atom indices pairs and the distance. 415 * If both iAtoms and jAtoms are defined then contacts are between iAtoms and jAtoms, 416 * if jAtoms is null, then contacts are within the iAtoms. 417 * @return 418 */ 419 public List<Contact> getIndicesContacts() { 420 421 List<Contact> list = new ArrayList<>(); 422 423 // if the 2 sets of atoms are not overlapping they are too far away and no need to calculate anything 424 // this won't apply if there's only one set of atoms (iAtoms), where we would want all-to-all contacts 425 if (noOverlap) return list; 426 427 428 for (int xind=0;xind<cells.length;xind++) { 429 for (int yind=0;yind<cells[xind].length;yind++) { 430 for (int zind=0;zind<cells[xind][yind].length;zind++) { 431 // distances of points within this cell 432 GridCell thisCell = cells[xind][yind][zind]; 433 if (thisCell==null) continue; 434 435 list.addAll(thisCell.getContactsWithinCell()); 436 437 // distances of points from this box to all neighbouring boxes: 26 iterations (26 neighbouring boxes) 438 for (int x=xind-1;x<=xind+1;x++) { 439 for (int y=yind-1;y<=yind+1;y++) { 440 for (int z=zind-1;z<=zind+1;z++) { 441 if (x==xind && y==yind && z==zind) continue; 442 443 if (x>=0 && x<cells.length && y>=0 && y<cells[x].length && z>=0 && z<cells[x][y].length) { 444 if (cells[x][y][z] == null) continue; 445 446 list.addAll(thisCell.getContactsToOtherCell(cells[x][y][z])); 447 } 448 } 449 } 450 } 451 } 452 } 453 } 454 455 return list; 456 } 457 458 /** 459 * Fast determination of whether any atoms from a given set fall within 460 * the cutoff of iAtoms. If {@link #addAtoms(Atom[], Atom[])} was called 461 * with two sets of atoms, contacts to either set are considered. 462 * @param atoms 463 * @return 464 */ 465 public boolean hasAnyContact(Point3d[] atoms) { 466 return hasAnyContact(Arrays.asList(atoms)); 467 } 468 public boolean hasAnyContact(Collection<Point3d> atoms) { 469 for(Point3d atom : atoms) { 470 // Calculate Grid cell for the atom 471 int xind = xintgrid2xgridindex(getFloor(atom.x)); 472 int yind = yintgrid2ygridindex(getFloor(atom.y)); 473 int zind = zintgrid2zgridindex(getFloor(atom.z)); 474 475 // Consider 3x3x3 grid of cells around point 476 for (int x=xind-1;x<=xind+1;x++) { 477 if( x<0 || cells.length<=x) continue; 478 for (int y=yind-1;y<=yind+1;y++) { 479 if( y<0 || cells[x].length<=y ) continue; 480 for (int z=zind-1;z<=zind+1;z++) { 481 if( z<0 || cells[x][y].length<=z ) continue; 482 483 GridCell cell = cells[x][y][z]; 484 // Check for contacts in this cell 485 if(cell != null && cell.hasContactToAtom(iAtoms, jAtoms, atom, cutoff)) { 486 return true; 487 } 488 } 489 } 490 } 491 } 492 return false; 493 } 494 495 public double getCutoff() { 496 return cutoff; 497 } 498 499 /** 500 * Tells whether (after having added atoms to grid) the i and j grids are not overlapping. 501 * Overlap is defined as enclosing bounds of the 2 grids being no more than one cell size apart. 502 * @return true if the 2 grids don't overlap, false if they do 503 */ 504 public boolean isNoOverlap() { 505 return noOverlap; 506 } 507 508 protected Point3d[] getIAtoms() { 509 return iAtoms; 510 } 511 512 protected Point3d[] getJAtoms() { 513 return jAtoms; 514 } 515 516}