tut03_transform.cpp

/// @file tut03_transform.cpp
/// @brief Guided tutorial for GSTransform, the rigid 4x4 homogeneous transform
///        that places a canonical local frame into the world.
///
/// @page tutorial_03 Tutorial 03 — GSTransform: rigid transforms
///
/// This tutorial walks through GSTransform, the L0-math object GOST uses to map
/// a component's LOCAL canonical frame (a z=0 plane, an origin-centred sphere)
/// into the shared WORLD frame and back. It is a runnable, self-contained C++11
/// program that links only the GOST core (no ROOT); every claim it makes is
/// checked with an assert helper, so the tutorial doubles as a smoke test and
/// returns a non-zero exit code on any failure.
///
/// @section tut03_contracts Domain contracts
/// The tutorial demonstrates the contracts that are easy to get wrong:
///  - A GSTransform is read as LOCAL -> PARENT/WORLD: TransformPoint() takes a
///    point in the local frame and returns it in the world frame.
///  - Composition uses operator* with the PARENT ON THE LEFT: @c (A*B) applies
///    @c B first then @c A, so it is non-commutative and order matters.
///  - Inverse() is the analytic rigid inverse (world -> local) and round-trips.
///  - TransformPoint applies rotation + translation; TransformDirection applies
///    rotation only and is ALSO how a surface normal transforms — valid only
///    because the transform is rigid (orthonormal rotation).
///  - GSVector is passed BY VALUE throughout (its accessors are non-const).
///
/// @section tut03_sections Walkthrough
///  1. the identity transform; 2. pure translation; 3. axis rotations X/Y/Z;
///  4. axis-angle rotation; 5. composition and parent-on-the-left ordering;
///  6. the analytic Inverse(); 7. TransformPoint vs TransformDirection.
///
/// @section tut03_build Build and run
/// Built as the @c tut03_transform target at the build root (see
/// @c tutorials/CMakeLists.txt). It runs standalone with no arguments and exits
/// 0 when all checks pass. Manual build from the repo root:
/// @code
/// g++ -std=c++11 -Iinclude tutorials/tut03_transform.cpp src/*.cpp -o tut03
/// @endcode

// tut03_transform.cpp — GOST tutorial: GSTransform (L0 math layer).
//
// GSTransform is the rigid-body 4x4 homogeneous transform that GOST uses to
// PLACE a canonical geometry (a z=0 plane, an origin-centred sphere) somewhere
// in the world. Every optical element in GOST is defined in its own LOCAL
// canonical frame (see DESIGN.md §3); GSTransform is the bridge that maps that
// local frame to the WORLD frame and back.
//
// This tutorial is a guided walkthrough. Each numbered section explains WHY a
// piece of the API exists and the domain contract it obeys, then demonstrates
// it with a few lines of code whose results are checked with a tiny Assert()
// helper. Because every claim is asserted, this tutorial doubles as a smoke
// test: main() returns non-zero if any check fails.
//
// Build (from the repo root, no ROOT needed):
//   g++ -std=c++11 -Iinclude tutorials/tut03_transform.cpp src/*.cpp -o tut03
//
// KEY DOMAIN CONTRACTS demonstrated here:
//   * A GSTransform is interpreted as LOCAL -> PARENT/WORLD. TransformPoint()
//     takes a point expressed in the local frame and returns it in the world.
//   * Composition uses operator*, and the PARENT GOES ON THE LEFT:
//     (A * B) applies B first, then A. Order matters (non-commutative).
//   * GSVector is passed BY VALUE everywhere (project idiom): its X()/Y()/Z()
//     accessors are non-const, so a const& would not compile. Every helper
//     below therefore takes GSVector by value.

#include "GSVector.hpp"
#include "GSTransform.hpp"
#include "GSConstants.hpp"

#include "cmath"
#include "iostream"
#include "iomanip"
#include "string"

// ------------------------------------------------------------------ utilities
static int gFail = 0;

// Minimal check helper: prints [OK]/[FAIL], counts failures. Same spirit as the
// Assert() in exe/demo_procedure.cpp so the tutorial is also a self-test.
static void Assert(bool cond, const std::string& what){
    std::cout << (cond ? "  [OK]   " : "  [FAIL] ") << what << std::endl;
    if(!cond) ++gFail;
}

// Approximate scalar equality (floating-point round-off tolerance).
static bool Approx(double a, double b, double tol = 1e-9){
    return std::fabs(a - b) <= tol;
}

// Approximate vector equality. GSVector is taken BY VALUE (project idiom: the
// X()/Y()/Z() accessors are non-const).
static bool VecEq(GSVector a, GSVector b, double tol = 1e-9){
    return Approx(a.X(), b.X(), tol) && Approx(a.Y(), b.Y(), tol) &&
           Approx(a.Z(), b.Z(), tol);
}

int main(){
    std::cout << "======================================================\n";
    std::cout << " GOST tutorial 03 — GSTransform (rigid 4x4 transform)\n";
    std::cout << "======================================================\n";
    std::cout << std::fixed << std::setprecision(4);

    // ==== 1. The identity transform ====
    // A default-constructed GSTransform is the identity: it leaves every point
    // and direction untouched. Physically it means "the local frame IS the
    // world frame" — no displacement, no rotation.
    {
        GSTransform I;                       // identity
        GSVector p(1, 2, 3);
        Assert(VecEq(I.TransformPoint(p), GSVector(1, 2, 3)),
               "1. identity leaves a point unchanged");
        Assert(VecEq(I.TransformDirection(GSVector(0, 0, 1)), GSVector(0, 0, 1)),
               "1. identity leaves a direction unchanged");
    }

    // ==== 2. Pure translation ====
    // Translation(t) builds a transform whose rotation is identity and whose
    // translation is t (t is given in the parent/world frame). It slides the
    // local origin to t. This is how you position a lens surface along the
    // optical axis, e.g. Translation((0,0,z)) to move a sphere centre to z.
    {
        GSVector t(10, 20, 30);
        GSTransform T = GSTransform::Translation(t);
        GSVector localPoint(1, 2, 3);
        // A point sitting at (1,2,3) in the LOCAL frame lands at (11,22,33) in
        // the WORLD frame once the frame is translated by (10,20,30).
        GSVector world = T.TransformPoint(localPoint);
        std::cout << "    local (1,2,3) -> world " << world.Print(true) << "\n";
        Assert(VecEq(world, GSVector(11, 22, 33)),
               "2. Translation maps local point into the world frame");
    }

    // ==== 3. Axis rotations RotationX / RotationY / RotationZ ====
    // The three factory rotations are right-handed about the named world axis
    // and carry ZERO translation. Angles are in RADIANS. Handy conversion:
    // GSConstant::deg() is the radians-per-degree factor, and
    // GSConstant::Pi()/2 is a clean 90 degrees.
    //
    // Right-hand rule sanity (90 degrees):
    //   RotationZ: +x -> +y      RotationX: +y -> +z      RotationY: +z -> +x
    {
        const double half = GSConstant::Pi() / 2.0;   // 90 degrees
        GSTransform Rz = GSTransform::RotationZ(half);
        GSTransform Rx = GSTransform::RotationX(half);
        GSTransform Ry = GSTransform::RotationY(half);

        Assert(VecEq(Rz.TransformPoint(GSVector(1, 0, 0)), GSVector(0, 1, 0)),
               "3. RotationZ(90) sends +x to +y");
        Assert(VecEq(Rx.TransformPoint(GSVector(0, 1, 0)), GSVector(0, 0, 1)),
               "3. RotationX(90) sends +y to +z");
        Assert(VecEq(Ry.TransformPoint(GSVector(0, 0, 1)), GSVector(1, 0, 0)),
               "3. RotationY(90) sends +z to +x");

        // A rotation leaves the origin fixed (no translation component).
        Assert(VecEq(Rz.TransformPoint(GSVector(0, 0, 0)), GSVector(0, 0, 0)),
               "3. a pure rotation fixes the origin");
    }

    // ==== 4. Axis-angle rotation (Rodrigues) ====
    // Rotation(axis, rad) rotates right-handedly about an arbitrary axis. The
    // axis need NOT be a unit vector — it is normalised internally — so you can
    // pass any convenient direction. Rotating about +Z must reproduce
    // RotationZ exactly, which is the check below.
    {
        const double half = GSConstant::Pi() / 2.0;
        GSTransform Ra = GSTransform::Rotation(GSVector(0, 0, 5), half); // non-unit axis
        GSTransform Rz = GSTransform::RotationZ(half);
        bool same = true;
        for(int i = 0; i < 4; ++i)
            for(int j = 0; j < 4; ++j)
                if(!Approx(Ra.At(i, j), Rz.At(i, j))) same = false;
        Assert(same, "4. Rotation(axis=+Z, 90) equals RotationZ(90)");
    }

    // ==== 5. Composition with operator*  —  PARENT ON THE LEFT ====
    // operator* composes transforms. The crucial rule: the LEFT operand is the
    // OUTER (parent) transform, applied AFTER the right one. Formally
    //     (A * B).TransformPoint(p) == A.TransformPoint( B.TransformPoint(p) ).
    // So to "rotate first, then translate" you write Translation * Rotation.
    // This matches assembly flattening in DESIGN.md §5: T_world = T_parent * T_child.
    {
        GSVector p(1, 2, 3);
        GSVector t(10, 20, 30);
        GSTransform T  = GSTransform::Translation(t);
        GSTransform Rz = GSTransform::RotationZ(GSConstant::Pi() / 2.0);

        GSTransform rotThenTrans = T * Rz;     // rotate FIRST, then translate
        GSVector viaComposed = rotThenTrans.TransformPoint(p);
        GSVector viaManual   = T.TransformPoint(Rz.TransformPoint(p));
        Assert(VecEq(viaComposed, viaManual),
               "5. (A*B).apply(p) == A.apply(B.apply(p))  (parent on the left)");

        // Concretely: Rz(1,2,3) = (-2,1,3); then +t = (8,21,33).
        std::cout << "    T*Rz applied to (1,2,3) = " << viaComposed.Print(true) << "\n";
        Assert(VecEq(viaComposed, GSVector(8, 21, 33)),
               "5. rotate-then-translate gives (8,21,33)");

        // ORDER MATTERS: swapping the operands changes the result. Translating
        // first and THEN rotating moves the point somewhere else entirely.
        GSTransform transThenRot = Rz * T;     // translate first, then rotate
        GSVector swapped = transThenRot.TransformPoint(p);
        std::cout << "    Rz*T applied to (1,2,3) = " << swapped.Print(true)
                  << "  (different!)\n";
        Assert(!VecEq(swapped, viaComposed),
               "5. composition is NON-commutative: T*Rz != Rz*T");
    }

    // ==== 6. Analytic Inverse() — world -> local ====
    // Because a GSTransform is rigid (orthonormal rotation + translation), its
    // inverse is computed analytically (R^-1 = R^T, t^-1 = -R^T t) rather than
    // by general matrix inversion. Inverse() maps WORLD -> LOCAL, undoing the
    // placement. Round-trip: T^-1 . T applied to any point returns that point.
    {
        // A non-trivial composite placement.
        GSTransform T = GSTransform::Translation(GSVector(10, 20, 30))
                      * GSTransform::RotationX(0.7)
                      * GSTransform::RotationY(-1.3);
        GSVector p(1, 2, 3);
        GSVector roundTrip = T.Inverse().TransformPoint(T.TransformPoint(p));
        Assert(VecEq(roundTrip, p),
               "6. Inverse round-trip on a point: T^-1 (T p) == p");

        // And T * T^-1 is the 4x4 identity matrix.
        GSMatrix I = (T * T.Inverse()).Matrix();
        bool ident = true;
        for(int i = 0; i < 4; ++i)
            for(int j = 0; j < 4; ++j)
                if(!Approx(I.At(i, j), (i == j ? 1.0 : 0.0))) ident = false;
        Assert(ident, "6. T * T^-1 == identity matrix");
    }

    // ==== 7. TransformPoint vs TransformDirection ====
    // These differ in ONE thing: whether the translation part is applied.
    //   * TransformPoint    — rotation + translation (homogeneous w = 1). Use it
    //                         for positions (a ray origin, a hit point).
    //   * TransformDirection— rotation ONLY (homogeneous w = 0, translation
    //                         ignored). Use it for directions (a ray direction)
    //                         and for SURFACE NORMALS.
    //
    // Why does the SAME direction routine also transform a normal correctly?
    // Because the transform is RIGID: its rotation block is orthonormal, so the
    // inverse-transpose that a general (non-rigid) transform would need for
    // normals collapses back to the rotation itself. This shortcut is valid
    // ONLY for rigid transforms — if GOST ever allowed shear/scale, normals
    // would need the inverse-transpose instead.
    {
        GSVector t(10, 20, 30);
        GSTransform T = GSTransform::Translation(t);
        GSVector dir(1, 0, 0);

        // A pure translation must NOT move a direction (only points shift).
        Assert(VecEq(T.TransformDirection(dir), GSVector(1, 0, 0)),
               "7. TransformDirection ignores translation");
        // ...while the same transform DOES move a point.
        Assert(VecEq(T.TransformPoint(dir), GSVector(11, 20, 30)),
               "7. TransformPoint applies translation");

        // Rigid transforms preserve unit length: a unit direction in local maps
        // to a unit direction in world (no renormalisation needed).
        GSTransform R = GSTransform::RotationX(0.7) * GSTransform::RotationY(-1.3);
        GSVector n = R.TransformDirection(GSVector(0, 0, 1)); // a unit normal
        Assert(Approx(n.Norm(), 1.0),
               "7. rigid transform preserves unit length of a direction/normal");

        // The normal round-trips through Inverse() just like a direction.
        GSVector back = R.Inverse().TransformDirection(n);
        Assert(VecEq(back, GSVector(0, 0, 1)),
               "7. normal maps world->local via Inverse().TransformDirection");
    }

    // ------------------------------------------------------------------ summary
    std::cout << "\n------------------------------------------------------\n";
    if(gFail == 0){
        std::cout << "tut03: all checks passed.\n";
        return 0;
    }
    std::cout << "tut03: " << gFail << " check(s) FAILED.\n";
    return 1;
}