0
$\begingroup$

Is there a good strategy to speed up RegionIntersection computations like the one below:

ClearAll[reg1, reg2, res];
reg1 = Parallelepiped[
    {-223.812, -231.195, -3449.08},
    {
        {117.905, -236.219, -163.835},
        {-37.0885, -477.182, 661.318},
        {1056.68, 324.114, 293.13}
    }
];
reg2 = Polyhedron[{
    {
        {1231.48, 357.95, -4211.93},
        {-357.48, 135.948, -4812.15},
        {-702.666, 566.55, -4057.59}
    },
    {
        {-1159.38, -1559.96, -2061.98},
        {846.6, -756.866, -337.978},
        {853.676, -702.691, -322.517},
        {-1420.85, 64.8043, -1970.76}
    },
    {
        {1999.35, -2007.18, -2510.94},
        {3018.78, -199.677, -3076.07},
        {853.676, -702.691, -322.517},
        {846.6, -756.866, -337.978}
    },
    {
        {846.6, -756.866, -337.978},
        {-1159.38, -1559.96, -2061.98},
        {-518.838, -2359.01, -3462.16},
        {1999.35, -2007.18, -2510.94}
    },
    {
        {1231.48, 357.95, -4211.93},
        {-702.666, 566.55, -4057.59},
        {-1420.85, 64.8043, -1970.76},
        {853.676, -702.691, -322.517},
        {3018.78, -199.677, -3076.07}
    },
    {
        {-702.666, 566.55, -4057.59},
        {-357.48, 135.948, -4812.15},
        {-518.838, -2359.01, -3462.16},
        {-1159.38, -1559.96, -2061.98},
        {-1420.85, 64.8043, -1970.76}
    },
    {
        {3018.78, -199.677, -3076.07},
        {1999.35, -2007.18, -2510.94},
        {-518.838, -2359.01, -3462.16},
        {-357.48, 135.948, -4812.15},
        {1231.48, 357.95, -4211.93}
    }
}];
res = RegionIntersection[reg1, reg2];
(* Computes forever, before aborting *)

In my problem, reg2 is constant, and I need to compute its intersection with many smaller regions of the kind in reg1.

So far, the only strategy that I have found to speed up the computation--and, indeed, to return a result at all--is to compute axis-aligned BoundingRegions. But then the errors are much larger.

Windows 11, Mathematica 13.2, 12th generation i9, 64 GB of 4800 MHz DDR5 RAM (so the hardware should not really be the problem here)

How reg2 was obtained

I first obtained two PDB files from RCSB PDB. Here are the PDB files, uploaded to Google Drive for your convenience:

PD1

PDL1

Then I obtained the BoundingRegion for each of these and computed the intersection as follows:

ClearAll[pd1coords, pdl1coords, reg2];
pd1coords = Import[pd1File, {"PDB", "VertexCoordinates"}];
pdl1coords = Import[pdl1File, {"PDB", "VertexCoordinates"}];
reg2 = {pd1coords, pdl1coords} // (
    Map[BoundingRegion[#, "FastOrientedCuboid"] &] /*
    Apply[RegionIntersection]
);
$\endgroup$
5
  • $\begingroup$ I think "Polyhedron" is not a primitive for regions, Only "UniformPolyhedron" is. $\endgroup$ Jan 11, 2023 at 19:46
  • $\begingroup$ @DanielHuber, reg2 was actually obtained by intersecting two BoundingRegion[coords, "FastOrientedCuboid"] expressions. Is RegionIntersection not a closed operator? $\endgroup$ Jan 11, 2023 at 21:01
  • $\begingroup$ BoundaryDiscretizeGraphics[Graphics3D[reg2]] shows that the Polyhedron is not valid. Can you show how you obtained reg2? $\endgroup$
    – user21
    Jan 11, 2023 at 21:04
  • $\begingroup$ @user21 please see the edit above. $\endgroup$ Jan 11, 2023 at 21:25
  • $\begingroup$ I think this is a question Wolfram should answer. Report it to: [email protected] and let us know what they say. $\endgroup$ Jan 12, 2023 at 9:14

1 Answer 1

6
$\begingroup$

Edit

$Version
13.2.0 for Microsoft Windows (64-bit) (November 18, 2022)
reg3 = ConvexHullRegion[Flatten[reg2[[1]], 1]];
Manipulate[
 Show[Region[
   Style[RegionIntersection[
     TransformedRegion[reg1, TranslationTransform[t {1, 1, 1}]], 
     reg3], Red]], Graphics3D[Style[reg3, Opacity[.2]]]], {t, -1000, 
  600}]

enter image description here

Original

It seems that the problem come from reg2.

reg3 = ConvexHullRegion[Flatten[reg2[[1]], 1]]
intersection = RegionIntersection[reg1, reg3]
RegionMeasure[intersection]

enter image description here

$\endgroup$
7
  • $\begingroup$ Unfortunately, in my TimeConstrained run, this fails to return a result in 30 seconds. In fact, the kernel still aborts. How long did it take your computer to return a result? $\endgroup$ Jan 11, 2023 at 20:46
  • $\begingroup$ @Shredderroy Only 1 seconds. $\endgroup$
    – cvgmt
    Jan 11, 2023 at 21:13
  • $\begingroup$ Hmm...fascinating! I have restarted my kernel, but something seems to be off... $\endgroup$ Jan 11, 2023 at 21:25
  • $\begingroup$ @Shredderroy My code also work fine in v13.2 on Linux Mint 21.1. $\endgroup$
    – cvgmt
    Jan 11, 2023 at 21:45
  • $\begingroup$ @Shredderroy with 13.2.0 on macOS I managed to run the suggested solution by cvgmt very quickly $\endgroup$
    – bmf
    Jan 12, 2023 at 1:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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