Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a set of 3-space coordinates for the atoms of a molecule (I could also transform them into spheres with radii corresponding to the atoms they represent). I would like to place this molecule into the tightest possible 3D rectangular bounding box and determine the coordinates for the box vertices. Is there a graphics processing routine in Mathematica that will do this automatically for me, or do I need to implement an algorithm from the literature?

To be more specific - I would like to obtain the minimum bounding box for my molecule allowing rotation etc. I was able to find a solution for the special case of the molecule I'm interested in using a sort of lame trick with a rotational symmetry. I'm asking this question because I'd like to have a general-case solution.

One way to proceed would be for me to perform a set of random rotations, calculate the volume of the bounding box Mathematica generates, and keep going until I obtain a sufficiently tight fit. Is there a way to get the coordinates for the Graphics3D bounding box vertices?

share|improve this question
1  
Do you happen to have a representative example (coordinates and radii)? –  Daniel Lichtblau Mar 6 '13 at 23:46
    
@DanielLichtblau The number of coordinates is too large to post here, but basically I'm extracting atomistic coordinates from a PDB file like the following: rcsb.org/pdb/explore/explore.do?structureId=2UZ3 –  Roger Harris Mar 7 '13 at 5:30

2 Answers 2

up vote 1 down vote accepted

Here is an idea for an approximate method. Center the data, compute the singular value decomposition, and use the right factor rotation matrix to align the singular values axes with the coordinate axes (I'm not saying that very well but it gives the rotation matrix we want). Now take min and max values along coordinate axes. This gives the dimensions of what could be a fairly tight bounding box. Also the mean and rotation matrix provide a way to go back to the original coordinates via geometric transformation.

Here is an example. It may work better than your actual ones though, because it uses "nice" random values.

SeedRandom[11111];
dmat = DiagonalMatrix[{4, 7, 13}];
rmat = RandomReal[{-1, 1}, {3, 3}];
mat = dmat + rmat + Transpose[rmat]

(* Out[115]= {{3.05494151932, 
  1.22535650049, -0.103665150251}, {1.22535650049, 5.67393092991, 
  1.38468670897}, {-0.103665150251, 1.38468670897, 12.2854081966}} *)

One easily checks that the eigenvalues are all positive, hence this can be used as correlation matrix for generating random normals in 3D.

vals = RandomVariate[MultinormalDistribution[{0, 0, 0}, mat], 10^4];

We'll have a look.

ListPointPlot3D[vals]

enter image description here

Now find the "best" rotation matrix via SVD. I remark that this should be done differently if one is on a 32 bit or otherwise memory-limited machine, since it will take substantial memory if done as below.

means = Mean[vals];
newvals = Map[# - means &, vals];
{uu, ww, vv} = SingularValueDecomposition[newvals];

Now we can rotate these centered points.

rotatedvals = newvals.vv;
ListPointPlot3D[rotatedvals]

enter image description here

Finally we find max and min values for a plausible bounding box.

mins = Map[Min[rotatedvals[[All, #]]] &, Range[3]]
maxes = Map[Max[rotatedvals[[All, #]]] &, Range[3]]

(* Out[125]= {-12.6346155672, -9.21373247283, -6.11807437861}

Out[126]= {13.2151960928, 9.81786657885, 5.87492886205} *)
share|improve this answer
    
I like this answer a lot. I will also eventually post an implementation of O'Rourke's bounding box algorithm. –  Roger Harris Mar 7 '13 at 22:45

Assuming you're not trying to rotate the molecule to minimize some parameter (volume? diagonal?) of the bounding box, take a set of positions and radii:

pos = RandomReal[{0, 10}, {10, 3}];
rad = RandomReal[{1, 2}, {10}];

Add and subtract a radius from each position and dimension, then take the bounds:

box = {Min[#], Max[#]}& /@ Transpose[Join @@ MapThread[{#1 + #2, #1 - #2}&, {pos, rad}]];

Here it is:

Graphics3D[{Gray, Sphere[pos, rad]}, 
 PlotRange -> box, PlotRangePadding -> None, 
 Lighting -> "Neutral"]

Molecule in a box

share|improve this answer
    
Rotating the molecule to obtain the minimum bounding box is precisely what I'm trying to do... unfortunately. –  Roger Harris Mar 7 '13 at 5:30
    
That is a much harder problem. The solution in C++ is at Containment; one in Python is at Minimum Bounding Box. I don't believe anything is available in Mathematica. I'm afraid I haven't the time to cook it up just now. Note that these are actually for points; you'd probably want to make a sampling of the sphere surface for points belonging to the hull and then re-hull to use the algorithm for spheres. –  Xerxes Mar 7 '13 at 6:25
    
Joseph O'Rourke (of MathOverflow fame) has a 1985 paper on an $O(n^3)$ algorithm: "Finding Minimal Enclosing Boxes" <cs.smith.edu/~orourke/Papers/MinVolBox.pdf>;. –  Roger Harris Mar 7 '13 at 6:51
1  
This may also be of interest to you: Efficiently Approximating the Minimum-Volume Bounding Box of a Point Set in Three Dimensions. The standard algorithm runs in $O(n^3)$, but this heuristic runs in $O(n+\epsilon^{-3})$, where $1+\epsilon$ is how much worse the resulting box is than the true minimum. Might be good for very large sets of points. –  Xerxes Mar 7 '13 at 6:51
    
Ok, well, I guess there have been improvements. =) –  Roger Harris Mar 7 '13 at 6:52

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.