John Collins

John Collins

Developer

© 2024

The Perfect Loft...

Doesn’t exist. Any seasoned veteran in solid modeling world will tell you loft is hard. Be careful with loft. This is probably true but I’ll tell you even dirt simple loft algorithms are incredibly useful.

Rather than try to implement the perfect Do-What-I-Mean loft algorithm, I’ve decided to start simple and capture the 80% use case. The simplest algorithm simply connects vertices one-to-one between positioned cross sections. This is the full implementation:


manifold::Manifold IsomorphicLoft(const std::vector<manifold::Polygons>& sections, const std::vector<glm::mat4x3>& transforms) {
    std::vector<glm::vec3> vertPos;
    std::vector<glm::ivec3> triVerts;

    if (sections.size() != transforms.size()) {
        throw std::runtime_error("Mismatched number of sections and transforms");
    }

    std::size_t offset = 0;
    std::size_t nVerticesInEachSection = 0;

    for (std::size_t i = 0; i < sections.size(); ++i) {
        const manifold::Polygons polygons = sections[i];
        glm::mat4x3 transform = transforms[i];

        for (const auto& polygon : polygons) {
            for (const glm::vec2& vertex : polygon) {
                glm::vec3 translatedVertex = MatrixTransforms::Translate(transform, glm::vec3(vertex.x, vertex.y, 0))[3];
                vertPos.push_back(translatedVertex);
            }
        }

        if (i == 0) {
            nVerticesInEachSection = vertPos.size();
        } else if ((vertPos.size() % nVerticesInEachSection) != 0)  {
            throw std::runtime_error("Recieved CrossSection with different number of vertices");
        }

        if (i < sections.size() - 1) {
            std::size_t currentOffset = offset;
            std::size_t nextOffset = offset + nVerticesInEachSection;

            for (std::size_t j = 0; j < polygons.size(); ++j) {
                const auto& polygon = polygons[j];

                for (std::size_t k = 0; k < polygon.size(); ++k) {
                    std::size_t nextIndex = (k + 1) % polygon.size();

                    glm::ivec3 triangle1(currentOffset + k, currentOffset + nextIndex, nextOffset + k);
                    glm::ivec3 triangle2(currentOffset + nextIndex, nextOffset + nextIndex, nextOffset + k);

                    triVerts.push_back(triangle1);
                    triVerts.push_back(triangle2);
                }
                currentOffset += polygon.size();
                nextOffset += polygon.size();
            }
        }

        offset += nVerticesInEachSection;
    }

    auto frontPolygons = sections.front();
    auto frontTriangles = manifold::Triangulate(frontPolygons, -1.0);
    for (auto& tri : frontTriangles) {
        triVerts.push_back({tri.z, tri.y, tri.x});
    }

    auto backPolygons = sections.back();
    auto backTriangles = manifold::Triangulate(backPolygons, -1.0);
    for (auto& triangle : backTriangles) {
        triangle.x += offset - nVerticesInEachSection;
        triangle.y += offset - nVerticesInEachSection;
        triangle.z += offset - nVerticesInEachSection;
        triVerts.push_back(triangle);
    }

    manifold::Mesh mesh;
    mesh.triVerts = triVerts;
    mesh.vertPos = vertPos;
    return manifold::Manifold(mesh);
}

Even this simple implementation is so useful I sometimes winder how I used to model without it.

Many “f(z)” type constructions become much simpler. For example, this main body segment for a hydroponic tower:

(ns spiralized-hydroponic-tower.core
  (:refer-clojure :exclude [set])
  (:require
   [plexus.utils :as pu]
   [spiralized-hydroponic-tower.utils :as u]
   [spiralized-hydroponic-tower.math :refer [pi|2 pi|3 pi|4 pi|5 pi|6 two-pi
                             cos acos sin asin tan atan pi
                             TT T T|2 T|3 T|4 T|5 T|6 sqrt
                             vec-sub vec-normalize vec-scale]]
   [plexus.core :as p]
   [clj-manifold3d.core :as m]))

(def cup-holder-inner-radius 53/2)
(def base-extra-radius 25)
(def tower-wall-thickness 0.85)
(def tower-segment-radius 57.8)
(def tower-segment-amplitude 12.7)
(def spacer-segment-height 208)
(def spacer-segment-offset (+ 2 1.2 base-extra-radius))
(def spacer-segment-angle (atan (/ spacer-segment-height spacer-segment-offset)))
(def spacer-segment-length (/ (abs spacer-segment-offset) (cos spacer-segment-angle)))
(def top-joint--height 10)
(def top-joint-offset 0.5)
(def top-joint-angle (atan (/ 0.5 10)))

(defn housing-shape
  "Sum a circle function with four revolutions of a sin function."
  [circle-radius sin-wave-amplitude wall-thickness ratio]
  (let [section
        (m/cross-section
         (let [angle (* 2 pi)
               n-steps 120
               step-size (/ angle n-steps)]
           (for [step (range (inc n-steps))
                 :let [x (* step step-size)
                       a (+ (+ circle-radius wall-thickness)
                            (* sin-wave-amplitude ratio (+ 1 (Math/sin (- (* 3 (- x (/ pi 2))) (* 3/3 pi)))))
                            (* sin-wave-amplitude (- 1 ratio) (+ 1 (Math/sin (* 3 (- x (/ pi 2)))))))]]
             [(* a (cos x))
              (* a (sin x))])))]
    (m/difference section (m/offset section -1))))

(def main-body-vertical-contour
  (p/points
   :axes [:x :y]
   :meta-props [:segment]
   (p/frame :name :origin :meta {:segment :bottom})
   (p/rotate :x (- pi|2))
   (p/translate :x tower-segment-radius)
   (p/forward :length 160 :n-steps 55)
   (p/set-meta :segment :top)
   (u/curve-segment :delta 0.8
                    :height 3
                    :cs 4)
   (p/forward :length 5)
   (p/rotate :y top-joint-angle)
   (p/forward :length 10)))

(let [[bottom top] (partition-by (fn [x] (-> x meta :segment))
                                  main-body-vertical-contour)]
  (def tower-segment-bottom-contour bottom)
  (def tower-segment-top-contour top))

(def hydro-tower-main-body
  (m/loft
   (concat
    (let [[radius amplitude] [(- tower-segment-radius 1) 12.5]]
      [{:cross-section (housing-shape radius amplitude 0 0)
        :frame (m/frame 1)}])
    (for [[i [radius z]] (take-nth 1 (map-indexed list tower-segment-bottom-contour))]
      (let [ratio (/ i (dec (count tower-segment-bottom-contour)))]
        {:cross-section (housing-shape radius tower-segment-amplitude 0 ratio)
         :frame (m/translate (m/frame 1) [0 0 (+ 5 z)])}))
    (for [[_ [radius z]] (map-indexed list tower-segment-top-contour)]
      (let [ratio 1]
        {:cross-section (housing-shape radius tower-segment-amplitude 0 ratio)
         :frame (m/translate (m/frame 1) [0 0 (+ 5 z)])})))))

Tower Body

Lofting isomorphic cross sections is great, but of course you end up wanting a bit more. The next simplest loft algorithm simply constructs edges by eagerly adding the edge of minimum distance as it zips around consecutive polygons. This one handles one-to-many and many-to-one cross section transitions no problem. It’s actually pretty hard to trip the algo up and generate non-sense meshes.

manifold::Manifold EagerNearestNeighborLoft(const std::vector<manifold::Polygons>& sections, const std::vector<glm::mat4x3>& transforms) {
    if (sections.size() != transforms.size()) {
      throw std::runtime_error("Mismatched number of sections and transforms");
    }
    if (sections.size() < 2) {
      throw std::runtime_error("Loft requires at least two sections.");
    }

    size_t nVerts = 0;
    std::vector<size_t> sectionSizes;
    sectionSizes.reserve(sections.size());
    for (auto& section: sections) {
        size_t sectionSize = 0;
        for (auto& poly: section) {
            sectionSize += poly.size();
        }
        nVerts + sectionSize;
        sectionSizes.push_back(sectionSize);
    }

    std::vector<glm::vec3> vertPos;
    vertPos.reserve(nVerts);
    std::vector<glm::ivec3> triVerts;
    triVerts.reserve(2*nVerts);

    size_t botSectionOffset = 0;
    for (std::size_t i = 0; i < sections.size() - 1; ++i) {
        const manifold::Polygons& botPolygons = sections[i];
        const manifold::Polygons& topPolygons = sections[i + 1];
        const glm::mat4x3& botTransform = transforms[i];

        if (botPolygons.size() != topPolygons.size()) {
          throw std::runtime_error("Cross sections must be composed of euqal number of polygons.");
        }

        size_t botSectionSize = sectionSizes[i];
        size_t topSectionOffset = botSectionOffset + botSectionSize;

        size_t botPolyOffset = 0;
        size_t topPolyOffset = 0;
        auto currPolyIt = botPolygons.begin();
        auto nextPolyIt = topPolygons.begin();
        for (int idx = 0; currPolyIt != botPolygons.end(); idx++, currPolyIt++, nextPolyIt++) {
          auto botPolygon = *currPolyIt;
          auto topPolygon = *nextPolyIt;

          glm::vec2 botCentroid = calculatePolygonCentroid(botPolygon);
          glm::vec2 topCentroid = calculatePolygonCentroid(topPolygon);

          glm::vec2 centroidOffset = topCentroid - botCentroid;

          for (const auto& vertex : botPolygon) {
              vertPos.push_back(MatrixTransforms::Translate(botTransform, glm::vec3(vertex.x, vertex.y, 0))[3]);
          }

          float minDistance = std::numeric_limits<float>::max();
          size_t botStartVertOffset = 0,
            topStartVertOffset = 0;
          for (size_t j = 0; j < topPolygon.size(); ++j) {
            float dist = glm::distance(botPolygon[0], topPolygon[j] - centroidOffset);
            if (dist < minDistance) {
              minDistance = dist;
              topStartVertOffset = j;
            }
          }

          bool botHasMoved = false,
            topHasMoved = false;
          size_t botVertOffset = botStartVertOffset,
            topVertOffset = topStartVertOffset;
          do {
              size_t botNextVertOffset = (botVertOffset + 1) % botPolygon.size();
              size_t topNextVertOffset = (topVertOffset + 1) % topPolygon.size();

              float distBotNextToTop = glm::distance(botPolygon[botNextVertOffset], topPolygon[topVertOffset] - centroidOffset);
              float distBotToTopNext = glm::distance(botPolygon[botVertOffset], topPolygon[topNextVertOffset] - centroidOffset);
              float distBotNextToTopNext = glm::distance(botPolygon[botNextVertOffset], topPolygon[topNextVertOffset] - centroidOffset);

              bool botHasNext = botNextVertOffset != (botStartVertOffset + 1) % botPolygon.size() || !botHasMoved;
              bool topHasNext = topNextVertOffset != (topStartVertOffset + 1) % topPolygon.size() || !topHasMoved;

              if (distBotNextToTopNext < distBotNextToTop && distBotNextToTopNext <= distBotToTopNext && botHasNext && topHasNext) {
                  triVerts.emplace_back(botSectionOffset + botPolyOffset + botVertOffset,
                                        topSectionOffset + topPolyOffset + topNextVertOffset,
                                        topSectionOffset + topPolyOffset + topVertOffset);
                  triVerts.emplace_back(botSectionOffset + botPolyOffset + botVertOffset,
                                        botSectionOffset + botPolyOffset + botNextVertOffset,
                                        topSectionOffset + topPolyOffset + topNextVertOffset);
                  botVertOffset = botNextVertOffset;
                  topVertOffset = topNextVertOffset;
                  botHasMoved = true;
                  topHasMoved = true;
              } else if (distBotNextToTop < distBotToTopNext && botHasNext) {
                  triVerts.emplace_back(botSectionOffset + botPolyOffset + botVertOffset,
                                        botSectionOffset + botPolyOffset + botNextVertOffset,
                                        topSectionOffset + topPolyOffset + topVertOffset);
                  botVertOffset = botNextVertOffset;
                  botHasMoved = true;
              } else {
                  triVerts.emplace_back(botSectionOffset + botPolyOffset + botVertOffset,
                                        topSectionOffset + topPolyOffset + topNextVertOffset,
                                        topSectionOffset + topPolyOffset + topVertOffset);
                  topVertOffset = topNextVertOffset;
                  topHasMoved = true;
              }

          } while (botVertOffset != botStartVertOffset || topVertOffset != topStartVertOffset);
          botPolyOffset += botPolygon.size();
          topPolyOffset += topPolygon.size();
        }
        botSectionOffset += botSectionSize;
    }

    auto frontPolygons = sections.front();
    auto frontTriangles = manifold::Triangulate(frontPolygons, -1.0);
    for (auto& tri : frontTriangles) {
      triVerts.push_back({tri[2], tri[1], tri[0]});
    }

    auto backPolygons = sections.back();
    auto backTransform = transforms.back();
    for (const auto& poly: backPolygons) {
      for (const auto& vertex : poly) {
        vertPos.push_back(MatrixTransforms::Translate(backTransform, glm::vec3(vertex.x, vertex.y, 0))[3]);
      }
    }
    auto backTriangles = manifold::Triangulate(backPolygons, -1.0);

    for (auto& triangle : backTriangles) {
        triangle[0] += botSectionOffset;
        triangle[1] += botSectionOffset;
        triangle[2] += botSectionOffset;
        triVerts.push_back(triangle);
    }

    manifold::Mesh mesh;
    mesh.triVerts = triVerts;
    mesh.vertPos = vertPos;
    auto man = manifold::Manifold(mesh);
    return man;
}


Not much more to it than the isomorphic variant to be honest, except this loft algorithm has performed satisfactory for just everything I’ve thrown at it. It’s damn useful. I’d go so far as to say any CSG modeler should have loft in their toolkit. Stop handling every complex shape transition with hull. Use loft.

Try it out with Manifold for Clojure.