#include "nmPebiRegionPressureCalculator.h"

#include "nmPebiResultSnapshot.h"

#include <QRectF>

#include <vtkCell.h>
#include <vtkCellType.h>
#include <vtkDoubleArray.h>
#include <vtkPoints.h>
#include <vtkUnstructuredGrid.h>

#include <algorithm>
#include <float.h>
#include <math.h>

#if !defined(_MSC_VER)
#include <cmath>
#endif

namespace {

// 所有区域统计几何均投影到 XY 平面，三角形是剪裁和求交的最小单元。
struct Triangle2D
{
    QPointF m_a;
    QPointF m_b;
    QPointF m_c;
};

struct PreparedPolygon
{
    QVector<QPointF> m_vecPoints;
    QVector<Triangle2D> m_vecTriangles;
    QRectF m_oBounds;
    double m_dArea;
};

struct CachedCell
{
    CachedCell()
        : m_nCellId(-1),
          m_dArea(0.0),
          m_dPoreVolumeWeight(0.0),
          m_bWeightValid(false)
    {
    }

    vtkIdType m_nCellId;
    double m_dArea;
    double m_dPoreVolumeWeight;
    bool m_bWeightValid;
};

bool isFiniteValue(double dValue)
{
#if defined(_MSC_VER)
    return _finite(dValue) != 0;
#else
    return std::isfinite(dValue);
#endif
}

double crossProduct(const QPointF& a,
                    const QPointF& b,
                    const QPointF& c)
{
    return (b.x() - a.x()) * (c.y() - a.y()) -
           (b.y() - a.y()) * (c.x() - a.x());
}

double squaredDistance(const QPointF& a, const QPointF& b)
{
    const double dX = a.x() - b.x();
    const double dY = a.y() - b.y();
    return dX * dX + dY * dY;
}

double signedDoubleArea(const QVector<QPointF>& vecPoints)
{
    if(vecPoints.size() < 3) {
        return 0.0;
    }

    // 先平移到首点附近再做鞋带求和，避免大地坐标中的大数相减吞掉小区域面积。
    const QPointF oOrigin = vecPoints.first();
    double dArea = 0.0;
    for(int nIndex = 0; nIndex < vecPoints.size(); ++nIndex) {
        const QPointF& a = vecPoints[nIndex];
        const QPointF& b = vecPoints[(nIndex + 1) % vecPoints.size()];
        const double dAx = a.x() - oOrigin.x();
        const double dAy = a.y() - oOrigin.y();
        const double dBx = b.x() - oOrigin.x();
        const double dBy = b.y() - oOrigin.y();
        dArea += dAx * dBy - dBx * dAy;
    }
    return dArea;
}

QRectF polygonBounds(const QVector<QPointF>& vecPoints)
{
    double dMinimumX = vecPoints[0].x();
    double dMaximumX = dMinimumX;
    double dMinimumY = vecPoints[0].y();
    double dMaximumY = dMinimumY;
    for(int nIndex = 1; nIndex < vecPoints.size(); ++nIndex) {
        dMinimumX = qMin(dMinimumX, vecPoints[nIndex].x());
        dMaximumX = qMax(dMaximumX, vecPoints[nIndex].x());
        dMinimumY = qMin(dMinimumY, vecPoints[nIndex].y());
        dMaximumY = qMax(dMaximumY, vecPoints[nIndex].y());
    }
    return QRectF(QPointF(dMinimumX, dMinimumY),
                  QPointF(dMaximumX, dMaximumY));
}

bool pointOnSegment(const QPointF& p,
                    const QPointF& a,
                    const QPointF& b,
                    double dLinearTolerance,
                    double dAreaTolerance)
{
    if(qAbs(crossProduct(a, b, p)) > dAreaTolerance) {
        return false;
    }
    return p.x() >= qMin(a.x(), b.x()) - dLinearTolerance &&
           p.x() <= qMax(a.x(), b.x()) + dLinearTolerance &&
           p.y() >= qMin(a.y(), b.y()) - dLinearTolerance &&
           p.y() <= qMax(a.y(), b.y()) + dLinearTolerance;
}

int orientationSign(const QPointF& a,
                    const QPointF& b,
                    const QPointF& c,
                    double dAreaTolerance)
{
    const double dCross = crossProduct(a, b, c);
    if(dCross > dAreaTolerance) {
        return 1;
    }
    if(dCross < -dAreaTolerance) {
        return -1;
    }
    return 0;
}

bool segmentsIntersect(const QPointF& a,
                       const QPointF& b,
                       const QPointF& c,
                       const QPointF& d,
                       double dLinearTolerance,
                       double dAreaTolerance)
{
    const int nAbC = orientationSign(a, b, c, dAreaTolerance);
    const int nAbD = orientationSign(a, b, d, dAreaTolerance);
    const int nCdA = orientationSign(c, d, a, dAreaTolerance);
    const int nCdB = orientationSign(c, d, b, dAreaTolerance);
    if(nAbC * nAbD < 0 && nCdA * nCdB < 0) {
        return true;
    }
    return (nAbC == 0 && pointOnSegment(c, a, b,
                                        dLinearTolerance, dAreaTolerance)) ||
           (nAbD == 0 && pointOnSegment(d, a, b,
                                        dLinearTolerance, dAreaTolerance)) ||
           (nCdA == 0 && pointOnSegment(a, c, d,
                                        dLinearTolerance, dAreaTolerance)) ||
           (nCdB == 0 && pointOnSegment(b, c, d,
                                        dLinearTolerance, dAreaTolerance));
}

bool isSimplePolygon(const QVector<QPointF>& vecPoints,
                     double dLinearTolerance,
                     double dAreaTolerance)
{
    // 只检查非相邻边；相邻边共享端点是合法多边形的必要条件。
    const int nPointCount = vecPoints.size();
    for(int nFirst = 0; nFirst < nPointCount; ++nFirst) {
        const int nFirstNext = (nFirst + 1) % nPointCount;
        for(int nSecond = nFirst + 1; nSecond < nPointCount; ++nSecond) {
            const int nSecondNext = (nSecond + 1) % nPointCount;
            if(nFirst == nSecond || nFirstNext == nSecond ||
               nSecondNext == nFirst) {
                continue;
            }
            if(segmentsIntersect(vecPoints[nFirst], vecPoints[nFirstNext],
                                 vecPoints[nSecond], vecPoints[nSecondNext],
                                 dLinearTolerance, dAreaTolerance)) {
                return false;
            }
        }
    }
    return true;
}

bool pointInsideTriangle(const QPointF& p,
                         const QPointF& a,
                         const QPointF& b,
                         const QPointF& c,
                         double dAreaTolerance)
{
    return crossProduct(a, b, p) >= -dAreaTolerance &&
           crossProduct(b, c, p) >= -dAreaTolerance &&
           crossProduct(c, a, p) >= -dAreaTolerance;
}

bool triangulatePolygon(const QVector<QPointF>& vecPoints,
                        double dAreaTolerance,
                        QVector<Triangle2D>& vecTriangles)
{
    // 耳切法同时支持凸、凹简单多边形，无需依赖 VTK 的三维三角化管线。
    vecTriangles.clear();
    QVector<int> vecIndices;
    for(int nIndex = 0; nIndex < vecPoints.size(); ++nIndex) {
        vecIndices.append(nIndex);
    }

    while(vecIndices.size() > 3) {
        bool bEarFound = false;
        for(int nPosition = 0; nPosition < vecIndices.size(); ++nPosition) {
            const int nPrevious = vecIndices[
                    (nPosition + vecIndices.size() - 1) % vecIndices.size()];
            const int nCurrent = vecIndices[nPosition];
            const int nNext = vecIndices[(nPosition + 1) % vecIndices.size()];
            const QPointF& a = vecPoints[nPrevious];
            const QPointF& b = vecPoints[nCurrent];
            const QPointF& c = vecPoints[nNext];
            if(crossProduct(a, b, c) <= dAreaTolerance) {
                continue;
            }

            bool bContainsOtherPoint = false;
            for(int nOtherPosition = 0;
                nOtherPosition < vecIndices.size(); ++nOtherPosition) {
                const int nOther = vecIndices[nOtherPosition];
                if(nOther == nPrevious || nOther == nCurrent ||
                   nOther == nNext) {
                    continue;
                }
                if(pointInsideTriangle(vecPoints[nOther], a, b, c,
                                       dAreaTolerance)) {
                    bContainsOtherPoint = true;
                    break;
                }
            }
            if(bContainsOtherPoint) {
                continue;
            }

            Triangle2D oTriangle;
            oTriangle.m_a = a;
            oTriangle.m_b = b;
            oTriangle.m_c = c;
            vecTriangles.append(oTriangle);
            vecIndices.remove(nPosition);
            bEarFound = true;
            break;
        }
        if(!bEarFound) {
            vecTriangles.clear();
            return false;
        }
    }

    if(vecIndices.size() != 3) {
        vecTriangles.clear();
        return false;
    }
    Triangle2D oTriangle;
    oTriangle.m_a = vecPoints[vecIndices[0]];
    oTriangle.m_b = vecPoints[vecIndices[1]];
    oTriangle.m_c = vecPoints[vecIndices[2]];
    if(crossProduct(oTriangle.m_a, oTriangle.m_b, oTriangle.m_c) <=
       dAreaTolerance) {
        vecTriangles.clear();
        return false;
    }
    vecTriangles.append(oTriangle);
    return true;
}

bool preparePolygon(const QVector<QPointF>& vecInput,
                    PreparedPolygon& oPolygon)
{
    oPolygon = PreparedPolygon();
    if(vecInput.size() < 3) {
        return false;
    }
    for(int nIndex = 0; nIndex < vecInput.size(); ++nIndex) {
        if(!isFiniteValue(vecInput[nIndex].x()) ||
           !isFiniteValue(vecInput[nIndex].y())) {
            return false;
        }
    }

    // 容差跟随当前多边形尺度，兼容小尺寸网格和大地坐标。
    const QRectF oInputBounds = polygonBounds(vecInput);
    const double dScale = qMax(1.0, qMax(oInputBounds.width(),
                                        oInputBounds.height()));
    const double dLinearTolerance = dScale * 1.0e-10;
    const double dAreaTolerance = dScale * dScale * 1.0e-12;
    const double dSquaredTolerance =
            dLinearTolerance * dLinearTolerance;

    // 先去掉相邻重复点及显式闭合点，避免产生零长边。
    QVector<QPointF> vecPoints;
    for(int nIndex = 0; nIndex < vecInput.size(); ++nIndex) {
        if(vecPoints.isEmpty() ||
           squaredDistance(vecPoints.last(), vecInput[nIndex]) >
                dSquaredTolerance) {
            vecPoints.append(vecInput[nIndex]);
        }
    }
    if(vecPoints.size() > 1 &&
       squaredDistance(vecPoints.first(), vecPoints.last()) <=
            dSquaredTolerance) {
        vecPoints.remove(vecPoints.size() - 1);
    }

    // 反复清理共线中间点，直到点集稳定，降低耳切法退化风险。
    bool bRemovedPoint = true;
    while(bRemovedPoint && vecPoints.size() >= 3) {
        bRemovedPoint = false;
        for(int nIndex = 0; nIndex < vecPoints.size(); ++nIndex) {
            const QPointF& a = vecPoints[
                    (nIndex + vecPoints.size() - 1) % vecPoints.size()];
            const QPointF& p = vecPoints[nIndex];
            const QPointF& b = vecPoints[(nIndex + 1) % vecPoints.size()];
            if(pointOnSegment(p, a, b,
                              dLinearTolerance, dAreaTolerance)) {
                vecPoints.remove(nIndex);
                bRemovedPoint = true;
                break;
            }
        }
    }
    if(vecPoints.size() < 3 ||
       !isSimplePolygon(vecPoints,
                        dLinearTolerance, dAreaTolerance)) {
        return false;
    }

    double dDoubleArea = signedDoubleArea(vecPoints);
    if(!isFiniteValue(dDoubleArea) ||
       qAbs(dDoubleArea) <= dAreaTolerance) {
        return false;
    }
    // 后续半平面剪裁约定顶点为逆时针顺序。
    if(dDoubleArea < 0.0) {
        std::reverse(vecPoints.begin(), vecPoints.end());
        dDoubleArea = -dDoubleArea;
    }

    QVector<Triangle2D> vecTriangles;
    if(!triangulatePolygon(vecPoints,
                           dAreaTolerance, vecTriangles)) {
        return false;
    }
    oPolygon.m_vecPoints = vecPoints;
    oPolygon.m_vecTriangles = vecTriangles;
    oPolygon.m_oBounds = polygonBounds(vecPoints);
    oPolygon.m_dArea = dDoubleArea * 0.5;
    return isFiniteValue(oPolygon.m_dArea) && oPolygon.m_dArea > 0.0;
}

QPointF segmentLineIntersection(const QPointF& s,
                                const QPointF& e,
                                double dStartDistance,
                                double dEndDistance)
{
    const double dDenominator = dStartDistance - dEndDistance;
    const double dT = dDenominator == 0.0
            ? 0.0 : dStartDistance / dDenominator;
    return QPointF(s.x() + dT * (e.x() - s.x()),
                   s.y() + dT * (e.y() - s.y()));
}

QVector<QPointF> clipByEdge(const QVector<QPointF>& vecSubject,
                            const QPointF& a,
                            const QPointF& b,
                            double dAreaTolerance)
{
    // Sutherland-Hodgman 半平面剪裁，裁剪边左侧视为三角形内部。
    QVector<QPointF> vecOutput;
    if(vecSubject.isEmpty()) {
        return vecOutput;
    }
    QPointF s = vecSubject.last();
    double dStartDistance = crossProduct(a, b, s);
    bool bStartInside = dStartDistance >= -dAreaTolerance;
    for(int nIndex = 0; nIndex < vecSubject.size(); ++nIndex) {
        const QPointF e = vecSubject[nIndex];
        const double dEndDistance = crossProduct(a, b, e);
        const bool bEndInside = dEndDistance >= -dAreaTolerance;
        if(bEndInside) {
            if(!bStartInside) {
                vecOutput.append(segmentLineIntersection(
                    s, e, dStartDistance, dEndDistance));
            }
            vecOutput.append(e);
        } else if(bStartInside) {
            vecOutput.append(segmentLineIntersection(
                s, e, dStartDistance, dEndDistance));
        }
        s = e;
        dStartDistance = dEndDistance;
        bStartInside = bEndInside;
    }
    return vecOutput;
}

double triangleIntersectionArea(const Triangle2D& oSubject,
                                const Triangle2D& oClip,
                                double dAreaTolerance)
{
    QVector<QPointF> vecClipped;
    vecClipped.append(oSubject.m_a);
    vecClipped.append(oSubject.m_b);
    vecClipped.append(oSubject.m_c);
    vecClipped = clipByEdge(vecClipped,
                            oClip.m_a, oClip.m_b, dAreaTolerance);
    vecClipped = clipByEdge(vecClipped,
                            oClip.m_b, oClip.m_c, dAreaTolerance);
    vecClipped = clipByEdge(vecClipped,
                            oClip.m_c, oClip.m_a, dAreaTolerance);
    return vecClipped.size() >= 3
            ? qAbs(signedDoubleArea(vecClipped)) * 0.5 : 0.0;
}

bool rectanglesOverlap(const QRectF& a,
                       const QRectF& b,
                       double dTolerance)
{
    return a.right() >= b.left() - dTolerance &&
           b.right() >= a.left() - dTolerance &&
           a.bottom() >= b.top() - dTolerance &&
           b.bottom() >= a.top() - dTolerance;
}

bool polygonsHavePositiveAreaIntersection(
        const PreparedPolygon& a,
        const PreparedPolygon& b)
{
    // 业务语义是“面积相交”：仅边界或顶点接触不把单元纳入统计。
    const double dFirstScale =
            qMax(a.m_oBounds.width(), a.m_oBounds.height());
    const double dSecondScale =
            qMax(b.m_oBounds.width(), b.m_oBounds.height());
    // 交叠容差跟随较小多边形，避免大选区吞掉小单元的真实窄交叠。
    const double dScale = qMax(1.0, qMin(dFirstScale, dSecondScale));
    const double dLinearTolerance = dScale * 1.0e-10;
    if(!rectanglesOverlap(a.m_oBounds, b.m_oBounds,
                          dLinearTolerance)) {
        return false;
    }
    const double dAreaTolerance = dScale * dScale * 1.0e-13;
    const double dPositiveAreaThreshold = qMax(
                dAreaTolerance,
                qMin(a.m_dArea, b.m_dArea) * 1.0e-12);
    for(int nA = 0; nA < a.m_vecTriangles.size(); ++nA) {
        for(int nB = 0; nB < b.m_vecTriangles.size(); ++nB) {
            if(triangleIntersectionArea(
                   a.m_vecTriangles[nA], b.m_vecTriangles[nB],
                   dAreaTolerance) > dPositiveAreaThreshold) {
                return true;
            }
        }
    }
    return false;
}

bool readCellPolygon(vtkCell* pCell,
                     QVector<QPointF>& vecPoints)
{
    // 当前 PEBI 结果只接受可直接解释为平面多边形的单元类型。
    vecPoints.clear();
    if(pCell == NULL ||
       (pCell->GetCellType() != VTK_TRIANGLE &&
        pCell->GetCellType() != VTK_QUAD &&
        pCell->GetCellType() != VTK_POLYGON) ||
       pCell->GetPoints() == NULL ||
       pCell->GetNumberOfPoints() < 3) {
        return false;
    }
    double aPoint[3] = { 0.0, 0.0, 0.0 };
    for(vtkIdType nPointIndex = 0;
        nPointIndex < pCell->GetNumberOfPoints(); ++nPointIndex) {
        pCell->GetPoints()->GetPoint(nPointIndex, aPoint);
        if(!isFiniteValue(aPoint[0]) || !isFiniteValue(aPoint[1])) {
            vecPoints.clear();
            return false;
        }
        vecPoints.append(QPointF(aPoint[0], aPoint[1]));
    }
    return true;
}

void addCompensated(double dValue,
                    double& dSum,
                    double& dCorrection)
{
    // Kahan 补偿求和用于缓解单元数量较多或权重差异较大时的累计误差。
    const double dAdjusted = dValue - dCorrection;
    const double dCandidate = dSum + dAdjusted;
    dCorrection = (dCandidate - dSum) - dAdjusted;
    dSum = dCandidate;
}

}

class nmPebiRegionPressureCalculatorPrivate
{
public:
    nmPebiRegionPressureCalculatorPrivate()
        : m_bPrepared(false),
          m_dSelectedCellArea(0.0)
    {
    }

    bool m_bPrepared;
    QSharedPointer<const nmPebiResultSnapshot> m_pSnapshot;
    QVector<QPointF> m_vecPolygonPoints;
    QVector<vtkIdType> m_vecIntersectingCellIds;
    QVector<CachedCell> m_vecCells;
    double m_dSelectedCellArea;
};

nmPebiRegionPressureStatistics::nmPebiRegionPressureStatistics()
{
    clear();
}

void nmPebiRegionPressureStatistics::clear()
{
    m_bAvailable = false;
    m_nPressureFrameIndex = -1;
    m_dPressureTimeHours = 0.0;
    m_nIntersectingCellCount = 0;
    m_nValidCellCount = 0;
    m_nIgnoredCellCount = 0;
    m_dSelectedCellAreaSquareMeters = 0.0;
    m_dEffectivePoreVolumeCubicMeters = 0.0;
    m_dMinimumPressureMpa = 0.0;
    m_dMaximumPressureMpa = 0.0;
    m_dPoreVolumeWeightedPressureMpa = 0.0;
}

nmPebiRegionPressureCalculator::nmPebiRegionPressureCalculator()
    : m_pPrivate(new nmPebiRegionPressureCalculatorPrivate)
{
}

nmPebiRegionPressureCalculator::~nmPebiRegionPressureCalculator()
{
    delete m_pPrivate;
    m_pPrivate = NULL;
}

bool nmPebiRegionPressureCalculator::prepare(
        const QSharedPointer<const nmPebiResultSnapshot>& pSnapshot,
        const QVector<QPointF>& vecPolygonPoints,
        QString* pError)
{
    if(pError != NULL) {
        pError->clear();
    }
    const vtkIdType nCellCount = pSnapshot.isNull()
            ? 0 : pSnapshot->getResultCellCount();
    if(m_pPrivate == NULL || pSnapshot.isNull() ||
       !pSnapshot->m_bComplete || pSnapshot->m_pResultGrid == NULL ||
       nCellCount <= 0 || pSnapshot->m_pPorosity == NULL ||
       pSnapshot->m_pThickness == NULL ||
       pSnapshot->m_pPorosity->GetNumberOfComponents() != 1 ||
       pSnapshot->m_pThickness->GetNumberOfComponents() != 1 ||
       pSnapshot->m_pPorosity->GetNumberOfTuples() != nCellCount ||
       pSnapshot->m_pThickness->GetNumberOfTuples() != nCellCount ||
       pSnapshot->m_vecPressureFrames.isEmpty()) {
        if(pError != NULL) {
            *pError = QString::fromLatin1(
                        "The PEBI result snapshot is invalid.");
        }
        return false;
    }

    PreparedPolygon oSelection;
    if(!preparePolygon(vecPolygonPoints, oSelection)) {
        if(pError != NULL) {
            *pError = QString::fromLatin1(
                        "The region polygon is invalid or self-intersecting.");
        }
        return false;
    }

    // 先在候选状态中完成全部校验和缓存，失败时不破坏上一个可用区域。
    nmPebiRegionPressureCalculatorPrivate oCandidate;
    oCandidate.m_pSnapshot = pSnapshot;
    oCandidate.m_vecPolygonPoints = oSelection.m_vecPoints;

    vtkUnstructuredGrid* pGrid = pSnapshot->m_pResultGrid;
    QVector<QPointF> vecCellPoints;
    double dAreaCorrection = 0.0;
    for(vtkIdType nCellId = 0;
        nCellId < pGrid->GetNumberOfCells(); ++nCellId) {
        vtkCell* pCell = pGrid->GetCell(nCellId);
        if(!readCellPolygon(pCell, vecCellPoints)) {
            continue;
        }
        PreparedPolygon oCellPolygon;
        if(!preparePolygon(vecCellPoints, oCellPolygon) ||
           !polygonsHavePositiveAreaIntersection(
                oSelection, oCellPolygon)) {
            continue;
        }

        CachedCell oCell;
        oCell.m_nCellId = nCellId;
        // 选中规则是“有正面积相交即纳入”，纳入后使用整个单元面积计算孔隙体积权重。
        oCell.m_dArea = oCellPolygon.m_dArea;
        const double dPorosity =
                pSnapshot->m_pPorosity->GetValue(nCellId);
        const double dThickness =
                pSnapshot->m_pThickness->GetValue(nCellId);
        const double dWeight = oCell.m_dArea * dThickness * dPorosity;
        // 几何相交单元保留在计数中；非法物性只在压力统计时被忽略。
        oCell.m_bWeightValid = isFiniteValue(dPorosity) &&
                isFiniteValue(dThickness) && isFiniteValue(dWeight) &&
                dPorosity > 0.0 && dThickness > 0.0 && dWeight > 0.0;
        oCell.m_dPoreVolumeWeight = oCell.m_bWeightValid ? dWeight : 0.0;
        oCandidate.m_vecCells.append(oCell);
        oCandidate.m_vecIntersectingCellIds.append(nCellId);
        addCompensated(oCell.m_dArea,
                       oCandidate.m_dSelectedCellArea,
                       dAreaCorrection);
    }
    if(!isFiniteValue(oCandidate.m_dSelectedCellArea)) {
        if(pError != NULL) {
            *pError = QString::fromLatin1(
                        "The selected cell area is invalid.");
        }
        return false;
    }

    oCandidate.m_bPrepared = true;
    *m_pPrivate = oCandidate;
    return true;
}

void nmPebiRegionPressureCalculator::clear()
{
    if(m_pPrivate != NULL) {
        *m_pPrivate = nmPebiRegionPressureCalculatorPrivate();
    }
}

bool nmPebiRegionPressureCalculator::isPrepared() const
{
    return m_pPrivate != NULL && m_pPrivate->m_bPrepared &&
            !m_pPrivate->m_pSnapshot.isNull();
}

QString nmPebiRegionPressureCalculator::getSnapshotId() const
{
    return isPrepared()
            ? m_pPrivate->m_pSnapshot->getSnapshotId() : QString();
}

const QVector<QPointF>&
nmPebiRegionPressureCalculator::getPolygonPoints() const
{
    return m_pPrivate->m_vecPolygonPoints;
}

const QVector<vtkIdType>&
nmPebiRegionPressureCalculator::getIntersectingCellIds() const
{
    return m_pPrivate->m_vecIntersectingCellIds;
}

bool nmPebiRegionPressureCalculator::calculateFrame(
        int nPressureFrameIndex,
        nmPebiRegionPressureStatistics& oStatistics,
        QString* pError) const
{
    oStatistics.clear();
    if(pError != NULL) {
        pError->clear();
    }
    if(!isPrepared() || nPressureFrameIndex < 0 ||
       nPressureFrameIndex >=
            m_pPrivate->m_pSnapshot->m_vecPressureFrames.size()) {
        if(pError != NULL) {
            *pError = QString::fromLatin1(
                        "The pressure frame is invalid.");
        }
        return false;
    }

    const nmPebiResultSnapshot::PressureFrame& oFrame =
            m_pPrivate->m_pSnapshot->m_vecPressureFrames[
                nPressureFrameIndex];
    vtkDoubleArray* pPressure = oFrame.m_pPressure;
    if(pPressure == NULL || pPressure->GetNumberOfComponents() != 1 ||
       pPressure->GetNumberOfTuples() !=
            m_pPrivate->m_pSnapshot->getResultCellCount() ||
       !isFiniteValue(oFrame.m_dTime)) {
        if(pError != NULL) {
            *pError = QString::fromLatin1(
                        "The pressure frame data is invalid.");
        }
        return false;
    }

    oStatistics.m_nPressureFrameIndex = nPressureFrameIndex;
    oStatistics.m_dPressureTimeHours = oFrame.m_dTime;
    oStatistics.m_nIntersectingCellCount =
            m_pPrivate->m_vecCells.size();
    oStatistics.m_dSelectedCellAreaSquareMeters =
            m_pPrivate->m_dSelectedCellArea;

    // 交集及静态权重已在 prepare() 缓存，逐帧只读压力并做线性汇总。
    double dPoreVolumeSum = 0.0;
    double dPoreVolumeCorrection = 0.0;
    double dWeightedPressureSum = 0.0;
    double dWeightedPressureCorrection = 0.0;
    double dMinimumPressure = DBL_MAX;
    double dMaximumPressure = -DBL_MAX;
    for(int nIndex = 0;
        nIndex < m_pPrivate->m_vecCells.size(); ++nIndex) {
        const CachedCell& oCell = m_pPrivate->m_vecCells[nIndex];
        const double dPressure = pPressure->GetValue(oCell.m_nCellId);
        if(!oCell.m_bWeightValid || !isFiniteValue(dPressure)) {
            continue;
        }
        addCompensated(oCell.m_dPoreVolumeWeight,
                       dPoreVolumeSum, dPoreVolumeCorrection);
        addCompensated(oCell.m_dPoreVolumeWeight * dPressure,
                       dWeightedPressureSum,
                       dWeightedPressureCorrection);
        dMinimumPressure = qMin(dMinimumPressure, dPressure);
        dMaximumPressure = qMax(dMaximumPressure, dPressure);
        ++oStatistics.m_nValidCellCount;
    }
    oStatistics.m_nIgnoredCellCount =
            oStatistics.m_nIntersectingCellCount -
            oStatistics.m_nValidCellCount;
    // 空选区或全部权重无效都是可展示结果，不把它们当成计算错误。
    if(oStatistics.m_nValidCellCount == 0) {
        return true;
    }
    if(!isFiniteValue(dPoreVolumeSum) || dPoreVolumeSum <= 0.0 ||
       !isFiniteValue(dWeightedPressureSum)) {
        if(pError != NULL) {
            *pError = QString::fromLatin1(
                        "The pore-volume weighted pressure is invalid.");
        }
        return false;
    }

    const double dAveragePressure =
            dWeightedPressureSum / dPoreVolumeSum;
    if(!isFiniteValue(dMinimumPressure) ||
       !isFiniteValue(dMaximumPressure) ||
       !isFiniteValue(dAveragePressure)) {
        if(pError != NULL) {
            *pError = QString::fromLatin1(
                        "The regional pressure statistics are invalid.");
        }
        return false;
    }
    oStatistics.m_bAvailable = true;
    oStatistics.m_dEffectivePoreVolumeCubicMeters = dPoreVolumeSum;
    oStatistics.m_dMinimumPressureMpa = dMinimumPressure;
    oStatistics.m_dMaximumPressureMpa = dMaximumPressure;
    oStatistics.m_dPoreVolumeWeightedPressureMpa = dAveragePressure;
    return true;
}