Efficiently check if points are inside a surface.
| 740 | |
| 741 | |
| 742 | class _CheckInside: |
| 743 | """Efficiently check if points are inside a surface.""" |
| 744 | |
| 745 | @verbose |
| 746 | def __init__(self, surf, *, mode="old", verbose=None): |
| 747 | assert mode in ("pyvista", "old") |
| 748 | self.mode = mode |
| 749 | t0 = time.time() |
| 750 | self.surf = surf |
| 751 | if self.mode == "pyvista": |
| 752 | self._init_pyvista() |
| 753 | else: |
| 754 | self._init_old() |
| 755 | logger.debug( |
| 756 | f"Setting up {mode} interior check for {len(self.surf['rr'])} " |
| 757 | f"points took {(time.time() - t0) * 1000:0.1f} ms" |
| 758 | ) |
| 759 | |
| 760 | def _init_old(self): |
| 761 | self.inner_r = None |
| 762 | self.center = self.surf["rr"].mean(0) |
| 763 | # We could use Delaunay or ConvexHull here, Delaunay is slightly slower |
| 764 | # to construct but faster to evaluate |
| 765 | # See https://stackoverflow.com/questions/16750618/whats-an-efficient-way-to-find-if-a-point-lies-in-the-convex-hull-of-a-point-cl # noqa |
| 766 | self.del_tri = Delaunay(self.surf["rr"]) |
| 767 | if self.del_tri.find_simplex(self.center) >= 0: |
| 768 | # Immediately cull some points from the checks |
| 769 | dists = np.linalg.norm(self.surf["rr"] - self.center, axis=-1) |
| 770 | self.inner_r = dists.min() |
| 771 | self.outer_r = dists.max() |
| 772 | |
| 773 | def _init_pyvista(self): |
| 774 | if not isinstance(self.surf, dict): |
| 775 | self.pdata = self.surf |
| 776 | self.surf = _polydata_to_surface(self.pdata) |
| 777 | else: |
| 778 | self.pdata = _surface_to_polydata(self.surf).clean() |
| 779 | |
| 780 | @verbose |
| 781 | def __call__(self, rr, *, n_jobs=None, verbose=None): |
| 782 | n_orig = len(rr) |
| 783 | logger.info( |
| 784 | f"Checking surface interior status for {n_orig} point{_pl(n_orig, ' ')}..." |
| 785 | ) |
| 786 | t0 = time.time() |
| 787 | if self.mode == "pyvista": |
| 788 | inside = self._call_pyvista(rr) |
| 789 | else: |
| 790 | inside = self._call_old(rr, n_jobs) |
| 791 | n = inside.sum() |
| 792 | logger.info(f" Total {n}/{n_orig} point{_pl(n, ' ')} inside the surface") |
| 793 | logger.info(f"Interior check completed in {(time.time() - t0) * 1000:0.1f} ms") |
| 794 | return inside |
| 795 | |
| 796 | def query(self, rr): |
| 797 | """Get the distance to the nearest point.""" |
| 798 | if not hasattr(self, "_tree"): # compute on the fly only when needed |
| 799 | from scipy.spatial import KDTree |
no outgoing calls
no test coverage detected