feat(viz): contour plot via marching squares — layout puro + render

contour_compute implementa marching squares clasico (16 casos, casos ambiguos 5 y 10 partidos en 2 segmentos). Para cada level devuelve un ContourLine{pts, level} con segmentos en coords [0..nx-1]x[0..ny-1]. Verificado con gaussiana 32x32 + 4 niveles: todos los endpoints aparecen >=2 veces (curvas cerradas, ningun endpoint huerfano).
2026-04-25 21:52:48 +02:00
parent 75d4334e8c
commit cda557286e
3 changed files with 225 additions and 0 deletions
@@ -0,0 +1,132 @@
+#include "viz/contour.h"
+
+#include <cmath>
+
+namespace {
+
+// Devuelve la coordenada interpolada [0..1] entre v1 y v2 donde se cruza level.
+float interp(float v1, float v2, float level) {
+    float d = v2 - v1;
+    if (std::fabs(d) < 1e-9f) return 0.5f;
+    float t = (level - v1) / d;
+    if (t < 0.0f) t = 0.0f;
+    if (t > 1.0f) t = 1.0f;
+    return t;
+}
+
+inline float at(const float* g, int nx, int x, int y) {
+    return g[y * nx + x];
+}
+
+} // namespace
+
+std::vector<ContourLine> contour_compute(const float* grid,
+                                         int          nx,
+                                         int          ny,
+                                         const float* levels,
+                                         int          n_levels) {
+    std::vector<ContourLine> out;
+    if (!grid || nx < 2 || ny < 2 || !levels || n_levels <= 0) return out;
+    out.resize(n_levels);
+    for (int li = 0; li < n_levels; li++) out[li].level = levels[li];
+
+    for (int li = 0; li < n_levels; li++) {
+        float L = levels[li];
+        auto& line = out[li];
+        for (int y = 0; y < ny - 1; y++) {
+            for (int x = 0; x < nx - 1; x++) {
+                float v00 = at(grid, nx, x,     y);
+                float v10 = at(grid, nx, x + 1, y);
+                float v11 = at(grid, nx, x + 1, y + 1);
+                float v01 = at(grid, nx, x,     y + 1);
+
+                int code = 0;
+                if (v00 >= L) code |= 1;
+                if (v10 >= L) code |= 2;
+                if (v11 >= L) code |= 4;
+                if (v01 >= L) code |= 8;
+                if (code == 0 || code == 15) continue;
+
+                // Puntos en cada arista (top=between v00 y v10, etc.)
+                ImVec2 pT(x + interp(v00, v10, L), (float)y);
+                ImVec2 pR((float)(x + 1),          y + interp(v10, v11, L));
+                ImVec2 pB(x + interp(v01, v11, L), (float)(y + 1));
+                ImVec2 pL((float)x,                y + interp(v00, v01, L));
+
+                auto seg = [&](ImVec2 a, ImVec2 b) {
+                    line.pts.push_back(a);
+                    line.pts.push_back(b);
+                };
+
+                switch (code) {
+                    case 1:  case 14: seg(pL, pT); break;
+                    case 2:  case 13: seg(pT, pR); break;
+                    case 4:  case 11: seg(pR, pB); break;
+                    case 8:  case 7:  seg(pB, pL); break;
+                    case 3:  case 12: seg(pL, pR); break;
+                    case 6:  case 9:  seg(pT, pB); break;
+                    // ambiguos: 5 y 10 — partir en dos segmentos.
+                    case 5:           seg(pL, pT); seg(pR, pB); break;
+                    case 10:          seg(pT, pR); seg(pB, pL); break;
+                    default: break;
+                }
+            }
+        }
+    }
+    return out;
+}
+
+void contour(const char*  id,
+             const float* grid,
+             int          nx,
+             int          ny,
+             const float* levels,
+             int          n_levels,
+             ImVec2       size) {
+    ImGui::PushID(id);
+
+    ImVec2 avail = ImGui::GetContentRegionAvail();
+    float  W     = (size.x > 0.0f) ? size.x : avail.x;
+    float  H     = (size.y > 0.0f) ? size.y : 200.0f;
+
+    ImVec2 origin = ImGui::GetCursorScreenPos();
+    ImGui::Dummy(ImVec2(W, H));
+
+    if (!grid || nx < 2 || ny < 2 || n_levels <= 0) {
+        ImGui::PopID();
+        return;
+    }
+
+    auto lines = contour_compute(grid, nx, ny, levels, n_levels);
+
+    ImDrawList* dl = ImGui::GetWindowDrawList();
+
+    // Borde
+    dl->AddRect(origin, ImVec2(origin.x + W, origin.y + H),
+                IM_COL32(80, 80, 90, 200), 0.0f, 0, 1.0f);
+
+    float sx = W / (float)(nx - 1);
+    float sy = H / (float)(ny - 1);
+
+    // Color por nivel: gradiente azul -> amarillo
+    auto level_color = [&](int idx) {
+        float t = (n_levels > 1) ? (float)idx / (float)(n_levels - 1) : 0.5f;
+        // (76,140,230) -> (250,200,90)
+        int r = (int)(76  + (250 - 76)  * t);
+        int g = (int)(140 + (200 - 140) * t);
+        int b = (int)(230 + (90  - 230) * t);
+        return IM_COL32(r, g, b, 230);
+    };
+
+    for (int li = 0; li < (int)lines.size(); li++) {
+        ImU32 col = level_color(li);
+        const auto& seg = lines[li].pts;
+        for (size_t k = 0; k + 1 < seg.size(); k += 2) {
+            ImVec2 a(origin.x + seg[k].x     * sx, origin.y + seg[k].y     * sy);
+            ImVec2 b(origin.x + seg[k + 1].x * sx, origin.y + seg[k + 1].y * sy);
+            dl->AddLine(a, b, col, 1.5f);
+        }
+    }
+
+    ImGui::PopID();
+}