/* Compile (Linux):
 *   gcc loft.c -o loft -DUSE_SDL2 `sdl2-config --cflags --libs` -lm
 * or without SDL2:
 *   gcc loft.c -o loft -lm
 *
 * Run and redirect to OpenSCAD:
 *   ./loft > loft.scad
 *   (you may need "sudo ./loft > loft.scad" if running SDL2 programs requires root access on your system)
 *
 * The generated OpenSCAD file, loft.scad, describes the lofted solid as a polyhedron()
 * built from multiple interpolated cross-sections (“slices”). [web:70][web:72][web:77]
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#ifdef USE_SDL2
#include <SDL2/SDL.h>
#endif

/* ----------------- Basic math types ----------------- */

typedef struct { double x, y; } Vec2;
typedef struct { double x, y, z; } Vec3;

/* Plane frame: P = origin + u*x_axis + v*y_axis */
typedef struct {
    Vec3 origin;
    Vec3 x_axis;
    Vec3 y_axis;
    Vec3 normal;
} PlaneFrame;

/* 4x4 column-major matrix */
typedef struct {
    double m[16];
} Mat4;

/* ----------------- 2D/3D helpers ----------------- */

static Vec2 v2_add(Vec2 a, Vec2 b){ Vec2 r = {a.x + b.x, a.y + b.y}; return r; }
static Vec2 v2_sub(Vec2 a, Vec2 b){ Vec2 r = {a.x - b.x, a.y - b.y}; return r; }
static Vec2 v2_scale(Vec2 a, double s){ Vec2 r = {a.x * s, a.y * s}; return r; }
static double v2_dot(Vec2 a, Vec2 b){ return a.x*b.x + a.y*b.y; }
static double v2_len(Vec2 a){ return sqrt(v2_dot(a,a)); }

static Vec3 v3_add(Vec3 a, Vec3 b){ Vec3 r = {a.x + b.x, a.y + b.y, a.z + b.z}; return r; }
static Vec3 v3_sub(Vec3 a, Vec3 b){ Vec3 r = {a.x - b.x, a.y - b.y, a.z - b.z}; return r; }
static Vec3 v3_scale(Vec3 a, double s){ Vec3 r = {a.x * s, a.y * s, a.z * s}; return r; }
static double v3_dot(Vec3 a, Vec3 b){ return a.x*b.x + a.y*b.y + a.z*b.z; }
static Vec3 v3_cross(Vec3 a, Vec3 b){
    Vec3 r = {a.y*b.z - a.z*b.y,
              a.z*b.x - a.x*b.z,
              a.x*b.y - a.y*b.x};
    return r;
}
static double v3_len(Vec3 a){ return sqrt(v3_dot(a,a)); }
static Vec3 v3_norm(Vec3 a){
    double L = v3_len(a);
    if (L == 0.0) return a;
    return v3_scale(a, 1.0 / L);
}

/* ----------------- 4x4 matrices ----------------- */

static Mat4 mat4_identity(void)
{
    Mat4 M;
    for (int i = 0; i < 16; ++i) M.m[i] = 0.0;
    M.m[ 0] = 1.0;
    M.m[ 5] = 1.0;
    M.m[10] = 1.0;
    M.m[15] = 1.0;
    return M;
}

static Mat4 mat4_mul(Mat4 A, Mat4 B)
{
    Mat4 R;
    for (int c = 0; c < 4; ++c) {
        for (int r = 0; r < 4; ++r) {
            double v = 0.0;
            for (int k = 0; k < 4; ++k) {
                v += A.m[k*4 + r] * B.m[c*4 + k];
            }
            R.m[c*4 + r] = v;
        }
    }
    return R;
}

static Vec3 mat4_mul_point(Mat4 M, Vec3 p)
{
    double x = p.x, y = p.y, z = p.z;
    double rx = M.m[0]*x + M.m[4]*y + M.m[8]*z  + M.m[12];
    double ry = M.m[1]*x + M.m[5]*y + M.m[9]*z  + M.m[13];
    double rz = M.m[2]*x + M.m[6]*y + M.m[10]*z + M.m[14];
    double rw = M.m[3]*x + M.m[7]*y + M.m[11]*z + M.m[15];
    if (rw != 0.0) {
        rx /= rw; ry /= rw; rz /= rw;
    }
    Vec3 r = {rx, ry, rz};
    return r;
}

/* Build transform from 2D plane coords (u,v) to 3D space. */
static Mat4 plane_to_world(const PlaneFrame *pf)
{
    Mat4 M = mat4_identity();

    M.m[0]  = pf->x_axis.x;
    M.m[1]  = pf->x_axis.y;
    M.m[2]  = pf->x_axis.z;

    M.m[4]  = pf->y_axis.x;
    M.m[5]  = pf->y_axis.y;
    M.m[6]  = pf->y_axis.z;

    M.m[8]  = pf->normal.x;
    M.m[9]  = pf->normal.y;
    M.m[10] = pf->normal.z;

    M.m[12] = pf->origin.x;
    M.m[13] = pf->origin.y;
    M.m[14] = pf->origin.z;

    return M;
}

/* ----------------- 2D Bézier path ----------------- */

typedef struct {
    Vec2 p0, p1, p2, p3;
} CubicBezier2D;

typedef struct {
    int           count;
    CubicBezier2D *seg;
} BezierPath2D;

typedef struct {
    int   count;
    Vec2 *pt;
} Polyline2D;

/* Cubic Bézier 2D eval/flatness/subdivide */

static Vec2 cubic2_eval(const CubicBezier2D *c, double t)
{
    double u  = 1.0 - t;
    double tt = t * t;
    double uu = u * u;
    double uuu = uu * u;
    double ttt = tt * t;

    Vec2 p;
    p.x = uuu * c->p0.x;
    p.y = uuu * c->p0.y;

    p.x += 3.0 * uu * t * c->p1.x;
    p.y += 3.0 * uu * t * c->p1.y;

    p.x += 3.0 * u * tt * c->p2.x;
    p.y += 3.0 * u * tt * c->p2.y;

    p.x += ttt * c->p3.x;
    p.y += ttt * c->p3.y;

    return p;
}

static double cubic2_flatness_sq(const CubicBezier2D *c)
{
    Vec2 a = c->p0;
    Vec2 b = c->p3;
    Vec2 ab = v2_sub(b, a);

    double ab_len = v2_len(ab);
    if (ab_len == 0.0) {
        double d1 = v2_len(v2_sub(c->p1, a));
        double d2 = v2_len(v2_sub(c->p2, a));
        double d = d1 > d2 ? d1 : d2;
        return d * d;
    }

    double abx = ab.x;
    double aby = ab.y;

    double dist_sq(Vec2 p){
        Vec2 ap = v2_sub(p, a);
        double cross = ap.x * aby - ap.y * abx;
        double d = cross / ab_len;
        return d * d;
    }

    double d1 = dist_sq(c->p1);
    double d2 = dist_sq(c->p2);
    return d1 > d2 ? d1 : d2;
}

static void cubic2_subdivide(const CubicBezier2D *c,
                             CubicBezier2D *left,
                             CubicBezier2D *right)
{
    Vec2 p01  = v2_scale(v2_add(c->p0, c->p1), 0.5);
    Vec2 p12  = v2_scale(v2_add(c->p1, c->p2), 0.5);
    Vec2 p23  = v2_scale(v2_add(c->p2, c->p3), 0.5);

    Vec2 p012 = v2_scale(v2_add(p01, p12), 0.5);
    Vec2 p123 = v2_scale(v2_add(p12, p23), 0.5);

    Vec2 p0123 = v2_scale(v2_add(p012, p123), 0.5);

    left->p0 = c->p0;
    left->p1 = p01;
    left->p2 = p012;
    left->p3 = p0123;

    right->p0 = p0123;
    right->p1 = p123;
    right->p2 = p23;
    right->p3 = c->p3;
}

/* ----------------- Vec2 dynamic buffer and flattening ----------------- */

typedef struct {
    Vec2 *data;
    int   count;
    int   cap;
} Vec2Buffer;

static void buf2_init(Vec2Buffer *b){ b->data = NULL; b->count = b->cap = 0; }
static void buf2_free(Vec2Buffer *b){ free(b->data); b->data = NULL; b->count = b->cap = 0; }

static int buf2_push(Vec2Buffer *b, Vec2 p)
{
    if (b->count >= b->cap) {
        int ncap = b->cap ? b->cap * 2 : 64;
        Vec2 *nd = (Vec2*)realloc(b->data, (size_t)ncap * sizeof(Vec2));
        if (!nd) return 0;
        b->data = nd;
        b->cap  = ncap;
    }
    b->data[b->count++] = p;
    return 1;
}

static int flatten_cubic2_rec(const CubicBezier2D *c,
                              double tol_sq,
                              Vec2Buffer *out)
{
    if (cubic2_flatness_sq(c) <= tol_sq) {
        return buf2_push(out, c->p3);
    } else {
        CubicBezier2D l, r;
        cubic2_subdivide(c, &l, &r);
        if (!flatten_cubic2_rec(&l, tol_sq, out)) return 0;
        if (!flatten_cubic2_rec(&r, tol_sq, out)) return 0;
        return 1;
    }
}

static int flatten_path2(const BezierPath2D *path,
                         double tol,
                         Polyline2D *poly_out)
{
    Vec2Buffer buf;
    buf2_init(&buf);
    double tol_sq = tol * tol;

    if (path->count <= 0) {
        poly_out->count = 0;
        poly_out->pt = NULL;
        return 1;
    }

    if (!buf2_push(&buf, path->seg[0].p0)) {
        buf2_free(&buf);
        return 0;
    }

    for (int i = 0; i < path->count; ++i) {
        CubicBezier2D c = path->seg[i];
        if (!flatten_cubic2_rec(&c, tol_sq, &buf)) {
            buf2_free(&buf);
            return 0;
        }
    }

    if (buf.count > 0) {
        Vec2 first = buf.data[0];
        Vec2 last  = buf.data[buf.count - 1];
        if (fabs(first.x - last.x) > 1e-9 || fabs(first.y - last.y) > 1e-9) {
            if (!buf2_push(&buf, first)) {
                buf2_free(&buf);
                return 0;
            }
        }
    }

    poly_out->count = buf.count;
    poly_out->pt    = buf.data;
    return 1;
}

/* ----------------- Polyline arclength/resampling ----------------- */

static double *poly2_arclen(const Polyline2D *pl)
{
    if (pl->count <= 0) return NULL;
    double *s = (double*)malloc((size_t)pl->count * sizeof(double));
    if (!s) return NULL;
    s[0] = 0.0;
    for (int i = 1; i < pl->count; ++i) {
        Vec2 d = v2_sub(pl->pt[i], pl->pt[i-1]);
        s[i] = s[i-1] + v2_len(d);
    }
    return s;
}

static Vec2 poly2_sample_at(const Polyline2D *pl,
                            const double *s,
                            double d)
{
    int n = pl->count;
    double total = s[n-1];
    if (total <= 0.0) return pl->pt[0];

    d = fmod(d, total);
    if (d < 0.0) d += total;

    int i = 0;
    while (i+1 < n && s[i+1] < d) ++i;

    int j = (i+1 < n) ? i+1 : i;
    double si = s[i];
    double sj = s[j];
    double t;
    if (sj <= si) t = 0.0;
    else t = (d - si) / (sj - si);

    Vec2 a = pl->pt[i];
    Vec2 b = pl->pt[j];
    return v2_add(a, v2_scale(v2_sub(b, a), t));
}

static int resample_polyline2(const Polyline2D *pl,
                              int N,
                              Polyline2D *out)
{
    if (N <= 0 || pl->count <= 1) {
        out->count = 0;
        out->pt = NULL;
        return 1;
    }

    double *s = poly2_arclen(pl);
    if (!s) return 0;

    Vec2 *pts = (Vec2*)malloc((size_t)N * sizeof(Vec2));
    if (!pts) { free(s); return 0; }

    double total = s[pl->count - 1];
    for (int i = 0; i < N; ++i) {
        double d = total * ((double)i / (double)N);
        pts[i] = poly2_sample_at(pl, s, d);
    }

    free(s);
    out->count = N;
    out->pt    = pts;
    return 1;
}

/* ----------------- Alignment and mapping ----------------- */

static double cyclic_cost2(const Vec2 *A, const Vec2 *B, int n, int shift)
{
    double c = 0.0;
    for (int i = 0; i < n; ++i) {
        int j = i + shift;
        if (j >= n) j -= n;
        double dx = A[i].x - B[j].x;
        double dy = A[i].y - B[j].y;
        c += dx*dx + dy*dy;
    }
    return c;
}

static double cyclic_cost2_rev(const Vec2 *A, const Vec2 *B,
                               int n, int shift)
{
    double c = 0.0;
    for (int i = 0; i < n; ++i) {
        int bi = i + shift;
        if (bi >= n) bi -= n;
        int j = n - 1 - bi;
        double dx = A[i].x - B[j].x;
        double dy = A[i].y - B[j].y;
        c += dx*dx + dy*dy;
    }
    return c;
}

static void align_closed2(const Vec2 *A, const Vec2 *B, int n,
                          int *out_shift, int *out_reverse)
{
    double best_cost = 1e300;
    int best_shift = 0;
    int best_rev = 0;

    for (int shift = 0; shift < n; ++shift) {
        double c0 = cyclic_cost2(A, B, n, shift);
        if (c0 < best_cost) {
            best_cost = c0;
            best_shift = shift;
            best_rev   = 0;
        }
        double c1 = cyclic_cost2_rev(A, B, n, shift);
        if (c1 < best_cost) {
            best_cost = c1;
            best_shift = shift;
            best_rev   = 1;
        }
    }
    *out_shift   = best_shift;
    *out_reverse = best_rev;
}

static Vec2 *build_aligned2(const Vec2 *B, int n,
                            int shift, int reverse)
{
    Vec2 *B2 = (Vec2*)malloc((size_t)n * sizeof(Vec2));
    if (!B2) return NULL;

    if (!reverse) {
        for (int i = 0; i < n; ++i) {
            int j = i + shift;
            if (j >= n) j -= n;
            B2[i] = B[j];
        }
    } else {
        for (int i = 0; i < n; ++i) {
            int bi = i + shift;
            if (bi >= n) bi -= n;
            int j = n - 1 - bi;
            B2[i] = B[j];
        }
    }
    return B2;
}

typedef struct {
    int   N;
    Vec2 *A;
    Vec2 *B;
} ShapeMapping2D;

static void free_shape_mapping2(ShapeMapping2D *m)
{
    if (!m) return;
    free(m->A);
    free(m->B);
    m->A = m->B = NULL;
    m->N = 0;
}

static int choose_sample_count2(int ctrl_src, int ctrl_dst, double tol)
{
    (void)tol;
    int base = (ctrl_src > ctrl_dst) ? ctrl_src : ctrl_dst;
    int k = 4;
    int N = base * k;
    if (N < 32) N = 32;
    return N;
}

static int build_shape_mapping2(const BezierPath2D *src,
                                const BezierPath2D *dst,
                                double flat_tol,
                                ShapeMapping2D *out)
{
    Polyline2D pl_src, pl_dst;

    if (!flatten_path2(src, flat_tol, &pl_src)) return 0;
    if (!flatten_path2(dst, flat_tol, &pl_dst)) {
        free(pl_src.pt);
        return 0;
    }

    int ctrl_src = src->count * 4;
    int ctrl_dst = dst->count * 4;

    int N = choose_sample_count2(ctrl_src, ctrl_dst, flat_tol);

    Polyline2D rs_src, rs_dst;
    if (!resample_polyline2(&pl_src, N, &rs_src)) {
        free(pl_src.pt);
        free(pl_dst.pt);
        return 0;
    }
    if (!resample_polyline2(&pl_dst, N, &rs_dst)) {
        free(pl_src.pt);
        free(pl_dst.pt);
        free(rs_src.pt);
        return 0;
    }

    free(pl_src.pt);
    free(pl_dst.pt);

    int shift, rev;
    align_closed2(rs_src.pt, rs_dst.pt, N, &shift, &rev);

    Vec2 *A = rs_src.pt;
    Vec2 *B = build_aligned2(rs_dst.pt, N, shift, rev);
    free(rs_dst.pt);
    if (!B) {
        free(A);
        return 0;
    }

    out->N = N;
    out->A = A;
    out->B = B;
    return 1;
}

static Polyline2D tween_polyline2(const ShapeMapping2D *m, double t)
{
    Polyline2D pl;
    pl.count = m->N + 1;
    pl.pt = (Vec2*)malloc((size_t)pl.count * sizeof(Vec2));
    if (!pl.pt) {
        pl.count = 0;
        return pl;
    }
    double omt = 1.0 - t;
    for (int i = 0; i < m->N; ++i) {
        Vec2 a = m->A[i];
        Vec2 b = m->B[i];
        Vec2 p;
        p.x = omt * a.x + t * b.x;
        p.y = omt * a.y + t * b.y;
        pl.pt[i] = p;
    }
    pl.pt[m->N] = pl.pt[0];
    return pl;
}

/* ----------------- Loft: lifting 2D slices to 3D ----------------- */

typedef struct {
    int   count;
    Vec3 *pt;
} Polyline3D;

static void free_polyline2(Polyline2D *pl)
{
    if (!pl) return;
    free(pl->pt);
    pl->pt = NULL;
    pl->count = 0;
}

static void free_polyline3(Polyline3D *pl)
{
    if (!pl) return;
    free(pl->pt);
    pl->pt = NULL;
    pl->count = 0;
}

static Polyline3D lift_polyline_to_plane(const Polyline2D *pl,
                                         const PlaneFrame *plane)
{
    Polyline3D out;
    out.count = pl->count;
    out.pt = (Vec3*)malloc((size_t)out.count * sizeof(Vec3));
    if (!out.pt) {
        out.count = 0;
        return out;
    }

    Mat4 T = plane_to_world(plane);
    for (int i = 0; i < pl->count; ++i) {
        Vec2 q = pl->pt[i];
        Vec3 p = {q.x, q.y, 0.0};
        out.pt[i] = mat4_mul_point(T, p);
    }
    return out;
}

/* ----------------- Camera/projection (for SDL view) ----------------- */

typedef struct {
    Vec3 pos;
    Vec3 forward;
    Vec3 up;
} Camera;

static Mat4 camera_view_matrix(const Camera *cam)
{
    Vec3 f = v3_norm(cam->forward);
    Vec3 r = v3_norm(v3_cross(f, cam->up));
    Vec3 u = v3_cross(r, f);

    Mat4 V = mat4_identity();

    V.m[0] = r.x; V.m[4] = r.y; V.m[8]  = r.z;
    V.m[1] = u.x; V.m[5] = u.y; V.m[9]  = u.z;
    V.m[2] = -f.x;V.m[6] = -f.y;V.m[10] = -f.z;

    Vec3 t;
    t.x = -v3_dot(r, cam->pos);
    t.y = -v3_dot(u, cam->pos);
    t.z =  v3_dot(f, cam->pos);

    V.m[12] = t.x;
    V.m[13] = t.y;
    V.m[14] = t.z;

    return V;
}

static Mat4 perspective_matrix(double fov_y_deg,
                               double aspect,
                               double z_near,
                               double z_far)
{
    double fov_rad = fov_y_deg * (M_PI / 180.0);
    double f = 1.0 / tan(fov_rad * 0.5);
    Mat4 P;
    for (int i = 0; i < 16; ++i) P.m[i] = 0.0;

    P.m[0]  = f / aspect;
    P.m[5]  = f;
    P.m[10] = (z_far + z_near) / (z_near - z_far);
    P.m[11] = -1.0;
    P.m[14] = (2.0 * z_far * z_near) / (z_near - z_far);
    return P;
}

static int project_point(const Mat4 *VP, Vec3 p,
                         double *x_ndc, double *y_ndc, double *z_ndc)
{
    Vec3 q = mat4_mul_point(*VP, p);
    if (q.z <= 0.0) return 0;
    *x_ndc = q.x;
    *y_ndc = q.y;
    *z_ndc = q.z;
    (void)z_ndc;
    return 1;
}

#ifdef USE_SDL2
static void draw_polyline3d(SDL_Renderer *ren,
                            const Mat4 *VP,
                            const Polyline3D *pl,
                            int w, int h)
{
    if (pl->count < 2) return;

    for (int i = 0; i < pl->count - 1; ++i) {
        Vec3 a = pl->pt[i];
        Vec3 b = pl->pt[i+1];
        double ax, ay, az;
        double bx, by, bz;
        if (!project_point(VP, a, &ax, &ay, &az)) continue;
        if (!project_point(VP, b, &bx, &by, &bz)) continue;
        int sx1 = (int)((ax * 0.5 + 0.5) * (double)w);
        int sy1 = (int)((-ay * 0.5 + 0.5) * (double)h);
        int sx2 = (int)((bx * 0.5 + 0.5) * (double)w);
        int sy2 = (int)((-by * 0.5 + 0.5) * (double)h);
        SDL_RenderDrawLine(ren, sx1, sy1, sx2, sy2);
    }
}
#endif

/* ----------------- Shape construction: square & circle ----------------- */

static BezierPath2D make_square2(void)
{
    BezierPath2D path;
    path.count = 4;
    path.seg = (CubicBezier2D*)malloc(4 * sizeof(CubicBezier2D));
    double s = 1.0;

    path.seg[0].p0.x = -s; path.seg[0].p0.y = -s;
    path.seg[0].p1.x = -s; path.seg[0].p1.y = -s;
    path.seg[0].p2.x =  s; path.seg[0].p2.y = -s;
    path.seg[0].p3.x =  s; path.seg[0].p3.y = -s;

    path.seg[1].p0.x =  s; path.seg[1].p0.y = -s;
    path.seg[1].p1.x =  s; path.seg[1].p1.y = -s;
    path.seg[1].p2.x =  s; path.seg[1].p2.y =  s;
    path.seg[1].p3.x =  s; path.seg[1].p3.y =  s;

    path.seg[2].p0.x =  s; path.seg[2].p0.y =  s;
    path.seg[2].p1.x =  s; path.seg[2].p1.y =  s;
    path.seg[2].p2.x = -s; path.seg[2].p2.y =  s;
    path.seg[2].p3.x = -s; path.seg[2].p3.y =  s;

    path.seg[3].p0.x = -s; path.seg[3].p0.y =  s;
    path.seg[3].p1.x = -s; path.seg[3].p1.y =  s;
    path.seg[3].p2.x = -s; path.seg[3].p2.y = -s;
    path.seg[3].p3.x = -s; path.seg[3].p3.y = -s;

    return path;
}

/* Circle via 4 cubics. [web:36][web:39] */
static BezierPath2D make_circle2(void)
{
    BezierPath2D path;
    path.count = 4;
    path.seg = (CubicBezier2D*)malloc(4 * sizeof(CubicBezier2D));

    double r = 1.0;
    double k = 4.0 * (sqrt(2.0) - 1.0) / 3.0;

    path.seg[0].p0.x =  r; path.seg[0].p0.y = 0.0;
    path.seg[0].p1.x =  r; path.seg[0].p1.y =  r * k;
    path.seg[0].p2.x =  r * k; path.seg[0].p2.y =  r;
    path.seg[0].p3.x =  0.0;   path.seg[0].p3.y =  r;

    path.seg[1].p0 = path.seg[0].p3;
    path.seg[1].p1.x = -r * k; path.seg[1].p1.y =  r;
    path.seg[1].p2.x = -r;     path.seg[1].p2.y =  r * k;
    path.seg[1].p3.x = -r;     path.seg[1].p3.y =  0.0;

    path.seg[2].p0 = path.seg[1].p3;
    path.seg[2].p1.x = -r;     path.seg[2].p1.y = -r * k;
    path.seg[2].p2.x = -r * k; path.seg[2].p2.y = -r;
    path.seg[2].p3.x =  0.0;   path.seg[2].p3.y = -r;

    path.seg[3].p0 = path.seg[2].p3;
    path.seg[3].p1.x =  r * k; path.seg[3].p1.y = -r;
    path.seg[3].p2.x =  r;     path.seg[3].p2.y = -r * k;
    path.seg[3].p3.x =  r;     path.seg[3].p3.y =  0.0;

    return path;
}

/* ----------------- Plane frames ----------------- */

static PlaneFrame make_bottom_plane(void)
{
    PlaneFrame pf;
    pf.origin.x = 0.0; pf.origin.y = 0.0; pf.origin.z = 0.0;
    pf.x_axis.x = 1.0; pf.x_axis.y = 0.0; pf.x_axis.z = 0.0;
    pf.y_axis.x = 0.0; pf.y_axis.y = 1.0; pf.y_axis.z = 0.0;
    pf.normal.x = 0.0; pf.normal.y = 0.0; pf.normal.z = 1.0;
    return pf;
}

static PlaneFrame make_top_plane(void)
{
    PlaneFrame pf;
    pf.origin.x = 0.0;
    pf.origin.y = 0.0;
    pf.origin.z = 4.0;

    Vec3 n = {0.2, 0.4, 0.87};
    pf.normal = v3_norm(n);

    Vec3 tmp = {1.0, 0.0, 0.0};
    Vec3 x = v3_sub(tmp, v3_scale(pf.normal, v3_dot(tmp, pf.normal)));
    if (v3_len(x) < 1e-6) {
        tmp.x = 0.0; tmp.y = 1.0; tmp.z = 0.0;
        x = v3_sub(tmp, v3_scale(pf.normal, v3_dot(tmp, pf.normal)));
    }
    pf.x_axis = v3_norm(x);
    pf.y_axis = v3_norm(v3_cross(pf.normal, pf.x_axis));

    return pf;
}

static PlaneFrame lerp_plane(const PlaneFrame *a,
                             const PlaneFrame *b,
                             double t)
{
    double omt = 1.0 - t;
    PlaneFrame pf;
    pf.origin.x = omt * a->origin.x + t * b->origin.x;
    pf.origin.y = omt * a->origin.y + t * b->origin.y;
    pf.origin.z = omt * a->origin.z + t * b->origin.z;

    pf.x_axis.x = omt * a->x_axis.x + t * b->x_axis.x;
    pf.x_axis.y = omt * a->x_axis.y + t * b->x_axis.y;
    pf.x_axis.z = omt * a->x_axis.z + t * b->x_axis.z;

    pf.y_axis.x = omt * a->y_axis.x + t * b->y_axis.x;
    pf.y_axis.y = omt * a->y_axis.y + t * b->y_axis.y;
    pf.y_axis.z = omt * a->y_axis.z + t * b->y_axis.z;

    pf.normal.x = omt * a->normal.x + t * b->normal.x;
    pf.normal.y = omt * a->normal.y + t * b->normal.y;
    pf.normal.z = omt * a->normal.z + t * b->normal.z;

    pf.x_axis = v3_norm(pf.x_axis);
    pf.y_axis = v3_norm(pf.y_axis);
    pf.normal = v3_norm(pf.normal);

    return pf;
}

/* ----------------- OpenSCAD polyhedron output ----------------- */

/* Build point list and face list for polyhedron based on slices.
 * - S slices (steps), each slice has N_pts points.
 * - Points indexed as idx = slice * N_pts + pt.
 * - Faces: bottom cap, top cap, and side quads split into triangles. [web:70][web:72][web:77]
 */
static void emit_openscad_polyhedron(const Polyline3D *sections,
                                     int steps,
                                     int N_pts)
{
    int total_points = steps * N_pts;

    /* Header and points */
    printf("// Auto-generated loft from C program\n");
    printf("// steps=%d, points_per_slice=%d\n\n", steps, N_pts);
    printf("points = [\n");
    for (int s = 0; s < steps; ++s) {
        for (int i = 0; i < N_pts; ++i) {
            int idx = s * N_pts + i;
            Vec3 p = sections[s].pt[i];
            printf("  [% .6f, % .6f, % .6f]%s\n",
                   p.x, p.y, p.z,
                   (idx == total_points-1) ? "" : ",");
        }
    }
    printf("];\n\n");

    printf("faces = [\n");

    /* Bottom cap (slice 0): fan around point 0. */
    for (int i = 1; i < N_pts-1; ++i) {
        printf("  [ %d, %d, %d ],\n", 0, i, i+1);
    }

    /* Top cap (slice steps-1): fan, but reverse winding. */
    int base_top = (steps-1) * N_pts;
    for (int i = 1; i < N_pts-1; ++i) {
        printf("  [ %d, %d, %d ],\n",
               base_top, base_top + i+1, base_top + i);
    }

    /* Side faces between slices s and s+1.
     * For each segment i..i+1, form quad: (s,i)->(s,i+1)->(s+1,i+1)->(s+1,i)
     * Then split quad into two triangles. [web:72][web:77]
     */
    for (int s = 0; s < steps-1; ++s) {
        int base0 = s * N_pts;
        int base1 = (s+1) * N_pts;
        for (int i = 0; i < N_pts-1; ++i) {
            int a = base0 + i;
            int b = base0 + i+1;
            int c = base1 + i+1;
            int d = base1 + i;
            printf("  [ %d, %d, %d ],\n", a, b, c);
            printf("  [ %d, %d, %d ],\n", a, c, d);
        }
    }

    printf("];\n\n");
    printf("polyhedron(points = points, faces = faces, convexity = 10);\n");
}

/* ----------------- Main: loft + OpenSCAD output + optional SDL ----------------- */

int main(int argc, char **argv)
{
    (void)argc; (void)argv;

    BezierPath2D path_bottom = make_square2();
    BezierPath2D path_top    = make_circle2();

    ShapeMapping2D map2;
    double flat_tol = 0.002;
    if (!build_shape_mapping2(&path_bottom, &path_top, flat_tol, &map2)) {
        fprintf(stderr, "Failed to build 2D shape mapping\n");
        free(path_bottom.seg);
        free(path_top.seg);
        return 1;
    }

    PlaneFrame plane_bottom = make_bottom_plane();
    PlaneFrame plane_top    = make_top_plane();

    int steps = 16;
    if (steps < 2) steps = 2;

    Polyline3D *sections = (Polyline3D*)calloc((size_t)steps,
                                               sizeof(Polyline3D));
    if (!sections) {
        free_shape_mapping2(&map2);
        free(path_bottom.seg);
        free(path_top.seg);
        return 1;
    }

    for (int k = 0; k < steps; ++k) {
        double t = (double)k / (double)(steps - 1);
        PlaneFrame pf = lerp_plane(&plane_bottom, &plane_top, t);
        Polyline2D pl2 = tween_polyline2(&map2, t);
        sections[k] = lift_polyline_to_plane(&pl2, &pf);
        free_polyline2(&pl2);
    }

    /* Emit OpenSCAD source describing the loft as a polyhedron. [web:70][web:72] */
    if (sections[0].count > 1) {
        emit_openscad_polyhedron(sections, steps, sections[0].count);
    }

#ifdef USE_SDL2
    /* Optional SDL2 visualization: perspective wireframe of all slices. */
    if (SDL_Init(SDL_INIT_VIDEO) == 0) {
        int width = 800, height = 800;
        SDL_Window *win = SDL_CreateWindow("Loft wireframe",
                                           SDL_WINDOWPOS_CENTERED,
                                           SDL_WINDOWPOS_CENTERED,
                                           width, height,
                                           SDL_WINDOW_SHOWN);
        if (win) {
            SDL_Renderer *ren = SDL_CreateRenderer(win, -1,
                                                   SDL_RENDERER_ACCELERATED |
                                                   SDL_RENDERER_PRESENTVSYNC);
            if (ren) {
                Camera cam;
                cam.pos.x = 6.0; cam.pos.y = 5.0; cam.pos.z = 6.0;
                Vec3 target = {0.0, 0.0, 2.0};
                cam.forward = v3_norm(v3_sub(target, cam.pos));
                cam.up.x = 0.0; cam.up.y = 0.0; cam.up.z = 1.0;

                Mat4 V = camera_view_matrix(&cam);
                Mat4 P = perspective_matrix(60.0,
                                            (double)width / (double)height,
                                            0.1, 50.0);
                Mat4 VP = mat4_mul(P, V);

                int running = 1;
                while (running) {
                    SDL_Event ev;
                    while (SDL_PollEvent(&ev)) {
                        if (ev.type == SDL_QUIT) running = 0;
                        if (ev.type == SDL_KEYDOWN &&
                            ev.key.keysym.sym == SDLK_ESCAPE) running = 0;
                    }

                    SDL_SetRenderDrawColor(ren, 0, 0, 0, 255);
                    SDL_RenderClear(ren);

                    SDL_SetRenderDrawColor(ren, 160, 160, 160, 255);
                    for (int k = 0; k < steps; ++k) {
                        draw_polyline3d(ren, &VP, &sections[k], width, height);
                    }

                    SDL_SetRenderDrawColor(ren, 80, 255, 80, 255);
                    int N_pts = sections[0].count;
                    for (int i = 0; i < N_pts; ++i) {
                        for (int k = 0; k < steps - 1; ++k) {
                            Vec3 a = sections[k].pt[i];
                            Vec3 b = sections[k+1].pt[i];
                            double ax, ay, az, bx, by, bz;
                            if (!project_point(&VP, a, &ax, &ay, &az)) continue;
                            if (!project_point(&VP, b, &bx, &by, &bz)) continue;
                            int sx1 = (int)((ax * 0.5 + 0.5) * (double)width);
                            int sy1 = (int)((-ay * 0.5 + 0.5) * (double)height);
                            int sx2 = (int)((bx * 0.5 + 0.5) * (double)width);
                            int sy2 = (int)((-by * 0.5 + 0.5) * (double)height);
                            SDL_RenderDrawLine(ren, sx1, sy1, sx2, sy2);
                        }
                    }

                    SDL_RenderPresent(ren);
                }

                SDL_DestroyRenderer(ren);
            }
            SDL_DestroyWindow(win);
        }
        SDL_Quit();
    }
#endif

    for (int k = 0; k < steps; ++k) free_polyline3(&sections[k]);
    free(sections);
    free_shape_mapping2(&map2);
    free(path_bottom.seg);
    free(path_top.seg);

    return 0;
}