diff --git a/openspec/changes/use-tangle-cutout-meshing/.openspec.yaml b/openspec/changes/use-tangle-cutout-meshing/.openspec.yaml new file mode 100644 index 0000000..c238415 --- /dev/null +++ b/openspec/changes/use-tangle-cutout-meshing/.openspec.yaml @@ -0,0 +1,2 @@ +schema: spec-driven +created: 2026-09-13 diff --git a/openspec/changes/use-tangle-cutout-meshing/design.md b/openspec/changes/use-tangle-cutout-meshing/design.md new file mode 100644 index 0000000..ba6716f --- /dev/null +++ b/openspec/changes/use-tangle-cutout-meshing/design.md @@ -0,0 +1,64 @@ +## Context + +当前 Fine Outline 使用轮廓近似、内部支撑、自由点松弛和 Blender CDT，再通过 `support_flat_regions` 拆分边界弦及扇形补面；普通模式在 `cutout_domain` 中采用另一套三角化及 `refine_balloon_triangles`。两条最终建网路径需要统一，减少质量细化后的再次碎分。 + +临时实验采用相同轮廓、延伸中线和 10°目标，Tangle 的鱼/甲虫分别为 2797/655 面，Triangle 为 2529/605 面；Tangle 的六张样例均通过面积、孔洞数、中线覆盖和非退化检查。这些是外部 Python 实验，尚不能替代 Blender 验收。实验保留全部有效孔洞，因此新增微孔策略需要独立验证。 + +现有世界 Edge Length 与相对下限已完成，所有几何尺度沿用其有效间距 h。早期变化中的 Mesh Detail 档位描述以该世界边长行为为准。 + +## Goals / Non-Goals + +**Goals:** 在统一基础网格中保持轮廓、有效孔洞和细分支内部支撑，控制普通区域密度，使用可随 Blender 分发的纯 Python 后端。 + +**Non-Goals:** 本次不新增算法选择器、纯平面专用算法或 C++ 二进制，保留现有高度算法。 + +## Decisions + +### 1. 随包分发固定版本的 Tangle Python 源码 + +第三方实现放在 `src/anyimage/operators/cutout_tool/tangle.py`，MIT 文本放在现有 `src/anyimage/licenses/`。实施时固定上游提交及源码哈希，记录来源与必要本地改动，避免跟随浮动 master。实验版本标识为 0.4.18，实际分发以核验的提交为准。 + +将数值调用整理为内存输入输出：点、带身份的约束边、孔内种子、质量参数，返回顶点、三角面和约束来源。保留上游算法结构，业务预处理和验证由 Cutout 网格模块负责。调用不依赖命令行参数、临时 `.poly` 文件或额外 Python 进程；隔离 Decimal 精度、模块状态和退出处理，重复调用互不影响。优先使用已有标准库功能，Tangle 本身不新增第三方包。 + +### 2. Fine Outline 决定建网区域与轮廓精度 + +从当前 Alpha Threshold 和 Selection 得到前景，沿用极小独立前景分量清理规则。先按世界平面度量确定有效 h，再做孔洞策略，之后重新提取组件与轮廓： + +| 模式 | 内部封闭孔洞 | 外轮廓 | 支撑依据 | +|---|---|---|---| +| 开启 | 保留有效孔洞，填充微孔 | 随 h 保留精细形状 | 微孔处理后的区域 | +| 关闭 | 全部填充 | 随 h 平滑、概括 | 填孔及轮廓概括后的区域 | + +微孔初始规则采用双条件：世界平面中的包围盒最长边 ≤ 0.25h，且面积 ≤ 0.04h²，才填充。长而窄的孔保留；0.9h 的孔不会因略小于 h 被移除。实施中先按此保守规则建立边界用例，真实图只用于校准与记录，不依赖其像素尺寸定义阈值。填孔后重建拓扑，处理孔中岛合并，避免重复面。 + +几何填孔不改写图片或材质 Alpha。背景连通的凹口不是封闭孔洞，仍参与外轮廓。简化后的最终轮廓必须通过有效性检查，内部支撑直接来自最终网格。 + +### 3. 质量网格按需补足内部支撑 + +输入 PSLG 仅包含最终轮廓约束，在与 h 关联的归一化坐标中合并重合点、分割点落边并保留约束身份。各组件以辅助三角化的最大面中心作为初始内部点，孔内种子必须位于对应孔内。像素角自接触处按相邻边长的 1/8 分离，消除零宽连接；简化结果重新验证，保留有效孔洞。 + +Tangle 使用 10°目标最小角和最大面积 √3/4 × h²。完成面积、孔洞、组件和边界约束校验后，给两端都在边界上的内部边添加共享中点，并一致拆分相邻面；只有孤立全边界三角形需要补重心。已有内部支撑的区域保持原网格。此操作单次完成，新增点数受组件累计顶点上限约束，不重新进行全局三角化。 + +局部拆分保持轮廓、面积和面朝向，形成细分支到主体的连续内部支撑。它可能降低部分三角形最小角，10°作为基础三角化目标，最终角度与凸起效果共同验收。删除不再使用的中线图、宽度距离场和路径约束处理；辅助 Blender CDT 用于初始点和孔种子。 + +低分辨率深度投影可能将相邻轮廓采样映射到同一点。深度 Cutout 在投影和平滑后，以平均边长的百万分之一合并已连接重合点，再求壳体厚度，避免数值退化面翻转法线。节点资产与实际求值测试同步验证。 + +### 4. 验收同时覆盖平面几何和高度 + +- 自动化用例使用可再生的合成多孔形状、弯曲细腿、尖角、孔中岛和点落边约束，避免测试依赖用户临时图片。 +- 最终网格验证面积、孔洞及组件数量、非退化面、轮廓覆盖；细分支到主体的内部网格边连续连接，内部支撑高度为有限正值。 +- 六张用户图用于本地视觉和性能复核，重点比较鱼的多孔区域、甲虫上角到根部及细腿。报告面数、低于 10°占比、小面面积占比、压平面积与位置、耗时，不要求三角面索引或面数与实验逐项相同。 +- 分别验证 Fine Outline 两种模式及世界缩放、非均匀缩放、分辨率和透明画布变化。关闭模式使用填孔后的期望拓扑。 +- 在项目支持的 Blender 环境中验证源码导入、重复数值调用、基础网格创建与同步/AI 返回路径；包括现有分发覆盖的 Python 3.11/3.13 环境。缺少其中环境时明确记录，不以外部解释器验证代替。 + +## Risks / Trade-offs + +- 纯 Python 在密集输入上变慢 → 用预热与重复测量分别记录三角化、完整建网与 Blender 创建；对大规模多孔输入实现确定的迭代/顶点资源上限和明确失败结果，不静默降级。 +- 上游命令行实现含进程级状态 → 在数值接口测试 Decimal 上下文、错误退出、重复调用与失败后恢复。 +- 孔洞密集或边界尖角仍要求小面 → 保守微孔策略和质量统计共同控制，避免为了禁止短边破坏主体。 +- 局部支撑拆分降低部分面角度 → 复核最终角度和粗网格凸起效果。 +- 六张图覆盖有限 → 合成回归、Blender 实测与高复杂度资源测试共同验收。 + +## Migration Plan + +先固定后端与数值接口，再接入统一区域、中线及三角化，最后清理旧调用。检查源码和许可被现有打包规则收录；仅在实际改变依赖元数据时同步相关清单。运行相关测试，不构建扩展包或同步其他说明文档。若验收失败，在当前实现阶段修复；已集成变更的回退使用版本控制恢复整个变更，不保留运行时旧算法分支。 diff --git a/openspec/changes/use-tangle-cutout-meshing/proposal.md b/openspec/changes/use-tangle-cutout-meshing/proposal.md new file mode 100644 index 0000000..1e71cff --- /dev/null +++ b/openspec/changes/use-tangle-cutout-meshing/proposal.md @@ -0,0 +1,28 @@ +## Why + +Cutout 在多孔区域和细长结构附近产生拉伸三角形，现有补面细分进一步增加碎面。六张透明图的实验表明，Tangle Python 版能够保留孔洞与延伸中线，效果接近已认可的 Triangle，且仅依赖 Python 标准库、采用 MIT 许可，适合随 Blender 插件分发。 + +## What Changes + +- 使用随插件分发的 Tangle Python 实现统一构建 Cutout 基础三角网格，质量目标为 10°，密度由有效世界 Edge Length 驱动。 +- Fine Outline 开启时保留精细轮廓和有效孔洞，仅过滤远小于目标尺度的微孔；关闭时填充全部封闭内部孔洞并概括外轮廓。 +- 两种模式均先生成质量网格，拆分两端位于边界的内部弦，为孤立的全边界三角形补一个内部点，保持细长分支到主体的连续支撑。 +- 用质量三角化替换最终建网中的旧布点与扇形补面流程，检查孔洞、实际约束链、退化面和高度结果。 +- 固定第三方源码版本并保留许可；验证 Blender 实际调用与资源边界。 + +## Capabilities + +### New Capabilities + +- `cutout-quality-meshing`: 定义统一三角化、Fine Outline 孔洞策略、内部支撑与 Blender 接入验收。当前 `openspec/specs/` 为空，本变更建立此能力规格；世界边长行为承接已完成的 `use-world-space-cutout-edge-length`。 + +### Modified Capabilities + +无已归档主规格需要修改。 + +## Impact + +- 涉及 `src/anyimage/operators/cutout_tool/geometry.py`、统一建网模块 `mesh.py`（替换 `fine_mesh.py`）及基础网格调用链，新增第三方 Tangle 源码和 `licenses/tangle.txt`。 +- 涉及 Cutout 几何、世界边长、Operator 与依赖分发相关测试。 +- 运行时新增依赖为随包携带的纯 Python 源码；原有 NumPy、SciPy 继续服务于轮廓和高度计算。 +- Fine Outline 的设置捕获继续覆盖同步与 AI 返回路径，基础网格供现有 Cutout Shape 使用。 diff --git a/openspec/changes/use-tangle-cutout-meshing/specs/cutout-quality-meshing/spec.md b/openspec/changes/use-tangle-cutout-meshing/specs/cutout-quality-meshing/spec.md new file mode 100644 index 0000000..dee1315 --- /dev/null +++ b/openspec/changes/use-tangle-cutout-meshing/specs/cutout-quality-meshing/spec.md @@ -0,0 +1,71 @@ +## ADDED Requirements + +### Requirement: Unified quality meshing +系统 SHALL 在 Fine Outline 两种状态下使用 Tangle Python 版构建最终基础三角网格，采用 10°质量目标及与有效世界边长 h 对应的最大面积 √3/4 × h²。系统 SHALL 保持现有世界 Edge Length、相对下限和非均匀缩放语义，允许细结构产生短于 h 的边。 + +#### Scenario: Shared density controls +- **WHEN** 相同物理主体更换分辨率、透明画布或对象旋转后创建 Cutout +- **THEN** 同一有效 h 驱动三角化；两种 Fine Outline 模式均使用 Tangle + +#### Scenario: Acute constrained feature +- **WHEN** 必须保留的细结构包含尖角或短约束 +- **THEN** 系统保留结构，报告实际质量统计，不以删除有效结构强行满足最短边或角度目标 + +### Requirement: Outline mode defines holes +系统 SHALL 在关闭 Fine Outline 时填充所有封闭内部孔洞并概括外轮廓；开启时仅在孔最长世界跨度 ≤ 0.25h 且面积 ≤ 0.04h² 时填孔，保留其他孔洞。系统 SHALL 从处理后区域重建组件，且不修改原图片和材质 Alpha。 + +#### Scenario: Fine mode keeps a moderately small hole +- **WHEN** 孔最长跨度为 0.9h，或细长孔的最长跨度大于 0.25h +- **THEN** 开启 Fine Outline 时孔被保留，即使孔面积很小 + +#### Scenario: Fine mode fills a micro hole +- **WHEN** 封闭孔满足两个微孔条件 +- **THEN** 开启 Fine Outline 也填充该孔，后续支撑和三角化使用填孔结果 + +#### Scenario: Coarse mode fills a nested hole +- **WHEN** 关闭 Fine Outline，内部孔包含独立前景岛 +- **THEN** 填孔后区域重新合并，不产生重叠面；与外部背景相通的凹口仍属于外轮廓 + +### Requirement: Continuous interior support +系统 SHALL 在两种模式下按实际质量网格补足内部支撑。系统 SHALL 拆分两端位于边界的内部弦，为孤立全边界三角形添加内部点，保持细分支经关节到主体的连续内部连接。 + +#### Scenario: Coarse branch requires support +- **WHEN** 低密度网格的细腿或根部缺少内部顶点 +- **THEN** 系统添加局部支撑，Poisson 求解后内部顶点高度为正，不产生全零高度三角形 + +#### Scenario: Dense region already has support +- **WHEN** 网格区域已经具有连续内部支撑 +- **THEN** 系统保持该区域，不生成额外全局中线约束，不重新三角化 + +### Requirement: Final constraints and topology are validated +系统 SHALL 验证面积、轮廓约束链、孔洞和组件数量、非退化性，且局部支撑拆分 SHALL 保持已验证的轮廓、拓扑与面朝向。像素角自接触 SHALL 做亚像素分离并重新验证。 + +#### Scenario: Pixel-corner hole contact +- **WHEN** 孔洞轮廓在像素角形成零宽自接触 +- **THEN** 系统分离接触点，生成拓扑有效且保留孔洞的网格 + +#### Scenario: Mesher loses a boundary constraint +- **WHEN** 后端遗漏边界约束或产生错误面积 +- **THEN** 系统明确失败，不创建残缺对象 + +### Requirement: Height support preserves mesh quality +系统 SHALL 在单次局部拆分中补足支撑，并限制累计顶点数。最终角度允许低于基础三角化的 10°目标，验收 SHALL 同时检查角度、连续正高度与压平面积。 + +#### Scenario: Connected depth samples collapse +- **WHEN** 深度投影将相邻顶点映射到数值上重合的位置 +- **THEN** 系统按当前平均边长的百万分之一合并已连接重合点，壳体厚度保持正确方向 + +### Requirement: Blender distribution and bounded execution +系统 SHALL 分发固定版本的 Tangle Python 源码和 MIT 许可，通过内存数值接口调用，仅使用标准库作为 Tangle 自身依赖。系统 SHALL 隔离调用状态，保持 Blender 数据创建在主线程，并对资源超限明确终止。 + +#### Scenario: Clean Blender installation +- **WHEN** 在项目支持的 Blender 环境加载插件并创建 Cutout +- **THEN** Tangle 无需用户 pip 安装或外部编译器即可调用，源码和许可被分发规则收录 + +#### Scenario: Repeated and failed calls +- **WHEN** 连续构建不同输入，包含一次无效输入或资源超限 +- **THEN** 后续调用不受先前状态影响，Decimal 上下文保持隔离，失败不创建残缺对象 + +#### Scenario: Local and AI return paths +- **WHEN** 同步创建或 AI 返回后创建 Cutout +- **THEN** 两者采用捕获的 Fine Outline 设置及同一有效世界边长与网格校验规则 diff --git a/openspec/changes/use-tangle-cutout-meshing/tasks.md b/openspec/changes/use-tangle-cutout-meshing/tasks.md new file mode 100644 index 0000000..3c117d4 --- /dev/null +++ b/openspec/changes/use-tangle-cutout-meshing/tasks.md @@ -0,0 +1,52 @@ +## 1. 固定数值后端 + +- [x] 1.1 核验并固定 Tangle Python 上游提交、源码哈希和 MIT 许可，加入随插件分发的源码与许可文件。 +- [x] 1.2 提供内存数值入口，返回三角面和约束来源；隔离 Decimal 上下文、命令行退出与可变状态，验证重复调用和失败后恢复。 +- [x] 1.3 验证现有分发规则收录源码与许可；检查依赖清单一致性，不构建扩展包。 + +## 2. 统一轮廓和孔洞策略 + +- [x] 2.1 在现有有效世界 Edge Length 度量下实现 Fine Outline 开启的保守微孔双阈值，以及关闭时的全部封闭孔填充。 +- [x] 2.2 填孔后重建组件与轮廓，处理孔中岛、背景连通凹口和空输入；保持原图片和材质 Alpha。 +- [x] 2.3 接入两种外轮廓精度策略，确保中线距离场、域查询和最终轮廓一致；覆盖分辨率、画布与非均匀缩放用例。 + +## 3. 延伸中线并接入 Tangle + +- [x] 3.1 将细分支支撑沿有效内部图延伸至所在主体，保留稀疏性并检查关节、根部和跨孔路径。 +- [x] 3.2 规范化 PSLG 重合点、点落边和重复约束，生成可靠孔内种子并保留中线身份。 +- [x] 3.3 两种模式共用 10°和最大面积参数的 Tangle 最终三角化，输出现有基础网格所需组件、UV、边界及高度输入。 +- [x] 3.4 验证最终约束链、面积、孔洞、组件与非退化性；增加有界细化和明确错误传播，失败不创建残缺对象。 +- [x] 3.5 删除最终建网中失去用途的旧布点、三角化及扇形补面路径，保留仍用于中线构图的必要几何函数并理顺命名。 + +## 4. 验证与复核 + +- [x] 4.1 补充合成多孔、细腿根部、微孔阈值两侧、细长孔、孔中岛与点落边回归；检查实际边链及正高度支撑。 +- [x] 4.2 用六张用户图复核两种模式，重点输出鱼与甲虫的网格和高度对照；记录面数、角度、小面和压平面积、位置及耗时。 +- [x] 4.3 在支持的 Blender/Python 环境执行导入、重复调用、基础网格创建及同步/AI 返回路径验证，记录实际覆盖版本。 +- [x] 4.4 对高复杂度多孔输入验证迭代/顶点限制、失败恢复及耗时，确认不会无界细分或静默切回旧算法。 +- [x] 4.5 运行相关 Cutout 几何、边长、对象、Operator 和分发测试，检查最终差异与旧引用；仅在发现关联故障时扩大测试范围。 + +## 验证记录 + +- 相关几何、边长、对象、Operator、交互、节点消费与分发规则测试：141 passed，7 subtests passed。分发测试使用内存归档替身，未构建扩展包。 +- Blender 4.5.10 / Python 3.11.11、Blender 5.2.0 / Python 3.13.13：各 81 passed、4 subtests passed，正常退出。使用 Blender 内置 NumPy 和对应随包 SciPy wheel；5.2 测试通过运行时 socket API 读取修改器输入。 +- 六张原图的两种模式均通过；各预热一次并测量两次。鱼精细/普通为 2266/972 面、25/0 孔，甲虫为 669/382 面。精细模式压平面积：鱼 0.056%，甲虫 0%；精细建网中位耗时分别约 0.67s、0.13s（不含图片读取及高度求解）。 +- 256 孔合成图：2K 与 4K 在同一物理间距下均得到 3448 顶点、6044 面，构建分别约 2.25s、6.23s；资源上限及失败后恢复另有自动化测试覆盖。 +- 临时对照结果：`C:/Users/icrdr/AppData/Local/Temp/anyimage-mesh-trial-e5fd3fe134d24686900f07f9308db158/integrated-results.html`，含六图网格、Poisson 高度投影与两种模式对照。 +- 最终源码使用 `mesh.py` 承担两种模式的统一建网职责，替换原 `fine_mesh.py`；第三方版本与原始哈希记录于 `tangle.py` 文件头。 + +## 5. 按实际网格补足支撑 + +- [x] 5.1 用单次局部补点替换中线图和距离场，保留连续内部支撑及累计顶点上限，删除旧中线调用与专用测试。 +- [x] 5.2 分离像素角孔洞自接触，保持面积、孔洞和边界约束校验，复核密鱼能正常完成。 +- [x] 5.3 修复低分辨率深度投影的相邻点重合，构建节点资产并验证壳体方向。 +- [x] 5.4 运行粗密网格、两种 Fine Outline、真实样本和 Blender 相关回归，记录结果。 + +### 局部支撑复核 + +- 删除中线图、宽度场和全局支撑约束，单次局部补点保持内部连通；新增粗密分支、孤立三角形、共享弦、资源上限、像素角接触及边界损坏测试。最终网格专项 45 passed。 +- 六张用户原图、两种 Fine Outline 共 20 组配置通过；鱼和甲虫包含最长边 / 5、/ 20、/ 100 三档。全部 Poisson 高度有限、内部顶点为正，无全零高度面。密鱼精细模式 27077 面、55 孔，约 2–3 秒；密甲虫 9272 面，约 0.67–1.03 秒（本机单次建网，区间含并行测试负载）。 +- 极粗精细模式：甲虫 156 面，鱼 585 面，压平面积均为 0%。局部拆分允许部分角度低于 10°。 +- 深度投影合并数值重合的已连接顶点，修复角部退化面导致的壳体厚度翻转。已运行 `uv run node-group build`：10 个保存节点组加载验证通过，构建流程全部 1024 测试通过并正常退出。 +- Blender 4.5.10 / Python 3.11.11、Blender 5.2.0 / Python 3.13.13：各 81 passed、4 subtests passed，正常退出；未构建扩展包。 +- 临时原图结果：`C:/Users/icrdr/AppData/Local/Temp/anyimage-support-integrated/results.json`。 diff --git a/src/anyimage/assets/O_AnyImage.blend b/src/anyimage/assets/O_AnyImage.blend index 0e85b93..44142f0 100644 --- a/src/anyimage/assets/O_AnyImage.blend +++ b/src/anyimage/assets/O_AnyImage.blend @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:ee7629693ef157c5a77ba6319637cfbab3f775a69f7e68831fde647a4ad23b8a -size 9733198 +oid sha256:dfb219aae1aaaa7031dad450b021a9204cf35b3524bd79a8aa57fa3267dc938b +size 9802561 diff --git a/src/anyimage/licenses/tangle.txt b/src/anyimage/licenses/tangle.txt new file mode 100644 index 0000000..379ef9b --- /dev/null +++ b/src/anyimage/licenses/tangle.txt @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2026 David Meeker + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/src/anyimage/operators/cutout_tool/fine_mesh.py b/src/anyimage/operators/cutout_tool/fine_mesh.py deleted file mode 100644 index 9994ca8..0000000 --- a/src/anyimage/operators/cutout_tool/fine_mesh.py +++ /dev/null @@ -1,636 +0,0 @@ -"""Fine Cutout contours, sparse constrained support and interior sampling.""" - -import heapq - -import numpy as np -from mathutils import Vector -from mathutils.geometry import delaunay_2d_cdt -from scipy import ndimage -from scipy.spatial import cKDTree -from .polygon import approximate_polygon - -from . import geometry as g - - -def valid_contours(contours): - """Reject crossing rings and holes moved outside the outer ring.""" - domain = ContourDomain(contours) - count = len(domain.starts) - following = np.concatenate( - [ - np.roll(np.arange(len(contour)) + sum(map(len, contours[:index])), -1) - for index, contour in enumerate(contours) - ] - ) - for first in range(0, count, 32): - ids = np.arange(first, min(first + 32, count))[:, None] - other = np.arange(count)[None, :] - a = domain.starts[first : first + 32, None] - v = domain.delta[first : first + 32, None] - w = domain.delta[None] - d = domain.starts[None] - a - determinant = v[..., 0] * w[..., 1] - v[..., 1] * w[..., 0] - divisor = np.where(np.abs(determinant) > 1e-12, determinant, 1.0) - t = (d[..., 0] * w[..., 1] - d[..., 1] * w[..., 0]) / divisor - u = (d[..., 0] * v[..., 1] - d[..., 1] * v[..., 0]) / divisor - adjacent = (other == following[ids]) | (ids == following[other]) - if np.any( - (other > ids) - & ~adjacent - & (np.abs(determinant) > 1e-12) - & (t > 1e-9) - & (t < 1 - 1e-9) - & (u > 1e-9) - & (u < 1 - 1e-9) - ): - return False - return all( - g.points_inside_contours(hole[:1], contours[:1])[0] for hole in contours[1:] - ) - - -def support_flat_regions(points, faces): - """Split interior boundary chords so inflated faces remain connected.""" - boundary = g.boundary_vertices_from_faces(len(points), faces) - edges = np.sort(np.concatenate((faces[:, :2], faces[:, 1:], faces[:, ::2])), axis=1) - unique, counts = np.unique(edges, axis=0, return_counts=True) - chords = unique[(counts == 2) & boundary[unique].all(axis=1)] - midpoints = {tuple(edge): len(points) + index for index, edge in enumerate(chords)} - vertices = points.tolist() + points[chords].mean(axis=1).tolist() - triangles = [] - for face in faces: - ring = [] - for index, start in enumerate(face): - end = face[(index + 1) % 3] - ring.append(int(start)) - midpoint = midpoints.get(tuple(sorted((int(start), int(end))))) - if midpoint is not None: - ring.append(midpoint) - if len(ring) > 3 or boundary[face].all(): - center = len(vertices) - vertices.append(points[face].mean(axis=0).tolist()) - triangles.extend( - (center, start, ring[(index + 1) % len(ring)]) - for index, start in enumerate(ring) - ) - else: - triangles.append(tuple(face)) - return np.asarray(vertices), np.asarray(triangles, dtype=np.int32) - - -def build_fine_mesh(alpha, spacing, threshold, *, pixel_transform=None): - """Return pixel-space triangles and component ranges for height solving.""" - alpha = np.asarray(alpha, dtype=float) - if alpha.ndim != 2 or not alpha.size or not np.isfinite(alpha).all(): - raise ValueError("The Cutout Alpha must be a finite two-dimensional array") - if not np.isfinite(spacing) or spacing <= 0: - raise ValueError("Cutout point spacing must be greater than zero") - filtered, _statistics = filter_alpha(alpha, threshold) - mask = (filtered > 0) & (filtered >= threshold) - if not mask.any(): - raise ValueError("The Cutout selection contains no visible image pixels") - reference = g.extract_selection_contours( - filtered, max_resolution=max(mask.shape), threshold=threshold - ) - # Normalize geometric tolerances and support sampling to the target length. - transform = ( - np.eye(2) if pixel_transform is None else np.asarray(pixel_transform) - ) * (32 / spacing) - inverse = np.linalg.inv(transform) - width_sizes = WidthSamples(mask, transform=transform) - spacing = 32.0 - points, faces, components = [], [], [] - offset = 0 - for component in reference.components: - contours = contours_for_component([np.asarray(c) @ transform for c in component], spacing) - spine, path_edges = interior_paths( - contours, spacing, width_sizes - ) - contours, spine, path_edges = subdivide_constraints( - contours, spine, path_edges, spacing - ) - free = relax_points(contours, spine, spacing) - vertices, triangles, constraints = triangulate( - contours, spine, path_edges, free, spacing - ) - if not len(triangles): - raise ValueError( - "Fine Outline could not triangulate a foreground component" - ) - vertices, triangles = support_flat_regions(vertices, triangles) - boundary = g.boundary_vertices_from_faces(len(vertices), triangles) - if boundary[constraints].any(): - raise ValueError( - "Fine Outline placed an interior constraint on the boundary" - ) - points.append(vertices @ inverse) - faces.append(triangles + offset) - components.append((offset, offset + len(vertices), triangles, contours[0] @ inverse)) - offset += len(vertices) - points, faces = np.concatenate(points), np.concatenate(faces) - return points, faces, g.boundary_vertices_from_faces(len(points), faces), components - - -def clean_mesh(points, faces): - faces = np.asarray(faces, dtype=np.int32).reshape(-1, 3) - triangles = points[faces] - area = np.abs( - np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0]) - ) - faces = faces[area > 1e-8] - used = np.unique(faces) - remap = np.full(len(points), -1, dtype=np.int32) - remap[used] = np.arange(len(used)) - return points[used], remap[faces] - - -def filter_alpha(alpha, threshold, min_area=16): - mask = (alpha > 0) & (alpha >= threshold) - labels, count = ndimage.label(mask) - areas = np.bincount(labels.ravel(), minlength=count + 1) - keep = areas >= min_area - keep[0] = False - if count: - keep[1 + np.argmax(areas[1:])] = True - removed = mask & ~keep[labels] - result = alpha.copy() - result[removed] = 0 - return ( - result, - { - "components_before": count, - "components_after": int(keep.sum()), - "removed_pixels": int(removed.sum()), - "foreground_pixels": int(mask.sum()), - }, - ) - - -def width_ridge_samples(image): - """Supply ridge disks for width estimation, not a connected skeleton.""" - distance = ndimage.distance_transform_edt(image) - ridge = np.zeros(image.shape, dtype=bool) - dominated = np.zeros(image.shape, dtype=bool) - for dy, dx in ((0, 1), (1, 0), (1, 1), (1, -1)): - before = np.roll(distance, (dy, dx), axis=(0, 1)) - after = np.roll(distance, (-dy, -dx), axis=(0, 1)) - ridge |= ( - (distance >= before) - & (distance >= after) - & ((distance > before) | (distance > after)) - ) - dominated |= np.maximum(before, after) >= distance + np.hypot(dx, dy) - 1e-9 - ridge &= image & ~dominated - return ridge, distance - - -class WidthSamples: - """Estimate local width from indexed distance-ridge disks.""" - - def __init__(self, mask, transform=None): - transform = np.eye(2) if transform is None else transform - self.inverse = np.linalg.inv(transform) - self.size_scale = np.linalg.svd(transform, compute_uv=False).min() - ridge, distance = width_ridge_samples(np.pad(mask, 1)) - ridge = ridge[1:-1, 1:-1] - self.distance = np.maximum(distance[1:-1, 1:-1] - 0.5, 0.3) - yy, xx = np.nonzero(ridge) - self.centers = np.column_stack((xx + 0.5, yy + 0.5)) - self.radius = self.distance[yy, xx] - self.reach_squared = (self.radius + 0.75) ** 2 - buckets = np.floor(np.log2(self.radius)).astype(int) - self.groups = [] - for bucket in np.unique(buckets)[::-1]: - ids = np.flatnonzero(buckets == bucket) - self.groups.append( - (ids, cKDTree(self.centers[ids]), float(self.radius[ids].max() + 0.75)) - ) - - def __call__(self, points): - points = points @ self.inverse - radius = np.full(len(points), 0.5) - pending = np.arange(len(points)) - for ids, tree, reach in self.groups: - if not len(pending): - break - neighbors = tree.query_ball_point(points[pending], reach) - solved = np.zeros(len(pending), dtype=bool) - for index, local in enumerate(neighbors): - if not len(local): - continue - candidates = ids[np.asarray(local)] - delta = self.centers[candidates] - points[pending[index]] - covered = ( - np.sum(delta * delta, axis=1) <= self.reach_squared[candidates] - ) - if covered.any(): - radius[pending[index]] = max( - 0.5, float(self.radius[candidates[covered]].max()) - ) - solved[index] = True - pending = pending[~solved] - distance = ndimage.map_coordinates( - self.distance, - [points[:, 1] - 0.5, points[:, 0] - 0.5], - order=1, - mode="nearest", - ) - return np.maximum( - 0.7, np.minimum(32, np.minimum(16, radius * self.size_scale * 0.9) + 0.4 * distance * self.size_scale) - ) - - -class ContourDomain: - """Reuse boundary segments for exact internal segment validation.""" - - def __init__(self, contours): - self.contours = contours - self.starts = np.concatenate(contours) - self.ends = np.concatenate([np.roll(c, -1, axis=0) for c in contours]) - self.delta = self.ends - self.starts - self.lengths = np.maximum(np.sum(self.delta**2, axis=1), 1e-20) - - def clearance(self, points): - result = np.empty(len(points)) - for first in range(0, len(points), 64): - d = points[first : first + 64, None, :] - self.starts[None, :, :] - t = np.clip(np.sum(d * self.delta, axis=2) / self.lengths, 0, 1) - result[first : first + 64] = np.sqrt( - np.min(np.sum((d - t[:, :, None] * self.delta) ** 2, axis=2), axis=1) - ) - return result - - def contains(self, points): - result = np.empty(len(points), dtype=bool) - a = self.starts - b = self.ends - dy = b[:, 1] - a[:, 1] - divisor = np.where(dy != 0, dy, 1.0) - for first in range(0, len(points), 64): - p = points[first : first + 64] - y = p[:, 1, None] - x = p[:, 0, None] - crosses = (a[:, 1] > y) != (b[:, 1] > y) - hit = a[:, 0] + (y - a[:, 1]) * (b[:, 0] - a[:, 0]) / divisor - result[first : first + 64] = ( - np.count_nonzero(crosses & (x < hit), axis=1) % 2 == 1 - ) - return result - - def valid_segments(self, starts, ends): - starts = np.asarray(starts) - ends = np.asarray(ends) - valid = self.contains(np.concatenate((starts, ends))).reshape(2, -1).all(axis=0) - valid &= np.minimum(self.clearance(starts), self.clearance(ends)) > 1e-05 - for first in range(0, len(starts), 32): - a = starts[first : first + 32, None, :] - v = (ends - starts)[first : first + 32, None, :] - d = self.starts[None, :, :] - a - w = self.delta[None, :, :] - det = v[:, :, 0] * w[:, :, 1] - v[:, :, 1] * w[:, :, 0] - nonparallel = np.abs(det) > 1e-12 - divisor = np.where(nonparallel, det, 1.0) - t = (d[:, :, 0] * w[:, :, 1] - d[:, :, 1] * w[:, :, 0]) / divisor - u = (d[:, :, 0] * v[:, :, 1] - d[:, :, 1] * v[:, :, 0]) / divisor - crossing = (nonparallel & (t >= 0) & (t <= 1) & (u >= 0) & (u <= 1)).any( - axis=1 - ) - factors = np.clip( - np.sum(d * v, axis=2) / np.maximum(np.sum(v * v, axis=2), 1e-20), 0, 1 - ) - distance = np.min( - np.sum((d - factors[:, :, None] * v) ** 2, axis=2), axis=1 - ) - valid[first : first + 32] &= ~crossing & (distance > 1e-10) - return valid - - -def contours_for_component(component, spacing, strength=1.0): - result = [] - tolerance = ( - 2.5 if spacing > 16 else 1.25 if spacing > 8 else 0.625 if spacing > 4 else 0.5 - ) - for contour in component: - original = np.asarray(contour, dtype=float) - if spacing <= 16 and strength > 0: - original = ndimage.gaussian_filter1d( - g.resample_closed_contour(original, 1.0), - 0.8 * strength, - axis=0, - mode="wrap", - ) - for error in (tolerance * strength, tolerance * strength * 0.5, 0): - simplified = approximate_polygon(np.vstack((original, original[0])), error)[ - :-1 - ] - if len(simplified) >= 3 and abs(g.contour_area(simplified)) > 0.05: - break - points = [] - for a, b in zip(simplified, np.roll(simplified, -1, axis=0)): - count = max(1, int(np.ceil(np.linalg.norm(b - a) / 32))) - points.extend(a + np.arange(count)[:, None] / count * (b - a)) - result.append(np.asarray(points)) - if strength > 0 and not valid_contours(result): - return contours_for_component( - component, spacing, strength * 0.5 if strength > 0.125 else 0 - ) - return result - - -def interior_paths(contours, spacing, width_sizes): - p, f, _ = g._triangulate_contour_component(contours, np.empty((0, 2)), spacing) - p, f = clean_mesh(p, f) - users = {} - for ti, face in enumerate(f): - for i in range(3): - edge = tuple(sorted((int(face[i]), int(face[(i + 1) % 3])))) - users.setdefault(edge, []).append(ti) - nodes = [] - connections = [] - midpoints = {} - for edge, owners in users.items(): - if len(owners) == 2: - midpoints[edge] = len(nodes) - nodes.append(p[list(edge)].mean(axis=0)) - for face in f: - mids = [ - midpoints[e] - for i in range(3) - if (e := tuple(sorted((int(face[i]), int(face[(i + 1) % 3]))))) in midpoints - ] - if len(mids) == 2: - connections.append(tuple(mids)) - else: - center = len(nodes) - nodes.append(p[face].mean(axis=0)) - connections.extend(((m, center) for m in mids)) - if not nodes: - return (np.empty((0, 2)), []) - nodes = np.asarray(nodes) - domain = ContourDomain(contours) - radius = domain.clearance(nodes) - local_sizes = width_sizes(nodes) - connections = [(a, b) for a, b in connections if min(radius[a], radius[b]) > 0.0001] - if connections: - pairs = np.asarray(connections) - valid = domain.valid_segments(nodes[pairs[:, 0]], nodes[pairs[:, 1]]) - connections = [edge for edge, ok in zip(connections, valid) if ok] - full_graph = {} - for a, b in connections: - length = float(np.linalg.norm(nodes[a] - nodes[b])) - full_graph.setdefault(a, []).append((b, length)) - full_graph.setdefault(b, []).append((a, length)) - selected = { - tuple(sorted((a, b))) - for a, b in connections - if min(radius[a], radius[b]) < spacing * 0.4 - and min(local_sizes[a], local_sizes[b]) < min(16, spacing * 0.6) - } - if selected: - parents = {i: i for edge in selected for i in edge} - - def root(i): - while parents[i] != i: - parents[i] = parents[parents[i]] - i = parents[i] - return i - - for a, b in selected: - parents[root(a)] = root(b) - owners = {i: root(i) for i in parents} - distances = {i: 0.0 for i in parents} - previous = {} - queue = [(0.0, i) for i in parents] - heapq.heapify(queue) - while queue: - distance, node = heapq.heappop(queue) - if distance != distances[node]: - continue - for neighbor, length in full_graph.get(node, []): - if max(radius[node], radius[neighbor]) > spacing * 0.75: - continue - value = distance + length - if value < distances.get(neighbor, np.inf): - distances[neighbor] = value - owners[neighbor] = owners[node] - previous[neighbor] = node - heapq.heappush(queue, (value, neighbor)) - candidates = [] - for a, b in connections: - if ( - a in owners - and b in owners - and (owners[a] != owners[b]) - and (max(radius[a], radius[b]) <= spacing * 0.75) - ): - candidates.append( - ( - distances[a] - + np.linalg.norm(nodes[a] - nodes[b]) - + distances[b], - a, - b, - ) - ) - for _, a, b in sorted(candidates): - ra, rb = (root(owners[a]), root(owners[b])) - if ra == rb: - continue - parents[ra] = rb - selected.add(tuple(sorted((a, b)))) - for cursor in (a, b): - while cursor in previous: - other = previous[cursor] - selected.add(tuple(sorted((cursor, other)))) - cursor = other - connections = sorted(selected) - graph = {} - for a, b in connections: - graph.setdefault(a, set()).add(b) - graph.setdefault(b, set()).add(a) - used = set() - paths = [] - starts = [i for i in graph if len(graph[i]) != 2] + list(graph) - for start in starts: - for neighbor in graph[start]: - edge = tuple(sorted((start, neighbor))) - if edge in used: - continue - used.add(edge) - ids = [start, neighbor] - while len(graph[ids[-1]]) == 2: - nxt = next((x for x in graph[ids[-1]] if x != ids[-2])) - edge = tuple(sorted((ids[-1], nxt))) - if edge in used: - break - used.add(edge) - ids.append(nxt) - path = nodes[ids] - length = np.linalg.norm(np.diff(path, axis=0), axis=1).sum() - if length < 5 and (len(graph[ids[0]]) == 1 or len(graph[ids[-1]]) == 1): - continue - simplified = path - for tol in (2.0, 1.0, 0.25, 0): - candidate = approximate_polygon(path, tol) - if domain.valid_segments(candidate[:-1], candidate[1:]).all(): - simplified = candidate - break - expanded = [] - for a, b in zip(simplified, simplified[1:]): - n = max(1, int(np.ceil(np.linalg.norm(b - a) / 24))) - expanded.extend(a + np.arange(n)[:, None] / n * (b - a)) - expanded.append(simplified[-1]) - paths.append(np.asarray(expanded)) - points = [] - edges = [] - lookup = {} - for path in paths: - ids = [] - for p in path: - key = tuple(np.round(p, 7)) - if key not in lookup: - lookup[key] = len(points) - points.append(p) - ids.append(lookup[key]) - edges.extend(zip(ids, ids[1:])) - return (np.asarray(points).reshape(-1, 2), edges) - - -def relax_points(contours, spine, spacing): - free = g.uniform_interior_points_for_contours(contours, spacing) - if len(free) and len(spine): - free = free[cKDTree(spine).query(free)[0] > spacing * 0.55] - step = max(1, int(spacing / 4)) - phase = step // 2 - points = np.concatenate(contours) - minimum = points.min(axis=0) - maximum = points.max(axis=0) - xx = np.arange(minimum[0] + phase, maximum[0], step) + 0.5 - yy = np.arange(minimum[1] + phase, maximum[1], step) + 0.5 - xx = xx[(xx >= minimum[0]) & (xx <= maximum[0])] - yy = yy[(yy >= minimum[1]) & (yy <= maximum[1])] - xx, yy = np.meshgrid(xx, yy) - samples = np.column_stack((xx.ravel(), yy.ravel())) - samples = samples[g.points_inside_contours(samples, contours)] - if not len(samples): - return free - fixed = np.concatenate((*contours, spine)) - for iteration in range(12): - if not len(free): - break - sites = np.concatenate((fixed, free)) - labels = cKDTree(sites).query(samples)[1] - len(fixed) - next_free = free.copy() - valid = labels >= 0 - owners = labels[valid] - owned = samples[valid] - counts = np.bincount(owners, minlength=len(free)) - used = counts > 0 - for axis in range(2): - sums = np.bincount(owners, weights=owned[:, axis], minlength=len(free)) - next_free[used, axis] = sums[used] / counts[used] - proposed = free * 0.3 + next_free * 0.7 - usable = g.points_inside_contours(proposed, contours) - movement = np.linalg.norm(proposed[usable] - free[usable], axis=1) - free[usable] = proposed[usable] - if iteration >= 3 and ( - not len(movement) or np.sqrt(np.mean(movement**2)) < spacing * 0.01 - ): - break - return free - - -def subdivide_constraints(contours, spine, path_edges, spacing): - dense = [] - for c in contours: - points = [] - for a, b in zip(c, np.roll(c, -1, axis=0)): - count = max(1, int(np.ceil(np.linalg.norm(b - a) / spacing))) - points.extend(a + np.arange(count)[:, None] / count * (b - a)) - dense.append(np.asarray(points)) - points = spine.tolist() - edges = [] - for a, b in path_edges: - count = max(1, int(np.ceil(np.linalg.norm(spine[b] - spine[a]) / spacing))) - ids = [a] - for i in range(1, count): - ids.append(len(points)) - points.append((spine[a] + i / count * (spine[b] - spine[a])).tolist()) - ids.append(b) - edges.extend(zip(ids, ids[1:])) - return (dense, np.asarray(points).reshape(-1, 2), edges) - - -def triangulate(contours, spine, path_edges, free, spacing): - boundary_count = sum(map(len, contours)) - p = np.concatenate((*contours, spine, free)) - edges = [] - offset = 0 - for c in contours: - edges.extend(((offset + i, offset + (i + 1) % len(c)) for i in range(len(c)))) - offset += len(c) - constraint_start = len(edges) - edges.extend(((boundary_count + a, boundary_count + b) for a, b in path_edges)) - result = delaunay_2d_cdt([Vector(v) for v in p], edges, [], 0, 1e-07, True) - points = np.asarray([tuple(v) for v in result[0]]) - # Restore input precision before testing area and containment. Float32 CDT - # coordinates can turn collinear outline samples into artificial slivers. - for vertex, originals in enumerate(result[3]): - if originals: - points[vertex] = p[originals].mean(axis=0) - faces = np.asarray(result[2], dtype=np.int32) - faces = faces[g.points_inside_contours(points[faces].mean(axis=1), contours)] - area = np.abs( - np.cross( - points[faces[:, 1]] - points[faces[:, 0]], - points[faces[:, 2]] - points[faces[:, 0]], - ) - ) - faces = faces[area > 1e-08] - used = np.unique(faces) - remap = np.full(len(points), -1, dtype=np.int32) - remap[used] = np.arange(len(used)) - face_edges = { - tuple(sorted((int(face[i]), int(face[(i + 1) % 3])))) - for face in faces - for i in range(3) - } - path_result = [] - constraint_graphs = [{} for _ in path_edges] - for edge, originals in zip(result[1], result[4]): - if tuple(sorted(edge)) not in face_edges: - continue - for original in originals: - if constraint_start <= original < len(edges): - graph = constraint_graphs[original - constraint_start] - a, b = edge - graph.setdefault(a, set()).add(b) - graph.setdefault(b, set()).add(a) - path_result.append(tuple(remap[list(edge)])) - # CDT coordinates are float32. Validate endpoint connectivity via provenance, - # so short constraints at large pixel coordinates survive rounding. - original_vertices = { - original: vertex - for vertex, originals in enumerate(result[3]) - for original in originals - } - for (a, b), graph in zip(edges[constraint_start:], constraint_graphs): - start = original_vertices.get(a) - end = original_vertices.get(b) - seen = set() - pending = [start] if start in graph else [] - while pending: - vertex = pending.pop() - if vertex in seen: - continue - seen.add(vertex) - pending.extend(graph[vertex] - seen) - if end not in seen: - raise ValueError( - "Fine Outline lost an interior constraint during triangulation" - ) - return ( - points[used], - remap[faces], - np.asarray(path_result, dtype=np.int32).reshape(-1, 2), - ) diff --git a/src/anyimage/operators/cutout_tool/geometry.py b/src/anyimage/operators/cutout_tool/geometry.py index 45fc0c0..aad0be9 100644 --- a/src/anyimage/operators/cutout_tool/geometry.py +++ b/src/anyimage/operators/cutout_tool/geometry.py @@ -380,88 +380,6 @@ def points_inside_contours(points, contours): return inside -def uniform_interior_points_for_contours(contours, spacing): - import numpy as np - - contours = [np.asarray(contour, dtype=np.float64) for contour in contours] - boundary_points = np.concatenate(contours, axis=0) - minimum = boundary_points.min(axis=0) - maximum = boundary_points.max(axis=0) - row_spacing = spacing * math.sqrt(3.0) * 0.5 - row_count = int(math.ceil((maximum[1] - minimum[1]) / row_spacing)) + 1 - column_count = int(math.ceil((maximum[0] - minimum[0]) / spacing)) + 1 - rows = np.arange(row_count, dtype=np.float64) - columns = np.arange(column_count, dtype=np.float64) - y = minimum[1] + rows * row_spacing - x = ( - minimum[0] - + columns[None, :] * spacing - + (rows.astype(np.int64) % 2)[:, None] * spacing * 0.5 - ) - x = np.broadcast_to(x, (row_count, column_count)) - points = np.column_stack((x.ravel(), np.repeat(y, column_count))) - inside = points_inside_contours(points, contours).reshape(row_count, column_count) - - clearance = spacing * 0.4 - near_boundary = np.zeros_like(inside) - for contour in contours: - for start, end in zip(contour, np.roll(contour, -1, axis=0)): - segment = end - start - length_squared = float(segment.dot(segment)) - if length_squared <= 1e-16: - continue - first_row = max( - 0, - int( - math.ceil( - (min(start[1], end[1]) - clearance - minimum[1]) / row_spacing - ) - ), - ) - last_row = min( - row_count - 1, - int( - math.floor( - (max(start[1], end[1]) + clearance - minimum[1]) / row_spacing - ) - ), - ) - for row in range(first_row, last_row + 1): - offset = spacing * 0.5 if row % 2 else 0.0 - first_column = max( - 0, - int( - math.ceil( - (min(start[0], end[0]) - clearance - minimum[0] - offset) - / spacing - ) - ), - ) - last_column = min( - column_count - 1, - int( - math.floor( - (max(start[0], end[0]) + clearance - minimum[0] - offset) - / spacing - ) - ), - ) - if last_column < first_column: - continue - selected = np.arange(first_column, last_column + 1) - row_points = np.column_stack( - (x[row, selected], np.full(len(selected), y[row])) - ) - factor = np.clip( - (row_points - start) @ segment / length_squared, 0.0, 1.0 - ) - nearest = start + factor[:, None] * segment - close = np.sum((row_points - nearest) ** 2, axis=1) < clearance**2 - near_boundary[row, selected[close]] = True - usable = inside & ~near_boundary - return np.column_stack((x[usable], np.broadcast_to(y[:, None], x.shape)[usable])) - - def metric_transform(world_metric=None): import numpy as np @@ -478,82 +396,13 @@ def metric_transform(world_metric=None): raise ValueError("The Image Empty has a degenerate world transform") from error -def prepare_planar_contours( - contours, - image_size, - local_bounds, - spacing=POINT_SPACING, - world_metric=None, -): - import numpy as np - - if spacing <= 0.0: - raise ValueError("Cutout point spacing must be greater than zero") - image_width, image_height = image_size - x_min, x_max, y_min, y_max = local_bounds - local_width = x_max - x_min - local_height = y_max - y_min - transform = metric_transform(world_metric) - prepared = [] - for contour in contours: - pixels = np.asarray(contour, dtype=np.float64) - if pixels.ndim != 2 or pixels.shape[0] < 3 or pixels.shape[1] != 2: - continue - local = np.empty_like(pixels) - local[:, 0] = x_min + local_width * pixels[:, 0] / image_width - local[:, 1] = y_max - local_height * pixels[:, 1] / image_height - planar = local @ transform - perimeter = float( - np.linalg.norm(np.roll(planar, -1, axis=0) - planar, axis=1).sum() - ) - if perimeter <= 1e-12: - continue - contour_spacing = min(spacing, perimeter / 8.0) - smoothed = smooth_selection_contour(planar, contour_spacing) - smoothed_local = smoothed @ np.linalg.inv(transform) - smoothed_local[:, 0] = np.clip(smoothed_local[:, 0], x_min, x_max) - smoothed_local[:, 1] = np.clip(smoothed_local[:, 1], y_min, y_max) - prepared.append(smoothed_local @ transform) - if not prepared: - raise ValueError("The Cutout has no selection large enough for this resolution") - interior = uniform_interior_points_for_contours(prepared, spacing) - return prepared, interior, spacing - - -def _prepare_selection_contours( - selection_contours, - local_bounds, - spacing=POINT_SPACING, - world_metric=None, -): - prepared = [] - errors = [] - for component in selection_contours.components: - try: - contours, interior, _spacing = prepare_planar_contours( - component, - selection_contours.image_size, - local_bounds, - spacing=spacing, - world_metric=world_metric, - ) - except ValueError as error: - errors.append(error) - continue - prepared.append((contours, interior)) - if not prepared: - if errors: - raise errors[0] - raise ValueError("The Cutout has no usable selection") - return prepared, spacing - - -def _triangulate_contour_component(contours, interior, spacing): +def contour_interior_point(contours, spacing): + """Choose an interior seed from the largest auxiliary triangle.""" import numpy as np from mathutils import Vector from mathutils.geometry import delaunay_2d_cdt - input_points = np.concatenate((*contours, interior), axis=0) + input_points = np.concatenate(contours, axis=0) boundary_edges = [] offset = 0 for contour in contours: @@ -583,15 +432,15 @@ def _triangulate_contour_component(contours, interior, spacing): faces = faces[points_inside_contours(centers, contours)] if not len(faces): raise ValueError("Unable to triangulate the Cutout selection components") - used = np.unique(faces) - if len(used) != len(planar_points): - remap = np.full(len(planar_points), -1, dtype=np.int32) - remap[used] = np.arange(len(used), dtype=np.int32) - planar_points = planar_points[used] - faces = remap[faces] - boundary = boundary_vertices_from_faces(len(planar_points), faces) - planar_faces = [tuple(int(value) for value in face) for face in faces] - return planar_points, planar_faces, boundary + triangles = planar_points[faces] + first = triangles[:, 1] - triangles[:, 0] + second = triangles[:, 2] - triangles[:, 0] + areas = abs(first[:, 0] * second[:, 1] - first[:, 1] * second[:, 0]) + seed = triangles[np.argmax(areas)].mean(axis=0) + if not points_inside_contours(seed[None], contours)[0]: + raise ValueError("Cutout seed is outside its contour") + return seed + def boundary_vertices_from_faces(vertex_count, faces): @@ -614,54 +463,6 @@ def boundary_vertices_from_faces(vertex_count, faces): return boundary -def triangulate_selection_contours( - selection_contours, - local_bounds, - spacing=POINT_SPACING, - world_metric=None, - return_components=False, -): - import numpy as np - - prepared, spacing = _prepare_selection_contours( - selection_contours, - local_bounds, - spacing=spacing, - world_metric=world_metric, - ) - points = [] - faces = [] - boundaries = [] - components = [] - offset = 0 - for contours, interior in prepared: - component_points, component_faces, component_boundary = _triangulate_contour_component( - contours, - interior, - spacing, - ) - components.append( - ( - offset, - offset + len(component_points), - component_faces, - contours[0], - ) - ) - points.append(component_points) - faces.extend(tuple(offset + value for value in face) for face in component_faces) - boundaries.append(component_boundary) - offset += len(component_points) - if not points: - raise ValueError("Unable to triangulate the Cutout selection components") - result = ( - np.concatenate(points, axis=0), - faces, - np.concatenate(boundaries), - ) - return (*result, components) if return_components else result - - CUTOUT_ALPHA_THRESHOLD = 0.9 CUTOUT_CLEANUP_UI_SCALE = 100.0 @@ -674,84 +475,26 @@ class BaseShape: balloon: object -def refine_balloon_triangles(points, faces, boundary): - """Add interior samples to boundary ears and interior boundary chords.""" - from collections import Counter - import numpy as np - - edges = Counter(tuple(sorted((face[i], face[(i + 1) % 3]))) for face in faces for i in range(3)) - values, midpoints = list(points), {} - for edge, count in edges.items(): - if count == 2 and all(boundary[index] for index in edge): - midpoints[edge] = len(values) - values.append(points[list(edge)].mean(axis=0)) - refined = [] - for face in faces: - ring = [] - for index, vertex in enumerate(face): - ring.append(vertex) - midpoint = midpoints.get(tuple(sorted((vertex, face[(index + 1) % 3])))) - if midpoint is not None: - ring.append(midpoint) - if len(ring) == 3 and not all(boundary[vertex] for vertex in face): - refined.append(tuple(face)) - continue - center = len(values) - values.append(points[list(face)].mean(axis=0)) - refined.extend((center, ring[i], ring[(i + 1) % len(ring)]) for i in range(len(ring))) - points = np.asarray(values) - return points, refined, boundary_vertices_from_faces(len(points), refined) - - -def cutout_domain( - selection_contours, - bounds, - spacing, - world_metric, - balloon_algorithm="POISSON", - *, - fine_mesh=None, -): +def cutout_domain(mesh, bounds, world_metric, balloon_algorithm="POISSON"): + """Map the quality mesh into local coordinates and solve its height.""" import numpy as np x_min, x_max, y_min, y_max = bounds + pixels, planar_faces, boundary, pixel_components, image_size = mesh + image_width, image_height = image_size + transform = metric_transform(world_metric) - if fine_mesh is None: - planar, planar_faces, boundary, components = triangulate_selection_contours( - selection_contours, - bounds, - spacing=spacing, - world_metric=world_metric, - return_components=True, - ) - else: - pixels, planar_faces, boundary, pixel_components, image_size = fine_mesh - image_width, image_height = image_size - transform = metric_transform(world_metric) + def pixel_to_planar(points): + local = np.asarray(points, dtype=np.float64).copy() + local[:, 0] = x_min + local[:, 0] * (x_max - x_min) / image_width + local[:, 1] = y_max - local[:, 1] * (y_max - y_min) / image_height + return local @ transform - def pixel_to_planar(points): - local = np.asarray(points, dtype=np.float64).copy() - local[:, 0] = x_min + local[:, 0] * (x_max - x_min) / image_width - local[:, 1] = y_max - local[:, 1] * (y_max - y_min) / image_height - return local @ transform - - planar = pixel_to_planar(pixels) - planar_faces = planar_faces.tolist() - components = [ - (start, end, faces.tolist(), pixel_to_planar(contour)) - for start, end, faces, contour in pixel_components - ] - refined_points, refined_faces, refined_components = [], [], [] - offset = 0 - for start, end, component_faces, contour in components: - points, faces, _ = refine_balloon_triangles(planar[start:end], component_faces, boundary[start:end]) - refined_points.append(points) - refined_faces.extend(tuple(offset + index for index in face) for face in faces) - refined_components.append((offset, offset + len(points), faces, contour)) - offset += len(points) - planar = np.concatenate(refined_points) - planar_faces, components = refined_faces, refined_components - boundary = boundary_vertices_from_faces(len(planar), planar_faces) + planar = pixel_to_planar(pixels) + components = [ + (start, end, faces.tolist(), pixel_to_planar(contour)) + for start, end, faces, contour in pixel_components + ] profile = balloon.balloon_profile( balloon_algorithm, planar, planar_faces, boundary, components ) @@ -836,7 +579,11 @@ def extract_cutout_contours( def effective_cutout_edge_length( - context, source_object, edge_length, content_values, content_bounds, + context, + source_object, + edge_length, + content_values, + content_bounds, threshold=CUTOUT_ALPHA_THRESHOLD, ): """Resolve target spacing in Blender world units from visible content bounds.""" @@ -850,16 +597,25 @@ def effective_cutout_edge_length( columns = np.flatnonzero(visible.any(axis=0)) if not len(rows): raise ValueError("The Cutout selection contains no visible image pixels") - bounds = crop_plane_bounds(source_object, content_bounds, tuple(source_object.data.size)) + bounds = crop_plane_bounds( + source_object, content_bounds, tuple(source_object.data.size) + ) metric, _ = world_plane_metric(source_object.matrix_world) metric_transform(metric) - width = (bounds[1] - bounds[0]) * (columns.max() + 1 - columns.min()) / values.shape[1] + width = ( + (bounds[1] - bounds[0]) * (columns.max() + 1 - columns.min()) / values.shape[1] + ) height = (bounds[3] - bounds[2]) * (rows.max() + 1 - rows.min()) / values.shape[0] longest = max(width * math.sqrt(metric[0][0]), height * math.sqrt(metric[1][1])) - return max(float(edge_length), longest * configured_min_cutout_relative_edge_length(context)) + return max( + float(edge_length), + longest * configured_min_cutout_relative_edge_length(context), + ) -def selection_calibration(content_values, target_shape, threshold=CUTOUT_ALPHA_THRESHOLD): +def selection_calibration( + content_values, target_shape, threshold=CUTOUT_ALPHA_THRESHOLD +): """Return Cutout validity and transparency ratio at geometry resolution.""" import numpy as np @@ -900,39 +656,27 @@ def build_base_shape( world_metric, world_z_scale = world_plane_metric(source_object.matrix_world) if world_z_scale <= 1e-10: raise ValueError("The Image Empty has a degenerate world transform") - pixel_transform = np.diag(( - (bounds[1] - bounds[0]) / content_values.shape[1], - -(bounds[3] - bounds[2]) / content_values.shape[0], - )) @ metric_transform(world_metric) - longest = max( - np.linalg.norm(pixel_transform[0]) * content_values.shape[1], - np.linalg.norm(pixel_transform[1]) * content_values.shape[0], - ) - fine_mesh = None - selection_contours = None - if fine_outline: - from .fine_mesh import build_fine_mesh + pixel_transform = np.diag( + ( + (bounds[1] - bounds[0]) / content_values.shape[1], + -(bounds[3] - bounds[2]) / content_values.shape[0], + ) + ) @ metric_transform(world_metric) + from .mesh import build_mesh - fine_mesh = ( - *build_fine_mesh(content_values, edge_length, alpha_threshold, pixel_transform=pixel_transform), - (content_values.shape[1], content_values.shape[0]), - ) - else: - selection_contours = extract_cutout_contours( + mesh = ( + *build_mesh( content_values, - longest / edge_length, - threshold=alpha_threshold, - ) + edge_length, + alpha_threshold, + pixel_transform=pixel_transform, + fine_outline=fine_outline, + ), + (content_values.shape[1], content_values.shape[0]), + ) preferences = addon_preferences(context) balloon_algorithm = getattr(preferences, "balloon_profile_algorithm", "POISSON") - domain = cutout_domain( - selection_contours, - bounds, - edge_length, - world_metric, - balloon_algorithm, - fine_mesh=fine_mesh, - ) + domain = cutout_domain(mesh, bounds, world_metric, balloon_algorithm) source_vertices, faces, uv, _boundary, planar, balloon = domain return BaseShape( vertices=np.asarray(source_vertices, dtype=np.float64), diff --git a/src/anyimage/operators/cutout_tool/mesh.py b/src/anyimage/operators/cutout_tool/mesh.py new file mode 100644 index 0000000..7851a96 --- /dev/null +++ b/src/anyimage/operators/cutout_tool/mesh.py @@ -0,0 +1,426 @@ +"""Cutout contours, quality triangulation and local height support.""" + +import numpy as np +from scipy import ndimage +from scipy.spatial import cKDTree +from .polygon import approximate_polygon + +from . import geometry as g + + +# Keep predicate tolerances independent of scene units and requested density. +NORMALIZED_EDGE_LENGTH = 32.0 + + +def separate_contour_contacts(contours): + """Open pixel-corner contacts with a subpixel bevel in each incident ring.""" + points = np.concatenate(contours) + pairs = cKDTree(points).query_pairs(1e-8) + if not pairs: + return contours + contacts = {i for pair in pairs for i in pair} + result = [] + offset = 0 + for contour in contours: + ring = contour.copy() + for i in range(len(ring)): + if offset + i in contacts: + before = contour[i-1] - contour[i] + after = contour[(i+1) % len(ring)] - contour[i] + length = min(np.linalg.norm(before), np.linalg.norm(after)) + direction = ( + before / np.linalg.norm(before) + after / np.linalg.norm(after) + ) + ring[i] += direction * length * 0.125 + result.append(ring) + offset += len(ring) + return result + + +def support_mesh(points, faces, *, max_vertices=50000): + """Split boundary chords once so every face has an interior height sample.""" + boundary = g.boundary_vertices_from_faces(len(points), faces) + edges, counts = np.unique( + np.sort(np.concatenate((faces[:, :2], faces[:, 1:], faces[:, ::2])), axis=1), + axis=0, return_counts=True, + ) + chords = edges[(counts == 2) & boundary[edges].all(axis=1)] + midpoints = {tuple(edge): len(points) + i for i, edge in enumerate(chords)} + isolated = sum( + boundary[face].all() and not any( + tuple(sorted((face[i], face[(i+1) % 3]))) in midpoints for i in range(3) + ) for face in faces + ) + if len(points) + len(chords) + isolated > max_vertices: + raise ValueError("Cutout support vertex limit exceeded") + if not len(chords) and not isolated: + return points, faces + vertices = points.tolist() + points[chords].mean(axis=1).tolist() + output = [] + for face in faces: + pieces = [list(face)] + for i in range(3): + a, b = face[i], face[(i+1) % 3] + midpoint = midpoints.get(tuple(sorted((a, b)))) + if midpoint is None: + continue + split = [] + for piece in pieces: + if a in piece and b in piece: + k = next( + j for j in range(3) + if piece[j] in (a, b) and piece[(j+1) % 3] in (a, b) + ) + x, y, z = piece[k], piece[(k+1) % 3], piece[(k+2) % 3] + split.extend(((x, midpoint, z), (midpoint, y, z))) + else: + split.append(piece) + pieces = split + if len(pieces) == 1 and boundary[face].all(): + midpoint = len(vertices) + vertices.append(points[face].mean(axis=0).tolist()) + pieces = [(face[i], face[(i+1) % 3], midpoint) for i in range(3)] + output.extend(pieces) + return np.asarray(vertices), np.asarray(output, dtype=np.int32) + + +def valid_contours(contours): + """Reject crossing rings and holes moved outside the outer ring.""" + starts = np.concatenate(contours) + if cKDTree(starts).query_pairs(1e-8): + return False + ends = np.concatenate([np.roll(c, -1, axis=0) for c in contours]) + delta = ends - starts + count = len(starts) + following = np.concatenate( + [ + np.roll(np.arange(len(contour)) + sum(map(len, contours[:index])), -1) + for index, contour in enumerate(contours) + ] + ) + for first in range(0, count, 32): + ids = np.arange(first, min(first + 32, count))[:, None] + other = np.arange(count)[None, :] + a = starts[first : first + 32, None] + v = delta[first : first + 32, None] + w = delta[None] + d = starts[None] - a + determinant = v[..., 0] * w[..., 1] - v[..., 1] * w[..., 0] + divisor = np.where(np.abs(determinant) > 1e-12, determinant, 1.0) + t = (d[..., 0] * w[..., 1] - d[..., 1] * w[..., 0]) / divisor + u = (d[..., 0] * v[..., 1] - d[..., 1] * v[..., 0]) / divisor + adjacent = (other == following[ids]) | (ids == following[other]) + if np.any( + (other > ids) + & ~adjacent + & (np.abs(determinant) > 1e-12) + & (t >= -1e-9) + & (t <= 1 + 1e-9) + & (u >= -1e-9) + & (u <= 1 + 1e-9) + ): + return False + return all( + g.points_inside_contours(hole[:1], contours[:1])[0] for hole in contours[1:] + ) + + +def prepare_mask(alpha, threshold, transform, fine_outline): + """Fill holes according to the normalized world-space target length.""" + filtered = filter_alpha(alpha, threshold) + mask = (filtered > 0) & (filtered >= threshold) + filled = ndimage.binary_fill_holes(mask) + if not fine_outline: + return filled + holes, _ = ndimage.label(filled & ~mask) + areas = np.bincount(holes.ravel()) * abs(np.linalg.det(transform)) + for index, slices in enumerate(ndimage.find_objects(holes), 1): + if slices is None or areas[index] > 0.04 * NORMALIZED_EDGE_LENGTH**2: + continue + y, x = slices + corners = ( + np.array( + ( + (x.start, y.start), + (x.stop, y.start), + (x.stop, y.stop), + (x.start, y.stop), + ) + ) + @ transform + ) + if np.ptp(corners, axis=0).max() <= 0.25 * NORMALIZED_EDGE_LENGTH: + mask[slices] |= holes[slices] == index + return mask + + +def build_mesh(alpha, spacing, threshold, *, pixel_transform=None, fine_outline=True): + """Return pixel-space quality triangles and component ranges.""" + alpha = np.asarray(alpha, dtype=float) + if alpha.ndim != 2 or not alpha.size or not np.isfinite(alpha).all(): + raise ValueError("The Cutout Alpha must be a finite two-dimensional array") + if not np.isfinite(spacing) or spacing <= 0: + raise ValueError("Cutout point spacing must be greater than zero") + transform = ( + np.eye(2) if pixel_transform is None else np.asarray(pixel_transform) + ) * (NORMALIZED_EDGE_LENGTH / spacing) + inverse = np.linalg.inv(transform) + mask = prepare_mask(alpha, threshold, transform, fine_outline) + if not mask.any(): + raise ValueError("The Cutout selection contains no visible image pixels") + reference = g.extract_selection_contours( + mask.astype(float), max_resolution=max(mask.shape), threshold=0.5 + ) + prepared = [] + for component in reference.components: + raw = separate_contour_contacts([np.asarray(c) @ transform for c in component]) + if fine_outline: + contours = contours_for_component(raw, NORMALIZED_EDGE_LENGTH, pixel_transform=transform) + else: + contours = [ + g.smooth_selection_contour( + c, + min( + NORMALIZED_EDGE_LENGTH, + np.linalg.norm(np.roll(c, -1, axis=0) - c, axis=1).sum() / 8, + ), + ) + for c in raw + ] + contours = [ + np.clip(c @ inverse, (0, 0), mask.shape[::-1]) @ transform + for c in contours + ] + if not valid_contours(contours): + raise ValueError("Cutout simplified outline is invalid") + prepared.append(contours) + points, faces, components = [], [], [] + offset = 0 + for contours in prepared: + vertices, triangles = triangulate(contours, NORMALIZED_EDGE_LENGTH) + vertices, triangles = support_mesh( + vertices, triangles, max_vertices=50000-offset + ) + points.append(vertices @ inverse) + faces.append(triangles + offset) + components.append( + (offset, offset + len(vertices), triangles, contours[0] @ inverse) + ) + offset += len(vertices) + if offset > 50000: + raise ValueError("Cutout mesh vertex limit exceeded") + points, faces = np.concatenate(points), np.concatenate(faces) + return points, faces, g.boundary_vertices_from_faces(len(points), faces), components + + +def filter_alpha(alpha, threshold, min_area=16): + mask = (alpha > 0) & (alpha >= threshold) + labels, count = ndimage.label(mask) + areas = np.bincount(labels.ravel(), minlength=count + 1) + keep = areas >= min_area + keep[0] = False + if count: + keep[1 + np.argmax(areas[1:])] = True + removed = mask & ~keep[labels] + result = alpha.copy() + result[removed] = 0 + return result + + +def contours_for_component(component, spacing, strength=1.0, *, pixel_transform=None): + """Soften pixel steps before simplifying in the target edge-length metric.""" + result = [] + transform = np.eye(2) if pixel_transform is None else pixel_transform + inverse = np.linalg.inv(transform) + tolerance = spacing * 0.078125 + for contour in component: + original = np.asarray(contour, dtype=float) + if strength > 0: + original = ndimage.gaussian_filter1d( + g.resample_closed_contour(original @ inverse, 0.5), + 1.6 * strength, + axis=0, + mode="wrap", + ) @ transform + for error in (tolerance * strength, tolerance * strength * 0.5, 0): + simplified = approximate_polygon(np.vstack((original, original[0])), error)[ + :-1 + ] + if len(simplified) >= 3 and abs(g.contour_area(simplified)) > 0.05: + break + points = [] + for a, b in zip(simplified, np.roll(simplified, -1, axis=0)): + count = max(1, int(np.ceil(np.linalg.norm(b - a) / spacing))) + points.extend(a + np.arange(count)[:, None] / count * (b - a)) + result.append(np.asarray(points)) + if not valid_contours(result): + if strength == 0: + raise ValueError("Cutout outline is invalid") + return contours_for_component( + component, spacing, strength * 0.5 if strength > 0.125 else 0, + pixel_transform=transform, + ) + return result + + +def normalize_constraints(points, edges): + """Merge coincident vertices and split constraints at existing vertices.""" + tree = cKDTree(points) + parents = np.arange(len(points)) + + def root(index): + while parents[index] != index: + parents[index] = parents[parents[index]] + index = parents[index] + return index + + for a, b in sorted(tree.query_pairs(1e-8)): + parents[root(b)] = root(a) + representatives = np.array([root(i) for i in range(len(points))]) + used, inverse = np.unique(representatives, return_inverse=True) + points = points[used] + tree = cKDTree(points) + segments = [] + seen = {} + sources = {} + for marker, (a, b) in enumerate(edges, 1): + a, b = inverse[a], inverse[b] + delta = points[b] - points[a] + length = np.linalg.norm(delta) + if length <= 1e-8: + raise ValueError("Cutout constraint collapsed during normalization") + candidates = np.array( + tree.query_ball_point((points[a] + points[b]) / 2, length / 2 + 1e-8), + dtype=int, + ) + fraction = (points[candidates] - points[a]) @ delta / length**2 + distance = np.linalg.norm( + points[candidates] - points[a] - fraction[:, None] * delta, axis=1 + ) + valid = (fraction >= -1e-10) & (fraction <= 1 + 1e-10) & (distance < 1e-8) + ids = candidates[valid][np.argsort(fraction[valid])] + for x, y in zip(ids, ids[1:]): + key = tuple(sorted((int(x), int(y)))) + if key in seen: + sources[seen[key]].add(marker) + continue + identity = len(segments) + 1 + seen[key] = identity + sources[identity] = {marker} + segments.append((int(x), int(y), identity)) + return points, segments, sources + + +def triangulate(contours, spacing): + """Build and validate a quality mesh with complete constraint chains.""" + from . import tangle + from scipy.sparse import coo_matrix + from scipy.sparse.csgraph import connected_components + + points = np.concatenate(contours) + if len(points) > 20000: + raise ValueError("Cutout constraint vertex limit exceeded") + edges = [] + offset = 0 + for contour in contours: + edges.extend( + (offset + i, offset + (i + 1) % len(contour)) for i in range(len(contour)) + ) + offset += len(contour) + input_points = points + points, segments, sources = normalize_constraints(points, edges) + points = np.vstack((points, g.contour_interior_point(contours, spacing))) + holes = [g.contour_interior_point([contour], spacing).tolist() for contour in contours[1:]] + raw_points, raw_faces, raw_segments = tangle.build_mesh( + points.tolist(), segments, holes, np.sqrt(3) / 4 * spacing**2 + ) + points = np.asarray(raw_points) + faces = np.asarray(raw_faces, dtype=np.int32).reshape(-1, 3) + if not len(faces) or not np.isfinite(points).all(): + raise ValueError("Cutout triangulation is empty or nonfinite") + used = np.unique(faces) + remap = np.full(len(points), -1, dtype=np.int32) + remap[used] = np.arange(len(used)) + points, faces = points[used], remap[faces] + triangles = points[faces] + areas = ( + abs( + np.cross( + triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0] + ) + ) + / 2 + ) + expected_area = abs(g.contour_area(contours[0])) - sum( + abs(g.contour_area(c)) for c in contours[1:] + ) + if np.any(areas <= 1e-10) or not np.isclose( + areas.sum(), expected_area, rtol=1e-8, atol=1e-7 + ): + raise ValueError("Cutout triangulation has degenerate faces or incorrect area") + face_edges, counts = np.unique( + np.sort(np.concatenate((faces[:, :2], faces[:, 1:], faces[:, ::2])), axis=1), + axis=0, + return_counts=True, + ) + if counts.max() > 2: + raise ValueError("Cutout triangulation has nonmanifold edges") + graph = coo_matrix( + (np.ones(len(face_edges)), (face_edges[:, 0], face_edges[:, 1])), + shape=(len(points), len(points)), + ).tocsr() + components, _ = connected_components(graph, directed=False) + if ( + components != 1 + or components - len(points) + len(face_edges) - len(faces) != len(contours) - 1 + ): + raise ValueError( + "Cutout triangulation has incorrect component or hole topology" + ) + actual_edges = {tuple(e) for e in face_edges} + chains = [[] for _ in edges] + for a, b, marker in raw_segments: + if marker in sources and min(remap[a], remap[b]) >= 0: + edge = tuple(sorted((int(remap[a]), int(remap[b])))) + if edge in actual_edges: + for original in sources[marker]: + chains[original - 1].append(edge) + tree = cKDTree(points) + for (a, b), chain in zip(edges, chains): + error = "Cutout lost a boundary constraint" + if not chain: + raise ValueError(error) + chain = np.asarray(chain) + start, end = input_points[[a, b]] + delta = end - start + length = np.linalg.norm(delta) + vertices = np.unique(chain) + distance = np.linalg.norm( + points[vertices] + - start + - ((points[vertices] - start) @ delta / length**2)[:, None] * delta, + axis=1, + ) + covered = np.linalg.norm( + points[chain[:, 1]] - points[chain[:, 0]], axis=1 + ).sum() + endpoint_distance, endpoints = tree.query([start, end]) + adjacency = {} + for x, y in chain: + adjacency.setdefault(x, set()).add(y) + adjacency.setdefault(y, set()).add(x) + visited = set() + pending = [endpoints[0]] + while pending: + vertex = pending.pop() + if vertex not in visited: + visited.add(vertex) + pending.extend(adjacency.get(vertex, ())) + if ( + distance.max() > 1e-6 + or endpoint_distance.max() > 1e-6 + or endpoints[1] not in visited + or not np.isclose(covered, length, rtol=1e-7, atol=1e-6) + ): + raise ValueError(error) + return points, faces diff --git a/src/anyimage/operators/cutout_tool/tangle.py b/src/anyimage/operators/cutout_tool/tangle.py new file mode 100644 index 0000000..51ab46c --- /dev/null +++ b/src/anyimage/operators/cutout_tool/tangle.py @@ -0,0 +1,2690 @@ +# Upstream: https://github.com/dcm3c/tangle/tree/66221c398cd90f06a7beb0c0abf637b269313af3 +# Original tangle.py SHA256: e53fe29bb9fe5174d1735f2d7060b1d7d14f1005666ef60bb93e27b4e0caa679 +# Local changes: memory entry point, isolated numeric context, resource limits; +# removed CLI and file IO. Numerical meshing core retained under MIT. +# SPDX-License-Identifier: MIT +# Copyright (c) 2026 David Meeker +""" +tangle +A 2D Delaunay triangulation tool compatible with Shewchuk's Triangle format. + +Author: David Meeker +Generated with the assistance of Claude Code +Translated from tangle.cpp + +Version 0.4.18 +22 Aug 2026 + +MIT License + +Copyright (c) 2026 David Meeker + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import sys +import math +import bisect +from collections import deque +from dataclasses import dataclass, field +from decimal import Decimal, localcontext +from threading import RLock +from heapq import heappush, heappop + +# Exact-predicate fallback precision. orient2d/in_circle fall back to Decimal +# when the double filter is inconclusive; the products there need ~34 (orient2d) +# to ~51 (in_circle) significant digits to be exact, well above Decimal's default +# 28 (~93 bits, which is BELOW the float128 tangle.cpp found insufficient for +# mixed-scale coords -- see f0854c1, which moved C++ to float256 ~245 bits). +# 80 digits comfortably exceeds float256's precision; it only costs in the rare +# fallback (the double filter handles >99.9% of calls). +_CALL_LOCK = RLock() + +# Base relative tolerance. EPS is scaled to the bounding-box diagonal in +# build_delaunay(); EPS_SCALE controls that ratio and is also used directly +# for dimensionless comparisons. +EPS_SCALE = 1e-10 +EPS = EPS_SCALE + +# Relative tolerance for "too close" geometry: used both by cleanPSLG +# (merge nodes within bboxDiag * CLOSE_ENOUGH) and by quality refinement +# (reduce angle threshold for triangles with all edges < bboxDiag * CLOSE_ENOUGH). +CLOSE_ENOUGH = 1e-6 + +# Reduced minimum angle (degrees) for tiny triangles in locally dense +# regions where the full minimum angle is unachievable. +TINY_TRI_ANGLE = 15.0 + +# Hard cap on the requested minimum angle (degrees). Delaunay refinement is +# only guaranteed to terminate below ~33.8 deg (Ruppert/Shewchuk); matches the +# ceiling FEMM enforces. Under shortest-edge ordering refinement stays tame +# right up to this bound. +MINANGLE_MAX_VAL = 33.8 + + +# ============================================================ +# Options +# ============================================================ + +@dataclass +class Options: + pslg: bool = False + no_poly_out: bool = False + jettison: bool = False + quality: bool = False + min_angle: float = 20.0 + edges: bool = False + regions: bool = False + area_limit: bool = False + max_area: float = -1.0 + zero_indexed: bool = False + quiet: bool = False + suppress_iter: bool = False + no_holes: bool = False + no_steiner: bool = False + neighbors: bool = False + clean_pslg: bool = False + clean_tol: float = -1.0 # -1 = auto (bboxDiag * 1e-6) + reorder: bool = False # reverse Cuthill-McKee bandwidth/profile reduction + first_index: int = 1 + + +# ============================================================ +# Geometry primitives +# ============================================================ + +@dataclass +class Point: + x: float + y: float + id: int = 0 + marker: int = 0 + attribs: list = field(default_factory=list) + lfs: float = -1.0 # local feature size constraint (-1 = none) + + +@dataclass +class Segment: + v0: int + v1: int + marker: int = 0 + lfs: float = -1.0 # per-segment LFS constraint (-1 = none) + pbc_type: int = -1 # 0=periodic, 1=anti-periodic (-1 = none) + no_split: bool = False # true for AGE chord segments (must stay uniformly spaced) + pbc_side: int = -1 # 0/1: declared side of an (anti)periodic pair; -1 = n/a + + +@dataclass +class Hole: + x: float + y: float + + +@dataclass +class Region: + x: float + y: float + attrib: float = 0.0 + max_area: float = -1.0 + + +@dataclass +class Triangle: + v: list = field(default_factory=lambda: [-1, -1, -1]) + neighbors: list = field(default_factory=lambda: [-1, -1, -1]) + region_attrib: float = 0.0 + region_max_area: float = -1.0 + marker: int = 0 + generation: int = 0 + + +def orient2d_xy(ax, ay, bx, by, cx, cy): + """orient2d on raw coordinates — avoids constructing a Point. + Same exact-arithmetic fallback as orient2d.""" + Ax, Ay = bx - ax, by - ay + Bx, By = cx - ax, cy - ay + det = Ax * By - Bx * Ay + M = max(abs(Ax), abs(Ay), abs(Bx), abs(By)) + if abs(det) > 6.662e-16 * M * M: + return det + fAx = Decimal(bx) - Decimal(ax) + fAy = Decimal(by) - Decimal(ay) + fBx = Decimal(cx) - Decimal(ax) + fBy = Decimal(cy) - Decimal(ay) + return float(fAx * fBy - fBx * fAy) + + +def orient2d(a, b, c): + Ax, Ay = b.x - a.x, b.y - a.y + Bx, By = c.x - a.x, c.y - a.y + det = Ax * By - Bx * Ay + + # Error bound: det terms are O(M²), check if result is well above fp noise + M = max(abs(Ax), abs(Ay), abs(Bx), abs(By)) + if abs(det) > 6.662e-16 * M * M: + return det + + # Fall back to exact arithmetic via Decimal + fAx = Decimal(b.x) - Decimal(a.x) + fAy = Decimal(b.y) - Decimal(a.y) + fBx = Decimal(c.x) - Decimal(a.x) + fBy = Decimal(c.y) - Decimal(a.y) + return float(fAx * fBy - fBx * fAy) + + +def in_circle(a, b, c, d): + # 3x3 formulation: subtract d to eliminate one row + Ax, Ay = a.x - d.x, a.y - d.y + Bx, By = b.x - d.x, b.y - d.y + Cx, Cy = c.x - d.x, c.y - d.y + As = Ax * Ax + Ay * Ay + Bs = Bx * Bx + By * By + Cs = Cx * Cx + Cy * Cy + + det = (Ax * (By * Cs - Cy * Bs) + - Ay * (Bx * Cs - Cx * Bs) + + As * (Bx * Cy - Cx * By)) + + # Error bound: det terms are O(M⁴), check if result is well above fp noise + M = max(abs(Ax), abs(Ay), abs(Bx), abs(By), + abs(Cx), abs(Cy)) + M4 = M * M * M * M + if abs(det) > 1.333e-14 * M4: + return det + + # Fall back to exact arithmetic via Decimal (standard 3x3 formulation) + ax = Decimal(a.x) - Decimal(d.x) + ay = Decimal(a.y) - Decimal(d.y) + bx = Decimal(b.x) - Decimal(d.x) + by = Decimal(b.y) - Decimal(d.y) + cx = Decimal(c.x) - Decimal(d.x) + cy = Decimal(c.y) - Decimal(d.y) + cs = cx * cx + cy * cy + bs = bx * bx + by * by + fdet = (ax * (by * cs - cy * bs) + - ay * (bx * cs - cx * bs) + + (ax * ax + ay * ay) * (bx * cy - cx * by)) + return float(fdet) + + +def tri_area(a, b, c): + return 0.5 * abs(orient2d(a, b, c)) + + +def vertex_angle(a, b, c): + abx, aby = b.x - a.x, b.y - a.y + acx, acy = c.x - a.x, c.y - a.y + dot = abx * acx + aby * acy + cross = abx * acy - aby * acx + return math.atan2(abs(cross), dot) * 180.0 / math.pi + + +def min_angle(a, b, c): + la = math.sqrt((b.x-c.x)*(b.x-c.x)+(b.y-c.y)*(b.y-c.y)) + lb = math.sqrt((a.x-c.x)*(a.x-c.x)+(a.y-c.y)*(a.y-c.y)) + lc = math.sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)) + if la < EPS or lb < EPS or lc < EPS: + return 0.0 + cos_a = max(-1.0, min(1.0, (lb * lb + lc * lc - la * la) / (2 * lb * lc))) + cos_b = max(-1.0, min(1.0, (la * la + lc * lc - lb * lb) / (2 * la * lc))) + cos_c = max(-1.0, min(1.0, (la * la + lb * lb - lc * lc) / (2 * la * lb))) + return min(math.acos(cos_a), math.acos(cos_b), math.acos(cos_c)) * 180.0 / math.pi + + +def circumcenter(a, b, c): + ax, ay = b.x - a.x, b.y - a.y + bx, by = c.x - a.x, c.y - a.y + D = 2.0 * (ax * by - ay * bx) + if abs(D) < EPS: + return (a.x + b.x + c.x) / 3.0, (a.y + b.y + c.y) / 3.0 + ux = (by * (ax * ax + ay * ay) - ay * (bx * bx + by * by)) / D + uy = (ax * (bx * bx + by * by) - bx * (ax * ax + ay * ay)) / D + return a.x + ux, a.y + uy + + +def point_on_segment(p, a, b): + if orient2d(a, b, p) != 0: + return False + dx, dy = b.x - a.x, b.y - a.y + if abs(dx) > EPS: + t = (p.x - a.x) / dx + elif abs(dy) > EPS: + t = (p.y - a.y) / dy + else: + return False + return -EPS_SCALE <= t <= 1 + EPS_SCALE # dimensionless: fixed tolerance + + +def edge_key(a, b): + return (min(a, b), max(a, b)) + + +# ============================================================ +# Mesh structure +# ============================================================ + +@dataclass +class PBCNodePair: + node_a: int + node_b: int + type: int # 0=periodic, 1=anti-periodic + + +class Mesh: + def __init__(self): + self.vertices = [] + self.triangles = [] + self.segments = [] + self.holes = [] + self.regions = [] + self.edges = [] + self.pbc_pairs = [] # list of PBCNodePair, filled after refinement + self.pbc_twin = {} # dict: node -> twin node (bidirectional) + self.pbc_node_type = {} # dict: node -> pbc type (0 or 1) + self.pbc_defs = [] # PBC declarations: (marker_a, marker_b, type) + + def rebuild_adjacency(self): + nt = len(self.triangles) + edge_to_tri_edge = {} + for t in self.triangles: + t.neighbors = [-1, -1, -1] + for i in range(nt): + t = self.triangles[i] + if t.v[0] < 0: + continue + for j in range(3): + a, b = t.v[j], t.v[(j + 1) % 3] + mn, mx = min(a, b), max(a, b) + key = (mn << 32) | mx + if key in edge_to_tri_edge: + packed = edge_to_tri_edge[key] + ot, oe = packed >> 2, packed & 3 + self.triangles[i].neighbors[j] = ot + self.triangles[ot].neighbors[oe] = i + del edge_to_tri_edge[key] + else: + edge_to_tri_edge[key] = (i << 2) | j + + def locate_triangle(self, px, py, hint=0): + tris = self.triangles + if not tris: + return -1 + n = len(tris) + if hint < 0 or hint >= n: + hint = 0 + cur = hint + verts = self.vertices + for _ in range(n * 2): + t = tris[cur] + if t.v[0] < 0: + cur = (cur + 1) % n + continue + v0, v1, v2 = verts[t.v[0]], verts[t.v[1]], verts[t.v[2]] + o0 = orient2d_xy(v0.x, v0.y, v1.x, v1.y, px, py) + o1 = orient2d_xy(v1.x, v1.y, v2.x, v2.y, px, py) + o2 = orient2d_xy(v2.x, v2.y, v0.x, v0.y, px, py) + if o0 >= -EPS and o1 >= -EPS and o2 >= -EPS: + return cur + worst, wv = -1, 0 + if o0 < wv: + wv, worst = o0, 0 + if o1 < wv: + wv, worst = o1, 1 + if o2 < wv: + wv, worst = o2, 2 + if worst >= 0 and t.neighbors[worst] >= 0: + cur = t.neighbors[worst] + else: + break + # Fallback: linear scan + for i in range(n): + t = tris[i] + if t.v[0] < 0: + continue + v0, v1, v2 = verts[t.v[0]], verts[t.v[1]], verts[t.v[2]] + if (orient2d_xy(v0.x, v0.y, v1.x, v1.y, px, py) >= -EPS and + orient2d_xy(v1.x, v1.y, v2.x, v2.y, px, py) >= -EPS and + orient2d_xy(v2.x, v2.y, v0.x, v0.y, px, py) >= -EPS): + return i + return -1 + + def has_edge(self, u, v, v2t): + for ti in v2t[u]: + t = self.triangles[ti] + if t.v[0] < 0: + continue + for j in range(3): + if (t.v[j] == u and t.v[(j + 1) % 3] == v) or \ + (t.v[j] == v and t.v[(j + 1) % 3] == u): + return True + return False + + +# ============================================================ +# PSLG cleanup (-C flag) +# ============================================================ + +def clean_pslg(mesh, tol, quiet): + pts = mesh.vertices + segs = mesh.segments + nv = len(pts) + + # Compute auto tolerance if not specified. tol == -1 is the -C default + # (bb_diag * CLOSE_ENOUGH); tol <= -2 is the always-on hygiene pass run for + # every PSLG, with an essentially-exact tolerance (bb_diag * 1e-12, ~4 decimal + # digits above double ULP) that merges only true duplicates and geometry no + # refinement could ever resolve, while still splitting crossing segments (the + # intersection test is tolerance-independent). + if tol < 0: + min_x = min(p.x for p in pts) + max_x = max(p.x for p in pts) + min_y = min(p.y for p in pts) + max_y = max(p.y for p in pts) + bb_diag = math.sqrt((max_x - min_x)**2 + (max_y - min_y)**2) + tol = bb_diag * (1e-12 if tol <= -2.0 else CLOSE_ENOUGH) + tol2 = tol * tol + + merged_nodes = 0 + split_segs = 0 + removed_segs = 0 + + # 1. Merge near-duplicate nodes (x-sorted sweep to avoid O(n²)) + remap = list(range(nv)) + xorder = sorted(range(nv), key=lambda k: pts[k].x) + for ii in range(nv): + i = xorder[ii] + if remap[i] != i: + continue + for jj in range(ii + 1, nv): + j = xorder[jj] + if pts[j].x - pts[i].x > tol: + break # x-sorted: no more candidates + if remap[j] != j: + continue + dx = pts[j].x - pts[i].x + dy = pts[j].y - pts[i].y + if dx * dx + dy * dy < tol2: + remap[j] = i + merged_nodes += 1 + + # Apply remap to segments + for s in segs: + s.v0 = remap[s.v0] + s.v1 = remap[s.v1] + + # 2. Remove degenerate segments + good = [s for s in segs if s.v0 != s.v1] + removed_segs += len(segs) - len(good) + mesh.segments = good + segs = mesh.segments + + # 3. Remove duplicate segments + seen = set() + good = [] + for s in segs: + k = edge_key(s.v0, s.v1) + if k not in seen: + seen.add(k) + good.append(s) + else: + removed_segs += 1 + mesh.segments = good + segs = mesh.segments + + # Spatial grid shared by sections 4 and 5, replacing their O(S*N) / O(S^2) + # brute-force scans. A segment's candidate points/segments are gathered from + # the grid cells its supporting line passes through (sampled at <=1-cell + # steps, each expanded to a 3x3 halo so anything within tol of the line is + # covered) and processed lowest-index-first with the same tests as before, + # so the output is identical -- just without the quadratic scan. Grid side + # is capped at 1024, so cells >> tol always. + def grid_dims(): + xmn = min(p.x for p in pts); xmx = max(p.x for p in pts) + ymn = min(p.y for p in pts); ymx = max(p.y for p in pts) + side = min(1024, max(1, int(math.sqrt(len(pts))))) + cell = max((xmx - xmn) / side, (ymx - ymn) / side, 1e-30) + return xmn, ymn, cell, int((xmx - xmn) / cell) + 1, int((ymx - ymn) / cell) + 1 + + def seg_cells(x0, y0, x1, y1, gx0, gy0, cell, gnx, gny): + # cells the segment's line passes through, each haloed 3x3 (deduped) + cells = set() + ns = max(1, int(math.hypot(x1 - x0, y1 - y0) / cell) + 1) + for s in range(ns + 1): + tt = s / ns + cx = min(max(int((x0 + (x1 - x0) * tt - gx0) / cell), 0), gnx - 1) + cy = min(max(int((y0 + (y1 - y0) * tt - gy0) / cell), 0), gny - 1) + for ddx in (-1, 0, 1): + for ddy in (-1, 0, 1): + mx, my = cx + ddx, cy + ddy + if 0 <= mx < gnx and 0 <= my < gny: + cells.add(mx * gny + my) + return cells + + # 4. Split segments at near-coincident nodes (node lies on segment). + # Points are not added in this section, so the point grid is built once. + gx0, gy0, cell, gnx, gny = grid_dims() + pgrid = {} + for i in range(len(pts)): + if i < nv and remap[i] != i: + continue # merged away + cx = min(max(int((pts[i].x - gx0) / cell), 0), gnx - 1) + cy = min(max(int((pts[i].y - gy0) / cell), 0), gny - 1) + pgrid.setdefault(cx * gny + cy, []).append(i) + changed = True + while changed: + changed = False + for si in range(len(segs)): + v0, v1 = segs[si].v0, segs[si].v1 + ax, ay = pts[v0].x, pts[v0].y + bx, by = pts[v1].x, pts[v1].y + dx, dy = bx - ax, by - ay + len2 = dx * dx + dy * dy + if len2 < tol2: + continue + sqrt_len = math.sqrt(len2) + # segment bounding box (expanded by tol) + sx0 = min(ax, bx) - tol + sx1 = max(ax, bx) + tol + sy0 = min(ay, by) - tol + sy1 = max(ay, by) + tol + # lowest-index point lying on this segment, from the cells it crosses + best_i = -1 + for c in seg_cells(ax, ay, bx, by, gx0, gy0, cell, gnx, gny): + for i in pgrid.get(c, ()): + if i == v0 or i == v1: + continue + if best_i >= 0 and i >= best_i: + continue + if i < nv and remap[i] != i: + continue + if pts[i].x < sx0 or pts[i].x > sx1 or pts[i].y < sy0 or pts[i].y > sy1: + continue + px, py = pts[i].x - ax, pts[i].y - ay + t = (dx * px + dy * py) / len2 + if t <= tol / sqrt_len or t >= 1.0 - tol / sqrt_len: + continue + cross = dx * py - dy * px + if cross * cross > tol2 * len2: + continue + best_i = i + if best_i >= 0: + # Point best_i is on segment si — split it + s2 = Segment(best_i, segs[si].v1, segs[si].marker, segs[si].lfs, segs[si].pbc_type, segs[si].no_split, segs[si].pbc_side) + segs[si] = Segment(segs[si].v0, best_i, segs[si].marker, segs[si].lfs, segs[si].pbc_type, segs[si].no_split, segs[si].pbc_side) + segs.append(s2) + split_segs += 1 + changed = True + + # 5. Split segments at segment-segment intersections. Segments and points + # change on every split, so the segment grid is rebuilt each pass; two + # intersecting segments share the crossing cell, so each is a candidate for + # the other. + changed = True + while changed: + changed = False + gx0, gy0, cell, gnx, gny = grid_dims() + sgrid = {} + for s in range(len(segs)): + for c in seg_cells(pts[segs[s].v0].x, pts[segs[s].v0].y, + pts[segs[s].v1].x, pts[segs[s].v1].y, + gx0, gy0, cell, gnx, gny): + sgrid.setdefault(c, []).append(s) + for i in range(len(segs)): + if changed: + break + ax, ay = pts[segs[i].v0].x, pts[segs[i].v0].y + bx, by = pts[segs[i].v1].x, pts[segs[i].v1].y + ax0, ax1 = min(ax, bx), max(ax, bx) + ay0, ay1 = min(ay, by), max(ay, by) + d1x, d1y = bx - ax, by - ay + # lowest j>i whose segment intersects i, from the cells i crosses + best_j = -1 + ix_best = iy_best = 0.0 + for c in seg_cells(ax, ay, bx, by, gx0, gy0, cell, gnx, gny): + for j in sgrid.get(c, ()): + if j <= i: + continue + if best_j >= 0 and j >= best_j: + continue + # Skip if they share an endpoint + if (segs[i].v0 in (segs[j].v0, segs[j].v1) or + segs[i].v1 in (segs[j].v0, segs[j].v1)): + continue + # Bounding box pre-filter + bx0 = min(pts[segs[j].v0].x, pts[segs[j].v1].x) + bx1 = max(pts[segs[j].v0].x, pts[segs[j].v1].x) + if ax1 < bx0 or ax0 > bx1: + continue + by0 = min(pts[segs[j].v0].y, pts[segs[j].v1].y) + by1 = max(pts[segs[j].v0].y, pts[segs[j].v1].y) + if ay1 < by0 or ay0 > by1: + continue + cx, cy = pts[segs[j].v0].x, pts[segs[j].v0].y + d2x, d2y = pts[segs[j].v1].x, pts[segs[j].v1].y + ddx, ddy = d2x - cx, d2y - cy + denom = d1x * ddy - d1y * ddx + # Parallel-rejection must be SCALE-FREE: denom scales with the + # product of the segment lengths, so a fixed 1e-15 cutoff + # silently skipped genuine crossings between short near-parallel + # segments (and was needlessly far from underflow for long ones). + dmag = abs(d1x * ddy) + abs(d1y * ddx) + if dmag == 0.0 or abs(denom) < 1e-12 * dmag: + continue + t = ((cx - ax) * ddy - (cy - ay) * ddx) / denom + u = ((cx - ax) * d1y - (cy - ay) * d1x) / denom + if t <= 0 or t >= 1 or u <= 0 or u >= 1: + continue + best_j = j + ix_best, iy_best = ax + t * d1x, ay + t * d1y + if best_j >= 0: + j = best_j + ix, iy = ix_best, iy_best + # Check if intersection is near an existing node + near_node = -1 + for k in range(len(pts)): + ex, ey = pts[k].x - ix, pts[k].y - iy + if ex * ex + ey * ey < tol2: + near_node = k + break + if near_node < 0: + near_node = len(pts) + pts.append(Point(ix, iy, near_node, 0)) + # Split both segments + s_i2 = Segment(near_node, segs[i].v1, segs[i].marker, segs[i].lfs, segs[i].pbc_type, segs[i].no_split, segs[i].pbc_side) + segs[i] = Segment(segs[i].v0, near_node, segs[i].marker, segs[i].lfs, segs[i].pbc_type, segs[i].no_split, segs[i].pbc_side) + segs.append(s_i2) + s_j2 = Segment(near_node, segs[j].v1, segs[j].marker, segs[j].lfs, segs[j].pbc_type, segs[j].no_split, segs[j].pbc_side) + segs[j] = Segment(segs[j].v0, near_node, segs[j].marker, segs[j].lfs, segs[j].pbc_type, segs[j].no_split, segs[j].pbc_side) + segs.append(s_j2) + split_segs += 2 + changed = True + + # Compact: remove merged-away vertices and re-index + if merged_nodes > 0: + new_idx = [-1] * len(pts) + new_pts = [] + for i in range(len(pts)): + target = remap[i] if i < nv else i + if target != i: + continue + new_idx[i] = len(new_pts) + pts[i].id = len(new_pts) + new_pts.append(pts[i]) + for i in range(nv): + if remap[i] != i: + new_idx[i] = new_idx[remap[i]] + for i in range(nv, len(pts)): + if new_idx[i] < 0: + new_idx[i] = len(new_pts) + pts[i].id = len(new_pts) + new_pts.append(pts[i]) + for s in segs: + s.v0 = new_idx[s.v0] + s.v1 = new_idx[s.v1] + mesh.vertices = new_pts + + if not quiet and (merged_nodes > 0 or split_segs > 0 or removed_segs > 0): + print(f" PSLG cleanup: merged {merged_nodes} nodes, split " + f"{split_segs} segments, removed {removed_segs} duplicates", + file=sys.stderr) + + +# ============================================================ +# Bowyer-Watson Delaunay triangulation +# ============================================================ + +def build_delaunay(mesh): + pts = mesh.vertices + n = len(pts) + + min_x = min(p.x for p in pts) + max_x = max(p.x for p in pts) + min_y = min(p.y for p in pts) + max_y = max(p.y for p in pts) + dx, dy = max_x - min_x, max_y - min_y + d = max(dx, dy) * 10.0 + + # Scale EPS to the bounding-box diagonal so tolerances stay meaningful + # regardless of coordinate magnitude. + global EPS + bbox_diag = math.sqrt(dx*dx + dy*dy) + EPS = bbox_diag * EPS_SCALE + pts.append(Point(min_x - d, min_y - 3 * d, n, 0)) + pts.append(Point(min_x + 3 * d, min_y - d, n + 1, 0)) + pts.append(Point(min_x - d, min_y + 3 * d, n + 2, 0)) + + tris = mesh.triangles + tris.clear() + tris.append(Triangle(v=[n, n + 1, n + 2], neighbors=[-1, -1, -1])) + + order = sorted(range(n), key=lambda i: (pts[i].x, pts[i].y)) + + last_tri = 0 + visit_epoch = [] + epoch = 0 + skipped_pts = 0 + + for ii in range(n): + pidx = order[ii] + p = pts[pidx] + + # Locate containing triangle by walking + start_tri = -1 + cur = min(last_tri, len(tris) - 1) + for _ in range(len(tris) * 2): + t = tris[cur] + if t.v[0] < 0: + cur = 0 + continue + o0 = orient2d(pts[t.v[0]], pts[t.v[1]], p) + o1 = orient2d(pts[t.v[1]], pts[t.v[2]], p) + o2 = orient2d(pts[t.v[2]], pts[t.v[0]], p) + if o0 >= -EPS and o1 >= -EPS and o2 >= -EPS: + start_tri = cur + break + w, wv = -1, 0 + if o0 < wv: + wv, w = o0, 0 + if o1 < wv: + wv, w = o1, 1 + if o2 < wv: + wv, w = o2, 2 + if w >= 0 and t.neighbors[w] >= 0: + cur = t.neighbors[w] + else: + break + if start_tri < 0: + for i in range(len(tris)): + t = tris[i] + if t.v[0] < 0: + continue + if (orient2d(pts[t.v[0]], pts[t.v[1]], p) >= -EPS and + orient2d(pts[t.v[1]], pts[t.v[2]], p) >= -EPS and + orient2d(pts[t.v[2]], pts[t.v[0]], p) >= -EPS): + start_tri = i + break + if start_tri < 0: + skipped_pts += 1 + continue + + # Precise re-walk with exact orientation tests (the Bowyer-Watson + # analogue of insert_point_lawson's exact refinement): the EPS-tolerant + # walk above can stop at a triangle the point is actually just OUTSIDE + # of, and near a vertex the circumcircle of that triangle can exclude the + # point -- yielding a falsely empty cavity below. + for _ in range(len(tris)): + t = tris[start_tri] + neg_edge = -1 + for j in range(3): + if orient2d(pts[t.v[j]], pts[t.v[(j + 1) % 3]], p) < 0.0: + neg_edge = j + break + if neg_edge < 0: + break # exact containment (or on edge) + nb = t.neighbors[neg_edge] + if nb < 0: + break # hull edge -- accept current + start_tri = nb + + # BFS cavity (epoch-based visited to avoid per-point allocation) + epoch += 1 + while len(visit_epoch) < len(tris): + visit_epoch.append(0) + bad = [] + stk = [start_tri] + while stk: + cur = stk.pop() + if cur < 0 or cur >= len(tris) or visit_epoch[cur] == epoch: + continue + visit_epoch[cur] = epoch + t = tris[cur] + if t.v[0] < 0: + continue + if in_circle(pts[t.v[0]], pts[t.v[1]], pts[t.v[2]], p) > 0: + bad.append(cur) + for j in range(3): + if t.neighbors[j] >= 0 and visit_epoch[t.neighbors[j]] != epoch: + stk.append(t.neighbors[j]) + + # Empty cavity: with exact predicates, a point strictly inside any + # triangle is strictly inside its circumcircle, so an empty cavity means + # the point coincides with an existing vertex (or the EPS walk mislocated + # it). Leave it unconnected, like Triangle's duplicate-vertex handling. + # (A previous revision force-fanned the located triangle here, which + # emitted zero-area triangles for exact duplicates and overlapping + # triangles for mislocations -- a corrupted triangulation either way.) + if not bad: + skipped_pts += 1 + continue + + # Collect boundary polygon + bad_set = set(bad) + poly = [] # list of (v0, v1, outer_tri, outer_local_edge) + for bi in bad: + t = tris[bi] + for j in range(3): + nb = t.neighbors[j] + if nb < 0 or nb not in bad_set: + oe = -1 + if nb >= 0: + for k in range(3): + if tris[nb].neighbors[k] == bi: + oe = k + break + poly.append((t.v[j], t.v[(j + 1) % 3], nb, oe)) + + # (A strict-CCW fan-refusal guard used to live here -- one orient2d per + # cavity-boundary edge, refusing the insertion if any fan triangle came + # out non-CCW. With exact predicates the precise re-walk above guarantees + # a correctly located point produces an all-CCW fan, so it refused 0 + # insertions suite-wide and was removed; outputs byte-identical without + # it. The empty-cavity skip above still catches exact duplicates. + # tangle.cpp 0.4.3.) + + n_new = len(poly) + bad.sort() + slots = [] + for i in range(min(n_new, len(bad))): + slots.append(bad[i]) + while len(slots) < n_new: + slots.append(len(tris)) + tris.append(Triangle()) + for i in range(n_new, len(bad)): + tris[bad[i]].v = [-1, -1, -1] + tris[bad[i]].neighbors = [-1, -1, -1] + + for i in range(n_new): + s = slots[i] + tris[s].v = [poly[i][0], poly[i][1], pidx] + tris[s].neighbors = [-1, -1, -1] + tris[s].region_attrib = 0 + tris[s].region_max_area = -1 + tris[s].marker = 0 + # Strictly CCW by the fan guard above -- no orientation fix-up. + + # Wire outer adjacency + for i in range(n_new): + s = slots[i] + if poly[i][2] < 0: + continue + for j in range(3): + a, b = tris[s].v[j], tris[s].v[(j + 1) % 3] + if (a == poly[i][0] and b == poly[i][1]) or \ + (a == poly[i][1] and b == poly[i][0]): + tris[s].neighbors[j] = poly[i][2] + if poly[i][3] >= 0: + tris[poly[i][2]].neighbors[poly[i][3]] = s + break + + # Wire internal adjacency + pidx_edges = {} + for i in range(n_new): + s = slots[i] + for j in range(3): + a, b = tris[s].v[j], tris[s].v[(j + 1) % 3] + if a != pidx and b != pidx: + continue + other = b if a == pidx else a + if other in pidx_edges: + os, oe = pidx_edges[other] + tris[s].neighbors[j] = os + tris[os].neighbors[oe] = s + del pidx_edges[other] + else: + pidx_edges[other] = (s, j) + last_tri = slots[0] + + if skipped_pts > 0: + sys.stderr.write( + f"Warning: {skipped_pts} point{'' if skipped_pts == 1 else 's'} " + "could not be inserted (duplicate or unlocatable); left unconnected.\n") + + # Remove dead triangles and super-triangle vertices + keep = [False] * len(tris) + for i in range(len(tris)): + if tris[i].v[0] < 0: + continue + if tris[i].v[0] >= n or tris[i].v[1] >= n or tris[i].v[2] >= n: + continue + keep[i] = True + remap = [-1] * len(tris) + new_tris = [] + for i in range(len(tris)): + if keep[i]: + remap[i] = len(new_tris) + new_tris.append(tris[i]) + for t in new_tris: + for j in range(3): + t.neighbors[j] = remap[t.neighbors[j]] if t.neighbors[j] >= 0 else -1 + mesh.triangles = new_tris + del pts[n:] + + # Detect and re-insert orphaned vertices + mesh.rebuild_adjacency() + tris = mesh.triangles + orphan_visit_epoch = [] + orphan_epoch = 0 + for retry in range(10): + in_mesh = [False] * n + for t in tris: + if t.v[0] >= 0: + for j in range(3): + if t.v[j] < n: + in_mesh[t.v[j]] = True + orphans = [i for i in range(n) if not in_mesh[i]] + if not orphans: + break + + for pidx in orphans: + p = pts[pidx] + contain_tri = mesh.locate_triangle(p.x, p.y) + if contain_tri < 0: + continue + + # Skip if nearly coincident + too_close = False + for j in range(3): + dx2 = pts[tris[contain_tri].v[j]].x - p.x + dy2 = pts[tris[contain_tri].v[j]].y - p.y + if dx2 * dx2 + dy2 * dy2 < EPS * EPS: + too_close = True + break + if too_close: + continue + + # Bowyer-Watson insertion with strict inCircle to avoid over-expanding + # the cavity and orphaning nearby vertices + orphan_epoch += 1 + while len(orphan_visit_epoch) < len(tris): + orphan_visit_epoch.append(0) + bad = [] + stk = [contain_tri] + while stk: + cur = stk.pop() + if cur < 0 or cur >= len(tris) or orphan_visit_epoch[cur] == orphan_epoch: + continue + orphan_visit_epoch[cur] = orphan_epoch + t = tris[cur] + if t.v[0] < 0: + continue + if in_circle(pts[t.v[0]], pts[t.v[1]], pts[t.v[2]], p) > 0: + bad.append(cur) + for j in range(3): + if t.neighbors[j] >= 0 and orphan_visit_epoch[t.neighbors[j]] != orphan_epoch: + stk.append(t.neighbors[j]) + if not bad: + bad.append(contain_tri) + + bad_set = set(bad) + poly = [] + for bi in bad: + t = tris[bi] + for j in range(3): + nb = t.neighbors[j] + if nb < 0 or nb not in bad_set: + oe = -1 + if nb >= 0: + for k in range(3): + if tris[nb].neighbors[k] == bi: + oe = k + break + poly.append((t.v[j], t.v[(j + 1) % 3], nb, oe)) + + n_new = len(poly) + bad.sort() + slots = [] + for i in range(min(n_new, len(bad))): + slots.append(bad[i]) + while len(slots) < n_new: + slots.append(len(tris)) + tris.append(Triangle()) + for i in range(n_new, len(bad)): + tris[bad[i]].v = [-1, -1, -1] + tris[bad[i]].neighbors = [-1, -1, -1] + + for i in range(n_new): + s = slots[i] + tris[s].v = [poly[i][0], poly[i][1], pidx] + tris[s].neighbors = [-1, -1, -1] + tris[s].region_attrib = 0 + tris[s].region_max_area = -1 + tris[s].marker = 0 + if orient2d(pts[tris[s].v[0]], pts[tris[s].v[1]], pts[tris[s].v[2]]) < 0: + tris[s].v[0], tris[s].v[1] = tris[s].v[1], tris[s].v[0] + + # Wire outer adjacency + for i in range(n_new): + s = slots[i] + if poly[i][2] < 0: + continue + for j in range(3): + a, b = tris[s].v[j], tris[s].v[(j + 1) % 3] + if (a == poly[i][0] and b == poly[i][1]) or \ + (a == poly[i][1] and b == poly[i][0]): + tris[s].neighbors[j] = poly[i][2] + if poly[i][3] >= 0: + tris[poly[i][2]].neighbors[poly[i][3]] = s + break + + # Wire internal adjacency + pidx_edges = {} + for i in range(n_new): + s = slots[i] + for j in range(3): + a, b = tris[s].v[j], tris[s].v[(j + 1) % 3] + if a != pidx and b != pidx: + continue + other = b if a == pidx else a + if other in pidx_edges: + os, oe = pidx_edges[other] + tris[s].neighbors[j] = os + tris[os].neighbors[oe] = s + del pidx_edges[other] + else: + pidx_edges[other] = (s, j) + + +# ============================================================ +# CDT: enforce PSLG segments via cavity retriangulation +# ============================================================ + +def enforce_constraints(mesh, quiet): + if not mesh.segments: + return 0 + + def build_v2t(): + v2t = [[] for _ in range(len(mesh.vertices))] + for i, t in enumerate(mesh.triangles): + if t.v[0] < 0: + continue + for j in range(3): + v2t[t.v[j]].append(i) + return v2t + + def ear_clip(poly_verts): + # Strict ear clipping (mirrors C++ earClip): clip only well-formed empty + # ears, and reject an ear when a vertex lies ON its new diagonal (the + # >= 0 term) -- clipping it would orphan that vertex as a hanging node + # (FEMM's near-collinear corner twins land exactly there). No relaxed + # pass / fan fallback. On a stall it returns an INCOMPLETE result; the + # caller checks the triangle count (a simple m-gon -> m-2 triangles) and + # declines, routing the segment to the midpoint-split + rebuild fallback. + result = [] + p = list(poly_verts) + while len(p) > 3: + clipped = False + nn = len(p) + for i in range(nn): + prev_v = p[(i + nn - 1) % nn] + cur_v = p[i] + next_v = p[(i + 1) % nn] + if orient2d(mesh.vertices[prev_v], mesh.vertices[cur_v], + mesh.vertices[next_v]) <= 0.0: + continue + blocked = False + for k in range(nn): + if k in ((i + nn - 1) % nn, i, (i + 1) % nn): + continue + vi = p[k] + if (orient2d(mesh.vertices[prev_v], mesh.vertices[cur_v], + mesh.vertices[vi]) > 0 and + orient2d(mesh.vertices[cur_v], mesh.vertices[next_v], + mesh.vertices[vi]) > 0 and + orient2d(mesh.vertices[next_v], mesh.vertices[prev_v], + mesh.vertices[vi]) >= 0): + blocked = True + break + if blocked: + continue + result.append([prev_v, cur_v, next_v]) + p.pop(i) + clipped = True + break + if not clipped: + break + if len(p) == 3: + o = orient2d(mesh.vertices[p[0]], mesh.vertices[p[1]], mesh.vertices[p[2]]) + if o > 0: + result.append([p[0], p[1], p[2]]) + elif o < 0: + result.append([p[0], p[2], p[1]]) + return result + + v2t = build_v2t() + + # Split segments at collinear intermediate vertices + # Sort vertex indices by x-coordinate for efficient range queries + xorder = sorted(range(len(mesh.vertices)), key=lambda i: mesh.vertices[i].x) + xsorted = [mesh.vertices[i].x for i in xorder] + + new_segs = [] + for seg in mesh.segments: + su, sv = seg.v0, seg.v1 + sx, sy = mesh.vertices[su].x, mesh.vertices[su].y + ddx, ddy = mesh.vertices[sv].x - sx, mesh.vertices[sv].y - sy + len2 = ddx * ddx + ddy * ddy + if len2 < EPS * EPS: + new_segs.append(seg) + continue + eps2_len2 = EPS * EPS * len2 # squared collinearity threshold + bxlo = min(sx, sx + ddx) - EPS + bxhi = max(sx, sx + ddx) + EPS + bylo = min(sy, sy + ddy) - EPS + byhi = max(sy, sy + ddy) + EPS + # Binary search for x-range [bxlo, bxhi] + ilo = bisect.bisect_left(xsorted, bxlo) + ihi = bisect.bisect_right(xsorted, bxhi) + on_seg = [] + for ii in range(ilo, ihi): + i = xorder[ii] + if i == su or i == sv: + continue + vy = mesh.vertices[i].y + if vy < bylo or vy > byhi: + continue + px, py = mesh.vertices[i].x - sx, vy - sy + cross = ddx * py - ddy * px + if cross * cross > eps2_len2: + continue + t_param = (ddx * px + ddy * py) / len2 + if t_param <= EPS_SCALE or t_param >= 1.0 - EPS_SCALE: + continue + on_seg.append((t_param, i)) + if not on_seg: + new_segs.append(seg) + else: + on_seg.sort() + prev_v = su + for _, vi in on_seg: + new_segs.append(Segment(prev_v, vi, seg.marker, seg.lfs, seg.pbc_type, seg.no_split, seg.pbc_side)) + prev_v = vi + new_segs.append(Segment(prev_v, sv, seg.marker, seg.lfs, seg.pbc_type, seg.no_split, seg.pbc_side)) + if len(new_segs) != len(mesh.segments): + if not quiet: + print(f" Split {len(mesh.segments)} segments into {len(new_segs)} sub-segments (collinear vertices)", + file=sys.stderr) + mesh.segments = new_segs + + # Direct adjacency check: flip-based + mesh.rebuild_adjacency() + v2t = build_v2t() + # Constrained-segment edges: a flip must never destroy one of these to + # enforce another, or two crossing segments oscillate A<->B forever. + seg_set = set() + for s in mesh.segments: + seg_set.add(edge_key(s.v0, s.v1)) + any_fixed = True + # Hard cap: each segment needs only a handful of flips; the cap stops a + # cocircular/collinear configuration from spinning. Segments not enforced + # here fall through to the corridor + midpoint-split stages below. + flip_passes = 2 * len(mesh.segments) + 16 + while any_fixed and flip_passes > 0: + flip_passes -= 1 + any_fixed = False + for seg in mesh.segments: + su2, sv2 = seg.v0, seg.v1 + if mesh.has_edge(su2, sv2, v2t): + continue + found = False + for ti in v2t[su2]: + if mesh.has_edge(su2, sv2, v2t): + break + t = mesh.triangles[ti] + if t.v[0] < 0: + continue + for j in range(3): + nb = t.neighbors[j] + if nb < 0: + continue + tn = mesh.triangles[nb] + has_sv = any(tn.v[k] == sv2 for k in range(3)) + if not has_sv: + continue + a, b = t.v[j], t.v[(j + 1) % 3] + if a == su2 or b == su2 or a == sv2 or b == sv2: + continue + if edge_key(a, b) in seg_set: + continue # never flip out a constrained segment + p_v, q_v = su2, sv2 + o1 = orient2d(mesh.vertices[p_v], mesh.vertices[a], mesh.vertices[q_v]) + o2 = orient2d(mesh.vertices[p_v], mesh.vertices[q_v], mesh.vertices[b]) + if o1 < 0 or o2 < 0: + continue + # Do the flip + n_pa, n_bp = -1, -1 + for k in range(3): + e0, e1 = t.v[k], t.v[(k + 1) % 3] + if (e0 == p_v and e1 == a) or (e0 == a and e1 == p_v): + n_pa = t.neighbors[k] + if (e0 == b and e1 == p_v) or (e0 == p_v and e1 == b): + n_bp = t.neighbors[k] + n_aq, n_qb = -1, -1 + for k in range(3): + e0, e1 = tn.v[k], tn.v[(k + 1) % 3] + if (e0 == a and e1 == q_v) or (e0 == q_v and e1 == a): + n_aq = tn.neighbors[k] + if (e0 == q_v and e1 == b) or (e0 == b and e1 == q_v): + n_qb = tn.neighbors[k] + t.v = [p_v, a, q_v] + t.neighbors = [n_pa, n_aq, nb] + tn.v = [p_v, q_v, b] + tn.neighbors = [ti, n_qb, n_bp] + if n_aq >= 0: + for k in range(3): + if mesh.triangles[n_aq].neighbors[k] == nb: + mesh.triangles[n_aq].neighbors[k] = ti + break + if n_bp >= 0: + for k in range(3): + if mesh.triangles[n_bp].neighbors[k] == ti: + mesh.triangles[n_bp].neighbors[k] = nb + break + v2t = build_v2t() + any_fixed = True + found = True + break + if found: + break + + # Cavity fallback for the few segments the flip block couldn't enforce + def segs_cross(v1, v2, v3, v4): + d1 = orient2d(mesh.vertices[v1], mesh.vertices[v2], mesh.vertices[v3]) + d2 = orient2d(mesh.vertices[v1], mesh.vertices[v2], mesh.vertices[v4]) + if d1 * d2 >= 0: + return False + d3 = orient2d(mesh.vertices[v3], mesh.vertices[v4], mesh.vertices[v1]) + d4 = orient2d(mesh.vertices[v3], mesh.vertices[v4], mesh.vertices[v2]) + if d3 * d4 >= 0: + return False + return True + + for seg in mesh.segments: + su, sv = seg.v0, seg.v1 + if mesh.has_edge(su, sv, v2t): + continue + + # Cavity-based fallback + crossed_set = set() + for ti in range(len(mesh.triangles)): + t = mesh.triangles[ti] + if t.v[0] < 0: + continue + for j in range(3): + a, b = t.v[j], t.v[(j + 1) % 3] + if segs_cross(su, sv, a, b): + crossed_set.add(ti) + break + for vi in (su, sv): + for ti in v2t[vi]: + t = mesh.triangles[ti] + if t.v[0] < 0: + continue + for j in range(3): + if t.neighbors[j] in crossed_set: + for k in range(3): + a, b = t.v[k], t.v[(k + 1) % 3] + if segs_cross(su, sv, a, b): + crossed_set.add(ti) + break + break + + if crossed_set: + bnd_edges_cav = [] + for ti in crossed_set: + t = mesh.triangles[ti] + for j in range(3): + nb = t.neighbors[j] + if nb < 0 or nb not in crossed_set: + bnd_edges_cav.append((t.v[j], t.v[(j + 1) % 3], nb)) + + # Chain boundary edges + adj_map = {} + for i, be in enumerate(bnd_edges_cav): + adj_map.setdefault(be[0], []).append((be[1], i)) + for v_key in adj_map: + if len(adj_map[v_key]) > 1: + adj_map[v_key].sort( + key=lambda x: (math.atan2( + mesh.vertices[x[0]].y - mesh.vertices[v_key].y, + mesh.vertices[x[0]].x - mesh.vertices[v_key].x), + x[0])) # deterministic tie-break on collinear edges + + poly_chain = [] + used_edges = set() + start_v = bnd_edges_cav[0][0] + cur_v = start_v + prev_v = -1 + chain_ok = True + + for _ in range(len(bnd_edges_cav) + 1): + poly_chain.append(cur_v) + outs = adj_map.get(cur_v, []) + best_idx = -1 + if len(outs) == 1 and outs[0][1] not in used_edges: + best_idx = 0 + elif prev_v >= 0: + in_dx = mesh.vertices[cur_v].x - mesh.vertices[prev_v].x + in_dy = mesh.vertices[cur_v].y - mesh.vertices[prev_v].y + in_angle = math.atan2(in_dy, in_dx) + best_turn = 1e30 + for k in range(len(outs)): + if outs[k][1] in used_edges: + continue + ox = mesh.vertices[outs[k][0]].x - mesh.vertices[cur_v].x + oy = mesh.vertices[outs[k][0]].y - mesh.vertices[cur_v].y + out_angle = math.atan2(oy, ox) + turn = out_angle - in_angle + while turn <= 0: + turn += 2 * math.pi + # deterministic tie-break: on equal turn (collinear + # outgoing edges) prefer the lower target-vertex index + if turn < best_turn or \ + (turn == best_turn and best_idx >= 0 and + outs[k][0] < outs[best_idx][0]): + best_turn = turn + best_idx = k + else: + for k in range(len(outs)): + if outs[k][1] not in used_edges: + best_idx = k + break + if best_idx < 0: + chain_ok = False + break + used_edges.add(outs[best_idx][1]) + prev_v = cur_v + cur_v = outs[best_idx][0] + if cur_v == start_v and len(poly_chain) > 2: + break + + if len(poly_chain) >= 3 and cur_v != start_v: + chain_ok = False + + new_tris_cav = [] + if chain_ok and len(poly_chain) >= 3: + su_idx = sv_idx = -1 + for i in range(len(poly_chain)): + if poly_chain[i] == su: + su_idx = i + if poly_chain[i] == sv: + sv_idx = i + if su_idx >= 0 and sv_idx >= 0: + nn = len(poly_chain) + poly1, poly2 = [], [] + i = su_idx + while True: + poly1.append(poly_chain[i]) + if i == sv_idx: + break + i = (i + 1) % nn + i = sv_idx + while True: + poly2.append(poly_chain[i]) + if i == su_idx: + break + i = (i + 1) % nn + # Accept only COMPLETE triangulations (m-gon -> m-2 tris). A + # short count means ear_clip declined a degenerate cavity; + # leave new_tris_cav empty so the segment falls through to the + # midpoint-split + rebuild fallback. + t1, t2 = ear_clip(poly1), ear_clip(poly2) + if len(t1) == max(0, len(poly1) - 2) and len(t2) == max(0, len(poly2) - 2): + new_tris_cav = t1 + t2 + else: + new_tris_cav = [] + else: + tt = ear_clip(poly_chain) + new_tris_cav = tt if len(tt) == max(0, len(poly_chain) - 2) else [] + + if new_tris_cav: + crossed_vec = sorted(crossed_set) + total_new = len(new_tris_cav) + total_old = len(crossed_vec) + + slots = [] + for i in range(min(total_old, total_new)): + slots.append(crossed_vec[i]) + while len(slots) < total_new: + slots.append(len(mesh.triangles)) + mesh.triangles.append(Triangle()) + for i in range(total_new, total_old): + mesh.triangles[crossed_vec[i]].v = [-1, -1, -1] + mesh.triangles[crossed_vec[i]].neighbors = [-1, -1, -1] + + for i in range(total_new): + s = slots[i] + mesh.triangles[s].v = list(new_tris_cav[i]) + mesh.triangles[s].neighbors = [-1, -1, -1] + + edge_map2 = {} + for i in range(total_new): + s = slots[i] + t = mesh.triangles[s] + for j in range(3): + key = edge_key(t.v[j], t.v[(j + 1) % 3]) + if key in edge_map2: + os, oe = edge_map2[key] + t.neighbors[j] = os + mesh.triangles[os].neighbors[oe] = s + del edge_map2[key] + else: + edge_map2[key] = (s, j) + + for be in bnd_edges_cav: + key = edge_key(be[0], be[1]) + if key in edge_map2: + s, j = edge_map2[key] + mesh.triangles[s].neighbors[j] = be[2] + if be[2] >= 0: + ot = mesh.triangles[be[2]] + for k in range(3): + if ot.neighbors[k] in crossed_set: + e0, e1 = ot.v[k], ot.v[(k + 1) % 3] + if edge_key(e0, e1) == key: + ot.neighbors[k] = s + break + v2t = build_v2t() + + # Count and return remaining missing segments + miss = 0 + v2t = build_v2t() + for s in mesh.segments: + if not mesh.has_edge(s.v0, s.v1, v2t): + miss += 1 + if not quiet: + print(f" UNENFORCED: v{s.v0} ({mesh.vertices[s.v0].x},{mesh.vertices[s.v0].y})" + f"->v{s.v1} ({mesh.vertices[s.v1].x},{mesh.vertices[s.v1].y})" + f" len={math.sqrt((mesh.vertices[s.v1].x-mesh.vertices[s.v0].x)*(mesh.vertices[s.v1].x-mesh.vertices[s.v0].x)+(mesh.vertices[s.v1].y-mesh.vertices[s.v0].y)*(mesh.vertices[s.v1].y-mesh.vertices[s.v0].y))}", + file=sys.stderr) + if not quiet and miss > 0: + print(f"Warning: {miss} of {len(mesh.segments)} segments could not be enforced.", + file=sys.stderr) + return miss + + +# ============================================================ +# Flood-fill hole removal +# ============================================================ + +def remove_holes(mesh, opts): + if opts.no_holes: + return + mesh.rebuild_adjacency() + + seg_edges = set() + for s in mesh.segments: + seg_edges.add(edge_key(s.v0, s.v1)) + + remove = [False] * len(mesh.triangles) + + # Remove exterior: flood-fill from boundary triangles whose boundary + # edge is not a constrained segment. + q = deque() + for i in range(len(mesh.triangles)): + if remove[i]: + continue + t = mesh.triangles[i] + for j in range(3): + if t.neighbors[j] >= 0: + continue + key = edge_key(t.v[j], t.v[(j + 1) % 3]) + if key in seg_edges: + continue + if not remove[i]: + remove[i] = True + q.append(i) + break + while q: + ti = q.popleft() + t = mesh.triangles[ti] + for j in range(3): + nb = t.neighbors[j] + if nb < 0 or remove[nb]: + continue + key = edge_key(t.v[j], t.v[(j + 1) % 3]) + if key in seg_edges: + continue + remove[nb] = True + q.append(nb) + + # Remove explicit holes via flood fill from seed points + for hole in mesh.holes: + seed_tri = mesh.locate_triangle(hole.x, hole.y) + if seed_tri < 0: + continue + q = deque([seed_tri]) + remove[seed_tri] = True + while q: + ti = q.popleft() + t = mesh.triangles[ti] + for j in range(3): + nb = t.neighbors[j] + if nb < 0 or remove[nb]: + continue + key = edge_key(t.v[j], t.v[(j + 1) % 3]) + if key in seg_edges: + continue + remove[nb] = True + q.append(nb) + + kept = [] + remap = [-1] * len(mesh.triangles) + for i in range(len(mesh.triangles)): + if not remove[i]: + remap[i] = len(kept) + kept.append(mesh.triangles[i]) + for t in kept: + for j in range(3): + t.neighbors[j] = remap[t.neighbors[j]] if t.neighbors[j] >= 0 else -1 + mesh.triangles = kept + + +# ============================================================ +# Flood-fill region assignment +# ============================================================ + +def assign_regions(mesh, opts): + if not opts.regions or not mesh.regions: + return + mesh.rebuild_adjacency() + + seg_edges = set() + for s in mesh.segments: + seg_edges.add(edge_key(s.v0, s.v1)) + + for t in mesh.triangles: + t.region_attrib = 0 + t.region_max_area = -1 + + for reg in mesh.regions: + seed_tri = mesh.locate_triangle(reg.x, reg.y) + if seed_tri < 0: + continue + q = deque([seed_tri]) + visited = [False] * len(mesh.triangles) + visited[seed_tri] = True + while q: + ti = q.popleft() + t = mesh.triangles[ti] + t.region_attrib = reg.attrib + t.region_max_area = reg.max_area + for j in range(3): + nb = t.neighbors[j] + if nb < 0 or visited[nb]: + continue + key = edge_key(t.v[j], t.v[(j + 1) % 3]) + if key in seg_edges: + continue + visited[nb] = True + q.append(nb) + + +# ============================================================ +# Point insertion for quality refinement (Lawson flip-based) +# ============================================================ + +def insert_point_lawson(mesh, np, hint_tri, constrained_edges, affected, min_dist2=0.0): + pidx = len(mesh.vertices) + if len(mesh.vertices) >= mesh.max_vertices: + raise ValueError("Cutout mesh vertex limit exceeded") + mesh.vertices.append(np) + + ti = mesh.locate_triangle(np.x, np.y, hint_tri) + if ti < 0: + mesh.vertices.pop() + return -1 + + # Precise refinement: use exact orient2d_xy to walk from the approximate + # triangle to the true containing triangle, à la Triangle's preciselocate(). + on_edge = -1 + for _ in range(len(mesh.triangles)): + tr = mesh.triangles[ti] + if tr.v[0] < 0: + mesh.vertices.pop() + return -1 + on_edge = -1 + neg_edge = -1 + for j in range(3): + o = orient2d_xy(mesh.vertices[tr.v[j]].x, mesh.vertices[tr.v[j]].y, + mesh.vertices[tr.v[(j+1)%3]].x, mesh.vertices[tr.v[(j+1)%3]].y, + np.x, np.y) + if o == 0.0: + on_edge = j + elif o < 0.0: + neg_edge = j + if neg_edge < 0: + break # all orient2d >= 0: point is inside (or on edge) + nb = tr.neighbors[neg_edge] + if nb < 0: + break # boundary — accept current triangle + ti = nb + on_edge = -1 + + t = mesh.triangles[ti] + + # If precise refinement didn't detect on-edge (midpoint rounding), + # fall back to EPS-based on-edge detection for segment midpoints. + if on_edge < 0: + for j in range(3): + va, vb = t.v[j], t.v[(j + 1) % 3] + edx = mesh.vertices[vb].x - mesh.vertices[va].x + edy = mesh.vertices[vb].y - mesh.vertices[va].y + o = orient2d_xy(mesh.vertices[va].x, mesh.vertices[va].y, + mesh.vertices[vb].x, mesh.vertices[vb].y, np.x, np.y) + if o * o < EPS * EPS * (edx * edx + edy * edy): + on_edge = j + break + + # Check near-coincident with existing vertex. The reject radius is normally + # EPS, but the quality-refinement caller passes min_dist2 = (0.5*shortest + # edge)^2 as a too-close guard: an off-center/circumcenter landing within half + # the bad triangle's shortest edge of an existing vertex is rejected (the + # triangle is then tolerated). Stops the insertion radius collapsing on + # near-degenerate geometry. No-op on well-conditioned meshes. + rej_r2 = max(EPS * EPS, min_dist2) + for j in range(3): + dx = mesh.vertices[t.v[j]].x - np.x + dy = mesh.vertices[t.v[j]].y - np.y + if dx * dx + dy * dy < rej_r2: + mesh.vertices.pop() + return -1 + + # Output-orientation validation (parity with tangle.cpp 156f54e): every + # triangle a split creates must be strictly CCW (positive area). The locate + # walk can accept a boundary triangle for a point actually OUTSIDE the domain; + # splitting there manufactures an inverted triangle, and one invalid triangle + # corrupts every circumcenter / flip / adjacency built on it. Refuse any such + # split (the bad triangle is tolerated). Exact predicate => no inverted or + # degenerate triangle can ever enter the mesh. + def tri_ccw(a, b, c): + return orient2d_xy(mesh.vertices[a].x, mesh.vertices[a].y, + mesh.vertices[b].x, mesh.vertices[b].y, + mesh.vertices[c].x, mesh.vertices[c].y) > 0.0 + + # Inherit region from the containing triangle + reg_attr = t.region_attrib + reg_area = t.region_max_area + + flip_stack = [] + + if on_edge < 0: + # Interior: split triangle into 3 + v0, v1, v2 = t.v[0], t.v[1], t.v[2] + if not (tri_ccw(v0, v1, pidx) and tri_ccw(v1, v2, pidx) and tri_ccw(v2, v0, pidx)): + mesh.vertices.pop() + return -1 + n0, n1, n2 = t.neighbors[0], t.neighbors[1], t.neighbors[2] + + t0 = ti + t1 = len(mesh.triangles) + mesh.triangles.append(Triangle()) + t2 = len(mesh.triangles) + mesh.triangles.append(Triangle()) + + mesh.triangles[t0] = Triangle(v=[v0, v1, pidx], neighbors=[n0, t1, t2], + region_attrib=reg_attr, region_max_area=reg_area) + mesh.triangles[t1] = Triangle(v=[v1, v2, pidx], neighbors=[n1, t2, t0], + region_attrib=reg_attr, region_max_area=reg_area) + mesh.triangles[t2] = Triangle(v=[v2, v0, pidx], neighbors=[n2, t0, t1], + region_attrib=reg_attr, region_max_area=reg_area) + + if n1 >= 0: + for k in range(3): + if mesh.triangles[n1].neighbors[k] == ti: + mesh.triangles[n1].neighbors[k] = t1 + break + if n2 >= 0: + for k in range(3): + if mesh.triangles[n2].neighbors[k] == ti: + mesh.triangles[n2].neighbors[k] = t2 + break + + flip_stack.extend([(t0, 0), (t1, 0), (t2, 0)]) + affected.extend([t0, t1, t2]) + else: + # On edge: split 2 triangles into 4 + va = t.v[on_edge] + vb = t.v[(on_edge + 1) % 3] + vc = t.v[(on_edge + 2) % 3] + n_ab = t.neighbors[on_edge] + n_bvc = t.neighbors[(on_edge + 1) % 3] + n_cva = t.neighbors[(on_edge + 2) % 3] + + if n_ab < 0: + # Boundary edge: split 1 triangle into 2 + if not (tri_ccw(va, pidx, vc) and tri_ccw(pidx, vb, vc)): + mesh.vertices.pop() + return -1 + t0 = ti + t1 = len(mesh.triangles) + mesh.triangles.append(Triangle()) + + mesh.triangles[t0] = Triangle(v=[va, pidx, vc], neighbors=[-1, t1, n_cva], + region_attrib=reg_attr, region_max_area=reg_area) + mesh.triangles[t1] = Triangle(v=[pidx, vb, vc], neighbors=[-1, n_bvc, t0], + region_attrib=reg_attr, region_max_area=reg_area) + + if n_bvc >= 0: + for k in range(3): + if mesh.triangles[n_bvc].neighbors[k] == ti: + mesh.triangles[n_bvc].neighbors[k] = t1 + break + + flip_stack.extend([(t0, 2), (t1, 1)]) + affected.extend([t0, t1]) + else: + # Both sides: split 2 triangles into 4 + tn = mesh.triangles[n_ab] + je = -1 + for k in range(3): + if tn.neighbors[k] == ti: + je = k + break + if je < 0: + mesh.vertices.pop() + return -1 + + vd = tn.v[(je + 2) % 3] + # Coincidence check for the neighbor's apex: the on-edge split creates + # edges to vd, but the rejection loop above only tested the CONTAINING + # triangle's vertices. A point within the reject radius of vd passed + # that test, and the strict-CCW guard below accepts any sub-EPS-thin + # but positively-oriented sliver -- a legal-but-degenerate triangle at vd. + dxv = mesh.vertices[vd].x - np.x + dyv = mesh.vertices[vd].y - np.y + if dxv * dxv + dyv * dyv < rej_r2: + mesh.vertices.pop() + return -1 + if tn.v[je] == vb and tn.v[(je + 1) % 3] == va: + n_avd = tn.neighbors[(je + 1) % 3] + n_dvb = tn.neighbors[(je + 2) % 3] + else: + n_dvb = tn.neighbors[(je + 1) % 3] + n_avd = tn.neighbors[(je + 2) % 3] + + if not (tri_ccw(va, pidx, vc) and tri_ccw(pidx, vb, vc) and + tri_ccw(vb, pidx, vd) and tri_ccw(pidx, va, vd)): + mesh.vertices.pop() + return -1 + # Inherit region from the neighbor triangle + reg_attr2 = tn.region_attrib + reg_area2 = tn.region_max_area + + t0 = ti + t1 = len(mesh.triangles) + mesh.triangles.append(Triangle()) + t2 = n_ab + t3 = len(mesh.triangles) + mesh.triangles.append(Triangle()) + + mesh.triangles[t0] = Triangle(v=[va, pidx, vc], neighbors=[t3, t1, n_cva], + region_attrib=reg_attr, region_max_area=reg_area) + mesh.triangles[t1] = Triangle(v=[pidx, vb, vc], neighbors=[t2, n_bvc, t0], + region_attrib=reg_attr, region_max_area=reg_area) + mesh.triangles[t2] = Triangle(v=[vb, pidx, vd], neighbors=[t1, t3, n_dvb], + region_attrib=reg_attr2, region_max_area=reg_area2) + mesh.triangles[t3] = Triangle(v=[pidx, va, vd], neighbors=[t0, n_avd, t2], + region_attrib=reg_attr2, region_max_area=reg_area2) + + if n_bvc >= 0: + for k in range(3): + if mesh.triangles[n_bvc].neighbors[k] == ti: + mesh.triangles[n_bvc].neighbors[k] = t1 + break + if n_avd >= 0: + for k in range(3): + if mesh.triangles[n_avd].neighbors[k] == n_ab: + mesh.triangles[n_avd].neighbors[k] = t3 + break + + flip_stack.extend([(t0, 2), (t1, 1), (t2, 2), (t3, 1)]) + affected.extend([t0, t1, t2, t3]) + + # Lawson flip loop + while flip_stack: + ft, fj = flip_stack.pop() + tri = mesh.triangles[ft] + if tri.v[0] < 0: + continue + fn = tri.neighbors[fj] + if fn < 0: + continue + + fa, fb, fp = tri.v[fj], tri.v[(fj + 1) % 3], tri.v[(fj + 2) % 3] + + if edge_key(fa, fb) in constrained_edges: + continue + + tri_n = mesh.triangles[fn] + if tri_n.v[0] < 0: + continue + fj2 = -1 + for k in range(3): + if tri_n.neighbors[k] == ft: + fj2 = k + break + if fj2 < 0: + continue + + fq = tri_n.v[(fj2 + 2) % 3] + + if in_circle(mesh.vertices[fa], mesh.vertices[fb], mesh.vertices[fp], + mesh.vertices[fq]) <= 0: + continue + + # No convexity check: with exact predicates, in_circle(fa,fb,fp,fq)>0 + # implies the quad (fa,fq,fb,fp) is convex, so both flipped triangles + # come out CCW. The belt-and-suspenders orient2d pair that lived here + # rejected 0 of millions of flips suite-wide (byte-identical without it) + # and was removed (tangle.cpp 0.4.3). + + n_fb_fp = tri.neighbors[(fj + 1) % 3] + n_fp_fa = tri.neighbors[(fj + 2) % 3] + n_ext1 = tri_n.neighbors[(fj2 + 1) % 3] + n_ext2 = tri_n.neighbors[(fj2 + 2) % 3] + + if tri_n.v[fj2] != fb: + n_ext1, n_ext2 = n_ext2, n_ext1 + + ra, ra2 = tri.region_attrib, tri.region_max_area + rb, rb2 = tri_n.region_attrib, tri_n.region_max_area + + mesh.triangles[ft] = Triangle(v=[fa, fq, fp], neighbors=[n_ext1, fn, n_fp_fa], + region_attrib=ra, region_max_area=ra2) + mesh.triangles[fn] = Triangle(v=[fb, fp, fq], neighbors=[n_fb_fp, ft, n_ext2], + region_attrib=rb, region_max_area=rb2) + + if n_ext1 >= 0: + for k in range(3): + if mesh.triangles[n_ext1].neighbors[k] == fn: + mesh.triangles[n_ext1].neighbors[k] = ft + break + if n_fb_fp >= 0: + for k in range(3): + if mesh.triangles[n_fb_fp].neighbors[k] == ft: + mesh.triangles[n_fb_fp].neighbors[k] = fn + break + + affected.extend([ft, fn]) + flip_stack.extend([(ft, 0), (fn, 2)]) + + # Bump generation + for ai in affected: + if 0 <= ai < len(mesh.triangles): + mesh.triangles[ai].generation += 1 + + return pidx + + +# ============================================================ +# Quality refinement +# ============================================================ + +def tri_metric(mesh, ti, min_ang_threshold, global_max_a, use_region_area, + cos_min_ang=-1.0): + t = mesh.triangles[ti] + if t.v[0] < 0: + return (0, False, False, 0, -1, 0, 0, 0, 0) + a = mesh.vertices[t.v[0]] + b = mesh.vertices[t.v[1]] + c = mesh.vertices[t.v[2]] + + # Squared edge lengths (computed once, reused for both area and angle checks) + la2 = (b.x-c.x)*(b.x-c.x) + (b.y-c.y)*(b.y-c.y) + lb2 = (a.x-c.x)*(a.x-c.x) + (a.y-c.y)*(a.y-c.y) + lc2 = (a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y) + + # Area via cross product (no orient2d exact-arithmetic overhead) + cross = (b.x-a.x)*(c.y-a.y) - (b.y-a.y)*(c.x-a.x) + area = 0.5 * abs(cross) + + effective_max_area = -1 + if use_region_area and t.region_max_area > 0: + effective_max_area = t.region_max_area + if global_max_a > 0: + if effective_max_area > 0: + effective_max_area = min(effective_max_area, global_max_a) + else: + effective_max_area = global_max_a + # LFS constraint: if any vertex has lfs defined, impose area <= sqrt(3)/4 * lfs^2 + for j in range(3): + vl = mesh.vertices[t.v[j]].lfs + if vl > 0: + lfs_area = 0.43301270189221932 * vl * vl # sqrt(3)/4 + if effective_max_area > 0: + effective_max_area = min(effective_max_area, lfs_area) + else: + effective_max_area = lfs_area + area_viol = effective_max_area > 0 and area > effective_max_area + EPS + + # Fast angle check using cosine comparison (no hypot/acos) + # Minimum angle < threshold iff maximum cosine > cosMinAng. + # cos(angle_opposite_la) = (lb²+lc²-la²) / (2·sqrt(lb²·lc²)) + # Compare squared: (lb²+lc²-la²)² > cos²Threshold · 4·lb²·lc² + ang_viol = False + max_cos_ratio = 0.0 + if min_ang_threshold > 0 and cos_min_ang > -0.5: + cos2t = cos_min_ang * cos_min_ang + for num, prod in ((lb2 + lc2 - la2, lb2 * lc2), + (la2 + lc2 - lb2, la2 * lc2), + (la2 + lb2 - lc2, la2 * lb2)): + if num > 0: + denom = cos2t * 4.0 * prod + if denom > 0: + ratio = num * num / denom + if ratio > 1.0: + ang_viol = True + if ratio > max_cos_ratio: + max_cos_ratio = ratio + elif min_ang_threshold > 0: + ang = min_angle(a, b, c) + ang_viol = ang < min_ang_threshold - EPS + if ang_viol: + max_cos_ratio = (min_ang_threshold - ang) / min_ang_threshold + 1.0 + + if not area_viol and not ang_viol: + return (0, False, False, area, effective_max_area, la2, lb2, lc2, cross) + metric = 0 + if area_viol: + metric = area / (effective_max_area if effective_max_area > 0 else 1.0) + if ang_viol: + metric = max(metric, max_cos_ratio) + return (metric, ang_viol, area_viol, area, effective_max_area, la2, lb2, lc2, cross) + + +def refine_quality(mesh, opts, n_input_verts): + if not opts.quality and not opts.area_limit: + return + # Refine to exactly the requested angle. (An earlier +0.1 deg overshoot + # "for floating-point precision" only added needless triangles: the refiner + # never stops at a triangle below the threshold it tests, so the output + # minimum can fall short by at most rounding in the angle evaluation, orders + # of magnitude below 0.1 deg.) + min_ang_val = opts.min_angle if opts.quality else 0.0 + cos_min_ang = math.cos(min_ang_val * math.pi / 180.0) if min_ang_val > 0 else -1.0 + global_max_a = opts.max_area if opts.area_limit else -1.0 + use_region_area = opts.area_limit + + # Estimate final triangle count from CDT using LFS integral with + # circle-of-influence model for refinement seed vertices. + nv = len(mesh.vertices) + nt = len(mesh.triangles) + + # Build vertex-to-triangle adjacency + v2t = [[] for _ in range(nv)] + for i in range(nt): + t = mesh.triangles[i] + if t.v[0] < 0: + continue + for j in range(3): + v2t[t.v[j]].append(i) + + # Bounding box diagonal — used to cap degenerate circumradii + gxmin = min(v.x for v in mesh.vertices) + gxmax = max(v.x for v in mesh.vertices) + gymin = min(v.y for v in mesh.vertices) + gymax = max(v.y for v in mesh.vertices) + bb_diag = math.sqrt((gxmax-gxmin)*(gxmax-gxmin) + (gymax-gymin)*(gymax-gymin)) + + # Compute LFS at each vertex as min circumradius of incident CDT triangles + lfs = [1e30] * nv + for i in range(nt): + t = mesh.triangles[i] + if t.v[0] < 0: + continue + p0, p1, p2 = mesh.vertices[t.v[0]], mesh.vertices[t.v[1]], mesh.vertices[t.v[2]] + area = tri_area(p0, p1, p2) + e0 = math.sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)) + e1 = math.sqrt((p0.x-p2.x)*(p0.x-p2.x)+(p0.y-p2.y)*(p0.y-p2.y)) + e2 = math.sqrt((p0.x-p1.x)*(p0.x-p1.x)+(p0.y-p1.y)*(p0.y-p1.y)) + R = e0 * e1 * e2 / (4.0 * area) if area > 0 else 0 + if R > 0: + R = min(R, bb_diag) + for j in range(3): + lfs[t.v[j]] = min(lfs[t.v[j]], R) + + # Fold in user-assigned LFS constraints + for vi in range(nv): + if mesh.vertices[vi].lfs > 0: + lfs[vi] = min(lfs[vi], mesh.vertices[vi].lfs) + + # Compute background mesh size H at each vertex from area constraints + # and median circumradius of incident CDT triangles + bg_size = [1e30] * nv + for vi in range(nv): + # Area constraint contribution (only if area limiting is active) + amax = global_max_a + if use_region_area: + for ti in v2t[vi]: + t = mesh.triangles[ti] + if t.v[0] < 0: + continue + if t.region_max_area > 0: + amax = min(amax, t.region_max_area) if amax > 0 else t.region_max_area + if amax > 0: + bg_size[vi] = math.sqrt(4.0 * amax / math.sqrt(3.0)) + + # Median circumradius of incident CDT triangles + inc_r = [] + for ti in v2t[vi]: + t = mesh.triangles[ti] + if t.v[0] < 0: + continue + p0, p1, p2 = mesh.vertices[t.v[0]], mesh.vertices[t.v[1]], mesh.vertices[t.v[2]] + a2 = tri_area(p0, p1, p2) + ee0 = math.sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)) + ee1 = math.sqrt((p0.x-p2.x)*(p0.x-p2.x)+(p0.y-p2.y)*(p0.y-p2.y)) + ee2 = math.sqrt((p0.x-p1.x)*(p0.x-p1.x)+(p0.y-p1.y)*(p0.y-p1.y)) + R2 = ee0 * ee1 * ee2 / (4.0 * a2) if a2 > 0 else 0 + if R2 > 0: + inc_r.append(min(R2, bb_diag)) + if inc_r: + inc_r.sort() + med_r = inc_r[len(inc_r) // 2] + bg_size[vi] = min(bg_size[vi], med_r) + + # Pass 1: background mesh estimate from CDT triangles + alpha = math.sqrt(3.0) / 4.0 + pi_over_alpha = math.pi / alpha + cos_min_angle = math.cos(min_ang_val * math.pi / 180.0) + + est_background = 0.0 + n_cdt = 0 + n_bad_cdt = 0 + total_area = 0.0 + for i in range(nt): + t = mesh.triangles[i] + if t.v[0] < 0: + continue + n_cdt += 1 + p0, p1, p2 = mesh.vertices[t.v[0]], mesh.vertices[t.v[1]], mesh.vertices[t.v[2]] + area = tri_area(p0, p1, p2) + total_area += area + H = (bg_size[t.v[0]] + bg_size[t.v[1]] + bg_size[t.v[2]]) / 3.0 + bg_est = area / (alpha * H * H) + # Check for bad angle + e0 = math.sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)) + e1 = math.sqrt((p0.x-p2.x)*(p0.x-p2.x)+(p0.y-p2.y)*(p0.y-p2.y)) + e2 = math.sqrt((p0.x-p1.x)*(p0.x-p1.x)+(p0.y-p1.y)*(p0.y-p1.y)) + cos_a0 = (e1*e1 + e2*e2 - e0*e0) / (2*e1*e2 + 1e-30) + cos_a1 = (e0*e0 + e2*e2 - e1*e1) / (2*e0*e2 + 1e-30) + cos_a2 = (e0*e0 + e1*e1 - e2*e2) / (2*e0*e1 + 1e-30) + cos_max = max(cos_a0, cos_a1, cos_a2) + is_bad = cos_max > cos_min_angle + if is_bad: + n_bad_cdt += 1 + # Angle-area interaction: 1.7x when area dominates on bad triangles + if is_bad and bg_est > 6.0: + est_background += bg_est * 1.7 + else: + est_background += max(bg_est, 6.0) + + # Pass 2: refinement seed contributions (circle-of-influence model) + est_seeds = 0.0 + n_seeds = 0 + for vi in range(nv): + h = lfs[vi] + H = bg_size[vi] + if h >= 1e29 or h >= H * 0.5: + continue + n_seeds += 1 + ratio = H / h + log_r = math.log(ratio) + # Correction: Ruppert's geometric grading creates quadratically more + # triangles than the linear grading model at high H/h ratios. + seed_contrib = pi_over_alpha * (1.0 + 2.0 * log_r) * max(1.0, log_r) + seed_contrib -= pi_over_alpha # subtract background already counted + if seed_contrib > 0: + est_seeds += seed_contrib + + # Pass 3: segment LFS contributions + # A segment with LFS constraint h generates a grading band on both sides. + # Theoretical integral (linear grading) telescopes to 2*L/(alpha*h), but + # Ruppert's geometric grading is ~3x denser, giving 6*L/(alpha*h). + # NOTE: this assumes grading on both sides (interior segment). Boundary + # segments only grade on one side, so this overestimates by ~2x for them. + # Conservative is fine for a safety cap; if a tighter estimate is ever + # needed, halve the contribution for boundary segments (those with only + # one incident triangle on the constrained edge). + est_seg_lfs = 0 + for s in mesh.segments: + if s.lfs <= 0: + continue + Ldx = mesh.vertices[s.v1].x - mesh.vertices[s.v0].x + Ldy = mesh.vertices[s.v1].y - mesh.vertices[s.v0].y + L = math.sqrt(Ldx*Ldx + Ldy*Ldy) + est_seg_lfs += 6.0 * L / (alpha * s.lfs) + + est_total_tris = est_background + est_seeds + est_seg_lfs + max_steiner = max(10000, int(est_total_tris * 5.0)) + if not opts.quiet: + sys.stderr.write(f" CDT: {n_cdt} tris ({n_bad_cdt} bad), area={total_area:.6g}, seeds: {n_seeds}\n") + sys.stderr.write(f" Est triangles: {int(est_total_tris)}, Steiner limit: {max_steiner}\n") + steiner_count = [0] # use list for mutability in closures + + constrained_edges = set() + for s in mesh.segments: + constrained_edges.add(edge_key(s.v0, s.v1)) + + # Detect acute apices and initialize shell tracking (concentric shells) + acute_apices = set() + vert_segs = {} + for si in range(len(mesh.segments)): + for v in (mesh.segments[si].v0, mesh.segments[si].v1): + vert_segs.setdefault(v, []).append(si) + for vi, slist in vert_segs.items(): + if vi >= n_input_verts or len(slist) < 2: + continue + for i in range(len(slist)): + for j in range(i + 1, len(slist)): + oi = mesh.segments[slist[i]].v1 if mesh.segments[slist[i]].v0 == vi else mesh.segments[slist[i]].v0 + oj = mesh.segments[slist[j]].v1 if mesh.segments[slist[j]].v0 == vi else mesh.segments[slist[j]].v0 + dix = mesh.vertices[oi].x - mesh.vertices[vi].x + diy = mesh.vertices[oi].y - mesh.vertices[vi].y + djx = mesh.vertices[oj].x - mesh.vertices[vi].x + djy = mesh.vertices[oj].y - mesh.vertices[vi].y + dot = dix * djx + diy * djy + cross = dix * djy - diy * djx + if dot > 0 and abs(cross) < dot * 1.15: # angle < ~60° + acute_apices.add(vi) + shell_apex = [-1] * len(mesh.vertices) + + gxmin -= 1 + gymin -= 1 + gxmax += 1 + gymax += 1 + # Size the grid so each cell holds ~O(1) segments at the final mesh density, + # not a fixed 50x50. Estimate final segment count = (boundary length)/h, + # h = sqrt(domArea/estTotalTris); target ~sqrt(estSegs) cells/side. Floor 50 + # keeps low-segment meshes from over-building; cap 700 bounds memory. + g_target = 50 + b_len = 0.0 + for s in mesh.segments: + sa, sb = mesh.vertices[s.v0], mesh.vertices[s.v1] + b_len += math.sqrt((sb.x - sa.x) ** 2 + (sb.y - sa.y) ** 2) + dom_a = (gxmax - gxmin) * (gymax - gymin) + h = math.sqrt(dom_a / max(1.0, est_total_tris)) + est_segs = b_len / h if h > 0 else 0.0 + g = int(math.sqrt(max(1.0, est_segs))) + if g > g_target: + g_target = g + if g_target > 700: + g_target = 700 + g_nx, g_ny = g_target, g_target + gcell = max((gxmax - gxmin) / g_nx, (gymax - gymin) / g_ny) + g_nx = int((gxmax - gxmin) / gcell) + 1 + g_ny = int((gymax - gymin) / gcell) + 1 + + seg_grid = {} + + def grid_add_seg(si): + s = mesh.segments[si] + sa, sb = mesh.vertices[s.v0], mesh.vertices[s.v1] + mx, my = (sa.x + sb.x) * 0.5, (sa.y + sb.y) * 0.5 + r = math.sqrt((sb.x-sa.x)*(sb.x-sa.x)+(sb.y-sa.y)*(sb.y-sa.y)) * 0.5 + cx0 = max(0, int((mx - r - gxmin) / gcell)) + cy0 = max(0, int((my - r - gymin) / gcell)) + cx1 = min(g_nx - 1, int((mx + r - gxmin) / gcell)) + cy1 = min(g_ny - 1, int((my + r - gymin) / gcell)) + for cx in range(cx0, cx1 + 1): + for cy in range(cy0, cy1 + 1): + key = cx * g_ny + cy + seg_grid.setdefault(key, []).append(si) + + def grid_remove_seg(si): + s = mesh.segments[si] + sa, sb = mesh.vertices[s.v0], mesh.vertices[s.v1] + mx, my = (sa.x + sb.x) * 0.5, (sa.y + sb.y) * 0.5 + r = math.sqrt((sb.x-sa.x)*(sb.x-sa.x)+(sb.y-sa.y)*(sb.y-sa.y)) * 0.5 + cx0 = max(0, int((mx - r - gxmin) / gcell)) + cy0 = max(0, int((my - r - gymin) / gcell)) + cx1 = min(g_nx - 1, int((mx + r - gxmin) / gcell)) + cy1 = min(g_ny - 1, int((my + r - gymin) / gcell)) + for cx in range(cx0, cx1 + 1): + for cy in range(cy0, cy1 + 1): + key = cx * g_ny + cy + if key in seg_grid: + lst = seg_grid[key] + if si in lst: + lst.remove(si) + + # edge_to_seg map: edge_key -> segment index (for PBC partner lookup) + edge_to_seg = {} + for si in range(len(mesh.segments)): + grid_add_seg(si) + edge_to_seg[edge_key(mesh.segments[si].v0, mesh.segments[si].v1)] = si + + # Vertex -> constrained-subsegment neighbor vertices, maintained across + # splits. Used by the corner-aligned split rule (split_fraction) below. + seg_adj = [[] for _ in range(len(mesh.vertices))] + for s in mesh.segments: + seg_adj[s.v0].append(s.v1) + seg_adj[s.v1].append(s.v0) + + def split_fraction(si): + """Corner-aligned split position (anti-resonance rule). + + Pure midpoint splitting pins the innermost subsegment lengths of a + corner's two arms to L1/2^k and L2/2^j: their folded ratio is fixed by + the input, and the wedge triangle between the innermost rungs is + self-similar under halving with an invariant min angle f(theta, ratio). + If the quality bound exceeds that invariant and the corner angle is near + ~100 deg (where the encroach-split ping-pong sustains), refinement + cascades to the FP floor (MotorB tooth roots: theta=97.757, ratio 5.694 + -> invariant 32.77 deg; see corner_lab.py). + + Cure: when an endpoint is an input-vertex corner in the (88,108)-deg + window and the folded-invariant wedge violates the bound, snap the split + so this arm's corner-incident length equals the neighbor arm's innermost + length (times a power of two, whichever lands in [1/3, 2/3]). Aligned + arms make the wedge isoceles (min angle >= (180-theta)/2) so the + resonance cannot form. Everywhere else: exact midpoint. + """ + if min_ang_val <= 0.0: + return 0.5 # no angle bound -> no resonance possible + seg = mesh.segments[si] + sa, sb = mesh.vertices[seg.v0], mesh.vertices[seg.v1] + sx, sy = sb.x - sa.x, sb.y - sa.y + l2seg = sx * sx + sy * sy + if l2seg <= 0.0: + return 0.5 + t, best_d = 0.5, 1e9 + for side in range(2): + e = seg.v1 if side else seg.v0 + far = seg.v0 if side else seg.v1 + if e >= len(seg_adj): + continue + pe, pf = mesh.vertices[e], mesh.vertices[far] + dx, dy = pf.x - pe.x, pf.y - pe.y + for n in seg_adj[e]: + if n == far: + continue + pn = mesh.vertices[n] + ox, oy = pn.x - pe.x, pn.y - pe.y + r2 = ox * ox + oy * oy + if r2 <= 0.0: + continue + dot = dx * ox + dy * oy + co = dot / math.sqrt(l2seg * r2) # cos(corner angle theta) + # THE guard: snap iff it flips the corner's limit state from + # failing the bound to meeting it. Midpoint dynamics converge to + # the folded-ratio wedge; snap dynamics converge to the aligned + # isoceles wedge (min angle min(theta, (180-theta)/2)). + # (a) The ALIGNED attractor must MEET the bound: + # min_ang <= theta <= 180 - 2*min_ang. This rejects acute + # corners (snapping cannot help; shell/give-up territory) and + # near-collinear junctions (arc-chord chains and every + # Steiner chain vertex), where aligning is meaningless. + if co > cos_min_ang: # theta < min_ang + continue + if co < 1.0 - 2.0 * cos_min_ang * cos_min_ang: # theta > 180-2*min_ang + continue + # (b) The MIDPOINT attractor must FAIL the bound: fold the + # midpoint ratio (L/2)/r by octaves into [1/sqrt2, sqrt2) + # and test the invariant wedge (sides 1 and rho at theta). + # Transiently-bad-but-convergent states stay midpoint. + rho = 0.5 * math.sqrt(l2seg / r2) + while rho < 0.70710678118654752: + rho *= 2.0 + while rho >= 1.4142135623730951: + rho *= 0.5 + c2 = 1.0 + rho * rho - 2.0 * rho * co + if c2 > 0.0: + c = math.sqrt(c2) + cos_a = (1.0 + c2 - rho * rho) / (2.0 * c) # opposite rho + cos_b = (rho * rho + c2 - 1.0) / (2.0 * rho * c) # opposite 1 + if co <= cos_min_ang and cos_a <= cos_min_ang and cos_b <= cos_min_ang: + continue # invariant wedge meets the bound: no snap + f = math.sqrt(r2 / l2seg) # neighbor arm length / this length + for cand in (0.5 * f, f, 2.0 * f): + if cand < 1.0 / 3.0 or cand > 2.0 / 3.0: + continue + d = abs(cand - 0.5) + if d < best_d: + best_d = d + t = 1.0 - cand if side else cand + return t + + def split_segment_core(si, out_affected, hint_tri=-1, t_frac=-1.0): + """Split segment (corner-aligned fraction, or forced t_frac for PBC + partner sync). Returns split-point index or -1 on failure.""" + seg = Segment(mesh.segments[si].v0, mesh.segments[si].v1, + mesh.segments[si].marker, mesh.segments[si].lfs, + mesh.segments[si].pbc_type, mesh.segments[si].no_split, + mesh.segments[si].pbc_side) + sa, sb = mesh.vertices[seg.v0], mesh.vertices[seg.v1] + t = t_frac if 0.0 < t_frac < 1.0 else split_fraction(si) + mx, my = sa.x + (sb.x - sa.x) * t, sa.y + (sb.y - sa.y) * t + mid = Point(mx, my, len(mesh.vertices), seg.marker, [], seg.lfs) + + ek = edge_key(seg.v0, seg.v1) + constrained_edges.discard(ek) + edge_to_seg.pop(ek, None) + grid_remove_seg(si) + + mid_idx = insert_point_lawson(mesh, mid, hint_tri, constrained_edges, out_affected) + if mid_idx < 0: + constrained_edges.add(ek) + edge_to_seg[ek] = si + grid_add_seg(si) + return -1 + steiner_count[0] += 1 + # Propagate shell membership for acute-angle wedge tracking + while len(shell_apex) <= mid_idx: + shell_apex.append(-1) + a0 = seg.v0 if seg.v0 < n_input_verts else shell_apex[seg.v0] + a1 = seg.v1 if seg.v1 < n_input_verts else shell_apex[seg.v1] + if a0 >= 0 and a0 in acute_apices: + shell_apex[mid_idx] = a0 + elif a1 >= 0 and a1 in acute_apices: + shell_apex[mid_idx] = a1 + # Maintain the constrained-neighbor adjacency for split_fraction. + while len(seg_adj) < len(mesh.vertices): + seg_adj.append([]) + seg_adj[seg.v0] = [mid_idx if nb == seg.v1 else nb for nb in seg_adj[seg.v0]] + seg_adj[seg.v1] = [mid_idx if nb == seg.v0 else nb for nb in seg_adj[seg.v1]] + seg_adj[mid_idx] = [seg.v0, seg.v1] + s1 = Segment(seg.v0, mid_idx, seg.marker, seg.lfs, seg.pbc_type, seg.no_split, seg.pbc_side) + s2 = Segment(mid_idx, seg.v1, seg.marker, seg.lfs, seg.pbc_type, seg.no_split, seg.pbc_side) + mesh.segments[si] = s1 + new_si = len(mesh.segments) + mesh.segments.append(s2) + ek1 = edge_key(s1.v0, s1.v1) + ek2 = edge_key(s2.v0, s2.v1) + constrained_edges.add(ek1) + constrained_edges.add(ek2) + edge_to_seg[ek1] = si + edge_to_seg[ek2] = new_si + grid_add_seg(si) + grid_add_seg(new_si) + return mid_idx + + def split_segment(si, out_affected, hint_tri=-1): + """Split segment; for PBC segments, also split the partner.""" + seg = mesh.segments[si] # read before split modifies it + pbc_type = seg.pbc_type + sv0, sv1 = seg.v0, seg.v1 + t_used = split_fraction(si) + mid_idx = split_segment_core(si, out_affected, hint_tri, t_used) + if mid_idx < 0: + return False + + # Synchronized PBC partner split + if pbc_type >= 0: + tv0 = mesh.pbc_twin.get(sv0, -1) + tv1 = mesh.pbc_twin.get(sv1, -1) + if tv0 >= 0 and tv1 >= 0: + pek = edge_key(tv0, tv1) + partner_si = edge_to_seg.get(pek, -1) + if partner_si >= 0: + partner_affected = [] + # Split the partner at the SAME parametric position (oriented: + # partner may be stored with endpoints in either order) so the + # periodic boundaries keep matching node positions. + pt = t_used if mesh.segments[partner_si].v0 == tv0 else 1.0 - t_used + partner_mid = split_segment_core(partner_si, partner_affected, -1, pt) + if partner_mid >= 0: + mesh.pbc_twin[mid_idx] = partner_mid + mesh.pbc_twin[partner_mid] = mid_idx + mesh.pbc_node_type[mid_idx] = pbc_type + mesh.pbc_node_type[partner_mid] = pbc_type + # Add partner-affected triangles so the caller can + # queue them for quality checking. + out_affected.extend(partner_affected) + return True + + def find_encroached(px, py): + gcx = max(0, min(g_nx - 1, int((px - gxmin) / gcell))) + gcy = max(0, min(g_ny - 1, int((py - gymin) / gcell))) + for ddx in range(-1, 2): + for ddy in range(-1, 2): + cx, cy = gcx + ddx, gcy + ddy + if cx < 0 or cx >= g_nx or cy < 0 or cy >= g_ny: + continue + key = cx * g_ny + cy + if key not in seg_grid: + continue + for si in seg_grid[key]: + s = mesh.segments[si] + if s.no_split: + continue # AGE chords must stay uniformly spaced + sa, sb = mesh.vertices[s.v0], mesh.vertices[s.v1] + mx2, my2 = (sa.x + sb.x) * 0.5, (sa.y + sb.y) * 0.5 + sr2 = ((sb.x-sa.x)*(sb.x-sa.x)+(sb.y-sa.y)*(sb.y-sa.y)) * 0.25 + dd2 = (px-mx2)*(px-mx2)+(py-my2)*(py-my2) + if dd2 < sr2: + return si + return -1 + + def walk_blocked(start_tri, px, py): + """Walk straight from start_tri toward (px,py); return the segment index + of the first CONSTRAINED edge the walk must cross (the candidate is + hidden behind it and encroaches it, even one outside the domain), -1 if + the point is reachable (visible), or -2 if the walk fails structurally + (drop the candidate). Without this, blocked candidates fail insertion + and their bad triangles are silently dropped (e.g. a vertex a few microns + off a constrained segment left a 0.107-degree sliver untouched).""" + ti = start_tri + for _ in range(len(mesh.triangles)): + if ti < 0 or ti >= len(mesh.triangles): + return -2 + t = mesh.triangles[ti] + if t.v[0] < 0: + return -2 + exit_edge = -1 + for j in range(3): + va = t.v[j] + vb = t.v[(j + 1) % 3] + if orient2d_xy(mesh.vertices[va].x, mesh.vertices[va].y, + mesh.vertices[vb].x, mesh.vertices[vb].y, + px, py) < 0.0: + exit_edge = j + break + if exit_edge < 0: + return -1 # point inside (or on boundary of) ti: visible + va = t.v[exit_edge] + vb = t.v[(exit_edge + 1) % 3] + if edge_key(va, vb) in constrained_edges: + return edge_to_seg.get(edge_key(va, vb), -2) + ti = t.neighbors[exit_edge] + if ti < 0: + return -2 # unconstrained hull edge (shouldn't happen in a PSLG) + return -2 + + # Priority queue (min-heap): order by SHORTEST EDGE first = smallest squared + # edge length, the order under which off-center insertion guarantees growth + # (Ungor). Tuple = (key, idx, generation, retry). (A scale-gated switch to + # worst-angle ordering for tiny triangles used to live here; with the + # circumcenter test below made scale-free it is unnecessary and actively + # harmful -- it broke the smallest-first growth guarantee and carpeted + # sub-tiny features at -q33.) + REQUEUE_CAP = 4 # capped badIdx re-queue after an encroachment split + GUARD_FACTOR = 0.5 # too-close-insertion guard: reject radius factor + + def prio_key(mr): + return min(mr[5], mr[6], mr[7]) # shortest squared edge + + pq = [] + for i in range(len(mesh.triangles)): + mr = tri_metric(mesh, i, min_ang_val, global_max_a, use_region_area, cos_min_ang) + if mr[0] > 0: + heappush(pq, (prio_key(mr), i, mesh.triangles[i].generation, 0)) + + # Precompute off-center constant (loop-invariant) + offconst = (0.475 * math.sqrt( + (1.0 + math.cos(min_ang_val * math.pi / 180.0)) / + (1.0 - math.cos(min_ang_val * math.pi / 180.0)))) if min_ang_val > 0 else 0 + + dbg_encroach = dbg_insert = dbg_skip = dbg_fail = 0 + + iterations = 0 + while pq and steiner_count[0] < max_steiner: + iterations += 1 + if iterations > mesh.max_vertices * 20: + raise ValueError("Cutout mesh refinement limit exceeded") + _prio, bad_idx, gen, retry = heappop(pq) + + if bad_idx >= len(mesh.triangles): + continue + bt = mesh.triangles[bad_idx] + if bt.v[0] < 0: + continue + if bt.generation != gen: + continue + tmr = tri_metric(mesh, bad_idx, min_ang_val, global_max_a, use_region_area, cos_min_ang) + if tmr[0] <= 0: + continue + ang_viol, area_viol, area, effective_max_area = tmr[1], tmr[2], tmr[3], tmr[4] + tmr_la2, tmr_lb2, tmr_lc2 = tmr[5], tmr[6], tmr[7] + + a = mesh.vertices[bt.v[0]] + b = mesh.vertices[bt.v[1]] + c = mesh.vertices[bt.v[2]] + + # Skip if both edges at smallest-angle vertex are segments + if ang_viol and not area_viol: + # Find vertex with smallest angle = vertex opposite the shortest edge + # Use cached squared edge lengths from tri_metric + la2, lb2, lc2 = tmr_la2, tmr_lb2, tmr_lc2 + min_v = 0 # opposite la (edge b-c) + if lb2 < la2 and lb2 < lc2: + min_v = 1 # opposite lb (edge a-c) + elif lc2 < la2: + min_v = 2 # opposite lc (edge a-b) + va = bt.v[min_v] + vb_v = bt.v[(min_v + 1) % 3] + vc_v = bt.v[(min_v + 2) % 3] + if (edge_key(va, vb_v) in constrained_edges and + edge_key(va, vc_v) in constrained_edges): + dbg_skip += 1 + continue + # Check if worst-angle vertex is inside a small-angle shell. + # Only skip if at least one edge at the vertex is constrained — + # this limits suppression to triangles directly on the shell boundary. + if va < len(shell_apex) and shell_apex[va] >= 0: + if (edge_key(va, vb_v) in constrained_edges or + edge_key(va, vc_v) in constrained_edges): + dbg_skip += 1 + continue + # For tiny triangles (all edges < bb_diag * 1e-9, i.e. at the floor + # where coordinates have ~7 ULPs of room left), accept a reduced + # minimum angle. This engages only for EXACTLY (or near-exactly) + # coincident input geometry that splitting cannot resolve -- anything + # with a representable gap is refined to full quality instead. + tiny_len2 = bb_diag * 1e-9 + tiny_len2 *= tiny_len2 + if la2 < tiny_len2 and lb2 < tiny_len2 and lc2 < tiny_len2: + ang = min_angle(a, b, c) + if ang >= min(min_ang_val, TINY_TRI_ANGLE): + dbg_skip += 1 + continue + + # Compute circumcenter or off-center + if ang_viol: + # Inline circumcenter using cached values from tri_metric + xdo, ydo = b.x - a.x, b.y - a.y + xao, yao = c.x - a.x, c.y - a.y + D = 2.0 * (xdo * yao - ydo * xao) + # Degenerate-circumcenter fallback must be SCALE-FREE: D scales like + # edge^2, so an absolute EPS classified every triangle below ~sqrt(EPS) + # as degenerate and inserted its centroid -> guaranteed insertion-radius + # descent (a carpet) once refinement reaches small scales. + d_mag = abs(xdo * yao) + abs(ydo * xao) + if abs(D) < 1e-12 * d_mag or d_mag == 0.0: + tcx, tcy = (a.x + b.x + c.x) / 3.0, (a.y + b.y + c.y) / 3.0 + else: + # dodist=xdo²+ydo²=lc2, aodist=xao²+yao²=lb2 (from tri_metric) + dodist_cc, aodist_cc = tmr_lc2, tmr_lb2 + tcx = a.x + (yao * dodist_cc - ydo * aodist_cc) / D + tcy = a.y + (xdo * aodist_cc - xao * dodist_cc) / D + dodist, aodist, dadist = tmr_lc2, tmr_lb2, tmr_la2 + ddx_cc, ddy_cc = tcx - a.x, tcy - a.y + if dodist < aodist and dodist < dadist: + dxoff = 0.5 * xdo - offconst * ydo + dyoff = 0.5 * ydo + offconst * xdo + if dxoff * dxoff + dyoff * dyoff < ddx_cc * ddx_cc + ddy_cc * ddy_cc: + tcx, tcy = a.x + dxoff, a.y + dyoff + elif aodist < dadist: + dxoff = 0.5 * xao + offconst * yao + dyoff = 0.5 * yao - offconst * xao + if dxoff * dxoff + dyoff * dyoff < ddx_cc * ddx_cc + ddy_cc * ddy_cc: + tcx, tcy = a.x + dxoff, a.y + dyoff + else: + dxoff = 0.5 * (c.x - b.x) - offconst * (c.y - b.y) + dyoff = 0.5 * (c.y - b.y) + offconst * (c.x - b.x) + if (dxoff * dxoff + dyoff * dyoff < + (ddx_cc - xdo) ** 2 + (ddy_cc - ydo) ** 2): + tcx, tcy = b.x + dxoff, b.y + dyoff + else: + tcx = (a.x + b.x + c.x) / 3.0 + tcy = (a.y + b.y + c.y) / 3.0 + + # Pre-insertion encroachment check: if the circumcenter/off-center + # encroaches a segment, or is hidden from its bad triangle behind a + # constrained edge (walk-blocked, which also covers candidates outside + # the domain), split that segment instead of inserting the point. + enc_seg = find_encroached(tcx, tcy) + if enc_seg < 0: + blk = walk_blocked(bad_idx, tcx, tcy) + if blk == -2: + dbg_fail += 1 + continue + if blk >= 0: + if mesh.segments[blk].no_split: + dbg_skip += 1 + continue + enc_seg = blk + if enc_seg >= 0 and opts.no_steiner and mesh.segments[enc_seg].marker != 0: + continue + if enc_seg >= 0: + seg_aff = [] + if split_segment(enc_seg, seg_aff, bad_idx): + dbg_encroach += 1 + # Triangles modified by the split are re-queued fresh (retry=0) + # via seg_aff below. If badIdx was NOT among them, the split left + # it unchanged but may have removed the segment its circumcenter + # encroached, making it fixable — so re-queue it too with an + # incremented retry. The cap stops the apex cascade. + if (bad_idx not in seg_aff and retry + 1 < REQUEUE_CAP and + bad_idx < len(mesh.triangles) and + mesh.triangles[bad_idx].v[0] >= 0): + mrb = tri_metric(mesh, bad_idx, min_ang_val, global_max_a, + use_region_area, cos_min_ang) + if mrb[0] > 0: + heappush(pq, (prio_key(mrb), bad_idx, + mesh.triangles[bad_idx].generation, retry + 1)) + for ti2 in seg_aff: + if 0 <= ti2 < len(mesh.triangles) and mesh.triangles[ti2].v[0] >= 0: + mr2 = tri_metric(mesh, ti2, min_ang_val, global_max_a, + use_region_area, cos_min_ang) + if mr2[0] > 0: + heappush(pq, (prio_key(mr2), ti2, + mesh.triangles[ti2].generation, 0)) + for k in range(3): + nb = mesh.triangles[ti2].neighbors[k] + if nb >= 0: + mr2 = tri_metric(mesh, nb, min_ang_val, global_max_a, + use_region_area, cos_min_ang) + if mr2[0] > 0: + heappush(pq, (prio_key(mr2), nb, + mesh.triangles[nb].generation, 0)) + continue + + # Insert, with the too-close guard: reject an off-center/circumcenter + # that would land within (GUARD_FACTOR * shortest edge of the bad + # triangle) of an existing vertex. + min_edge2 = min(tmr_la2, tmr_lb2, tmr_lc2) + guard_dist2 = GUARD_FACTOR * GUARD_FACTOR * min_edge2 + np_pt = Point(tcx, tcy, len(mesh.vertices), 0) + affected = [] + new_pt = insert_point_lawson(mesh, np_pt, bad_idx, constrained_edges, + affected, guard_dist2) + if new_pt < 0: + # No -2 vertex-coincident fallback: with the walk-blocked + # encroachment check above and the CDT empty-circumcircle property, + # the circumcenter is at least one circumradius from every visible + # vertex; a refused insertion just tolerates the bad triangle. + dbg_fail += 1 + continue + + steiner_count[0] += 1 + dbg_insert += 1 + + # Re-check only the affected triangles. insert_point_lawson records every + # modified triangle (split products + all flips) in `affected`, and a + # triangle's metric depends only on its own 3 vertices, so unmodified + # neighbors cannot have changed status — no neighbor sweep needed. + for ti in sorted(set(affected)): + if 0 <= ti < len(mesh.triangles): + mr = tri_metric(mesh, ti, min_ang_val, global_max_a, use_region_area, cos_min_ang) + if mr[0] > 0: + heappush(pq, (prio_key(mr), ti, mesh.triangles[ti].generation, 0)) + + if pq: + raise ValueError("Cutout mesh refinement limit exceeded") + if not opts.quiet: + print(f"Quality refinement: {dbg_insert} Steiner points, " + f"{dbg_encroach} segment splits, {dbg_skip} skipped, " + f"{dbg_fail} failed insertions", file=sys.stderr) + + +# ============================================================ +# Edge extraction +# ============================================================ + + +def build_mesh(points, segments, holes, max_area, *, max_vertices=50000): + """Build a constrained mesh in memory, preserving segment markers.""" + global EPS + if len(points) < 3 or len(points) > max_vertices: + raise ValueError("Cutout mesh input exceeds vertex limits") + if not math.isfinite(max_area) or max_area <= 0: + raise ValueError("Cutout mesh area must be positive") + if any(not math.isfinite(v) for p in points for v in p): + raise ValueError("Cutout mesh coordinates must be finite") + if any(a == b or min(a, b) < 0 or max(a, b) >= len(points) for a, b, _ in segments): + raise ValueError("Cutout mesh segment is invalid") + shift_x = (min(p[0] for p in points) + max(p[0] for p in points)) / 2 + shift_y = (min(p[1] for p in points) + max(p[1] for p in points)) / 2 + mesh = Mesh() + mesh.max_vertices = max_vertices + mesh.vertices = [Point(x-shift_x, y-shift_y, i) for i, (x, y) in enumerate(points)] + mesh.segments = [Segment(a, b, marker) for a, b, marker in segments] + mesh.holes = [Hole(x-shift_x, y-shift_y) for x, y in holes] + opts = Options(pslg=True, quality=True, min_angle=10., area_limit=True, max_area=max_area, quiet=True) + with _CALL_LOCK, localcontext() as context: + context.prec = 80 + previous_epsilon = EPS + try: + clean_pslg(mesh, -2.0, True) + for attempt in range(11): + if len(mesh.vertices) > max_vertices: + raise ValueError("Cutout mesh vertex limit exceeded") + mesh.triangles.clear() + build_delaunay(mesh) + mesh.rebuild_adjacency() + missing = enforce_constraints(mesh, True) + if not missing: + break + if attempt == 10: + raise ValueError("Cutout mesh could not enforce its constraints") + users = [[] for _ in mesh.vertices] + for i, triangle in enumerate(mesh.triangles): + if triangle.v[0] >= 0: + for v in triangle.v: + users[v].append(i) + for i in range(len(mesh.segments)-1, -1, -1): + segment = mesh.segments[i] + if mesh.has_edge(segment.v0, segment.v1, users): + continue + a, b = mesh.vertices[segment.v0], mesh.vertices[segment.v1] + index = len(mesh.vertices) + mesh.vertices.append(Point((a.x+b.x)/2, (a.y+b.y)/2, index)) + mesh.segments[i:i+1] = [Segment(segment.v0, index, segment.marker), Segment(index, segment.v1, segment.marker)] + remove_holes(mesh, opts) + if not mesh.triangles: + raise ValueError("Cutout mesh encloses no area") + refine_quality(mesh, opts, len(mesh.vertices)) + finally: + EPS = previous_epsilon + return ( + [(p.x+shift_x, p.y+shift_y) for p in mesh.vertices], + [t.v for t in mesh.triangles if t.v[0] >= 0], + [(s.v0, s.v1, s.marker) for s in mesh.segments], + ) diff --git a/src/anyimage/preferences.py b/src/anyimage/preferences.py index d8534eb..992d434 100644 --- a/src/anyimage/preferences.py +++ b/src/anyimage/preferences.py @@ -170,6 +170,8 @@ class AnyImagePreferences(bpy.types.AddonPreferences): default=1.0, min=0.001, max=100.0, + soft_min=0.5, + soft_max=10.0, subtype="PERCENTAGE", ) max_frame_resolution: bpy.props.IntProperty( diff --git a/src/anyimage/properties.py b/src/anyimage/properties.py index 9d52eec..7ec6d87 100644 --- a/src/anyimage/properties.py +++ b/src/anyimage/properties.py @@ -237,7 +237,7 @@ def cutout_alpha_threshold_property(): def cutout_fine_outline_property(): return bpy.props.BoolProperty( name="Fine Outline", - description="Fit detailed outlines and support thin parts; disable for simpler, faster meshes", + description="Preserve detailed outlines and meaningful holes; disable to fill holes and simplify outlines", default=False, ) @@ -265,6 +265,8 @@ def cutout_edge_length_property(*, store_meters=False, **options): description="Target mesh edge length in world space", default=0.1, min=1e-6, + soft_min=0.01, + soft_max=1.0, unit="LENGTH", subtype="DISTANCE", **options, diff --git a/tests/anyimage/operators/test_cutout_depth.py b/tests/anyimage/operators/test_cutout_depth.py index 71f5541..8669b25 100644 --- a/tests/anyimage/operators/test_cutout_depth.py +++ b/tests/anyimage/operators/test_cutout_depth.py @@ -3,7 +3,7 @@ import numpy as np import pytest from anyimage.operators.cutout_tool.object import load_cutout_node_group from anyimage.operators.cutout_tool.operators import normal_mode_for_shape -from anyimage.operators.cutout_tool.geometry import refine_balloon_triangles +from anyimage.operators.cutout_tool.mesh import build_mesh from anyimage.operators.cutout_tool.balloon import poisson_balloon_profile @@ -38,9 +38,8 @@ def test_depth_cutout_loads_uniform_thickness_revision_and_preserves_existing_gr assert load_cutout_node_group("DEPTH_SOLID") == current -def test_refined_boundary_ears_have_interior_balloon_samples(): - points = np.array([(0, 0), (1, 0), (1, 1), (0, 1)], dtype=float) - points, faces, boundary = refine_balloon_triangles(points, [(0, 1, 2), (0, 2, 3)], np.ones(4, bool)) +def test_small_quality_mesh_has_positive_balloon_support(): + points, faces, boundary, _ = build_mesh(np.ones((2, 2)), 32, .9) profile = poisson_balloon_profile(points, faces, boundary) assert np.all(profile[boundary] == 0) - assert all(any(profile[index] > 0 for index in face) for face in faces) + assert profile.max() > 0 diff --git a/tests/anyimage/operators/test_cutout_fine_mesh.py b/tests/anyimage/operators/test_cutout_fine_mesh.py deleted file mode 100644 index 3d73e6a..0000000 --- a/tests/anyimage/operators/test_cutout_fine_mesh.py +++ /dev/null @@ -1,300 +0,0 @@ -import numpy as np -import pytest -from scipy.sparse import coo_matrix -from scipy.sparse.csgraph import connected_components -from scipy.spatial import cKDTree - -from anyimage.operators.cutout_tool import geometry as g -from anyimage.operators.cutout_tool.balloon import poisson_balloon_profile -from anyimage.operators.cutout_tool import fine_mesh as fine - - -def branch_alpha(offset=0): - yy, xx = np.mgrid[:192, :192] - xx = xx - offset - yy = yy - offset - alpha = (xx - 46) ** 2 + (yy - 48) ** 2 < 30**2 - path = np.asarray(((48, 48), (91, 70), (130, 110), (149, 167))) - for a, b in zip(path[:-1], path[1:]): - d = np.stack((xx - a[0], yy - a[1]), axis=-1) - v = b - a - t = np.clip(d @ v / (v @ v), 0, 1) - alpha |= np.sum((d - t[..., None] * v) ** 2, axis=-1) < 2.5**2 - alpha |= (xx - 130) ** 2 + (yy - 110) ** 2 < 7**2 - return alpha.astype(float) - - -def assert_supported(points, faces, boundary, components): - triangles = points[faces] - area = np.abs( - np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0]) - ) - assert np.isfinite(points).all() and np.all(area > 1e-8) - assert not boundary[faces].all(axis=1).any() - edges = np.sort(np.concatenate((faces[:, :2], faces[:, 1:], faces[:, ::2])), axis=1) - edges, counts = np.unique(edges, axis=0, return_counts=True) - assert counts.max() <= 2 - assert not boundary[edges[counts == 2]].all(axis=1).any() - interior = np.flatnonzero(~boundary) - links = edges[~boundary[edges].any(axis=1)] - graph = coo_matrix( - (np.ones(len(links)), (links[:, 0], links[:, 1])), - shape=(len(points), len(points)), - ).tocsr() - count, _labels = connected_components(graph[interior][:, interior], directed=False) - assert count == len(components) - profile = poisson_balloon_profile(points, faces, boundary) - assert np.isfinite(profile).all() - assert np.all(profile[interior] > 0) - - -@pytest.mark.parametrize("spacing", (4, 8, 16, 28, 32, 36)) -@pytest.mark.parametrize("offset", (0, 1)) -def test_fine_bent_leg_and_joint_have_connected_positive_support(spacing, offset): - points, faces, boundary, components = fine.build_fine_mesh( - branch_alpha(offset), spacing, 0.9 - ) - assert_supported(points, faces, boundary, components) - interior = points[~boundary] - for checkpoint in ((90, 70), (130, 110), (146, 159)): - assert cKDTree(interior).query(np.asarray(checkpoint) + offset)[0] < 12 - - -@pytest.mark.parametrize("spacing", (4, 8, 16, 32)) -def test_fine_comb_hole_and_one_pixel_bridge(spacing): - alpha = np.zeros((160, 180)) - alpha[15:145, 15:42] = 1 - alpha[45:115, 125:165] = 1 - alpha[75, 40:126] = 1 - alpha[65:90, 140:150] = 0 - for row in range(20, 140, 12): - alpha[row : row + 2, 35:110] = 1 - result = fine.build_fine_mesh(alpha, spacing, 0.9) - assert_supported(*result) - points, faces, _boundary, components = result - assert len(components) == 1 - assert not np.any( - g.points_inside_contours( - points[faces].mean(axis=1), - [np.asarray(((140, 65), (150, 65), (150, 90), (140, 90)))], - ) - ) - - -def test_fine_cleanup_uses_threshold_and_preserves_input_and_largest_island(): - alpha = np.zeros((40, 60)) - alpha[5:25, 5:25] = 0.5 - alpha[12, 20:45] = 0.5 - alpha[32:34, 50:52] = 1 - original = alpha.copy() - filtered, stats = fine.filter_alpha(alpha, 0.3) - assert stats["components_after"] == 1 - assert filtered[12, 44] == 0.5 and filtered[32, 50] == 0 - np.testing.assert_array_equal(alpha, original) - assert_supported(*fine.build_fine_mesh(alpha, 32, 0.3)) - assert_supported(*fine.build_fine_mesh(alpha, 32, 0.9)) - with pytest.raises(ValueError, match="no visible"): - fine.build_fine_mesh(np.zeros((10, 10)), 32, 0.9) - - -def test_higher_detail_recovers_curve_instead_of_splitting_low_edges(): - yy, xx = np.mgrid[:240, :240] - alpha = ((xx - 120) ** 2 + (yy - 120) ** 2 < 92**2).astype(float) - reference = g.extract_selection_contours(alpha, max_resolution=240, threshold=0.9) - source = g.resample_closed_contour(np.asarray(reference.components[0][0]), 0.5) - errors = [] - for spacing in (32, 8, 4): - contours = fine.contours_for_component(reference.components[0], spacing) - sampled = g.resample_closed_contour(contours[0], 0.5) - forward = fine.ContourDomain(contours).clearance(source) - backward = fine.ContourDomain( - [np.asarray(reference.components[0][0])] - ).clearance(sampled) - errors.append(max(np.percentile(forward, 95), np.percentile(backward, 95))) - assert errors[1] < errors[0] - assert errors[2] <= errors[1] + 0.1 - - -def test_segment_validation_rejects_holes_tangency_and_outside(): - domain = fine.ContourDomain( - [ - np.asarray(((0, 0), (20, 0), (20, 20), (0, 20))), - np.asarray(((8, 8), (12, 8), (12, 12), (8, 12))), - ] - ) - starts = np.asarray(((2, 2), (2, 10), (2, 8), (-1, 2), (2, 2))) - ends = np.asarray(((18, 2), (18, 10), (18, 8), (2, 2), (2, 2))) - assert domain.valid_segments(starts, ends).tolist() == [ - True, - False, - False, - False, - True, - ] - - -def test_closed_mode_uses_original_geometry_without_fine_import(monkeypatch): - from types import SimpleNamespace - from mathutils import Matrix - from anyimage.common.image import image_empty_bounds - - source = SimpleNamespace( - data=SimpleNamespace(size=(80, 60)), - empty_display_size=2.0, - empty_image_offset=(-0.5, -0.5), - matrix_world=Matrix.Identity(4), - ) - monkeypatch.setattr( - g, - "addon_preferences", - lambda _context: SimpleNamespace(balloon_profile_algorithm="POISSON"), - ) - monkeypatch.setattr( - fine, - "build_fine_mesh", - lambda *_args: pytest.fail("Fine algorithm called while disabled"), - ) - alpha = np.ones((60, 80)) - actual = g.build_base_shape( - None, source, 0.4, alpha, (0, 0, 80, 60), fine_outline=False - ) - contours = g.extract_cutout_contours(alpha, 5, threshold=0.9) - expected = g.cutout_domain( - contours, - image_empty_bounds(source), - 0.4, - ((1.0, 0.0), (0.0, 1.0)), - ) - np.testing.assert_array_equal(actual.vertices, expected[0]) - np.testing.assert_array_equal(actual.uv, expected[2]) - np.testing.assert_array_equal(actual.balloon, expected[-1]) - - -def test_width_index_matches_exhaustive_covering_disks(): - from scipy import ndimage - - width = fine.WidthSamples(branch_alpha().astype(bool)) - points = np.asarray(((48, 48), (85, 70), (130, 110), (147, 160), (2, 2))) - radii = [] - for point in points: - covered = np.sum((width.centers - point) ** 2, axis=1) <= width.reach_squared - radii.append(max(0.5, width.radius[covered].max() if covered.any() else 0.5)) - distances = ndimage.map_coordinates( - width.distance, - [points[:, 1] - 0.5, points[:, 0] - 0.5], - order=1, - mode="nearest", - ) - expected = np.maximum( - 0.7, np.minimum(32, np.minimum(16, np.asarray(radii) * 0.9) + 0.4 * distances) - ) - np.testing.assert_allclose(width(points), expected) - - -def test_relaxation_reduces_sampling_error_without_adding_points(): - contours = [np.asarray(((4.0, 4.0), (194.0, 4.0), (194.0, 104.0), (4.0, 104.0)))] - spine = np.empty((0, 2)) - original = g.uniform_interior_points_for_contours(contours, 16) - relaxed = fine.relax_points(contours, spine, 16) - yy, xx = np.mgrid[6:104:4, 6:194:4] - samples = np.column_stack((xx.ravel() + 0.5, yy.ravel() + 0.5)) - - def energy(free): - return np.mean( - cKDTree(np.concatenate((*contours, free))).query(samples)[0] ** 2 - ) - - assert len(relaxed) == len(original) - assert energy(relaxed) < energy(original) - assert fine.ContourDomain(contours).contains(relaxed).all() - - -def test_cdt_preserves_constraint_split_by_an_interior_vertex(): - contours = [np.asarray(((0.0, 0.0), (20.0, 0.0), (20.0, 10.0), (0.0, 10.0)))] - spine = np.asarray(((2.0, 5.0), (18.0, 5.0))) - points, faces, constraints = fine.triangulate( - contours, spine, [(0, 1)], np.asarray(((10.0, 5.0),)), 8 - ) - assert len(constraints) == 2 - assert np.isclose( - np.linalg.norm( - points[constraints[:, 1]] - points[constraints[:, 0]], axis=1 - ).sum(), - 16, - ) - assert not g.boundary_vertices_from_faces(len(points), faces)[constraints].any() - - -@pytest.mark.parametrize("offset", [0, 4096, 16384]) -def test_cdt_collinear_outline_samples_do_not_become_sliver_faces(offset): - contour = np.asarray(((0., 0.), (240., 0.), (240., 160.), (0., 160.))) - contour = contour @ np.asarray(((.763, -.646), (.646, .763))) - contour += np.asarray((420.123456, 540.234567)) + offset - contours, spine, edges = fine.subdivide_constraints( - [contour], np.empty((0, 2)), [], 16 - ) - free = g.uniform_interior_points_for_contours(contours, 16) - points, faces, _ = fine.triangulate(contours, spine, edges, free, 16) - triangles = points[faces] - area = abs(np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0])) - longest_squared = np.max(np.sum((triangles - np.roll(triangles, 1, axis=1))**2, axis=2), axis=1) - assert np.min(area / longest_squared) > 1e-4 - assert area.sum() / 2 == pytest.approx(abs(g.contour_area(contour)), rel=1e-9) - boundary = g.boundary_vertices_from_faces(len(points), faces) - np.testing.assert_allclose(fine.ContourDomain([contour]).clearance(points[boundary]), 0, atol=1e-9) - # Every subdivided outline point survives; removing slivers must not skip it. - assert cKDTree(points[boundary]).query(contours[0])[0].max() < 1e-9 - - -def test_contour_validation_accepts_pixel_contacts_and_rejects_crossings(): - crossing = np.asarray(((0, 0), (4, 4), (0, 4), (4, 0))) - assert not fine.valid_contours([crossing]) - touching = np.asarray( - ((0, 0), (2, 0), (2, 2), (4, 2), (4, 4), (2, 4), (2, 2), (0, 2)) - ) - assert fine.valid_contours([touching]) - - -@pytest.mark.parametrize("offset", [0, 2048, 4096, 8192, 16384]) -@pytest.mark.parametrize("spacing", [4, 0.5]) -def test_cdt_preserves_short_constraints_at_large_pixel_coordinates(offset, spacing): - contour = np.asarray(((0.0, 0.0), (20.0, 0.0), (20.0, 20.0), (0.0, 20.0))) + offset - spine = np.asarray(((5.1234567, 5.654321), (6.456789, 6.234567))) + offset - free = g.uniform_interior_points_for_contours([contour], spacing) - points, faces, constraints = fine.triangulate( - [contour], spine, [(0, 1)], free, spacing - ) - assert len(constraints) >= 1 - assert not g.boundary_vertices_from_faces(len(points), faces)[constraints].any() - - -@pytest.mark.parametrize( - "damage", ["missing_edge", "missing_faces", "missing_endpoint"] -) -def test_cdt_rejects_broken_constraint_coverage(monkeypatch, damage): - real_cdt = fine.delaunay_2d_cdt - - def damaged_cdt(*args): - result = list(real_cdt(*args)) - result[2] = list(result[2]) - result[3] = [list(ids) for ids in result[3]] - result[4] = [list(ids) for ids in result[4]] - constraint = len(args[1]) - 1 - edge_id = next(i for i, ids in enumerate(result[4]) if constraint in ids) - if damage == "missing_edge": - result[4][edge_id].remove(constraint) - elif damage == "missing_faces": - a, b = result[1][edge_id] - result[2] = [face for face in result[2] if not (a in face and b in face)] - else: - endpoint = args[1][-1][1] - for originals in result[3]: - if endpoint in originals: - originals.remove(endpoint) - return result - - monkeypatch.setattr(fine, "delaunay_2d_cdt", damaged_cdt) - contour = np.asarray(((0.0, 0.0), (20.0, 0.0), (20.0, 10.0), (0.0, 10.0))) - spine = np.asarray(((2.0, 5.0), (18.0, 5.0))) - with pytest.raises(ValueError, match="lost an interior constraint"): - fine.triangulate([contour], spine, [(0, 1)], np.asarray(((10.0, 5.0),)), 8) diff --git a/tests/anyimage/operators/test_cutout_geometry.py b/tests/anyimage/operators/test_cutout_geometry.py index 553408d..0dc9a70 100644 --- a/tests/anyimage/operators/test_cutout_geometry.py +++ b/tests/anyimage/operators/test_cutout_geometry.py @@ -175,7 +175,7 @@ def test_cutout_cleanup_value_scales_the_ui_without_changing_node_area(): def test_cutout_domain_adapts_boundary_spacing_for_a_small_valid_selection(): components = ((((2046, 2046), (2050, 2046), (2050, 2050), (2046, 2050)),),) - domain = cutout_geometry.cutout_domain( + domain = quality_domain( SelectionContours(components, (4096, 4096)), (-1.0, 1.0, -1.0, 1.0), 0.5, @@ -189,7 +189,7 @@ def test_cutout_domain_adapts_boundary_spacing_for_a_small_valid_selection(): def test_cutout_domain_builds_a_full_image_at_minimum_detail(): components = ((((0, 0), (4096, 0), (4096, 4096), (0, 4096)),),) - domain = cutout_geometry.cutout_domain( + domain = quality_domain( SelectionContours(components, (4096, 4096)), (-1.0, 1.0, -1.0, 1.0), 0.5, @@ -203,7 +203,7 @@ def test_cutout_domain_builds_a_full_image_at_minimum_detail(): def test_cutout_domain_centers_vertices_on_crop_bounds(): components = ((((0, 0), (100, 0), (100, 50), (0, 50)),),) - vertices, _faces, uv, _boundary, _planar, _profile = cutout_geometry.cutout_domain( + vertices, _faces, uv, _boundary, _planar, _profile = quality_domain( SelectionContours(components, (100, 50)), (2.0, 6.0, -3.0, 1.0), 0.5, @@ -224,7 +224,7 @@ def test_cutout_domain_keeps_overlapping_selection_components_disconnected(): ), (((0, 20), (40, 20), (40, 50), (0, 50)),), ) - domain = cutout_geometry.cutout_domain( + domain = quality_domain( SelectionContours(components, (120, 80)), (0.0, 1.0, 0.0, 1.0), 0.025, @@ -286,3 +286,19 @@ class CutoutContourTest(unittest.TestCase): return float(np.std(points[selected, 1] - points[selected, 0])) self.assertLess(diagonal_noise(smooth), 0.05) + + +def quality_domain(selection, bounds, spacing, metric): + from anyimage.operators.cutout_tool import mesh + ps, fs, cs = [], [], [] + offset = 0 + width, height = selection.image_size + pixel_scale = max((bounds[1]-bounds[0])/width, (bounds[3]-bounds[2])/height) + for component in selection.components: + contours = [np.array(c, dtype=float) for c in component] + p, f = mesh.triangulate(contours, spacing/pixel_scale) + p, f = mesh.support_mesh(p, f) + ps.append(p);fs.append(f+offset) + cs.append((offset, offset+len(p), f, contours[0]));offset += len(p) + p, f = np.concatenate(ps), np.concatenate(fs) + return cutout_geometry.cutout_domain((p, f, cutout_geometry.boundary_vertices_from_faces(len(p), f), cs, selection.image_size), bounds, metric) diff --git a/tests/anyimage/operators/test_cutout_mesh.py b/tests/anyimage/operators/test_cutout_mesh.py new file mode 100644 index 0000000..ea7ab8b --- /dev/null +++ b/tests/anyimage/operators/test_cutout_mesh.py @@ -0,0 +1,334 @@ +import numpy as np +import pytest +from scipy.sparse import coo_matrix +from scipy.sparse.csgraph import connected_components +from scipy.spatial import cKDTree + +from anyimage.operators.cutout_tool import geometry as g +from anyimage.operators.cutout_tool.balloon import poisson_balloon_profile +from anyimage.operators.cutout_tool import mesh + + +def boundary_distance(points, contours): + starts = np.concatenate(contours) + delta = np.concatenate([np.roll(c, -1, axis=0) - c for c in contours]) + offset = points[:, None] - starts + fraction = np.clip(np.sum(offset * delta, axis=2) / np.sum(delta**2, axis=1), 0, 1) + return np.linalg.norm(offset - fraction[..., None] * delta, axis=2).min(axis=1) + + +@pytest.mark.parametrize("fine_outline", [False, True]) +@pytest.mark.parametrize("spacing", [2, 12, 48, 96]) +def test_local_support_connects_coarse_and_dense_branches(fine_outline, spacing): + result = mesh.build_mesh(branch_alpha(), spacing, 0.9, fine_outline=fine_outline) + assert_supported(*result) + + +def test_support_splits_shared_chords_and_preserves_boundary_and_area(): + points = np.array(((0., 0.), (20., 0.), (20., 2.), (0., 2.))) + faces = np.array(((0, 1, 2), (0, 2, 3))) + p, f = mesh.support_mesh(points, faces) + assert len(p) == 5 and len(f) == 4 + assert_supported(p, f, g.boundary_vertices_from_faces(len(p), f), [None]) + triangles = p[f] + assert np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0]).sum() == pytest.approx(80) + np.testing.assert_array_equal(p[:4], points) + with pytest.raises(ValueError, match="limit"): + mesh.support_mesh(points, faces, max_vertices=4) + + +def test_support_adds_one_sample_to_isolated_triangle_and_is_idempotent(): + points = np.array(((0., 0.), (20., 0.), (0., 2.))) + faces = np.array(((0, 1, 2),)) + p, f = mesh.support_mesh(points, faces) + assert len(p) == 4 and len(f) == 3 + assert_supported(p, f, g.boundary_vertices_from_faces(len(p), f), [None]) + again_p, again_f = mesh.support_mesh(p, f) + np.testing.assert_array_equal(p, again_p) + np.testing.assert_array_equal(f, again_f) + + +@pytest.mark.parametrize("spacing", [1, 4, 16]) +def test_pixel_corner_hole_contacts_produce_valid_supported_mesh(spacing): + alpha = np.ones((40, 40)) + alpha[10:17, 10:17] = 0 + alpha[17:24, 17:24] = 0 + original = alpha.copy() + result = mesh.build_mesh(alpha, spacing, .9) + assert_supported(*result) + np.testing.assert_array_equal(alpha, original) + p, f, b, _ = result + edges, counts = np.unique(np.sort(np.concatenate((f[:, :2], f[:, 1:], f[:, ::2])), axis=1), axis=0, return_counts=True) + assert np.all(np.bincount(edges[counts == 1].ravel(), minlength=len(p))[b] == 2) + centers = p[f].mean(axis=1) + # The pixel-scale smoothing may round the hole corners; its interior stays open. + assert not np.any((centers[:, 0] > 11) & (centers[:, 0] < 16) & (centers[:, 1] > 11) & (centers[:, 1] < 16)) + + +@pytest.mark.parametrize("damage", ["missing_edge", "missing_faces", "missing_endpoint"]) +def test_quality_mesh_rejects_broken_boundary_coverage(monkeypatch, damage): + from anyimage.operators.cutout_tool import tangle + + original = tangle.build_mesh + def damaged(*args, **kwargs): + points, faces, segments = original(*args, **kwargs) + if damage == "missing_edge": + segments = [segment for segment in segments if segment[2] != 1] + elif damage == "missing_faces": + faces = faces[1:] + else: + _, b, _ = next(segment for segment in segments if segment[2] == 1) + points[b] = (points[b][0] + .05, points[b][1] + .05) + return points, faces, segments + monkeypatch.setattr(tangle, "build_mesh", damaged) + contour = np.array(((0., 0.), (20., 0.), (20., 10.), (0., 10.))) + with pytest.raises(ValueError): + mesh.triangulate([contour], 8) + + +def branch_alpha(offset=0): + yy, xx = np.mgrid[:192, :192] + xx = xx - offset + yy = yy - offset + alpha = (xx - 46) ** 2 + (yy - 48) ** 2 < 30**2 + path = np.asarray(((48, 48), (91, 70), (130, 110), (149, 167))) + for a, b in zip(path[:-1], path[1:]): + d = np.stack((xx - a[0], yy - a[1]), axis=-1) + v = b - a + t = np.clip(d @ v / (v @ v), 0, 1) + alpha |= np.sum((d - t[..., None] * v) ** 2, axis=-1) < 2.5**2 + alpha |= (xx - 130) ** 2 + (yy - 110) ** 2 < 7**2 + return alpha.astype(float) + + +def assert_supported(points, faces, boundary, components): + triangles = points[faces] + area = np.abs( + np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0]) + ) + assert np.isfinite(points).all() and np.all(area > 1e-8) + edges = np.sort(np.concatenate((faces[:, :2], faces[:, 1:], faces[:, ::2])), axis=1) + edges, counts = np.unique(edges, axis=0, return_counts=True) + assert counts.max() <= 2 + interior = np.flatnonzero(~boundary) + links = edges[~boundary[edges].any(axis=1)] + graph = coo_matrix( + (np.ones(len(links)), (links[:, 0], links[:, 1])), + shape=(len(points), len(points)), + ).tocsr() + count, _labels = connected_components(graph[interior][:, interior], directed=False) + assert count == len(components) + assert not np.any(boundary[faces].all(axis=1)) + assert len(interior) >= len(components) + profile = poisson_balloon_profile(points, faces, boundary) + assert np.isfinite(profile).all() + assert np.all(profile[interior] > 0) + + +@pytest.mark.parametrize("spacing", (4, 8, 16, 28, 32, 36)) +@pytest.mark.parametrize("offset", (0, 1)) +def test_fine_bent_leg_and_joint_have_connected_positive_support(spacing, offset): + points, faces, boundary, components = mesh.build_mesh( + branch_alpha(offset), spacing, 0.9 + ) + assert_supported(points, faces, boundary, components) + interior = points[~boundary] + checkpoints = np.array(((48, 48), (90, 70), (130, 110), (146, 159))) + offset + distance, nearest = cKDTree(interior).query(checkpoints) + assert distance.max() < 12 + edges = np.unique( + np.sort(np.concatenate((faces[:, :2], faces[:, 1:], faces[:, ::2])), axis=1), + axis=0, + ) + links = edges[~boundary[edges].any(1)] + graph = coo_matrix( + (np.ones(len(links)), (links[:, 0], links[:, 1])), + shape=(len(points), len(points)), + ).tocsr() + _, labels = connected_components(graph, directed=False) + assert len(set(labels[np.flatnonzero(~boundary)[nearest]])) == 1 + + +@pytest.mark.parametrize("spacing", (4, 8, 16, 32)) +def test_fine_comb_hole_and_one_pixel_bridge(spacing): + alpha = np.zeros((160, 180)) + alpha[15:145, 15:42] = 1 + alpha[45:115, 125:165] = 1 + alpha[75, 40:126] = 1 + alpha[65:90, 140:150] = 0 + for row in range(20, 140, 12): + alpha[row : row + 2, 35:110] = 1 + result = mesh.build_mesh(alpha, spacing, 0.9) + assert_supported(*result) + points, faces, _boundary, components = result + assert len(components) == 1 + assert not np.any( + g.points_inside_contours( + points[faces].mean(axis=1), + [np.asarray(((140, 65), (150, 65), (150, 90), (140, 90)))], + ) + ) + + +def test_fine_cleanup_uses_threshold_and_preserves_input_and_largest_island(): + alpha = np.zeros((40, 60)) + alpha[5:25, 5:25] = 0.5 + alpha[12, 20:45] = 0.5 + alpha[32:34, 50:52] = 1 + original = alpha.copy() + filtered = mesh.filter_alpha(alpha, 0.3) + assert filtered[12, 44] == 0.5 and filtered[32, 50] == 0 + np.testing.assert_array_equal(alpha, original) + assert_supported(*mesh.build_mesh(alpha, 32, 0.3)) + assert_supported(*mesh.build_mesh(alpha, 32, 0.9)) + with pytest.raises(ValueError, match="no visible"): + mesh.build_mesh(np.zeros((10, 10)), 32, 0.9) + + +def test_higher_detail_recovers_curve_instead_of_splitting_low_edges(): + yy, xx = np.mgrid[:240, :240] + alpha = ((xx - 120) ** 2 + (yy - 120) ** 2 < 92**2).astype(float) + reference = g.extract_selection_contours(alpha, max_resolution=240, threshold=0.9) + source = g.resample_closed_contour(np.asarray(reference.components[0][0]), 0.5) + errors = [] + for spacing in (32, 8, 4): + contours = mesh.contours_for_component(reference.components[0], spacing) + sampled = g.resample_closed_contour(contours[0], 0.5) + forward = boundary_distance(source, contours) + backward = boundary_distance(sampled, [np.asarray(reference.components[0][0])]) + errors.append(max(np.percentile(forward, 95), np.percentile(backward, 95))) + assert errors[1] < errors[0] + assert errors[2] <= errors[1] + 0.1 + + +def test_fine_outline_softens_pixel_steps_in_normalized_coordinates(): + yy, xx = np.mgrid[:120, :120] + alpha = ((xx - 60)**2 + (yy - 60)**2 < 45**2).astype(float) + source = np.asarray(g.extract_selection_contours(alpha, max_resolution=120).components[0][0]) + transform = np.diag((8., 12.)) + contours = mesh.contours_for_component([source @ transform], 32, pixel_transform=transform) + smoothed = contours[0] @ np.linalg.inv(transform) + + def turning(points): + edges = np.roll(points, -1, axis=0) - points + angles = np.arctan2(edges[:, 1], edges[:, 0]) + change = (np.roll(angles, -1) - angles + np.pi) % (2*np.pi) - np.pi + return abs(change).sum() + + assert turning(smoothed) < turning(source) * .5 + assert boundary_distance(smoothed, [source]).max() < 1 + assert abs(g.contour_area(smoothed) / g.contour_area(source) - 1) < .02 + + +@pytest.mark.parametrize("offset", [0, 4096, 16384]) +def test_cdt_collinear_outline_samples_do_not_become_sliver_faces(offset): + contour = np.asarray(((0.0, 0.0), (240.0, 0.0), (240.0, 160.0), (0.0, 160.0))) + contour = contour @ np.asarray(((0.763, -0.646), (0.646, 0.763))) + contour += np.asarray((420.123456, 540.234567)) + offset + contours = [g.resample_closed_contour(contour, 16)] + points, faces = mesh.triangulate(contours, 16) + triangles = points[faces] + area = abs( + np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0]) + ) + longest_squared = np.max( + np.sum((triangles - np.roll(triangles, 1, axis=1)) ** 2, axis=2), axis=1 + ) + assert np.min(area / longest_squared) > 1e-4 + assert area.sum() / 2 == pytest.approx(abs(g.contour_area(contour)), rel=1e-9) + boundary = g.boundary_vertices_from_faces(len(points), faces) + np.testing.assert_allclose( + boundary_distance(points[boundary], [contour]), 0, atol=1e-9 + ) + # Every subdivided outline point survives; removing slivers must not skip it. + assert cKDTree(points[boundary]).query(contours[0])[0].max() < 1e-9 + + +def test_contour_validation_requires_separation_of_pixel_contacts(): + crossing = np.asarray(((0, 0), (4, 4), (0, 4), (4, 0))) + assert not mesh.valid_contours([crossing]) + touching = np.asarray( + ((0, 0), (2, 0), (2, 2), (4, 2), (4, 4), (2, 4), (2, 2), (0, 2)) + ) + assert not mesh.valid_contours([touching]) + assert mesh.valid_contours(mesh.separate_contour_contacts([touching.astype(float)])) + + +@pytest.mark.parametrize("fine_outline", [False, True]) +def test_hole_policy_preserves_alpha_and_external_notches(fine_outline): + alpha = np.ones((80, 80)) + alpha[10:12, 10:12] = 0 + alpha[20:29, 20:29] = 0 + alpha[40:41, 10:40] = 0 + alpha[:15, 60:65] = 0 + original = alpha.copy() + mask = mesh.prepare_mask(alpha, 0.9, np.eye(2) * 3.2, fine_outline) + assert mask[10, 10] + assert mask[24, 24] == (not fine_outline) + assert mask[40, 25] == (not fine_outline) + assert not mask[0, 60] + np.testing.assert_array_equal(alpha, original) + + +@pytest.mark.parametrize("size,filled", [(2, True), (3, False)]) +def test_micro_hole_threshold_respects_world_transform(size, filled): + alpha = np.ones((40, 40)) + alpha[10 : 10 + size, 10 : 10 + size] = 0 + mask = mesh.prepare_mask(alpha, 0.9, np.eye(2) * 3.2, True) + assert bool(mask[10, 10]) == filled + stretched = mesh.prepare_mask(alpha, 0.9, np.diag((6.4, 3.2)), True) + assert not stretched[10, 10] + + +def test_filling_nested_hole_merges_island_without_overlap(): + alpha = np.zeros((80, 80)) + alpha[5:75, 5:75] = 1 + alpha[20:60, 20:60] = 0 + alpha[30:50, 30:50] = 1 + result = mesh.build_mesh(alpha, 8, 0.9, fine_outline=False) + p, f, _, components = result + assert len(components) == 1 + assert_supported(*result) + edges = np.unique( + np.sort(np.concatenate((f[:, :2], f[:, 1:], f[:, ::2])), axis=1), axis=0 + ) + assert 1 - len(p) + len(edges) - len(f) == 0 + + +def test_tangle_numeric_context_and_failure_recovery(): + from decimal import getcontext + from anyimage.operators.cutout_tool import tangle + + points = [(0.0, 0.0), (20.0, 0.0), (20.0, 20.0), (0.0, 20.0)] + segments = [(0, 1, 1), (1, 2, 2), (2, 3, 3), (3, 0, 4)] + before = getcontext().copy() + epsilon = tangle.EPS + first = tangle.build_mesh(points, segments, [], 10.0) + with pytest.raises(ValueError, match="limit"): + tangle.build_mesh(points, segments, [], 0.01, max_vertices=12) + assert getcontext().prec == before.prec and getcontext().flags == before.flags + assert tangle.EPS == epsilon + assert first == tangle.build_mesh(points, segments, [], 10.0) + + +def test_large_hole_grid_stops_at_vertex_limit_and_recovers(): + from anyimage.operators.cutout_tool import tangle + + points = [(0.0, 0.0), (100.0, 0.0), (100.0, 100.0), (0.0, 100.0)] + segments = [(0, 1, 1), (1, 2, 2), (2, 3, 3), (3, 0, 4)] + holes = [] + for y in range(10, 90, 10): + for x in range(10, 90, 10): + start = len(points) + points.extend([(x, y), (x + 3, y), (x + 3, y + 3), (x, y + 3)]) + segments.extend( + (start + i, start + (i + 1) % 4, len(segments) + i + 1) + for i in range(4) + ) + holes.append((x + 1.5, y + 1.5)) + with pytest.raises(ValueError, match="limit"): + tangle.build_mesh(points, segments, holes, 0.001, max_vertices=300) + p, f, _ = tangle.build_mesh(points, segments, holes, 40.0, max_vertices=5000) + assert len(p) < 5000 and len(f) > 0 + + diff --git a/tests/anyimage/operators/test_cutout_objects.py b/tests/anyimage/operators/test_cutout_objects.py index a0c5c62..006cb57 100644 --- a/tests/anyimage/operators/test_cutout_objects.py +++ b/tests/anyimage/operators/test_cutout_objects.py @@ -39,6 +39,7 @@ def test_cutout_depth_axis_is_adjusted_only_in_modifier(shape, gesture): @ Matrix.Diagonal((1.5, 0.75, 2.0, 1.0)) ) bpy.context.collection.objects.link(source) + bpy.context.view_layer.update() positive_x = _create_cutout_shape( source, @@ -56,11 +57,11 @@ def test_cutout_depth_axis_is_adjusted_only_in_modifier(shape, gesture): for result in (positive_x, negative_z): modifier = result.modifiers[0] mode = object_data.modifier_input_identifier(modifier.node_group, "Mode") - assert modifier[mode] == (1 if gesture == "POLYLINE" else 0) + assert modifier_value(modifier, mode) == (1 if gesture == "POLYLINE" else 0) thickness = object_data.modifier_input_identifier( modifier.node_group, "Thickness", subtype="DISTANCE", ) - assert modifier[thickness] == (0.5 if shape == "SOLID" and gesture == "POLYLINE" else 0.0) + assert modifier_value(modifier, thickness) == (0.5 if shape == "SOLID" and gesture == "POLYLINE" else 0.0) local = np.asarray([tuple(vertex.co) for vertex in result.data.vertices]) assert np.allclose(local[:, 2], 0.0, atol=1e-8) assert np.ptp(local[:, 0]) > 0.0 @@ -70,7 +71,7 @@ def test_cutout_depth_axis_is_adjusted_only_in_modifier(shape, gesture): assert np.allclose(positive_x.matrix_world, source.matrix_world) modifier = positive_x.modifiers[0] axis = object_data.modifier_input_identifier(modifier.node_group, "Depth Axis") - assert modifier[axis] == 0 + assert modifier_value(modifier, axis) == 0 object_data.set_modifier_input(modifier, axis, 1) positive_x.update_tag(refresh={"DATA"}) bpy.context.view_layer.update() @@ -261,9 +262,9 @@ def test_depth_surface_protected_smooth_evaluates(gesture): assert cutout_shape.BALLOON_ATTRIBUTE_NAME in result.data.attributes modifier = result.modifiers[0] assert modifier.node_group.name == cutout_shape.CUTOUT_NODE_GROUP_NAMES["DEPTH_SOLID"] - assert modifier[object_data.modifier_input_identifier(modifier.node_group, "Mode")] == (1 if gesture == "POLYLINE" else 0) - assert modifier[object_data.modifier_input_identifier(modifier.node_group, "Thickness", subtype="DISTANCE")] == 0.5 - assert modifier[object_data.modifier_input_identifier(modifier.node_group, "Thickness", subtype="NONE")] == 1.0 + assert modifier_value(modifier, object_data.modifier_input_identifier(modifier.node_group, "Mode")) == (1 if gesture == "POLYLINE" else 0) + assert modifier_value(modifier, object_data.modifier_input_identifier(modifier.node_group, "Thickness", subtype="DISTANCE")) == 0.5 + assert modifier_value(modifier, object_data.modifier_input_identifier(modifier.node_group, "Thickness", subtype="NONE")) == 1.0 object_data.set_modifier_input( modifier, object_data.modifier_input_identifier(modifier.node_group, "Thickness", subtype="DISTANCE"), @@ -274,14 +275,12 @@ def test_depth_surface_protected_smooth_evaluates(gesture): cutout_geometry.world_plane_metric(source.matrix_world)[0], ) assert np.isclose( - modifier[object_data.modifier_input_identifier(modifier.node_group, "Cleanup Threshold")], + modifier_value(modifier, object_data.modifier_input_identifier(modifier.node_group, "Cleanup Threshold")), expected_cleanup, ) - assert np.isclose(modifier[object_data.modifier_input_identifier(modifier.node_group, "Smooth Weight")], 1.0) + assert np.isclose(modifier_value(modifier, object_data.modifier_input_identifier(modifier.node_group, "Smooth Weight")), 1.0) assert ( - modifier[ - object_data.modifier_input_identifier(modifier.node_group, "Smooth") - ] + modifier_value(modifier, object_data.modifier_input_identifier(modifier.node_group, "Smooth")) == 0 ) for name, value in (("Mode", 1), ("Front Inflation", 0.0)): @@ -290,7 +289,7 @@ def test_depth_surface_protected_smooth_evaluates(gesture): modifier.node_group, "Reference Depth", ) - assert modifier[reference_id] == 8.0 + assert modifier_value(modifier, reference_id) == 8.0 object_data.set_modifier_input( modifier, object_data.modifier_input_identifier(modifier.node_group, "Smooth Weight"), @@ -639,7 +638,7 @@ def test_depth_balloon_uses_metric_plane_inputs(): modifier.node_group, "Depth Image", ) - assert isinstance(modifier[depth_image_id], bpy.types.Image) + assert isinstance(modifier_value(modifier, depth_image_id), bpy.types.Image) expected_center = np.zeros(2, dtype=np.float64) crop_center = Vector( ( @@ -659,18 +658,18 @@ def test_depth_balloon_uses_metric_plane_inputs(): modifier.node_group, "Reference Depth", ) - assert not np.isclose(modifier[depth_reference_id], 8.0) + assert not np.isclose(modifier_value(modifier, depth_reference_id), 8.0) depth_direction_id = object_data.modifier_input_identifier( modifier.node_group, "Depth Direction", ) - direction = np.asarray(modifier[depth_direction_id], dtype=np.float64) + direction = np.asarray(modifier_value(modifier, depth_direction_id), dtype=np.float64) assert not np.allclose(direction, (0.0, 0.0, 1.0)) uniform_scale_id = object_data.modifier_input_identifier( modifier.node_group, "Uniform Scale", ) - assert modifier[uniform_scale_id] > 0.0 + assert modifier_value(modifier, uniform_scale_id) > 0.0 thickness_id = object_data.modifier_input_identifier( modifier.node_group, "Thickness", @@ -679,7 +678,7 @@ def test_depth_balloon_uses_metric_plane_inputs(): modifier.node_group, "Smooth Weight", ) - assert modifier[smooth_id] == 1.0 + assert modifier_value(modifier, smooth_id) == 1.0 object_data.set_modifier_input(modifier, smooth_id, 0.0) result.update_tag(refresh={"DATA"}) bpy.context.view_layer.update() @@ -707,7 +706,7 @@ def test_depth_balloon_uses_metric_plane_inputs(): object_data.set_modifier_input( modifier, depth_reference_id, - modifier[depth_reference_id] + 0.2, + modifier_value(modifier, depth_reference_id) + 0.2, ) result.update_tag(refresh={"DATA"}) bpy.context.view_layer.update() @@ -729,7 +728,7 @@ def test_depth_balloon_uses_metric_plane_inputs(): object_data.set_modifier_input( modifier, depth_reference_id, - modifier[depth_reference_id] + 0.15, + modifier_value(modifier, depth_reference_id) + 0.15, ) object_data.set_modifier_input(modifier, depth_direction_id, (0.2, 0.1, 1.0)) result.update_tag(refresh={"DATA"}) @@ -889,17 +888,15 @@ def test_local_cutout_builds_from_blender_image_pixels_without_a_job(): cutout_geometry.world_plane_metric(source.matrix_world)[0], ) assert np.isclose( - modifier[object_data.modifier_input_identifier(modifier.node_group, "Cleanup Threshold")], + modifier_value(modifier, object_data.modifier_input_identifier(modifier.node_group, "Cleanup Threshold")), expected_cleanup, ) assert ( - modifier[object_data.modifier_input_identifier(modifier.node_group, "Smooth Weight")] + modifier_value(modifier, object_data.modifier_input_identifier(modifier.node_group, "Smooth Weight")) == 1.0 ) assert ( - modifier[ - object_data.modifier_input_identifier(modifier.node_group, "Smooth") - ] + modifier_value(modifier, object_data.modifier_input_identifier(modifier.node_group, "Smooth")) == 0 ) assert len(result.data.vertices) > 4 @@ -931,3 +928,15 @@ def test_local_cutout_builds_from_blender_image_pixels_without_a_job(): assert np.allclose( copied[:, :, 3], np.flipud(rgba[top:bottom, left:right, 3]) ) + + +def modifier_value(modifier, identifier): + """Read modifier values through the runtime's socket API.""" + inputs = object_data._modifier_input_slots(modifier) + if inputs is None: + return modifier[identifier] + slot = getattr(inputs, identifier) + prop = slot.bl_rna.properties["value"] + if prop.type == "ENUM": + return list(prop.enum_items.keys()).index(slot.value) + return slot.value diff --git a/tests/tools/nodes/test_boundary_depth.py b/tests/tools/nodes/test_boundary_depth.py index dd66288..577d465 100644 --- a/tests/tools/nodes/test_boundary_depth.py +++ b/tests/tools/nodes/test_boundary_depth.py @@ -73,12 +73,12 @@ def test_camera_blur_preserves_rays_interior_protection_and_disconnected_surface def test_cutout_fine_outline_blurs_depth_without_moving_projection(): from tests.support.depth_surface import surface - from anyimage.operators.cutout_tool.fine_mesh import build_fine_mesh + from anyimage.operators.cutout_tool.mesh import build_mesh obj, set_value = surface(step=0, resolution=8) yy, xx = np.mgrid[:32, :64] alpha = ((xx - 32)**2 / 27**2 + (yy - 16)**2 / 13**2 < 1).astype(float) - points, faces, _, _ = build_fine_mesh(alpha, 3, .9) + points, faces, _, _ = build_mesh(alpha, 3, .9) obj.data.clear_geometry() obj.data.from_pydata([(x / 32, y / 32, 0) for x, y in points], [], faces.tolist()) uv = obj.data.uv_layers.get("UVMap") or obj.data.uv_layers.new(name="UVMap") diff --git a/tests/tools/test_pack.py b/tests/tools/test_pack.py index 58a691f..6c25215 100644 --- a/tests/tools/test_pack.py +++ b/tests/tools/test_pack.py @@ -94,7 +94,9 @@ class ExtensionPackagingTest(unittest.TestCase): "operators/cutout_tool/operators.py", "operators/cutout_tool/interaction.py", "operators/cutout_tool/geometry.py", - "operators/cutout_tool/fine_mesh.py", + "operators/cutout_tool/mesh.py", + "operators/cutout_tool/tangle.py", + "licenses/tangle.txt", "operators/cutout_tool/object.py", "operators/cutout_tool/shape.py", "operators/convert_to_plane/__init__.py", diff --git a/tools/nodes/groups/image_depth_cutout.py b/tools/nodes/groups/image_depth_cutout.py index 46f7947..7f1d5bd 100644 --- a/tools/nodes/groups/image_depth_cutout.py +++ b/tools/nodes/groups/image_depth_cutout.py @@ -278,6 +278,21 @@ def _build_depth_surface(group, geometry, controls): group, geometry, camera, scale, controls.outputs["Reference Depth"], controls.outputs["Depth Scale"], split, ) projection = smooth_cut_boundary(group, projection) + # Clamped depth pixels can collapse adjacent outline samples after projection. + edge = group.nodes.new("GeometryNodeInputMeshEdgeVertices") + length = group.nodes.new("ShaderNodeVectorMath") + length.operation = "DISTANCE" + links.new(edge.outputs["Position 1"], length.inputs[0]) + links.new(edge.outputs["Position 2"], length.inputs[1]) + typical = group.nodes.new("GeometryNodeAttributeStatistic") + typical.data_type, typical.domain = "FLOAT", "EDGE" + links.new(projection, typical.inputs["Geometry"]) + links.new(length.outputs["Value"], typical.inputs["Attribute"]) + merge = group.nodes.new("GeometryNodeMergeByDistance") + merge.mode = "CONNECTED" + links.new(projection, merge.inputs["Geometry"]) + links.new(_math(group, "MULTIPLY", typical.outputs["Mean"], 1e-6), merge.inputs["Distance"]) + projection = merge.outputs["Geometry"] projection = taper_split_profile(group, original, projection) # The mode field is shared by the geometry and thickness choices. mode = group.nodes.new("GeometryNodeMenuSwitch")