diff --git a/benchmarks/index/chain/MonotoneChainBuilderPerfTest.cpp b/benchmarks/index/chain/MonotoneChainBuilderPerfTest.cpp
index 2b7ecd8a27..4218305e83 100644
--- a/benchmarks/index/chain/MonotoneChainBuilderPerfTest.cpp
+++ b/benchmarks/index/chain/MonotoneChainBuilderPerfTest.cpp
@@ -60,7 +60,7 @@ static void BM_MonotoneChainBuilder(benchmark::State& state) {
 
     for (auto _ : state) {
         std::vector<MonotoneChain> chains;
-        MonotoneChainBuilder::getChains(&cs, nullptr, chains);
+        MonotoneChainBuilder::getChains(&cs, nullptr, 0, chains);
     }
 }
 
diff --git a/benchmarks/index/chain/MonotoneChainPerfTest.cpp b/benchmarks/index/chain/MonotoneChainPerfTest.cpp
index 2cf4a34c05..be6dfc99f2 100644
--- a/benchmarks/index/chain/MonotoneChainPerfTest.cpp
+++ b/benchmarks/index/chain/MonotoneChainPerfTest.cpp
@@ -48,8 +48,8 @@ static void BM_MonotoneChainOverlaps(benchmark::State& state) {
         prev = c;
     }
 
-    MonotoneChain mc1(cs1, 0, cs1.size(), nullptr);
-    MonotoneChain mc2(cs2, 0, cs1.size(), nullptr);
+    MonotoneChain mc1(cs1, 0, cs1.size(), nullptr, 0.0);
+    MonotoneChain mc2(cs2, 0, cs1.size(), nullptr, 0.0);
 
     struct EmptyOverlapAction : public MonotoneChainOverlapAction {
         virtual void
diff --git a/include/geos/geom/prep/PreparedLineString.h b/include/geos/geom/prep/PreparedLineString.h
index 017e06e72f..3a9d2c4279 100644
--- a/include/geos/geom/prep/PreparedLineString.h
+++ b/include/geos/geom/prep/PreparedLineString.h
@@ -40,10 +40,13 @@ namespace prep { // geos::geom::prep
  */
 class PreparedLineString : public BasicPreparedGeometry {
 private:
-    std::unique_ptr<noding::FastSegmentSetIntersectionFinder> segIntFinder;
+    mutable std::unique_ptr<noding::FastSegmentSetIntersectionFinder> segIntFinder;
     mutable noding::SegmentString::ConstVect segStrings;
     mutable std::unique_ptr<operation::distance::IndexedFacetDistance> indexedDistance;
 
+    mutable std::once_flag segIntFinderFlag;
+    mutable std::once_flag indexedDistanceFlag;
+
 protected:
 public:
     PreparedLineString(const Geometry* geom)
@@ -54,7 +57,7 @@ class PreparedLineString : public BasicPreparedGeometry {
 
     ~PreparedLineString() override;
 
-    noding::FastSegmentSetIntersectionFinder* getIntersectionFinder();
+    noding::FastSegmentSetIntersectionFinder* getIntersectionFinder() const;
 
     bool intersects(const geom::Geometry* g) const override;
     std::unique_ptr<geom::CoordinateSequence> nearestPoints(const geom::Geometry* g) const override;
diff --git a/include/geos/geom/prep/PreparedLineStringIntersects.h b/include/geos/geom/prep/PreparedLineStringIntersects.h
index 026455df5f..980b0b0ebc 100644
--- a/include/geos/geom/prep/PreparedLineStringIntersects.h
+++ b/include/geos/geom/prep/PreparedLineStringIntersects.h
@@ -58,10 +58,7 @@ class PreparedLineStringIntersects {
         return op.intersects(geom);
     }
 
-    /**
-     * \todo FIXME - mloskot: Why not taking linestring through const reference?
-     */
-    PreparedLineStringIntersects(PreparedLineString& prep)
+    PreparedLineStringIntersects(const PreparedLineString& prep)
         : prepLine(prep)
     { }
 
@@ -74,7 +71,7 @@ class PreparedLineStringIntersects {
     bool intersects(const geom::Geometry* g) const;
 
 protected:
-    PreparedLineString& prepLine;
+    const PreparedLineString& prepLine;
 
     /**
      * Tests whether any representative point of the test Geometry intersects
diff --git a/include/geos/geom/prep/PreparedPolygon.h b/include/geos/geom/prep/PreparedPolygon.h
index a2feee6c91..f8239f4c14 100644
--- a/include/geos/geom/prep/PreparedPolygon.h
+++ b/include/geos/geom/prep/PreparedPolygon.h
@@ -50,13 +50,19 @@ namespace prep { // geos::geom::prep
  */
 class PreparedPolygon : public BasicPreparedGeometry {
 private:
-    bool isRectangle;
     mutable std::unique_ptr<noding::FastSegmentSetIntersectionFinder> segIntFinder;
     mutable std::unique_ptr<algorithm::locate::PointOnGeometryLocator> ptOnGeomLoc;
     mutable std::unique_ptr<algorithm::locate::PointOnGeometryLocator> indexedPtOnGeomLoc;
-    mutable noding::SegmentString::ConstVect segStrings;
     mutable std::unique_ptr<operation::distance::IndexedFacetDistance> indexedDistance;
 
+    mutable noding::SegmentString::ConstVect segStrings;
+
+    mutable std::once_flag segIntFinderFlag;
+    mutable std::once_flag ptOnGeomLocFlag;
+    mutable std::once_flag indexedPtOnGeomLocFlag;
+    mutable std::once_flag indexedDistanceFlag;
+    bool isRectangle;
+
 protected:
 public:
     PreparedPolygon(const geom::Geometry* geom);
diff --git a/include/geos/index/chain/MonotoneChain.h b/include/geos/index/chain/MonotoneChain.h
index 4b0dbb4275..a8e3fc9806 100644
--- a/include/geos/index/chain/MonotoneChain.h
+++ b/include/geos/index/chain/MonotoneChain.h
@@ -95,14 +95,15 @@ class GEOS_DLL MonotoneChain {
     /// @param context
     ///   Ownership left to caller, this class holds a reference.
     ///
+    /// @param expansionDistance distance by which envelope should be expanded
+    ///
     MonotoneChain(const geom::CoordinateSequence& pts,
-                  std::size_t start, std::size_t end, void* context);
+                  std::size_t start, std::size_t end, void* context, double expansionDistance);
 
     ~MonotoneChain() = default;
 
     /// Returned envelope is owned by this class
     const geom::Envelope& getEnvelope() const;
-    const geom::Envelope& getEnvelope(double expansionDistance) const;
 
     size_t
     getStartIndex() const
diff --git a/include/geos/index/chain/MonotoneChainBuilder.h b/include/geos/index/chain/MonotoneChainBuilder.h
index 815c39b3d1..7db4bfd4ae 100644
--- a/include/geos/index/chain/MonotoneChainBuilder.h
+++ b/include/geos/index/chain/MonotoneChainBuilder.h
@@ -58,6 +58,7 @@ class GEOS_DLL MonotoneChainBuilder {
      */
     static void getChains(const geom::CoordinateSequence* pts,
                           void* context,
+                          double expansionDistance,
                           std::vector<MonotoneChain>& mcList);
 
     /**
diff --git a/include/geos/index/strtree/TemplateSTRtree.h b/include/geos/index/strtree/TemplateSTRtree.h
index 32263289b1..75aab54073 100644
--- a/include/geos/index/strtree/TemplateSTRtree.h
+++ b/include/geos/index/strtree/TemplateSTRtree.h
@@ -231,7 +231,7 @@ class TemplateSTRtreeImpl {
 
     template<typename ItemDistance>
     ItemType nearestNeighbour(const BoundsType& env, const ItemType& item, ItemDistance& itemDist) {
-        build();
+        buildIfNeeded();
 
         if (getRoot() == nullptr) {
             return nullptr;
@@ -273,9 +273,7 @@ class TemplateSTRtreeImpl {
     // false values will be taken as a signal to stop the query.
     template<typename Visitor>
     void query(const BoundsType& queryEnv, Visitor &&visitor) {
-        if (!built()) {
-            build();
-        }
+        buildIfNeeded();
 
         if (root && root->boundsIntersect(queryEnv)) {
             if (root->isLeaf()) {
@@ -294,9 +292,7 @@ class TemplateSTRtreeImpl {
     // false values will be taken as a signal to stop the query.
     template<typename Visitor>
     void queryPairs(Visitor&& visitor) {
-        if (!built()) {
-            build();
-        }
+        buildIfNeeded();
 
         if (numItems < 2) {
             return;
@@ -318,7 +314,7 @@ class TemplateSTRtreeImpl {
      * Returns a depth-first iterator over all items in the tree.
      */
     Items items() {
-        build();
+        buildIfNeeded();
         return Items(*this);
     }
 
@@ -341,7 +337,7 @@ class TemplateSTRtreeImpl {
     /// @{
 
     bool remove(const BoundsType& itemEnv, const ItemType& item) {
-        build();
+        buildIfNeeded();
 
         if (root == nullptr) {
             return false;
@@ -369,16 +365,18 @@ class TemplateSTRtreeImpl {
 
     /** Determine whether the tree has been built, and no more items may be added. */
     const Node* getRoot() {
-        build();
+        buildIfNeeded();
         return root;
     }
 
     /// @}
 
     /** Build the tree if it has not already been built. */
-    void build() {
-        std::lock_guard<std::mutex> lock(lock_);
+    void buildIfNeeded() {
+        std::call_once(builtFlag, [this]() { build(); });
+    }
 
+    void build() {
         if (built()) {
             return;
         }
@@ -410,7 +408,7 @@ class TemplateSTRtreeImpl {
     }
 
 protected:
-    std::mutex lock_;
+    std::once_flag builtFlag;
     NodeList nodes;      //**< a list of all leaf and branch nodes in the tree. */
     Node* root;          //**< a pointer to the root node, if the tree has been built. */
     size_t nodeCapacity; //*< maximum number of children of each node */
diff --git a/include/geos/math/DD.h b/include/geos/math/DD.h
index 9f2ea979c7..59787acd21 100644
--- a/include/geos/math/DD.h
+++ b/include/geos/math/DD.h
@@ -92,6 +92,8 @@
 
 #pragma once
 
+#include <geos/export.h>
+
 #include <cmath>
 
 namespace geos {
@@ -115,36 +117,36 @@ class GEOS_DLL DD {
 
 
     public:
-        DD(double p_hi, double p_lo) : hi(p_hi), lo(p_lo) {};
-        DD(double x) : hi(x), lo(0.0) {};
-        DD() : hi(0.0), lo(0.0) {};
+        constexpr DD(double p_hi, double p_lo) : hi(p_hi), lo(p_lo) {};
+        constexpr DD(double x) : hi(x), lo(0.0) {};
+        constexpr DD() : hi(0.0), lo(0.0) {};
 
-        bool operator==(const DD &rhs) const
+        constexpr bool operator==(const DD &rhs) const
         {
             return hi == rhs.hi && lo == rhs.lo;
         }
 
-        bool operator!=(const DD &rhs) const
+        constexpr bool operator!=(const DD &rhs) const
         {
             return hi != rhs.hi || lo != rhs.lo;
         }
 
-        bool operator<(const DD &rhs) const
+        constexpr bool operator<(const DD &rhs) const
         {
             return (hi < rhs.hi) || (hi == rhs.hi && lo < rhs.lo);
         }
 
-        bool operator<=(const DD &rhs) const
+        constexpr bool operator<=(const DD &rhs) const
         {
             return (hi < rhs.hi) || (hi == rhs.hi && lo <= rhs.lo);
         }
 
-        bool operator>(const DD &rhs) const
+        constexpr bool operator>(const DD &rhs) const
         {
             return (hi > rhs.hi) || (hi == rhs.hi && lo > rhs.lo);
         }
 
-        bool operator>=(const DD &rhs) const
+        constexpr bool operator>=(const DD &rhs) const
         {
             return (hi > rhs.hi) || (hi == rhs.hi && lo >= rhs.lo);
         }
diff --git a/include/geos/noding/FastSegmentSetIntersectionFinder.h b/include/geos/noding/FastSegmentSetIntersectionFinder.h
index c77ad86852..b6ba5df5bb 100644
--- a/include/geos/noding/FastSegmentSetIntersectionFinder.h
+++ b/include/geos/noding/FastSegmentSetIntersectionFinder.h
@@ -49,8 +49,7 @@ namespace noding { // geos::noding
  */
 class FastSegmentSetIntersectionFinder {
 private:
-    std::unique_ptr<MCIndexSegmentSetMutualIntersector> segSetMutInt;
-    std::unique_ptr<geos::algorithm::LineIntersector> lineIntersector;
+    MCIndexSegmentSetMutualIntersector segSetMutInt;
 
 protected:
 public:
@@ -67,7 +66,7 @@ class FastSegmentSetIntersectionFinder {
     const SegmentSetMutualIntersector*
     getSegmentSetIntersector() const
     {
-        return segSetMutInt.get();
+        return &segSetMutInt;
     }
 
     bool intersects(SegmentString::ConstVect* segStrings);
diff --git a/include/geos/noding/MCIndexSegmentSetMutualIntersector.h b/include/geos/noding/MCIndexSegmentSetMutualIntersector.h
index abc687919d..d3c1ae48f9 100644
--- a/include/geos/noding/MCIndexSegmentSetMutualIntersector.h
+++ b/include/geos/noding/MCIndexSegmentSetMutualIntersector.h
@@ -29,8 +29,6 @@ class SegmentIntersector;
 }
 }
 
-//using namespace geos::index::strtree;
-
 namespace geos {
 namespace noding { // geos::noding
 
@@ -45,12 +43,7 @@ class MCIndexSegmentSetMutualIntersector : public SegmentSetMutualIntersector {
 public:
 
     MCIndexSegmentSetMutualIntersector(double p_tolerance)
-        : monoChains()
-        , indexCounter(0)
-        , processCounter(0)
-        , nOverlaps(0)
-        , overlapTolerance(p_tolerance)
-        , indexBuilt(false)
+        : overlapTolerance(p_tolerance)
     {}
 
     MCIndexSegmentSetMutualIntersector()
@@ -68,9 +61,10 @@ class MCIndexSegmentSetMutualIntersector : public SegmentSetMutualIntersector {
 
     void setBaseSegments(SegmentString::ConstVect* segStrings) override;
 
-    // NOTE: re-populates the MonotoneChain vector with newly created chains
     void process(SegmentString::ConstVect* segStrings) override;
 
+    void process(SegmentString::ConstVect* segStrings, SegmentIntersector* si);
+
     class SegmentOverlapAction : public index::chain::MonotoneChainOverlapAction {
     private:
         SegmentIntersector& si;
@@ -98,7 +92,6 @@ class MCIndexSegmentSetMutualIntersector : public SegmentSetMutualIntersector {
 private:
 
     typedef std::vector<index::chain::MonotoneChain> MonoChains;
-    MonoChains monoChains;
 
     /*
      * The index::SpatialIndex used should be something that supports
@@ -106,23 +99,18 @@ class MCIndexSegmentSetMutualIntersector : public SegmentSetMutualIntersector {
      * or index::strtree::STRtree).
      */
     index::strtree::TemplateSTRtree<const index::chain::MonotoneChain*> index;
-    int indexCounter;
-    int processCounter;
-    // statistics
-    int nOverlaps;
-    double overlapTolerance;
+
+    const double overlapTolerance;
 
     /* memory management helper, holds MonotoneChain objects used
      * in the SpatialIndex. It's cleared when the SpatialIndex is
      */
-    bool indexBuilt;
+    std::once_flag indexBuilt;
     MonoChains indexChains;
 
-    void addToIndex(SegmentString* segStr);
-
-    void intersectChains();
+    void intersectChains(const MonoChains& chains, SegmentIntersector& segmentIntersector);
 
-    void addToMonoChains(SegmentString* segStr);
+    void addChains(const SegmentString* segStr, MonoChains& chains) const;
 
 };
 
diff --git a/include/geos/noding/SegmentIntersectionDetector.h b/include/geos/noding/SegmentIntersectionDetector.h
index 2165becb45..a99e9b0d5f 100644
--- a/include/geos/noding/SegmentIntersectionDetector.h
+++ b/include/geos/noding/SegmentIntersectionDetector.h
@@ -40,7 +40,7 @@ namespace noding { // geos::noding
  */
 class SegmentIntersectionDetector : public SegmentIntersector {
 private:
-    algorithm::LineIntersector* li;
+    algorithm::LineIntersector li;
 
     bool findProper;
     bool findAllTypes;
@@ -54,9 +54,8 @@ class SegmentIntersectionDetector : public SegmentIntersector {
 
 protected:
 public:
-    SegmentIntersectionDetector(algorithm::LineIntersector* p_li)
+    SegmentIntersectionDetector()
         :
-        li(p_li),
         findProper(false),
         findAllTypes(false),
         _hasIntersection(false),
@@ -68,7 +67,6 @@ class SegmentIntersectionDetector : public SegmentIntersector {
 
     ~SegmentIntersectionDetector() override
     {
-        //delete intPt;
         delete intSegments;
     }
 
diff --git a/src/algorithm/CGAlgorithmsDD.cpp b/src/algorithm/CGAlgorithmsDD.cpp
index b58acfe656..e4d2503d8f 100644
--- a/src/algorithm/CGAlgorithmsDD.cpp
+++ b/src/algorithm/CGAlgorithmsDD.cpp
@@ -31,7 +31,8 @@ namespace {
 inline int
 OrientationDD(const DD &dd)
 {
-    static DD const zero(0.0);
+    constexpr const DD zero(0.0);
+
     if(dd < zero) {
         return CGAlgorithmsDD::RIGHT;
     }
diff --git a/src/algorithm/locate/IndexedPointInAreaLocator.cpp b/src/algorithm/locate/IndexedPointInAreaLocator.cpp
index 0382c5deb8..f445fcb579 100644
--- a/src/algorithm/locate/IndexedPointInAreaLocator.cpp
+++ b/src/algorithm/locate/IndexedPointInAreaLocator.cpp
@@ -100,15 +100,12 @@ IndexedPointInAreaLocator::buildIndex(const geom::Geometry& g)
 IndexedPointInAreaLocator::IndexedPointInAreaLocator(const geom::Geometry& g)
     :	areaGeom(g)
 {
+    buildIndex(areaGeom);
 }
 
 geom::Location
 IndexedPointInAreaLocator::locate(const geom::CoordinateXY* /*const*/ p)
 {
-    if (index == nullptr) {
-        buildIndex(areaGeom);
-    }
-
     algorithm::RayCrossingCounter rcc(*p);
 
     index->query(p->y, p->y, [&rcc](const SegmentView& ls) {
diff --git a/src/geom/prep/AbstractPreparedPolygonContains.cpp b/src/geom/prep/AbstractPreparedPolygonContains.cpp
index cfc87436cd..b2821653bf 100644
--- a/src/geom/prep/AbstractPreparedPolygonContains.cpp
+++ b/src/geom/prep/AbstractPreparedPolygonContains.cpp
@@ -83,9 +83,7 @@ AbstractPreparedPolygonContains::findAndClassifyIntersections(const geom::Geomet
     noding::SegmentString::ConstVect lineSegStr;
     noding::SegmentStringUtil::extractSegmentStrings(geom, lineSegStr);
 
-    algorithm::LineIntersector li;
-
-    noding::SegmentIntersectionDetector intDetector(&li);
+    noding::SegmentIntersectionDetector intDetector;
 
     intDetector.setFindAllIntersectionTypes(true);
     prepPoly->getIntersectionFinder()->intersects(&lineSegStr, &intDetector);
diff --git a/src/geom/prep/PreparedLineString.cpp b/src/geom/prep/PreparedLineString.cpp
index 819a42d6cb..2da0586137 100644
--- a/src/geom/prep/PreparedLineString.cpp
+++ b/src/geom/prep/PreparedLineString.cpp
@@ -27,6 +27,8 @@
 #include <geos/noding/FastSegmentSetIntersectionFinder.h>
 #include <geos/operation/distance/IndexedFacetDistance.h>
 
+#include <mutex>
+
 namespace geos {
 namespace geom { // geos.geom
 namespace prep { // geos.geom.prep
@@ -44,12 +46,12 @@ PreparedLineString::~PreparedLineString()
 }
 
 noding::FastSegmentSetIntersectionFinder*
-PreparedLineString::getIntersectionFinder()
+PreparedLineString::getIntersectionFinder() const
 {
-    if(! segIntFinder) {
+    std::call_once(segIntFinderFlag, [this]() {
         noding::SegmentStringUtil::extractSegmentStrings(&getGeometry(), segStrings);
-        segIntFinder.reset(new noding::FastSegmentSetIntersectionFinder(&segStrings));
-    }
+        segIntFinder = detail::make_unique<noding::FastSegmentSetIntersectionFinder>(&segStrings);
+    });
 
     return segIntFinder.get();
 }
@@ -68,12 +70,12 @@ PreparedLineString::intersects(const geom::Geometry* g) const
 
 /* public */
 operation::distance::IndexedFacetDistance*
-PreparedLineString::
-getIndexedFacetDistance() const
+PreparedLineString::getIndexedFacetDistance() const
 {
-    if(! indexedDistance ) {
-        indexedDistance.reset(new operation::distance::IndexedFacetDistance(&getGeometry()));
-    }
+    std::call_once(indexedDistanceFlag, [this]() {
+        indexedDistance = detail::make_unique<operation::distance::IndexedFacetDistance>(&getGeometry());
+    });
+
     return indexedDistance.get();
 }
 
diff --git a/src/geom/prep/PreparedPolygon.cpp b/src/geom/prep/PreparedPolygon.cpp
index 332c4798e3..874b3f0468 100644
--- a/src/geom/prep/PreparedPolygon.cpp
+++ b/src/geom/prep/PreparedPolygon.cpp
@@ -58,10 +58,11 @@ noding::FastSegmentSetIntersectionFinder*
 PreparedPolygon::
 getIntersectionFinder() const
 {
-    if(! segIntFinder) {
+    std::call_once(segIntFinderFlag, [this]() {
         noding::SegmentStringUtil::extractSegmentStrings(&getGeometry(), segStrings);
-        segIntFinder.reset(new noding::FastSegmentSetIntersectionFinder(&segStrings));
-    }
+        segIntFinder = detail::make_unique<noding::FastSegmentSetIntersectionFinder>(&segStrings);
+    });
+
     return segIntFinder.get();
 }
 
@@ -73,13 +74,21 @@ getPointLocator() const
     // instead of an IndexedPointInAreaLocator. There's a reasonable chance we will only use this locator
     // once (for example, if we get here through Geometry::intersects). So we create a simple locator for the
     // first usage and switch to an indexed locator when it is clear we're in a multiple-use scenario.
-    if(! ptOnGeomLoc) {
+    bool first_call = false;
+
+    std::call_once(ptOnGeomLocFlag, [this, &first_call]() {
         ptOnGeomLoc = detail::make_unique<algorithm::locate::SimplePointInAreaLocator>(&getGeometry());
+        first_call = true;
+    });
+
+    if (first_call) {
         return ptOnGeomLoc.get();
-    } else if (!indexedPtOnGeomLoc) {
-        indexedPtOnGeomLoc = detail::make_unique<algorithm::locate::IndexedPointInAreaLocator>(getGeometry());
     }
 
+    std::call_once(indexedPtOnGeomLocFlag, [this]() {
+        indexedPtOnGeomLoc = detail::make_unique<algorithm::locate::IndexedPointInAreaLocator>(getGeometry());
+    });
+
     return indexedPtOnGeomLoc.get();
 }
 
@@ -157,9 +166,10 @@ operation::distance::IndexedFacetDistance*
 PreparedPolygon::
 getIndexedFacetDistance() const
 {
-    if(! indexedDistance ) {
-        indexedDistance.reset(new operation::distance::IndexedFacetDistance(&getGeometry()));
-    }
+    std::call_once(indexedDistanceFlag, [this]() {
+        indexedDistance = detail::make_unique<operation::distance::IndexedFacetDistance>(&getGeometry());
+    });
+
     return indexedDistance.get();
 }
 
diff --git a/src/index/chain/MonotoneChain.cpp b/src/index/chain/MonotoneChain.cpp
index e3d336a333..2e3703bf73 100644
--- a/src/index/chain/MonotoneChain.cpp
+++ b/src/index/chain/MonotoneChain.cpp
@@ -32,29 +32,22 @@ namespace index { // geos.index
 namespace chain { // geos.index.chain
 
 MonotoneChain::MonotoneChain(const geom::CoordinateSequence& newPts,
-                             std::size_t nstart, std::size_t nend, void* nContext)
+                             std::size_t nstart, std::size_t nend, void* nContext,
+                             double expansionDistance)
     : pts(&newPts)
     , context(nContext)
     , start(nstart)
     , end(nend)
-    , env()
-{}
-
-const Envelope&
-MonotoneChain::getEnvelope() const
+    , env(pts->getAt<CoordinateXY>(start), pts->getAt<CoordinateXY>(end))
 {
-    return getEnvelope(0.0);
+    if (expansionDistance > 0.0) {
+        env.expandBy(expansionDistance);
+    }
 }
 
 const Envelope&
-MonotoneChain::getEnvelope(double expansionDistance) const
+MonotoneChain::getEnvelope() const
 {
-    if (env.isNull()) {
-        env.init(pts->getAt<CoordinateXY>(start), pts->getAt<CoordinateXY>(end));
-        if (expansionDistance > 0.0) {
-            env.expandBy(expansionDistance);
-        }
-    }
     return env;
 }
 
diff --git a/src/index/chain/MonotoneChainBuilder.cpp b/src/index/chain/MonotoneChainBuilder.cpp
index 9185fe83ab..44bd9a579f 100644
--- a/src/index/chain/MonotoneChainBuilder.cpp
+++ b/src/index/chain/MonotoneChainBuilder.cpp
@@ -48,13 +48,14 @@ namespace chain { // geos.index.chain
  */
 class ChainBuilder : public CoordinateFilter {
 public:
-    ChainBuilder(const CoordinateSequence* pts, void* context, std::vector<MonotoneChain> & list) :
+    ChainBuilder(const CoordinateSequence* pts, void* context, double expansionDistance, std::vector<MonotoneChain> & list) :
      m_prev(nullptr),
      m_i(0),
      m_quadrant(-1),
      m_start(0),
      m_seq(pts),
      m_context(context),
+     m_distance(expansionDistance),
      m_list(list) {}
 
     void filter_ro(const CoordinateXY* c) override {
@@ -72,7 +73,7 @@ class ChainBuilder : public CoordinateFilter {
     void finishChain() {
         if ( m_i == 0 ) return;
         std::size_t chainEnd = m_i - 1;
-        m_list.emplace_back(*m_seq, m_start, chainEnd, m_context);
+        m_list.emplace_back(*m_seq, m_start, chainEnd, m_context, m_distance);
         m_start = chainEnd;
     }
 
@@ -99,6 +100,7 @@ class ChainBuilder : public CoordinateFilter {
     std::size_t m_start;
     const CoordinateSequence* m_seq;
     void* m_context;
+    double m_distance;
     std::vector<MonotoneChain>& m_list;
 };
 
@@ -106,8 +108,9 @@ class ChainBuilder : public CoordinateFilter {
 /* static public */
 void
 MonotoneChainBuilder::getChains(const CoordinateSequence* pts, void* context,
+                                double expansionDistance,
                                 std::vector<MonotoneChain>& mcList) {
-    ChainBuilder builder(pts, context, mcList);
+    ChainBuilder builder(pts, context, expansionDistance, mcList);
     pts->apply_ro(&builder);
     builder.finish();
 }
diff --git a/src/noding/FastSegmentSetIntersectionFinder.cpp b/src/noding/FastSegmentSetIntersectionFinder.cpp
index c9f528922b..8960b4e9bb 100644
--- a/src/noding/FastSegmentSetIntersectionFinder.cpp
+++ b/src/noding/FastSegmentSetIntersectionFinder.cpp
@@ -39,17 +39,15 @@ namespace noding { // geos::noding
  */
 FastSegmentSetIntersectionFinder::
 FastSegmentSetIntersectionFinder(noding::SegmentString::ConstVect* baseSegStrings)
-    :	segSetMutInt(new MCIndexSegmentSetMutualIntersector()),
-      lineIntersector(new algorithm::LineIntersector())
 {
-    segSetMutInt->setBaseSegments(baseSegStrings);
+    segSetMutInt.setBaseSegments(baseSegStrings);
 }
 
 bool
 FastSegmentSetIntersectionFinder::
 intersects(noding::SegmentString::ConstVect* segStrings)
 {
-    SegmentIntersectionDetector intFinder(lineIntersector.get());
+    SegmentIntersectionDetector intFinder;
 
     return this->intersects(segStrings, &intFinder);
 }
@@ -59,8 +57,7 @@ FastSegmentSetIntersectionFinder::
 intersects(noding::SegmentString::ConstVect* segStrings,
            SegmentIntersectionDetector* intDetector)
 {
-    segSetMutInt->setSegmentIntersector(intDetector);
-    segSetMutInt->process(segStrings);
+    segSetMutInt.process(segStrings, intDetector);
 
     return intDetector->hasIntersection();
 }
diff --git a/src/noding/MCIndexNoder.cpp b/src/noding/MCIndexNoder.cpp
index 47ffaa4b45..d414d2cd90 100644
--- a/src/noding/MCIndexNoder.cpp
+++ b/src/noding/MCIndexNoder.cpp
@@ -51,7 +51,7 @@ MCIndexNoder::computeNodes(SegmentString::NonConstVect* inputSegStrings)
 
     if (!indexBuilt) {
         for(const auto& mc : monoChains) {
-            index.insert(mc.getEnvelope(overlapTolerance), &mc);
+            index.insert(mc.getEnvelope(), &mc);
         }
         indexBuilt = true;
     }
@@ -84,7 +84,7 @@ MCIndexNoder::add(SegmentString* segStr)
 
     // segChains will contain newly allocated MonotoneChain objects
     MonotoneChainBuilder::getChains(segStr->getCoordinates(),
-                                    segStr, monoChains);
+                                    segStr, overlapTolerance, monoChains);
 
 }
 
diff --git a/src/noding/MCIndexSegmentSetMutualIntersector.cpp b/src/noding/MCIndexSegmentSetMutualIntersector.cpp
index 148ace410c..1f043b2508 100644
--- a/src/noding/MCIndexSegmentSetMutualIntersector.cpp
+++ b/src/noding/MCIndexSegmentSetMutualIntersector.cpp
@@ -13,14 +13,11 @@
  **********************************************************************/
 
 #include <geos/noding/MCIndexSegmentSetMutualIntersector.h>
-#include <geos/noding/SegmentSetMutualIntersector.h>
 #include <geos/noding/SegmentString.h>
 #include <geos/noding/SegmentIntersector.h>
-#include <geos/index/SpatialIndex.h>
 #include <geos/index/chain/MonotoneChain.h>
 #include <geos/index/chain/MonotoneChainBuilder.h>
 #include <geos/index/chain/MonotoneChainOverlapAction.h>
-#include <geos/index/strtree/SimpleSTRtree.h>
 
 // std
 #include <cstddef>
@@ -33,37 +30,26 @@ namespace noding { // geos::noding
 
 /*private*/
 void
-MCIndexSegmentSetMutualIntersector::addToIndex(SegmentString* segStr)
-{
-    MonotoneChainBuilder::getChains(segStr->getCoordinates(),
-                                    segStr, indexChains);
-
-}
-
-
-/*private*/
-void
-MCIndexSegmentSetMutualIntersector::addToMonoChains(SegmentString* segStr)
+MCIndexSegmentSetMutualIntersector::addChains(const SegmentString* segStr, MonoChains& chains) const
 {
     if (segStr->size() == 0)
         return;
     MonotoneChainBuilder::getChains(segStr->getCoordinates(),
-                                    segStr, monoChains);
+                                    (void*) segStr, overlapTolerance, chains);
 }
 
 
 /*private*/
 void
-MCIndexSegmentSetMutualIntersector::intersectChains()
+MCIndexSegmentSetMutualIntersector::intersectChains(const MonoChains& chains, SegmentIntersector& segmentIntersector)
 {
-    MCIndexSegmentSetMutualIntersector::SegmentOverlapAction overlapAction(*segInt);
+    MCIndexSegmentSetMutualIntersector::SegmentOverlapAction overlapAction(segmentIntersector);
 
-    for(auto& queryChain : monoChains) {
-        index.query(queryChain.getEnvelope(overlapTolerance), [&queryChain, &overlapAction, this](const MonotoneChain* testChain) -> bool {
+    for(auto& queryChain : chains) {
+        index.query(queryChain.getEnvelope(), [&queryChain, &overlapAction, &segmentIntersector, this](const MonotoneChain* testChain) -> bool {
             queryChain.computeOverlaps(testChain, overlapTolerance, &overlapAction);
-            nOverlaps++;
 
-            return !segInt->isDone(); // abort early if segInt->isDone()
+            return !segmentIntersector.isDone(); // abort early if segmentIntersector->isDone()
         });
     }
 }
@@ -73,13 +59,8 @@ MCIndexSegmentSetMutualIntersector::intersectChains()
 void
 MCIndexSegmentSetMutualIntersector::setBaseSegments(SegmentString::ConstVect* segStrings)
 {
-    // NOTE - mloskot: const qualifier is removed silently, dirty.
-
     for(const SegmentString* css: *segStrings) {
-        if (css->size() == 0)
-            continue;
-        SegmentString* ss = const_cast<SegmentString*>(css);
-        addToIndex(ss);
+        addChains(css, indexChains);
     }
 }
 
@@ -87,23 +68,29 @@ MCIndexSegmentSetMutualIntersector::setBaseSegments(SegmentString::ConstVect* se
 void
 MCIndexSegmentSetMutualIntersector::process(SegmentString::ConstVect* segStrings)
 {
-    if (!indexBuilt) {
+    process(segStrings, segInt);
+}
+
+/*public*/
+void
+MCIndexSegmentSetMutualIntersector::process(SegmentString::ConstVect* segStrings,
+                                            SegmentIntersector* segmentIntersector)
+{
+    std::call_once(indexBuilt, [this]() {
         for (auto& mc: indexChains) {
-            index.insert(&(mc.getEnvelope(overlapTolerance)), &mc);
+            index.insert(&(mc.getEnvelope()), &mc);
         }
-        indexBuilt = true;
-    }
-
-    // Reset counters for new inputs
-    monoChains.clear();
-    processCounter = indexCounter + 1;
-    nOverlaps = 0;
+    });
 
+    // TODO: Rework MonotoneChain extraction to take a callback, so we can pass the chains
+    // to intersectChains as they are identified.
+    MonoChains monoChains;
     for(const SegmentString* css: *segStrings) {
         SegmentString* ss = const_cast<SegmentString*>(css);
-        addToMonoChains(ss);
+        addChains(ss, monoChains);
     }
-    intersectChains();
+
+    intersectChains(monoChains, *segmentIntersector);
 }
 
 
diff --git a/src/noding/SegmentIntersectionDetector.cpp b/src/noding/SegmentIntersectionDetector.cpp
index 97924018b1..c25e4ffabe 100644
--- a/src/noding/SegmentIntersectionDetector.cpp
+++ b/src/noding/SegmentIntersectionDetector.cpp
@@ -43,13 +43,13 @@ processIntersections(
     const CoordinateXY& p10 = e1->getCoordinate<CoordinateXY>( segIndex1 );
     const CoordinateXY& p11 = e1->getCoordinate<CoordinateXY>( segIndex1 + 1 );
 
-    li->computeIntersection(p00, p01, p10, p11);
+    li.computeIntersection(p00, p01, p10, p11);
 
-    if(li->hasIntersection()) {
+    if(li.hasIntersection()) {
         // record intersection info
         _hasIntersection = true;
 
-        bool isProper = li->isProper();
+        bool isProper = li.isProper();
 
         if(isProper) {
             _hasProperIntersection = true;
@@ -69,7 +69,7 @@ processIntersections(
 
         if(!intPt || saveLocation) {
             // record intersection location (approximate)
-            intPt = &li->getIntersection(0);
+            intPt = &li.getIntersection(0);
 
             delete intSegments;
 
diff --git a/src/operation/buffer/SegmentMCIndex.cpp b/src/operation/buffer/SegmentMCIndex.cpp
index 9441f5a9e0..0c2b5f8124 100644
--- a/src/operation/buffer/SegmentMCIndex.cpp
+++ b/src/operation/buffer/SegmentMCIndex.cpp
@@ -37,7 +37,7 @@ SegmentMCIndex::SegmentMCIndex(const CoordinateSequence* segs)
 void
 SegmentMCIndex::buildIndex(const CoordinateSequence* segs)
 {
-    chain::MonotoneChainBuilder::getChains(segs, nullptr, segChains);
+    chain::MonotoneChainBuilder::getChains(segs, nullptr, 0.0, segChains);
     for (chain::MonotoneChain& mc : segChains) {
         index.insert(&(mc.getEnvelope()), &mc);
     }
diff --git a/tests/unit/geom/prep/PreparedGeometryTest.cpp b/tests/unit/geom/prep/PreparedGeometryTest.cpp
index efc5c4b821..1d0b318d1d 100644
--- a/tests/unit/geom/prep/PreparedGeometryTest.cpp
+++ b/tests/unit/geom/prep/PreparedGeometryTest.cpp
@@ -12,6 +12,7 @@
 #include <geos/io/WKTReader.h>
 // std
 #include <memory>
+#include <thread>
 
 using namespace geos::geom;
 using geos::geom::prep::PreparedGeometry;
@@ -68,4 +69,58 @@ void object::test<1>
     ensure( pg1->covers(g2.get()));
 }
 
+// Check prepared geometry can be used from multiple threads
+template<>
+template<>
+void object::test<2>
+()
+{
+    std::vector<std::unique_ptr<Geometry>> geoms;
+    std::vector<std::unique_ptr<PreparedGeometry>> ppolys;
+    std::vector<std::unique_ptr<PreparedGeometry>> plines;
+    std::vector<std::thread> threads;
+
+    constexpr std::size_t nrow = 10;
+    constexpr std::size_t ncol = 10;
+    constexpr std::size_t nthreads = 10;
+
+    for (std::size_t i = 0; i < ncol; i++) {
+        for (std::size_t j = 0; j < nrow; j++) {
+            CoordinateXY c(static_cast<double>(i), static_cast<double>(j));
+            auto pt = factory->createPoint(c);
+
+            geoms.emplace_back(pt->buffer(1.5));
+            ppolys.push_back(prep::PreparedGeometryFactory::prepare(geoms.back().get()));
+            plines.push_back(prep::PreparedGeometryFactory::prepare(static_cast<Polygon*>(geoms.back().get())->getExteriorRing()));
+        }
+    }
+
+    auto findIntersects = [&geoms](const PreparedGeometry* pg) {
+        for (const auto& geom : geoms) {
+            pg->intersects(geom.get());
+            pg->distance(geom.get());
+        }
+    };
+
+    // check PreparedPolygon
+    for (std::size_t i = 0; i < nthreads; i++) {
+        threads.emplace_back(findIntersects, ppolys[0].get());
+    }
+
+    for (auto& thread : threads) {
+        thread.join();
+    }
+
+    // check PreparedLineString
+    threads.clear();
+    for (std::size_t i = 0; i < nthreads; i++) {
+        threads.emplace_back(findIntersects, plines[0].get());
+    }
+
+    for (auto& thread : threads) {
+        thread.join();
+    }
+
+}
+
 } // namespace tut
diff --git a/tests/unit/index/chain/MonotoneChainBuilderTest.cpp b/tests/unit/index/chain/MonotoneChainBuilderTest.cpp
index 50b615ae63..4fd9265b20 100644
--- a/tests/unit/index/chain/MonotoneChainBuilderTest.cpp
+++ b/tests/unit/index/chain/MonotoneChainBuilderTest.cpp
@@ -40,7 +40,7 @@ void object::test<1>
     std::vector<geos::index::chain::MonotoneChain> chain;
     geos::geom::CoordinateSequence pts;
 
-    geos::index::chain::MonotoneChainBuilder::getChains(&pts, 0, chain);
+    geos::index::chain::MonotoneChainBuilder::getChains(&pts, nullptr, 0, chain);
 
     ensure_equals(chain.size(), 0u);
 }
diff --git a/tests/unit/operation/relate/RelateOpTest.cpp b/tests/unit/operation/relate/RelateOpTest.cpp
index f65cd9a6fb..5d93997cf2 100644
--- a/tests/unit/operation/relate/RelateOpTest.cpp
+++ b/tests/unit/operation/relate/RelateOpTest.cpp
@@ -111,7 +111,7 @@ void object::test<3>()
     // Launch some threads to check relationships between polygons
     std::vector<std::thread> threads;
     for (std::size_t i = 0; i < numThreads; i++) {
-        threads.emplace_back(runRelate, geoms[i].get(), std::ref(geoms));
+        threads.emplace_back(runRelate, geoms[0].get(), std::ref(geoms));
     }
 
     // Wait for threads to complete