# -*- coding: utf-8 -*-
#
# Isoliner - грид и изолинии (QGIS).
# © 2026 ООО «Информ++» (www.informpp.ru).
# SPDX-License-Identifier: GPL-2.0-or-later
#
"""Элементы залегания по трёхмерному следу выхода: ядро.

Если граница пласта или разлома оцифрована с отметками, плоскость по ней
восстанавливается однозначно, а из плоскости выпадают угол и азимут
падения. Правило трёх точек это минимальный случай; здесь считается по
всем вершинам сразу.

Метод по Альмендингеру («Structural Geology Algorithms: Vectors and
Tensors»): нормаль к плоскости это собственный вектор матрицы ориентации
при наименьшем собственном значении. Подгонка вида z = a*x + b*y + c
проще, но разваливается на крутых залеганиях - у вертикальной плоскости
коэффициенты уходят в бесконечность. Через собственные векторы крутизна
перестаёт быть особым случаем вовсе.

То же разложение бесплатно даёт меру доверия. Через одну прямую в
пространстве проходит бесконечно много плоскостей, поэтому прямой в плане
след залегания не определяет. Отношения собственных значений это видят:
два больших близки - плоскость определена уверенно, среднее и малое
близки - точки коллинеарны, и правильный ответ здесь отказ, а не уверенное
число из шума округления.
"""

import math

import numpy as np

# наименьшее отношение (s2 - s3) / s2, при котором след считается
# достаточно изогнутым, чтобы задать плоскость
MIN_PLANARITY = 0.15


def fit_plane(pts):
    """Плоскость по облаку точек: (нормаль, собственные значения, центр).

    pts - (N, 3). Нормаль единичная, направлена вверх (nz >= 0).
    Собственные значения по убыванию: s1 >= s2 >= s3. Нормаль отвечает s3.
    """
    p = np.asarray(pts, dtype=float)
    if p.shape[0] < 3:
        return None, None, None
    c = p.mean(axis=0)
    q = p - c
    # матрица ориентации (моментов) 3x3, собственные значения неотрицательны
    m = q.T @ q
    vals, vecs = np.linalg.eigh(m)
    order = np.argsort(vals)[::-1]
    vals = vals[order]
    vecs = vecs[:, order]
    n = vecs[:, 2]
    if n[2] < 0:
        n = -n
    return n, vals, c


def planarity(vals):
    """Мера того, задают ли точки плоскость: (s2 - s3) / s2 в [0, 1].

    Ноль означает коллинеарность: точки лежат на прямой, и плоскость через
    них не определена. Единица - идеально плоское облако с размахом в двух
    направлениях.
    """
    if vals is None or vals[1] <= 0.0:
        return 0.0
    return float(max(0.0, (vals[1] - vals[2]) / vals[1]))


def attitude_from_normal(n):
    """Угол и азимут падения из нормали: (dip, dip_az) в градусах.

    Азимут падения это направление наискорейшего спуска по плоскости,
    отсчёт от севера по часовой стрелке. Работа в направляющих косинусах
    снимает переходы через ноль и через 90 градусов, где расчёты в углах
    обычно и сбиваются.
    """
    nx, ny, nz = float(n[0]), float(n[1]), float(n[2])
    nz = min(1.0, max(-1.0, nz))
    dip = math.degrees(math.acos(abs(nz)))
    h = math.hypot(nx, ny)
    if h < 1e-12:
        return 0.0, 0.0            # горизонтальная плоскость: азимут не нужен
    # У плоскости z = a*x + b*y + c верхняя нормаль пропорциональна
    # (-a, -b, 1), а спуск идёт в сторону -(a, b). Значит горизонтальная
    # часть нормали смотрит ПО падению, а не против: знак здесь
    # положительный, и это проверяется правилом трёх точек в тестах.
    dx, dy = nx / h, ny / h
    az = (math.degrees(math.atan2(dx, dy)) + 360.0) % 360.0
    return dip, az


def attitude_of_trace(pts, min_planarity=MIN_PLANARITY):
    """Залегание по трёхмерному следу: dict или причина отказа.

    Возвращает dict(dip, dip_az, planarity, rms, n_pts) либо
    dict(reason=...) - отказ с причиной, а не число из шума.
    rms - средняя невязка точек от подогнанной плоскости, в единицах
    координат: она отвечает на другой вопрос, чем обусловленность, и
    ловит складку, где точки не лежат на одной плоскости вовсе.
    """
    p = np.asarray(pts, dtype=float)
    if p.shape[0] < 3:
        return {"reason": "меньше трёх вершин с отметками"}
    if np.ptp(p[:, 2]) <= 0.0:
        return {"reason": "все отметки одинаковы, плоскость горизонтальна"}
    n, vals, c = fit_plane(p)
    if n is None:
        return {"reason": "меньше трёх вершин с отметками"}
    pl = planarity(vals)
    if pl < min_planarity:
        return {"reason": "след почти прямой в плане, плоскость не определена"}
    dip, az = attitude_from_normal(n)
    rms = float(np.sqrt(np.mean(((p - c) @ n) ** 2)))
    return {"dip": dip, "dip_az": az, "planarity": pl, "rms": rms,
            "n_pts": int(p.shape[0])}


def attitude_windows(pts, window, min_planarity=MIN_PLANARITY):
    """Залегание по скользящему окну вдоль следа.

    Длинная граница редко лежит на одной плоскости: складка, флексура,
    смещение разломом. Одно число по всему следу дало бы среднее, физически
    не существующее. Окно даёт изменение залегания вдоль границы, что
    полезнее.

    window - число вершин в окне (не меньше трёх). Возвращает список
    (индекс середины окна, результат attitude_of_trace).
    """
    p = np.asarray(pts, dtype=float)
    w = max(3, int(window))
    out = []
    if p.shape[0] < w:
        return out
    for i in range(0, p.shape[0] - w + 1):
        res = attitude_of_trace(p[i:i + w], min_planarity=min_planarity)
        out.append((i + w // 2, res))
    return out
