Kuadrik Yüzeyler

Bu bölümün laboratuvarı, kitaptaki örnekleri yeniden üretmek ve sayısal sonuçları incelemek için hazırlanmıştır. Matematiksel tanımlar, ispatlar ve sorular kitapta yer alır.

Laboratuvar kodu

"""Bölüm 22: kuadrik yüzeyler laboratuvarı.

Uygulama, bilinen merkez, yarı eksen ve yönelime sahip sentetik bir üç
eksenli elipsoidi genel ikinci derece katsayılarına dönüştürür; ardından
merkezi, yarı eksenleri ve ana yönleri yalnız bu katsayılardan geri kazanır.
Veriler ölçülmüş değildir. Çalışma bir uyumlama, üretim toleransı veya kabul
modeli içermez.
"""

from __future__ import annotations

import math

import numpy as np

from agbook import (
    canonical_array_sha256,
    quadric_diagnostics,
    quadric_plane_trace,
    quadric_residuals,
    quadric_tangent_plane_diagnostics,
    quadratic_conic_diagnostics,
    transform_quadric,
)


def rotation_xyz(ax: float, ay: float, az: float) -> np.ndarray:
    """Sağ elli ``Rz(az) Ry(ay) Rx(ax)`` dönmesini kurar."""

    cx, sx = math.cos(ax), math.sin(ax)
    cy, sy = math.cos(ay), math.sin(ay)
    cz, sz = math.cos(az), math.sin(az)
    rx = np.array([[1, 0, 0], [0, cx, -sx], [0, sx, cx]], dtype=float)
    ry = np.array([[cy, 0, sy], [0, 1, 0], [-sy, 0, cy]], dtype=float)
    rz = np.array([[cz, -sz, 0], [sz, cz, 0], [0, 0, 1]], dtype=float)
    return rz @ ry @ rx


def world_coefficients(
    canonical_coefficients: np.ndarray,
    rotation: np.ndarray,
    center: np.ndarray,
) -> np.ndarray:
    """``X=center+rotation@Y`` modelini dünya katsayılarına çevirir."""

    return transform_quadric(
        canonical_coefficients,
        rotation.T,
        -rotation.T @ center,
    )


def main() -> None:
    length_unit = "mm"
    center_mm = np.array([120.0, -80.0, 45.0])
    semi_axes_mm = np.array([60.0, 35.0, 20.0])
    rotation = rotation_xyz(
        math.radians(15.0),
        math.radians(-20.0),
        math.radians(35.0),
    )
    canonical_coefficients = np.array(
        [
            1.0 / semi_axes_mm[0] ** 2,
            0.0,
            0.0,
            1.0 / semi_axes_mm[1] ** 2,
            0.0,
            1.0 / semi_axes_mm[2] ** 2,
            0.0,
            0.0,
            0.0,
            -1.0,
        ]
    )
    coefficients_mm = world_coefficients(
        canonical_coefficients,
        rotation,
        center_mm,
    )
    diagnostics_mm = quadric_diagnostics(
        coefficients_mm,
        coordinate_scale=100.0,
    )

    alignment = np.abs(rotation.T @ diagnostics_mm.principal_directions)
    principal_trace_types: list[str] = []
    principal_trace_axes: list[np.ndarray] = []
    for indices in ((0, 1), (0, 2), (1, 2)):
        basis = rotation[:, indices]
        trace = quadric_plane_trace(coefficients_mm, center_mm, basis)
        conic = quadratic_conic_diagnostics(
            trace.conic_coefficients,
            coordinate_scale=100.0,
        )
        principal_trace_types.append(conic.locus_type)
        principal_trace_axes.append(conic.semi_axes)

    surface_point_mm = center_mm + semi_axes_mm[0] * rotation[:, 0]
    tangent = quadric_tangent_plane_diagnostics(
        coefficients_mm,
        surface_point_mm,
        coordinate_scale=100.0,
        gradient_tolerance=1.0e-12,
        residual_tolerance=1.0e-12,
    )

    parameters = np.array(
        [
            [1.0, 0.0, 0.0],
            [-1.0, 0.0, 0.0],
            [0.0, 1.0, 0.0],
            [0.0, -1.0, 0.0],
            [0.0, 0.0, 1.0],
            [0.0, 0.0, -1.0],
        ]
    )
    surface_points_mm = center_mm + (parameters * semi_axes_mm) @ rotation.T
    residuals_mm = quadric_residuals(
        coefficients_mm,
        surface_points_mm,
        coordinate_scale=100.0,
    )

    motion = rotation_xyz(-0.19, 0.27, -0.41)
    translation_mm = np.array([-35.0, 70.0, 25.0])
    moved_coefficients = transform_quadric(
        coefficients_mm,
        motion.T,
        -motion.T @ translation_mm,
    )
    moved_diagnostics = quadric_diagnostics(
        moved_coefficients,
        coordinate_scale=100.0,
    )

    mm_per_m = 1000.0
    coefficients_m = coefficients_mm * np.array(
        [
            mm_per_m**2,
            mm_per_m**2,
            mm_per_m**2,
            mm_per_m**2,
            mm_per_m**2,
            mm_per_m**2,
            mm_per_m,
            mm_per_m,
            mm_per_m,
            1.0,
        ]
    )
    diagnostics_m = quadric_diagnostics(
        coefficients_m,
        coordinate_scale=0.1,
    )

    assert diagnostics_mm.locus_type == "ellipsoid"
    np.testing.assert_allclose(diagnostics_mm.center, center_mm, atol=3.0e-13)
    np.testing.assert_allclose(diagnostics_mm.semi_axes, semi_axes_mm, atol=3.0e-13)
    np.testing.assert_allclose(alignment, np.eye(3), atol=3.0e-14)
    assert principal_trace_types == ["ellipse", "ellipse", "ellipse"]
    np.testing.assert_allclose(principal_trace_axes[0], [60.0, 35.0])
    np.testing.assert_allclose(principal_trace_axes[1], [60.0, 20.0])
    np.testing.assert_allclose(principal_trace_axes[2], [35.0, 20.0])
    assert tangent.regular
    assert tangent.point_on_surface
    np.testing.assert_allclose(residuals_mm, np.zeros(6), atol=3.0e-15)
    assert moved_diagnostics.locus_type == "ellipsoid"
    np.testing.assert_allclose(moved_diagnostics.semi_axes, semi_axes_mm, atol=4.0e-13)
    np.testing.assert_allclose(
        moved_diagnostics.center,
        motion @ center_mm + translation_mm,
        atol=4.0e-13,
    )
    assert diagnostics_m.locus_type == "ellipsoid"
    np.testing.assert_allclose(diagnostics_m.center, center_mm / mm_per_m, atol=4.0e-16)
    np.testing.assert_allclose(
        diagnostics_m.semi_axes,
        semi_axes_mm / mm_per_m,
        atol=4.0e-16,
    )

    signature_values = np.concatenate(
        [
            center_mm,
            semi_axes_mm,
            rotation.ravel(),
            canonical_coefficients,
            coefficients_mm,
            diagnostics_mm.center,
            diagnostics_mm.semi_axes,
            diagnostics_mm.principal_directions.ravel(),
            np.concatenate(principal_trace_axes),
            surface_point_mm,
            tangent.gradient,
            tangent.plane_coefficients,
            surface_points_mm.ravel(),
            residuals_mm,
            moved_coefficients,
            diagnostics_m.center,
            diagnostics_m.semi_axes,
        ]
    )

    print("uzunluk_birimi:", length_unit)
    print("model: sentetik_uc_eksenli_elipsoid")
    print("olculmus_veri_mi: False")
    print("uyumlama_yapildi_mi: False")
    print("uretim_toleransi_modeli_mi: False")
    print("merkez_mm:", diagnostics_mm.center)
    print("yari_eksenler_mm:", diagnostics_mm.semi_axes)
    print("sinif:", diagnostics_mm.locus_type)
    print("kuadratik_rank:", diagnostics_mm.quadratic_rank)
    print("atalet:", diagnostics_mm.inertia)
    print("ana_yon_hizalama_mutlak:")
    print(alignment)
    print("asal_iz_turleri:", principal_trace_types)
    print("asal_iz_yari_eksenleri_mm:", principal_trace_axes)
    print("teget_noktasi_mm:", surface_point_mm)
    print("teget_duzlem_katsayilari:", tangent.plane_coefficients)
    print("en_buyuk_normalize_artik:", f"{np.max(np.abs(residuals_mm)):.12e}")
    print("rijit_harekette_sinif_korundu_mu:", moved_diagnostics.locus_type == "ellipsoid")
    print("metrede_yari_eksenler:", diagnostics_m.semi_axes)
    print("goreli_siniflandirma_esigi:", diagnostics_mm.relative_tolerance)
    print("sayisal_belirsizlik_var_mi:", diagnostics_mm.numerically_ambiguous)
    print("bilimsel_imza:", canonical_array_sha256(signature_values))


if __name__ == "__main__":
    main()

uzunluk_birimi: mm
model: sentetik_uc_eksenli_elipsoid
olculmus_veri_mi: False
uyumlama_yapildi_mi: False
uretim_toleransi_modeli_mi: False
merkez_mm: [120. -80.  45.]
yari_eksenler_mm: [60. 35. 20.]
sinif: ellipsoid
kuadratik_rank: 3
atalet: (3, 0, 0)
ana_yon_hizalama_mutlak:
[[1.00000000e+00 2.81065211e-16 7.03286187e-17]
 [2.93542648e-17 1.00000000e+00 1.10447058e-16]
 [0.00000000e+00 2.77555756e-17 1.00000000e+00]]
asal_iz_turleri: ['ellipse', 'ellipse', 'ellipse']
asal_iz_yari_eksenleri_mm: [array([60., 35.]), array([60., 20.]), array([35., 20.])]
teget_noktasi_mm: [166.18506788 -47.66086732  65.5212086 ]
teget_duzlem_katsayilari: [   0.76975113    0.53898554    0.34202014 -124.64219863]
en_buyuk_normalize_artik: 2.220446049250e-16
rijit_harekette_sinif_korundu_mu: True
metrede_yari_eksenler: [0.06  0.035 0.02 ]
goreli_siniflandirma_esigi: 0.0
sayisal_belirsizlik_var_mi: False
bilimsel_imza: 1bf892b47806c97e60aa819b4a5b927548f608c25f863e2d81178a2cc8d6be48

# Kitaptaki kod blokları; bu dosya içinde sırayla çalıştırılır.
# --- 1. Genel bir kuadriği sınıflandırma ve iz çıkarma ---
import numpy as np
from agbook import (
    quadric_diagnostics,
    quadric_plane_trace,
    quadratic_conic_diagnostics,
)

# x^2/4 + y^2/9 + z^2/16 = 1
coeff = np.array([1/4, 0, 0, 1/9, 0, 1/16, 0, 0, 0, -1])
d = quadric_diagnostics(coeff, coordinate_scale=4.0)
assert d.locus_type == "ellipsoid"
np.testing.assert_allclose(d.center, [0, 0, 0])
np.testing.assert_allclose(d.semi_axes, [4, 3, 2])

# z = 0 düzlemindeki iz
B = np.array([[1, 0], [0, 1], [0, 0]], dtype=float)
trace = quadric_plane_trace(coeff, [0, 0, 0], B)
conic = quadratic_conic_diagnostics(trace.conic_coefficients)
assert trace.orthonormal and conic.locus_type == "ellipse"

1. Genel bir kuadriği sınıflandırma ve iz çıkarma

# Kitaptaki kod blokları; bu dosya içinde sırayla çalıştırılır.
# --- 1. On katsayılı sınıfların ortak çarpan testi ---
import numpy as np
from agbook import quadric_diagnostics

cases = {
    "ellipsoid": [1,0,0,1,0,1,0,0,0,-1],
    "hyperboloid_one_sheet": [1,0,0,1,0,-1,0,0,0,-1],
    "hyperboloid_two_sheets": [-1,0,0,-1,0,1,0,0,0,-1],
    "elliptic_cone": [1,0,0,1,0,-1,0,0,0,0],
    "elliptic_paraboloid": [1,0,0,1,0,0,0,0,-1,0],
    "hyperbolic_paraboloid": [1,0,0,-1,0,0,0,0,-1,0],
    "elliptic_cylinder": [1,0,0,1,0,0,0,0,0,-1],
    "hyperbolic_cylinder": [1,0,0,-1,0,0,0,0,0,-1],
    "parabolic_cylinder": [1,0,0,0,0,0,0,-1,0,0],
    "intersecting_planes": [1,0,0,-1,0,0,0,0,0,0],
    "parallel_planes": [1,0,0,0,0,0,0,0,0,-1],
    "double_plane": [1,0,0,0,0,0,0,0,0,0],
    "line": [1,0,0,1,0,0,0,0,0,0],
    "point": [1,0,0,1,0,1,0,0,0,0],
    "plane": [0,0,0,0,0,0,1,0,0,-1],
    "all_space": [0,0,0,0,0,0,0,0,0,0],
    "empty": [1,0,0,1,0,1,0,0,0,1],
}
for factor in (1e-12, -3.0, 1.0, 1e8):
    for expected, coeff in cases.items():
        got = quadric_diagnostics(factor*np.array(coeff)).locus_type
        assert got == expected, (factor, expected, got)

# --- 2. Elipsoidin rijit hareket değişmezliği ---
import numpy as np
from agbook import quadric_diagnostics, transform_quadric

rng = np.random.default_rng(2201)
axes = np.array([6.0, 3.5, 2.0])
local = np.array([1/36,0,0,1/12.25,0,1/4,0,0,0,-1])
for _ in range(20):
    raw = rng.normal(size=(3,3))
    R, _ = np.linalg.qr(raw)
    if np.linalg.det(R) < 0:
        R[:, -1] *= -1
    C = rng.normal(size=3)
    world = transform_quadric(local, R.T, -R.T @ C)
    d = quadric_diagnostics(world, coordinate_scale=10.0)
    assert d.locus_type == "ellipsoid"
    np.testing.assert_allclose(d.center, C, atol=2e-13)
    np.testing.assert_allclose(d.semi_axes, axes, atol=2e-13)

# --- 3. 2B iz ile 3B polinomun eşdeğerliği ---
import numpy as np
from agbook import quadric_matrices, quadric_plane_trace

rng = np.random.default_rng(2202)
coeff = np.array([2,1,-2,3,1,4,-5,2,3,-7], dtype=float)
Q, g, _ = quadric_matrices(coeff)
normals = np.eye(3).tolist() + [[1,1,1]]
for normal in normals:
    n = np.asarray(normal, dtype=float)
    n /= np.linalg.norm(n)
    seed = np.eye(3)[np.argmin(np.abs(n))]
    u = seed - (seed @ n)*n
    u /= np.linalg.norm(u)
    v = np.cross(n, u)
    B = np.column_stack((u, v))
    P = 0.7*n
    trace = quadric_plane_trace(coeff, P, B)
    assert trace.orthonormal
    a,b,c,d,e,f = trace.conic_coefficients
    for st in rng.normal(size=(10,2)):
        s,t = st
        value2 = a*s*s + b*s*t + c*t*t + d*s + e*t + f
        X = P + B @ st
        value3 = X @ Q @ X + 2*g @ X + coeff[-1]
        np.testing.assert_allclose(value2, value3, atol=2e-13)

# --- 4. Sınıra yakın elipsoit--silindir ailesi ---
from agbook import quadric_diagnostics

for tolerance in (0.0, 1e-12, 1e-8):
    for k in range(15):
        epsilon = 10.0**(-k)
        coeff = [1,0,0,1,0,epsilon,0,0,0,-1]
        d = quadric_diagnostics(coeff, relative_tolerance=tolerance)
        print(tolerance, k, d.locus_type,
              d.quadratic_rank, d.numerically_ambiguous)

1. On katsayılı sınıfların ortak çarpan testi

2. Elipsoidin rijit hareket değişmezliği

3. 2B iz ile 3B polinomun eşdeğerliği

4. Sınıra yakın elipsoit--silindir ailesi
0.0 0 ellipsoid 3 False
0.0 1 ellipsoid 3 False
0.0 2 ellipsoid 3 False
0.0 3 ellipsoid 3 False
0.0 4 ellipsoid 3 False
0.0 5 ellipsoid 3 False
0.0 6 ellipsoid 3 False
0.0 7 ellipsoid 3 False
0.0 8 ellipsoid 3 False
0.0 9 ellipsoid 3 False
0.0 10 ellipsoid 3 False
0.0 11 ellipsoid 3 False
0.0 12 ellipsoid 3 False
0.0 13 ellipsoid 3 False
0.0 14 ellipsoid 3 False
1e-12 0 ellipsoid 3 False
1e-12 1 ellipsoid 3 False
1e-12 2 ellipsoid 3 False
1e-12 3 ellipsoid 3 False
1e-12 4 ellipsoid 3 False
1e-12 5 ellipsoid 3 False
1e-12 6 ellipsoid 3 False
1e-12 7 ellipsoid 3 False
1e-12 8 ellipsoid 3 False
1e-12 9 ellipsoid 3 False
1e-12 10 ellipsoid 3 False
1e-12 11 ellipsoid 3 False
1e-12 12 elliptic_cylinder 2 True
1e-12 13 elliptic_cylinder 2 True
1e-12 14 elliptic_cylinder 2 True
1e-08 0 ellipsoid 3 False
1e-08 1 ellipsoid 3 False
1e-08 2 ellipsoid 3 False
1e-08 3 ellipsoid 3 False
1e-08 4 ellipsoid 3 False
1e-08 5 ellipsoid 3 False
1e-08 6 ellipsoid 3 False
1e-08 7 ellipsoid 3 False
1e-08 8 elliptic_cylinder 2 True
1e-08 9 elliptic_cylinder 2 True
1e-08 10 elliptic_cylinder 2 True
1e-08 11 elliptic_cylinder 2 True
1e-08 12 elliptic_cylinder 2 True
1e-08 13 elliptic_cylinder 2 True
1e-08 14 elliptic_cylinder 2 True

cd analitik-geometri-python-v1.0
python3 -m venv .venv
source .venv/bin/activate
python -m pip install -e ".[dev]"
python labs/bolum_22_kuadrik_yuzeyler.py
python -m pytest

cd analitik-geometri-python-v1.0
py -3.13 -m venv .venv
.venv\Scripts\python.exe -m pip install -e ".[dev]"
.venv\Scripts\python.exe labs/bolum_22_kuadrik_yuzeyler.py
.venv\Scripts\python.exe -m pytest