fis_lib_measure/lib/process/calcuators/formulas/general.dart

// ignore_for_file:  constant_identifier_names

import 'dart:math' as math;

import 'package:fis_measure/configs/patient.dart';
import 'package:fis_measure/interfaces/date_types/point.dart';
import 'package:fis_measure/interfaces/enums/species.dart';

class GeneralFormulas {
  GeneralFormulas._();
  static const double VolumeCofficient = math.pi / 6.0;

  static final IGeneralFormulaStrategy _singleton =
      GlobalPatientConfig.speciesType == SpeciesType.mouse
          ? BaseGeneralFormulas()
          : AnimalsGeneralFormulas();

  static double volumeTwoLine(double d1, double d2,
          [bool useBigValue = true]) =>
      _singleton.volumeTwoLine(d1, d2, useBigValue);

  static double volumeWithCoefficient(
          double d1, double d2, double d3, double coefficient) =>
      _singleton.volumeWithCoefficient(d1, d2, d3, coefficient);

  static double countStenosis(double d1, double d2) =>
      _singleton.countStenosis(d1, d2);

  static bool doubleAlmostEquals(double num1, double num2,
          [double precision = 0.000001]) =>
      _singleton.doubleAlmostEquals(num1, num2, precision);

  static double distance2Line(DPoint pt1, DPoint pt2, DPoint pt) =>
      _singleton.distance2Line(pt1, pt2, pt);

  static double countSlope(DPoint p1, DPoint p2) =>
      _singleton.countSlope(p1, p2);

  static double countEnvelopeTime(List<DPoint> points) =>
      _singleton.countEnvelopeTime(points);

  static double countPressure(double v) => _singleton.countPressure(v);

  static double countRI(double ps, double mdEd) => _singleton.countRI(ps, mdEd);

  static double countMaxPG(double v1, double v2) =>
      _singleton.countMaxPG(v1, v2);

  static double countHR(double timeSpan) => _singleton.countHR(timeSpan);

  static List<double> countVTI(List<DPoint> tracePoints) =>
      _singleton.countVTI(tracePoints);

  static double svDiam(double diam, double vti) => _singleton.svDiam(diam, vti);

  static double csa(double diam) => _singleton.csa(diam);

  static double dpDtRatio(DPoint p1, DPoint p2) => _singleton.dpDtRatio(p1, p2);

  static double maxPG(double v1, double v2) => _singleton.maxPG(v1, v2);

  static double area(double d1, double d2) => _singleton.area(d1, d2);
  static double perimeter(double d1, double d2) => _singleton.perimeter(d1, d2);

  static double flowVolTAMAX(
          double flowArea, double tamean, double coefficient) =>
      _singleton.flowVolTAMAX(flowArea, tamean, coefficient);

  static double pi(double ps, double mdEd, double tamax) =>
      _singleton.pi(ps, mdEd, tamax);

  static double medianVelocity(double ps, double ed) =>
      _singleton.medianVelocity(ps, ed);
  static double volume(double d) => _singleton.volume(d);
}

abstract class IGeneralFormulaStrategy {
  double volumeTwoLine(double d1, double d2, [bool useBigValue = true]);
  double volumeWithCoefficient(
    double d1,
    double d2,
    double d3,
    double coefficient,
  );
  double countStenosis(double d1, double d2);

  bool doubleAlmostEquals(double num1, double num2,
      [double precision = 0.000001]);
  double distance2Line(DPoint pt1, DPoint pt2, DPoint pt);
  double countSlope(DPoint p1, DPoint p2);
  double countEnvelopeTime(List<DPoint> points);
  double countPressure(double v);
  double countRI(double ps, double mdEd);
  double countMaxPG(double v1, double v2);
  double countHR(double timeSpan);
  List<double> countVTI(
    List<DPoint> tracePoints,
  );
  double svDiam(double diam, double vti);
  double csa(double diam);
  double dpDtRatio(DPoint p1, DPoint p2);
  double maxPG(double v1, double v2);
  double area(double d1, double d2);
  double flowVolTAMAX(double flowArea, double tamean, double coefficient);
  double perimeter(double d1, double d2);
  double volume(double d);

  /// <summary>
  /// Pulastility Index
  /// <para>PI (ED):(PS - ED)/TAMAX</para>
  /// <para>PI (MD):(PS - MD)/TAMAX</para>
  /// </summary>
  /// <param name="ps">Unit cm/s</param>
  /// <param name="md_ed">Unit cm/s</param>
  /// <param name="tamax">Unit cm/s</param>
  /// <returns>Unit none</returns>
  /// 在TCCD下用Vmean替换tamax计算
  double pi(double ps, double mdEd, double tamax);

  /// <param name="ps">Unit cm/s</param>
  /// <param name="ed">Unit cm/s</param>
  /// <returns>Unit cm/s</returns>
  double medianVelocity(double ps, double ed);
}

class BaseGeneralFormulas implements IGeneralFormulaStrategy {
  /// Volume:1/6 x π
  static const double VolumeCofficient = math.pi / 6.0;

  /// 1/6 x π x D1^2 x D2 (D1 ＞ D2)
  ///
  /// 1/6 x π x D1 x D2^2 (D1 ≤ D2)
  ///
  /// return unit: cm³</returns>
  @override
  double volumeTwoLine(double d1, double d2, [bool useBigValue = true]) {
    double volume = 0;
    if (!doubleAlmostEquals(d1, 0) && !doubleAlmostEquals(d2, 0)) {
      if (useBigValue) {
        if (d1 > d2) {
          volume = d1 * d1 * d2 * math.pi / 6.0;
        } else {
          volume = d1 * d2 * d2 * math.pi / 6.0;
        }
      } else {
        if (d1 > d2) {
          volume = d2 * d2 * d1 * math.pi / 6.0;
        } else {
          volume = d1 * d1 * d2 * math.pi / 6.0;
        }
      }
    }
    return volume;
  }

  /// 带比率系数算体积
  ///
  /// [coefficient] 比率系数
  @override
  double volumeWithCoefficient(
    double d1,
    double d2,
    double d3,
    double coefficient,
  ) {
    double volume = 0;
    if (!doubleAlmostEquals(d1, 0) &&
        !doubleAlmostEquals(d2, 0) &&
        !doubleAlmostEquals(d3, 0)) {
      volume = d1 * d2 * d3 * coefficient;
    }
    return volume;
  }

  /// 狭窄率计算
  @override
  double countStenosis(double d1, double d2) {
    double residual;
    double lumen;
    if (d1 > d2) {
      residual = d2;
      lumen = d1;
    } else {
      residual = d1;
      lumen = d2;
    }
    return (1.0 - residual / lumen) * 100;
  }

  @override
  bool doubleAlmostEquals(double num1, double num2,
      [double precision = 0.000001]) {
    if (num1.isNaN && num2.isNaN) return true;
    return (num1 - num2).abs() <= precision;
  }

  /// 点到线的距离
  ///
  /// [pt1] 线端点1
  ///
  /// [pt2] 线端点2
  ///
  /// [pt] 点
  @override
  double distance2Line(DPoint pt1, DPoint pt2, DPoint pt) {
    double dis = 0;
    if (pt1.x == pt2.x) {
      dis = (pt.x - pt1.x).abs();
      return dis;
    }
    var lineK = (pt2.y - pt1.y) / (pt2.x - pt1.x);
    var lineC = (pt2.x * pt1.y - pt1.x * pt2.y) / (pt2.x - pt1.x);
    dis = (lineK * pt.x - pt.y + lineC).abs() / (math.sqrt(lineK * lineK + 1));
    return dis;
  }

  ///计算斜率
  @override
  double countSlope(DPoint p1, DPoint p2) {
    if (doubleAlmostEquals(p2.x, p1.x)) return 0;
    double vertical = p2.y - p1.y;
    double time = (p2.x - p1.x).abs();
    double slope = vertical / time;
    return slope;
  }

  ///计算包络时间
  ///Origin: Vinno.Modules.MeasureModule\Formulas\GeneralFormulas.cs [632:67]
  @override
  double countEnvelopeTime(List<DPoint> points) {
    if (points.length < 2) return 0;
    double startTime = double.maxFinite;
    double endTime = -double.maxFinite;
    for (var point in points) {
      if (point.x < startTime) {
        startTime = point.x;
      }
      if (point.x > endTime) {
        endTime = point.x;
      }
    }
    return (endTime - startTime).abs();
  }

  /// <summary>
  /// 4.0 x V^2
  /// </summary>
  /// <param name="v">Unit cm/s</param>
  /// <returns>Unit mmHg</returns>
  @override
  double countPressure(double v) {
    // The velocity is in cm/s, but it should be m/s in this formula, so divide it by 100.
    double pg = 4.0 * math.pow(v * 0.01, 2);
    return pg;
  }

  /// <summary>
  /// Resistivity Index 阻力指数
  /// <para>RI (ED)=(PS - ED)/PS  </para>
  /// <para>RI (MD)=(PS - MD)/PS  </para>
  /// </summary>
  /// <param name="ps">Unit cm/s</param>
  /// <param name="mdEd">Unit cm/s</param>
  @override
  double countRI(double ps, double mdEd) {
    ps = ps.abs();
    mdEd = mdEd.abs();
    double ri = (ps - mdEd) / ps;
    if (ri < 0) {
      return double.nan;
    }
    return ri;
  }

  /// <summary>
  /// 4.0 x |V1^2 - V2^2|
  /// </summary>
  /// <param name="v1">Unit cm/s</param>
  /// <param name="v2">Unit cm/s</param>
  /// <returns>Unit mmHg</returns>
  @override
  double countMaxPG(double v1, double v2) {
    // The velocity is in cm/s, but it should be m/s in this formula, so divide it by 100.
    double meanPG =
        4.0 * (math.pow(v1 * 0.01, 2) - math.pow(v2 * 0.01, 2)).abs();
    return meanPG;
  }

  // 计算心率
  // TODO: 心脏周期（可配置） DefaultHeartCycle = 2;
  // Origin: Vinno.Modules.MeasureModule\Primitives\DopplerTraceBase.cs
  @override
  double countHR(double timeSpan) {
    const defaultHeartCycle = 1;
    double hr = (defaultHeartCycle * 60 / timeSpan).roundToDouble();
    return hr;
  }

  // return [VTI,tiEnv,timeAveragedVelocity,pv,meanPG,hr];
  @override
  List<double> countVTI(
    List<DPoint> tracePoints,
  ) {
    var tiEnv = 0.0;
    var timeAveragedVelocity = 0.0;
    var pv = 0.0;
    var meanPG = 0.0;
    if (tracePoints.length < 2) return [0, 0, 0, 0, 0, 0];
    int n = tracePoints.length - 1;
    double dis = 0;
    double startTime = 0;
    double endTime = 0;
    double pgSum = 0;
    int positiveCount = 0;
    int negativeCount = 0;
    for (int i = 0; i < n; i++) {
      var p1 = tracePoints[i];
      var p2 = tracePoints[i + 1];
      if (pv.abs() < (p1.y).abs()) {
        pv = p1.y;
      }
      double meanVelocity = (p1.y + p2.y).abs() / 2;
      double timeSpan = (p1.x - p2.x).abs();
      dis += meanVelocity * timeSpan;
      pgSum += countPressure(meanVelocity) * timeSpan;

      if (i == 0) {
        startTime = p1.x;
        endTime = p2.x;
      } else {
        startTime = math.min(startTime, p2.x);
        endTime = math.max(endTime, p2.x);
      }

      if (p1.y < 0) {
        negativeCount++;
      } else {
        positiveCount++;
      }
    }
    if (tracePoints[n].y < 0) {
      negativeCount++;
    } else {
      positiveCount++;
    }
    if (positiveCount < negativeCount) {
      dis = -dis;
    }

    if (pv.abs() < (tracePoints[n].y).abs()) {
      pv = tracePoints[n].y;
    }

    tiEnv = (endTime - startTime).abs();
    meanPG = (pgSum / tiEnv);
    timeAveragedVelocity = (dis / tiEnv);
    var hr = countHR(tiEnv);

    return [dis, tiEnv, timeAveragedVelocity, pv, meanPG, hr];
  }

  /// SV Diam
  ///
  /// SV = 1/4 x π x (Diam)^2 x VTI
  ///
  /// [diam] `cm`
  ///
  /// [vti] `cm`
  ///
  /// return `cm3`
  @override
  double svDiam(double diam, double vti) {
    double sv = math.pi * math.pow(diam, 2) * vti / 4;
    return sv;
  }

  /// CSA
  ///
  /// CSA = 1/4 x π x Diam^2
  ///
  /// [diam] `cm`
  ///
  /// return `cm²`
  @override
  double csa(double diam) {
    double csa = math.pi * math.pow(diam, 2) / 4;
    return csa;
  }

  /// CW/PW Mode: dp/dt = (4* V2 * V2 - 4 * V1 * V1)/ (T2 - T1)
  ///
  /// return `mmHg/s`
  @override
  double dpDtRatio(DPoint p1, DPoint p2) {
    double dp = maxPG(p2.y, p1.y);
    double dt = (p2.x - p1.x).abs();
    double result = dp / dt;
    return result;
  }

  /// 4.0 x |V1^2 - V2^2|
  ///
  /// [v1] `cm/s`
  ///
  /// [v2] `cm/s`
  ///
  /// reutrn `mmHg`
  @override
  double maxPG(double v1, double v2) {
    // The velocity is in cm/s, but it should be m/s in this formula, so divide it by 100.
    double meanPG =
        4.0 * (math.pow(v1 * 0.01, 2) - math.pow(v2 * 0.01, 2)).abs();
    return meanPG;
  }

  /// 0.25 x π x D1 x D2
  ///
  /// [d1] `cm`
  ///
  /// [d2] `cm`
  ///
  /// return `cm²`
  @override
  double area(double d1, double d2) {
    double area = 0;
    if (!doubleAlmostEquals(d1, 0) && !doubleAlmostEquals(d2, 0)) {
      if ((d1 < 0.001 * d2) || (d2 < 0.001 * d1)) {
        // Stright line
        area = 0.0;
      } else {
        // Cirle
        area = 0.25 * d1 * d2 * math.pi;
      }
    }
    return area;
  }

  /// 0.25 x π x Diam^2 x TAMAX x Compensation Coefficient, where 0.5 ≤ Compensation Coefficient ≤ 1.0
  ///
  /// [flowArea] `cm²`
  ///
  /// [tamean] `cm/s`
  ///
  /// [coefficient] `[0.5,1.0]`
  ///
  /// return `cm³/s`
  @override
  double flowVolTAMAX(double flowArea, double tamean, double coefficient) {
    if (coefficient < 0.5 || coefficient > 1.0) {
      return double.nan;
    }
    double flowVol = flowArea * tamean.abs() * coefficient;
    return flowVol;
  }

  @override
  double pi(double ps, double mdEd, double tamax) {
    double pi = (ps - mdEd) / tamax;
    return pi.abs();
  }

  @override
  double medianVelocity(double ps, double ed) {
    ps = ps.abs();
    ed = ed.abs();
    return (ps + 2 * ed) / 3;
  }

  @override
  double perimeter(double d1, double d2) {
    double perimeter = 0.0;
    if (!almostEquals(d1, 0.0) && !almostEquals(d2, 0.0)) {
      if (almostEquals(d1, d2)) {
        // Circle
        perimeter = d1 * math.pi;
      } else if (d1 < 0.001 * d2) {
        // Straight line
        perimeter = 2.0 * d2;
      } else if (d2 < 0.001 * d1) {
        // Straight line
        perimeter = 2.0 * d1;
      } else {
        // Ellipse
        double a = d1 > d2 ? d1 / 2.0 : d2 / 2.0; // longest half axis
        double b = d1 < d2 ? d1 / 2.0 : d2 / 2.0; // shortest half axis

        if (a < 0.0001) {
          // Long axis too short to make computations. Return zero.
          perimeter = 0.0;
        } else {
          // Calculate k, kc
          double k = math.sqrt(a * a - b * b) / a;
          double kc2 = 1 - k * k;
          perimeter = 4.0 * a * ellipseCircum(math.sqrt(kc2), kc2);
        }
      }
    }
    return perimeter;
  }

  bool almostEquals(double a, double b, [double epsilon = 1e-10]) {
    return (a - b).abs() < epsilon;
  }

  double ellipseCircum(double k, double kc2) {
    // Placeholder for the EllipseCircum method equivalent.
    // You can implement the complete elliptic integral of the second kind here
    // or use an approximation method for ellipse circumference.
    // For the sake of this example, we will use Ramanujan's approximation.
    double h = math.pow((k - 1), 2) / math.pow((k + 1), 2);
    return math.pi * (1 + 3 * h / (10 + math.sqrt(4 - 3 * h)));
  }

  @override
  double volume(double d) {
    double volume = 0.0;
    if (!almostEquals(d, 0.0)) {
      volume = math.pow(d, 3) * math.pi / 6.0;
    }
    return volume;
  }
}

class AnimalsGeneralFormulas extends BaseGeneralFormulas {}