Aislo/B05_Profile/logic/B05_Profile_Sections_Ground.ts

/* =============================================================================
 * B05_Profile_Sections_Ground.ts
 * 노선 한 벌의 **지반선** — 종단 1 m 표본 · 측점 지반고 · 횡단 ±20 m 0.5 m · 미교차 덧붙임.
 * 노선을 고치는 동안 브라우저가 지형 격자(`common_util_surface_sampler.ts`)로 바로 셈한다(PLAN 344-1).
 *
 * ⚠ 파이썬·TS 짝 — `B05_Profile_Engine_Sections_Core.py` `generate_sections` 의 지반 부분
 * (종단 표본 · 측점 중심 · 횡단 표본 · `uphill_side`) ↔ `sectionsGround`.
 * ⚠ 파이썬·TS 짝 — `B05_Profile_Engine_Grade.py` `ground_profile` + 측점 `np.interp` ↔
 * `groundProfile` · `groundAt`.
 * ⚠ 파이썬·TS 짝 — `B06_Section_Engine_SampleExtend.py` `extend_unclosed` ↔ `extendUnclosed`.
 * 한쪽을 고치면 다른 쪽도 고치고 거울 시험 `resources/tester/test_344_1_sections_ground_mirror.py`.
 *
 * 측점 자리(규칙 · 스냅 · 구조물 측점 병합)는 부르는 쪽이 정해 넘김 — 틀은
 * `B05_Profile_UI_RouteEdit_Frame.ts`(이미 짝). 자릿수는 서버 출력과 같게 소수 6자리.
 * ========================================================================== */

import { type DtmGrid, heightAt } from "../../common_util/common_util_surface_sampler";
import type { CrossDesignEdge } from "../../common_util/common_util_cross_design";
import { stationFrameAt, stationsAlong } from "./B05_Profile_UI_RouteEdit_Frame";

type Vertex = [number, number];

/** 짝: `config_system_design.py` `SECTION_CROSS_HALF_WIDTH_M` · `_SAMPLE_INTERVAL_M` · `SECTION_LONG_…`. */
export const SECTION_CROSS_HALF_WIDTH_M = 20.0;
export const SECTION_CROSS_SAMPLE_INTERVAL_M = 0.5;
export const SECTION_LONG_SAMPLE_INTERVAL_M = 1.0;
/** 짝: `B06_Section_Engine_SampleExtend.EXTEND_STEP_M`. */
export const EXTEND_STEP_M = 5.0;
/** 짝: `B06_Section_Engine_Design._SLOPE_CLOSE_TOLERANCE_M`. */
const SLOPE_CLOSE_TOLERANCE_M = 0.01;

export interface GroundSample {
  x: number;
  y: number;
  z: number | null;
  elevation_m: number | null;
  valid: boolean;
}
export interface LongGroundSample extends GroundSample {
  chainage_m: number;
}
export interface CrossGroundPoint extends GroundSample {
  offset_m: number;
}

export interface StationGround {
  chainage_m: number;
  center_x: number;
  center_y: number;
  center_z: number | null;
  tangent_xy: Vertex;
  /** 좌측(+offset) 단위 벡터. */
  left_xy: Vertex;
  /** 횡단 기준 등고가 높은 쪽 — 좌우 유효 표본 평균 비교 · 모르면 null. */
  uphill_side: "left" | "right" | null;
  samples: CrossGroundPoint[];
}

export interface SectionsGround {
  length_m: number;
  /** 종단 지반 표본(1 m · 끝점 포함). */
  samples: LongGroundSample[];
  stations: StationGround[];
}

export interface GroundOptions {
  crossHalfWidthM?: number;
  crossSampleIntervalM?: number;
  longSampleIntervalM?: number;
}

/** 파이썬 `round(x, 6)` 자리. */
export const round6 = (value: number): number => Math.round(value * 1e6) / 1e6;
const round9 = (value: number): number => Math.round(value * 1e9) / 1e9;

function sample(grid: DtmGrid, x: number, y: number): GroundSample {
  const z = heightAt(grid, x, y);
  const valid = !Number.isNaN(z);
  const value = valid ? round6(z) : null;
  return { x: round6(x), y: round6(y), z: value, elevation_m: value, valid };
}

/** `_clean_polyline` — 앞 점과 길이 0(≤ 1e-8)인 점을 버리고 누가거리. */
function cleanLine(line: ReadonlyArray<Vertex>): {
  xs: number[];
  ys: number[];
  chainage: number[];
} {
  const kept = line.filter(
    (vertex, i) =>
      i === 0 || Math.hypot(vertex[0] - line[i - 1][0], vertex[1] - line[i - 1][1]) > 1e-8,
  );
  const chainage = [0];
  for (let i = 1; i < kept.length; i += 1) {
    chainage.push(
      chainage[i - 1] + Math.hypot(kept[i][0] - kept[i - 1][0], kept[i][1] - kept[i - 1][1]),
    );
  }
  return { xs: kept.map((v) => v[0]), ys: kept.map((v) => v[1]), chainage };
}

/** 횡단 offset 열 — `np.arange(-hw, hw + si/2, si)` · hw 넘는 것 버림 · 0 없으면 끼움. */
export function crossOffsets(halfWidth: number, interval: number): number[] {
  const count = Math.ceil((halfWidth + interval * 0.5 + halfWidth) / interval);
  const offsets: number[] = [];
  for (let i = 0; i < count; i += 1) {
    const value = -halfWidth + i * interval;
    if (value <= halfWidth + 1e-9) offsets.push(value);
  }
  if (!offsets.some((value) => Math.abs(value) <= 1e-9)) {
    offsets.push(0);
    offsets.sort((a, b) => a - b);
  }
  return offsets;
}

/** 종단 표본 자리 — 누가거리(간격 · 끝점 포함) · 자리 = `_interpolate_xy`(축마다 `np.interp`). */
export function longitudinalPoints(
  line: ReadonlyArray<Vertex>,
  intervalM: number = SECTION_LONG_SAMPLE_INTERVAL_M,
): { chainage: number[]; xy: Vertex[]; lengthM: number } {
  const { xs, ys, chainage: along } = cleanLine(line);
  const chainage = stationsAlong(line, intervalM, true);
  const xy = chainage.map((c): Vertex => [npInterp(c, along, xs), npInterp(c, along, ys)]);
  return { chainage, xy, lengthM: along[along.length - 1] ?? 0 };
}

/** 측점 자리 · 좌측 축 — 자리 = `_interpolate_xy` · 축 = 측점 틀 짝(`stationFrameAt`). 못 재면 뺌. */
export function placeStations(
  line: ReadonlyArray<Vertex>,
  stationsM: ReadonlyArray<number>,
  longIntervalM: number = SECTION_LONG_SAMPLE_INTERVAL_M,
): Array<{ chainage: number; origin: Vertex; left: Vertex }> {
  const { xs, ys, chainage: along } = cleanLine(line);
  const probe = Math.max(longIntervalM, 0.5);
  return stationsM.flatMap((chainage) => {
    const frame = stationFrameAt(line, chainage, probe);
    if (!frame) return [];
    const origin: Vertex = [npInterp(chainage, along, xs), npInterp(chainage, along, ys)];
    return [{ chainage, origin, left: frame.left }];
  });
}

/** 노선 지반선 — 종단 표본 · 측점마다 중심 · 좌측 축 · 횡단 표본. `stationsM` = 측점 누가거리. */
export function sectionsGround(
  line: ReadonlyArray<Vertex>,
  stationsM: ReadonlyArray<number>,
  grid: DtmGrid,
  options: GroundOptions = {},
): SectionsGround {
  const halfWidth = options.crossHalfWidthM ?? SECTION_CROSS_HALF_WIDTH_M;
  const interval = options.crossSampleIntervalM ?? SECTION_CROSS_SAMPLE_INTERVAL_M;
  const longInterval = options.longSampleIntervalM ?? SECTION_LONG_SAMPLE_INTERVAL_M;

  const long = longitudinalPoints(line, longInterval);
  const samples: LongGroundSample[] = long.chainage.map((chainage, i) => ({
    chainage_m: round6(chainage),
    ...sample(grid, long.xy[i][0], long.xy[i][1]),
  }));

  const offsets = crossOffsets(halfWidth, interval);
  let center = 0;
  offsets.forEach((value, index) => {
    if (Math.abs(value) < Math.abs(offsets[center])) center = index;
  });

  const stations: StationGround[] = [];
  for (const { chainage, origin, left } of placeStations(line, stationsM, longInterval)) {
    const [ox, oy] = origin;
    const [lx, ly] = left;
    const cross: CrossGroundPoint[] = [];
    const raw: number[] = [];
    for (const offset of offsets) {
      const x = ox + lx * offset;
      const y = oy + ly * offset;
      raw.push(heightAt(grid, x, y));
      cross.push({ offset_m: round6(offset), ...sample(grid, x, y) });
    }
    const mean = (pick: (offset: number) => boolean): number | null => {
      const values = raw.filter((z, i) => pick(offsets[i]) && !Number.isNaN(z));
      return values.length ? values.reduce((a, b) => a + b, 0) / values.length : null;
    };
    const leftMean = mean((offset) => offset > 1e-9);
    const rightMean = mean((offset) => offset < -1e-9);
    let uphill: StationGround["uphill_side"] = null;
    if (leftMean !== null && rightMean !== null && Math.abs(leftMean - rightMean) > 1e-6) {
      uphill = leftMean > rightMean ? "left" : "right";
    }
    const centerZ = raw[center];
    stations.push({
      chainage_m: round6(chainage),
      center_x: round6(ox),
      center_y: round6(oy),
      center_z: Number.isNaN(centerZ) ? null : round6(centerZ),
      tangent_xy: [round9(ly), round9(-lx)],
      left_xy: [round9(lx), round9(ly)],
      uphill_side: uphill,
      samples: cross,
    });
  }
  return { length_m: round6(long.lengthM), samples, stations };
}

/** `np.interp(x, xp, fp)` — xp 오름차순 · 범위 밖은 끝값 · 식 차례도 numpy 와 같음. */
function npInterp(x: number, xp: ReadonlyArray<number>, fp: ReadonlyArray<number>): number {
  const last = xp.length - 1;
  if (x <= xp[0]) return fp[0];
  if (x >= xp[last]) return fp[last];
  let low = 0;
  let high = last;
  while (high - low > 1) {
    const mid = (low + high) >> 1;
    if (xp[mid] <= x) low = mid;
    else high = mid;
  }
  const slope = (fp[low + 1] - fp[low]) / (xp[low + 1] - xp[low]);
  return slope * (x - xp[low]) + fp[low];
}

/** 종단 표본 → (누가거리, 지반고) · 결측은 이웃 유효 값으로 선형 보간. 유효 둘 미만이면 throw. */
export function groundProfile(samples: ReadonlyArray<Partial<LongGroundSample>>): {
  chainage: number[];
  ground: number[];
} {
  const rows = samples
    .map((s) => ({
      chainage: Number(s.chainage_m ?? 0),
      z: s.valid !== false && s.elevation_m != null ? Number(s.elevation_m) : NaN,
    }))
    .sort((a, b) => a.chainage - b.chainage);
  const known = rows.filter((row) => Number.isFinite(row.z));
  if (known.length < 2 || rows.length < 2) {
    throw new Error("계획선을 만들기에 유효한 종단 표고 샘플이 부족합니다.");
  }
  const xp = known.map((row) => row.chainage);
  const fp = known.map((row) => row.z);
  return {
    chainage: rows.map((row) => row.chainage),
    ground: rows.map((row) => npInterp(row.chainage, xp, fp)),
  };
}

/** 측점 지반고 — `groundProfile` 위 선형 보간(파이썬 `station_ground`). */
export function groundAt(profile: { chainage: number[]; ground: number[] }, chainageM: number) {
  return npInterp(chainageM, profile.chainage, profile.ground);
}

/** 횡단 설계가 덧붙임에 돌려줘야 할 칸. */
export interface UnclosedDesign {
  slope_unclosed: boolean;
  design_line?: ReadonlyArray<CrossDesignEdge> | null;
}

export interface ExtendInfo {
  left_m: number;
  right_m: number;
  surface_end: boolean;
}

function validSorted(samples: ReadonlyArray<Partial<CrossGroundPoint>>): Array<[number, number]> {
  return samples
    .filter((s) => s.valid !== false && s.elevation_m != null)
    .map((s): [number, number] => [Number(s.offset_m), Number(s.elevation_m)])
    .sort((a, b) => a[0] - b[0] || a[1] - b[1]);
}

/** 미교차 측점 지반 표본을 **열린 쪽만** +`stepM` 씩 지표면 끝까지 넓힘 — 닫혔으면 null.
 *  새 걸음에 무효 표본이 섞이면 그 쪽은 안 붙이고 `surface_end`. 상한 없음(파이썬과 같음). */
export function extendUnclosed<T extends Partial<CrossGroundPoint>>(
  samples: ReadonlyArray<T>,
  origin: Vertex,
  leftAxis: Vertex,
  compute: (samples: Array<T | CrossGroundPoint>) => UnclosedDesign,
  grid: DtmGrid,
  stepM: number = EXTEND_STEP_M,
): { samples: Array<T | CrossGroundPoint>; info: ExtendInfo } | null {
  let current: Array<T | CrossGroundPoint> = [...samples];
  let design = compute(current);
  if (!design.slope_unclosed) return null;
  let offsets = validSorted(current).map(([offset]) => offset);
  let spacing = Infinity;
  offsets.slice(1).forEach((value, index) => {
    const gap = value - offsets[index];
    if (gap > 1e-9 && gap < spacing) spacing = gap;
  });
  if (!Number.isFinite(spacing)) throw new Error("넓힐 표본 간격을 잴 수 없습니다.");
  const count = Math.max(Math.round(stepM / spacing), 1);
  const info: ExtendInfo = { left_m: 0, right_m: 0, surface_end: false };
  while (design.slope_unclosed) {
    const ground = validSorted(current);
    const lineEnds = (design.design_line ?? [])
      .map((p): [number, number] => [Number(p.offset_m), Number(p.elevation_m)])
      .sort((a, b) => a[0] - b[0] || a[1] - b[1]);
    const open =
      ground.length && lineEnds.length
        ? {
            right: Math.abs(ground[0][1] - lineEnds[0][1]) > SLOPE_CLOSE_TOLERANCE_M,
            left:
              Math.abs(ground[ground.length - 1][1] - lineEnds[lineEnds.length - 1][1]) >
              SLOPE_CLOSE_TOLERANCE_M,
          }
        : { right: false, left: false };
    const added: CrossGroundPoint[] = [];
    for (const [side, sign] of [
      ["right", -1],
      ["left", 1],
    ] as const) {
      if (!open[side]) continue;
      const edge = sign < 0 ? Math.min(...offsets) : Math.max(...offsets);
      const step: CrossGroundPoint[] = [];
      for (let k = 1; k <= count; k += 1) {
        const offset = edge + sign * spacing * k;
        const point = sample(
          grid,
          origin[0] + leftAxis[0] * offset,
          origin[1] + leftAxis[1] * offset,
        );
        step.push({ offset_m: round6(offset), ...point });
      }
      if (!step.every((point) => point.valid)) {
        info.surface_end = true;
        continue;
      }
      added.push(...step);
      info[`${side}_m`] += stepM;
    }
    if (!added.length) break;
    current = [...current, ...added].sort((a, b) => Number(a.offset_m) - Number(b.offset_m));
    offsets = validSorted(current).map(([offset]) => offset);
    design = compute(current);
  }
  return { samples: current, info };
}