Increased mesh size/triangle count on macOS with Python 3.13+ or numpy 2.X

[shred86]

Pulling this into its own issue since the root cause is still unknown. For now, I'll revert the next version to Python 3.12 to ensure consistency in the mesh output across macOS, Windows and Linux builds.

After a lot of variations of Python versions and dependencies, I think I'm convinced there's two cases that cause the increased triangle count.

I installed Python 3.12 with all of the original requirements in 1.40.07 and before. You can see the exact versions in the 1.40.8 pull request:

• With numpy==1.26.4, triangle count is 2,340,408
• With numpy==2.3.4, triangle count is 2,942,980

I then installed Python 3.13 with the numpy in version 1.40.07 and before (version 1.26.4) and as much of the original dependencies as I could (some were too old for Python 3.13 so I had to use a newer version):

• With numpy==1.26.4, triangle count is 2,942,980
• With numpy==2.3.4, triangle count is 2,942,980

So what I'm seeing is the increased triangle count is caused by either having Python 3.12 (or older I'm assuming) with numpy 2.x or Python 3.13 (or newer), regardless of the numpy version.

Edit: I even tried some more variations:

• Python 3.12 with latest requirements except numpy==1.26.4, triangle count is 2,340,408

Doing some research, the only thing I can gather is there were changes in numpy or some of these newer dependencies required by Python 3.13, such as "preserving higher precision in elevation data or coordinates compared to NumPy 1.x’s potential downcasting" that ultimately result in a higher resolution mesh with the same settings. So it's not necessarily an "issue", just that Ortho4XP is likely producing a higher quality mesh given the same Ortho4XP settings as before. The weird part is, it seems specific to macOS (I'm on an M4 so maybe specific to M-Series). The actual issue though is if the .mesh file size is too large, Ortho4XP will get stuck on step 3 which is a known issue... but for the tiles that don't get stuck, you're likely just ending up with a higher quality mesh.

Still would like to know why it's specifically happening on macOS. I suspect the only good option will be to revert the next build for macOS to Python 3.12.

Originally posted by @shred86 in #69

[w8sl]

Python 3.12 + Numpy 1.26.4 is the best solution for now. It may depend on the the way python wheels for numpy are build for latest macOS. There have been minor issues with accuracy. The problem is even worse with the default elevation data. Mesh build on macOS, Windows and Linux should be similar. I have no idea what to do with this. Never before was such a big discrepancy reported.

[w8sl]

https://files.pythonhosted.org/packages/3e/46/bdd3370dcea2f95ef14af79dbf81e6927102ddf1cc54adc0024d61252fd9/numpy-2.3.4-cp313-cp313-macosx_11_0_arm64.whl

There are two variants of NumPy wheels available on PyPI: macosx_14 and macosx_11. The latter are built in the same way as Linux wheels, using OpenBLAS, and will produce the same results. Only one tile, and only O4XP 1.40, is affected by this problem, which is weird. There is no difference with any other tile I have tested, including 30 -85.

pip download numpy --platform macosx_11_0_arm64 --only-binary=:all:
pip install numpy-2.3.4-cp313-cp313-macosx_11_0_arm64.whl

[shred86]

Oh interesting. Just so I understand what you're seeing, the macosx_11 version does or does not have the issue? Or you're saying it doesn't have the issue except for +30-085?

Edit: With NED 1", I'm seeing on +34-118:

• macosx_11: Mesh triangles: ---> 2,340,304
• macosx_14: Mesh triangles: ---> 2,942,888

For +30-085:

• macosx_11: Mesh triangles: ---> 3,148,325
• macosx_14: Mesh triangles: ---> 3,148,325

For +48-122 (just picked a random one with a lot of mountains):

• macosx_11: Mesh triangles: ---> 2,617,282
• macosx_14: Mesh triangles: ---> 2,617,294

For +34-084 (another random one):

• macosx_11: Mesh triangles: ---> 2,725,448
• macosx_14: Mesh triangles: ---> 2,725,448

So it looks like I'm seeing what you're seeing - this issue only seems to exist on +34-118 (perhaps others, but the only one we know of so far). With the macosx_11 version of numpy, this issue doesn't seem to exist but with macosx_14 version of numpy, it does exist for that one specific tile. This is very odd indeed.

I like the idea of keeping up with the latest dependencies and luckily we've been able to keep everything up-to-date except for numpy. However, I think it seems like now we can update numpy to the latest version but just make sure to manually install the macosx_11 version to avoid this odd issue. Just one more quirk of having to configure Ortho4XP. 🤣

Great job tracking this down! I'm still wondering why that one tile in particular (+34-118) seems to have issues but the others don't. I'm sure there's some other tiles out there that will have a difference too. I added this to the development wiki page as well since I use that page as sort of a checklist when working on Ortho4XP. I'll include the correct .whl file as well in the Utils/mac folder like I do with some of the Windows .whl files just so it's in the repo for folks doing the manual install.

[w8sl]

I have briefly looked into known issues related to NumPy wheels built with Accelerate, but I haven't seen anything similar. By the way, NumPy 1.26.4 is the last default version built with openBLAS. Since only one tile is affected, there might be something specific about the O4XP or Triangle4XP code.

[shred86]

Yeah this is certainly interesting and it seems whatever the "accelerate" build dependency is for numpy seems to be the difference. What I can tell from running the following with the different versions of numpy:

import numpy as np
from numpy.__config__ import show
show()

It's that, as you mentioned, numpy 1.26.4 but also numpy 2.3.4 (macosx_11_0_arm64) use a dependency for "blas" called "openblas64". Slightly different versions but still "openblas64" where numpy 2.3.4 (macosx_14_0_arm64) uses "accelerate".

numpy 1.26.4 (default version installed with pip) on python 3.12:

{
  "Compilers": {
    "c": {
      "name": "clang",
      "linker": "ld64",
      "version": "14.0.0",
      "commands": "cc",
      "args": "-fno-strict-aliasing, -DBLAS_SYMBOL_SUFFIX=64_, -DHAVE_BLAS_ILP64",
      "linker args": "-fno-strict-aliasing, -DBLAS_SYMBOL_SUFFIX=64_, -DHAVE_BLAS_ILP64"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "3.0.8",
      "commands": "cython"
    },
    "c++": {
      "name": "clang",
      "linker": "ld64",
      "version": "14.0.0",
      "commands": "c++",
      "args": "-DBLAS_SYMBOL_SUFFIX=64_, -DHAVE_BLAS_ILP64",
      "linker args": "-DBLAS_SYMBOL_SUFFIX=64_, -DHAVE_BLAS_ILP64"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "aarch64",
      "family": "aarch64",
      "endian": "little",
      "system": "darwin"
    },
    "build": {
      "cpu": "aarch64",
      "family": "aarch64",
      "endian": "little",
      "system": "darwin"
    }
  },
  "Build Dependencies": {
    "blas": {
      "name": "openblas64",
      "found": true,
      "version": "0.3.23.dev",
      "detection method": "pkgconfig",
      "include directory": "/opt/arm64-builds/include",
      "lib directory": "/opt/arm64-builds/lib",
      "openblas configuration": "USE_64BITINT=1 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= SANDYBRIDGE MAX_THREADS=3",
      "pc file directory": "/usr/local/lib/pkgconfig"
    },
    "lapack": {
      "name": "dep4335021056",
      "found": true,
      "version": "1.26.4",
      "detection method": "internal",
      "include directory": "unknown",
      "lib directory": "unknown",
      "openblas configuration": "unknown",
      "pc file directory": "unknown"
    }
  },
  "Python Information": {
    "path": "/private/var/folders/76/zy5ktkns50v6gt5g8r0sf6sc0000gn/T/cibw-run-q69bfk1p/cp312-macosx_arm64/build/venv/bin/python",
    "version": "3.12"
  },
  "SIMD Extensions": {
    "baseline": [
      "NEON",
      "NEON_FP16",
      "NEON_VFPV4",
      "ASIMD"
    ],
    "found": [
      "ASIMDHP"
    ],
    "not found": [
      "ASIMDFHM"
    ]
  }
}

numpy 2.3.4 (macosx_11_0_arm64) on python 3.13:

{
  "Compilers": {
    "c": {
      "name": "clang",
      "linker": "ld64",
      "version": "15.0.0",
      "commands": "cc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "3.1.4",
      "commands": "cython"
    },
    "c++": {
      "name": "clang",
      "linker": "ld64",
      "version": "15.0.0",
      "commands": "c++"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "aarch64",
      "family": "aarch64",
      "endian": "little",
      "system": "darwin"
    },
    "build": {
      "cpu": "aarch64",
      "family": "aarch64",
      "endian": "little",
      "system": "darwin"
    }
  },
  "Build Dependencies": {
    "blas": {
      "name": "scipy-openblas",
      "found": true,
      "version": "0.3.30",
      "detection method": "pkgconfig",
      "include directory": "/private/var/folders/y6/nj790rtn62lfktb1sh__79hc0000gn/T/cibw-run-vi2_0tkd/cp313-macosx_arm64/build/venv/lib/python3.13/site-packages/scipy_openblas64/include",
      "lib directory": "/private/var/folders/y6/nj790rtn62lfktb1sh__79hc0000gn/T/cibw-run-vi2_0tkd/cp313-macosx_arm64/build/venv/lib/python3.13/site-packages/scipy_openblas64/lib",
      "openblas configuration": "OpenBLAS 0.3.30  USE64BITINT DYNAMIC_ARCH NO_AFFINITY neoversen1 MAX_THREADS=64",
      "pc file directory": "/Users/runner/work/numpy-release/numpy-release/.openblas"
    },
    "lapack": {
      "name": "scipy-openblas",
      "found": true,
      "version": "0.3.30",
      "detection method": "pkgconfig",
      "include directory": "/private/var/folders/y6/nj790rtn62lfktb1sh__79hc0000gn/T/cibw-run-vi2_0tkd/cp313-macosx_arm64/build/venv/lib/python3.13/site-packages/scipy_openblas64/include",
      "lib directory": "/private/var/folders/y6/nj790rtn62lfktb1sh__79hc0000gn/T/cibw-run-vi2_0tkd/cp313-macosx_arm64/build/venv/lib/python3.13/site-packages/scipy_openblas64/lib",
      "openblas configuration": "OpenBLAS 0.3.30  USE64BITINT DYNAMIC_ARCH NO_AFFINITY neoversen1 MAX_THREADS=64",
      "pc file directory": "/Users/runner/work/numpy-release/numpy-release/.openblas"
    }
  },
  "Python Information": {
    "path": "/private/var/folders/y6/nj790rtn62lfktb1sh__79hc0000gn/T/build-env-tp996pz9/bin/python",
    "version": "3.13"
  },
  "SIMD Extensions": {
    "baseline": [
      "NEON",
      "NEON_FP16",
      "NEON_VFPV4",
      "ASIMD"
    ],
    "found": [
      "ASIMDHP"
    ],
    "not found": [
      "ASIMDFHM"
    ]
  }
}

numpy 2.3.4 (macosx_14_0_arm64) on python 3.13:

{
  "Compilers": {
    "c": {
      "name": "clang",
      "linker": "ld64",
      "version": "15.0.0",
      "commands": "cc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "3.1.4",
      "commands": "cython"
    },
    "c++": {
      "name": "clang",
      "linker": "ld64",
      "version": "15.0.0",
      "commands": "c++"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "aarch64",
      "family": "aarch64",
      "endian": "little",
      "system": "darwin"
    },
    "build": {
      "cpu": "aarch64",
      "family": "aarch64",
      "endian": "little",
      "system": "darwin"
    }
  },
  "Build Dependencies": {
    "blas": {
      "name": "accelerate",
      "found": true,
      "version": "unknown",
      "detection method": "system",
      "include directory": "unknown",
      "lib directory": "unknown",
      "openblas configuration": "unknown",
      "pc file directory": "unknown"
    },
    "lapack": {
      "name": "accelerate",
      "found": true,
      "version": "unknown",
      "detection method": "system",
      "include directory": "unknown",
      "lib directory": "unknown",
      "openblas configuration": "unknown",
      "pc file directory": "unknown"
    }
  },
  "Python Information": {
    "path": "/private/var/folders/y6/nj790rtn62lfktb1sh__79hc0000gn/T/build-env-hj_u2xtv/bin/python",
    "version": "3.13"
  },
  "SIMD Extensions": {
    "baseline": [
      "NEON",
      "NEON_FP16",
      "NEON_VFPV4",
      "ASIMD"
    ],
    "found": [
      "ASIMDHP"
    ],
    "not found": [
      "ASIMDFHM"
    ]
  }
}

Based on some questions I've asked using ChatGPT and Grok (I find LLMs pretty helpful for these type of questions), it seems the only difference between OpenBLAS and Apple's Accelerate is with linear algebra operations. Specifically, "openblas64 uses the ILP64 interface (64-bit int), while standard OpenBLAS and Accelerate use LP64 (32-bit int). This affects maximum array size and can change computation paths for large matrices."

Searching for some numpy linear algebra operations, they are being used in O4_Vector_Utils.py, such as numpy.linalg and numpy.dot. O4_Vectors_Utils.py is used by O4_Mesh_Utils.py which creates the mesh files. The weird part is the openblas64 actually supports larger arrays/datasets but we're talking in the 2+ billion range.

Anyways, there's obviously some differences between the two. I suspect it's maybe something to do with the linear algebra function/methods in numpy but not entirely sure.

[w8sl]

The problem with LLMs is, if the answer is unknown, they start to hallucinate. Anyway, the best answer i could get:

"If your meshing algorithm uses thresholds (e.g., max edge length, curvature), tiny differences in node spacing or alignment can tip elements over the threshold, triggering more refinement."

[w8sl]

Great catch: 5f0f7df
Is this related to this issue ?

I noticed that Oscar reduced a value in version 1.40 compared to 1.30 to prevent crashes: vector_map.snap_to_grid

I further decreased it to 8, which seems to address the problem of excessive triangles on some tiles while maintaining mesh accuracy well above available elevation and OSM data. However, I haven't had time for extensive testing over larger areas.

[shred86]

Great catch: 5f0f7df Is this related to this issue ?

I noticed that Oscar reduced a value in version 1.40 compared to 1.30 to prevent crashes: vector_map.snap_to_grid

I further decreased it to 8, which seems to address the problem of excessive triangles on some tiles while maintaining mesh accuracy well above available elevation and OSM data. However, I haven't had time for extensive testing over larger areas.

I was thinking the same thing but I don’t think it’s related. I tested the tile with the issue above:

Edit: With NED 1", I'm seeing on +34-118:

macosx_11: Mesh triangles: ---> 2,340,304
macosx_14: Mesh triangles: ---> 2,942,888

With macosx_11: triangle count was slightly higher (several thousand).

With macosx_14: triangle count matched what it was before, so still significantly higher and it was exact (2,942,888).

The odd part is, I tested a couple tiles across Windows and Mac, and just by recompiling Triangle4XP without the fix in Triangle4XP.c, the triangle count always seems slightly higher. With the fix in place, it didn’t seem to change (stayed the same as after recompiling without the fix).

I recall I reverted the Triangle4XP.exe on Windows based on an issue and it appeared a different version was accidentally included when oscarpilote uploaded the Mac native arm version, so it’s tough to tell exactly which source it’s built off of. It’s obviously clear quite a few changes were made in Triangle4XP.c from version 1.30.

All that to say, as of now, I don’t think it’s related to this particular issue but more testing is certainly needed.

[shred86]

I decided to have Claude (Opus 4.6) investigate this issue and figured I'd post the results here for reference. I may mess around with this at some point to compare the difference in triangle count with the suggested changes:

Root Cause: numpy.linalg operations on 2x2 matrices during vector map construction
The core problem is not in the DEM data itself or the Triangle4XP C binary — it's in the Python-side vector map building that happens before Triangle4XP is invoked. This is where different BLAS backends produce slightly different floating-point results, which cascade into different mesh topologies.
The smoking gun: O4_Vector_Utils.py
The are_encroached() method (line 252) uses several BLAS-dependent numpy operations on small 2D vectors:
numpy.linalg.det(A) (line 281) — determinant of a 2x2 matrix
numpy.linalg.solve(A, ac) (line 283) — solve a 2x2 linear system
numpy.linalg.norm() (line 271) — vector norms
numpy.dot() (line 273) — dot products
These results are compared against tight tolerances like eps = 1e-8 to decide whether two edge segments intersect. When they do, new intersection nodes are inserted into the polygon mesh, splitting existing edges.
How this causes different triangle counts
The insert_edge() method (line ~115) calls are_encroached() for every potential edge intersection when building the constraint polygon (the .poly file fed to Triangle4XP). A tiny floating-point difference in a determinant or solve result can flip a comparison like:
if abs(numpy.linalg.det(A)) > eps * pnorm: # line 281

This means:
With Accelerate (LP64): a det might compute to 1.00000001e-8 → intersection detected → node inserted → edge split
With OpenBLAS (ILP64): the same det computes to 9.99999999e-9 → no intersection → no split
Each extra/missing constraint node in the .poly file changes the Delaunay triangulation, and then the curvature-based refinement in Triangle4XP amplifies this — a slightly different initial triangulation leads to different triangles being flagged for subdivision, cascading into ~600K extra triangles.
Why it's tile-dependent
The tile +34-118 (Los Angeles area) has complex coastline and airport geometry, producing many edge intersections with values near the tolerance thresholds. Simpler tiles with fewer constraint edges are less likely to hit borderline cases.

The suggested fix:

Use pure Python math for the 2x2 operations (det, solve) instead of numpy.linalg — for 2x2 matrices this is trivial and eliminates BLAS entirely

Specifically:

Added import math at line 0
Replaced in are_encroached():
numpy.linalg.norm(ab) → math.sqrt(ab[0] * ab[0] + ab[1] * ab[1])
numpy.dot(ab, dc) → ab[0] * dc[0] + ab[1] * dc[1]
numpy.linalg.det(A) → ab[0] * dc[1] - ab[1] * dc[0] (explicit 2x2 determinant)
numpy.linalg.solve(A, ac) → Cramer's rule with explicit formulas
numpy.linalg.norm(ac) → math.sqrt(ac[0] * ac[0] + ac[1] * ac[1])
numpy.argmax(abs(ab)) → 0 if abs(ab[0]) >= abs(ab[1]) else 1

The specific change to the are_encroached method in O4_Vectors_Utils:

        ab = b - a
        dc = c - d
        norm_ab = math.sqrt(ab[0] * ab[0] + ab[1] * ab[1])
        norm_dc = math.sqrt(dc[0] * dc[0] + dc[1] * dc[1])
        pnorm = norm_ab * norm_dc
        ac = c - a
        ab_dot_dc = ab[0] * dc[0] + ab[1] * dc[1]
        if ((a == d).all() or (b == c).all()) and ab_dot_dc < 0.9999 * pnorm:
            return False
        if ((a == c).all() or (b == d).all()) and ab_dot_dc > -0.9999 * pnorm:
            return False
        eps = 1e-8
        oneminuseps = 1.0 - eps
        # det of 2x2 matrix [ab, dc]
        det_A = ab[0] * dc[1] - ab[1] * dc[0]
        if abs(det_A) > eps * pnorm:
            # ad and bc are not considered parallel
            # Cramer's rule for 2x2 solve: [ab, dc] * [alpha, beta] = ac
            alpha = (ac[0] * dc[1] - ac[1] * dc[0]) / det_A
            beta = (ab[0] * ac[1] - ab[1] * ac[0]) / det_A
            return (
                (alpha >= 0 and alpha <= 1)
                and (beta >= 0 and beta <= 1)
                and (
                    (alpha > eps and alpha < oneminuseps)
                    or (beta > eps and beta < oneminuseps)
                )
                and (alpha, beta)
            )
        elif abs(ab[0] * ac[1] - ab[1] * ac[0]) > eps * norm_ab * math.sqrt(ac[0] * ac[0] + ac[1] * ac[1]):
            # ad and bc are parallel but not colinear
            return False
        else:
            # ad and bc are parallel and colinear
            g_idx = 0 if abs(ab[0]) >= abs(ab[1]) else 1
            d_idx = 0 if abs(dc[0]) >= abs(dc[1]) else 1
            alpha0, alpha1 = ac[g_idx] / ab[g_idx], (d - a)[g_idx] / ab[g_idx]
            beta0, beta1 = ac[d_idx] / dc[d_idx], (c - b)[d_idx] / dc[d_idx]
            return (
                (alpha0 > eps or alpha1 > eps)
                and (alpha0 < oneminuseps or alpha1 < oneminuseps)
                and (alpha0, alpha1, beta0, beta1)
            )

It makes sense to me but I'm just not sure how different the end result (triangle count) would end up being using Python math versus numpy. Also, triangle count is one thing but I also wonder if there would be any difference in visual quality between the two which would be hard to test.

[w8sl]

Doesn't seem to change the number of triangles.

Numpy:

Math:

Accelerate with 8 digits precision gives result in the same range as OpenBLAS (8 or 9 digits):