// Quick smoke test: load MHD and sample three physical slices (no Qt GUI)
#include "DoseImage.h"

#include <cstdio>
#include <cmath>
#include <algorithm>
#include <vector>

static void sampleSlice(const DoseImage::Pointer& dose, int axisCode)
{
  double bmin[3], bmax[3];
  dose->worldBounds(bmin, bmax);
  const double cx = 0.5 * (bmin[0] + bmax[0]);
  const double cy = 0.5 * (bmin[1] + bmax[1]);
  const double cz = 0.5 * (bmin[2] + bmax[2]);

  double u[3] = {0, 0, 0}, v[3] = {0, 0, 0};
  if (axisCode == 0) { // axial
    u[0] = -1; v[1] = -1;
  } else if (axisCode == 1) { // sag
    u[1] = -1; v[2] = 1;
  } else {
    u[0] = 1; v[2] = 1;
  }

  const auto sp = dose->spacing();
  const double su = std::max(1e-3, std::abs(u[0]) * std::abs(sp[0]) + std::abs(u[1]) * std::abs(sp[1]) + std::abs(u[2]) * std::abs(sp[2]));
  const double sv = std::max(1e-3, std::abs(v[0]) * std::abs(sp[0]) + std::abs(v[1]) * std::abs(sp[1]) + std::abs(v[2]) * std::abs(sp[2]));

  double halfU = 0, halfV = 0;
  for (int iz = 0; iz < 2; ++iz)
    for (int iy = 0; iy < 2; ++iy)
      for (int ix = 0; ix < 2; ++ix) {
        const double px = (ix ? bmax[0] : bmin[0]) - cx;
        const double py = (iy ? bmax[1] : bmin[1]) - cy;
        const double pz = (iz ? bmax[2] : bmin[2]) - cz;
        halfU = std::max(halfU, std::abs(px * u[0] + py * u[1] + pz * u[2]));
        halfV = std::max(halfV, std::abs(px * v[0] + py * v[1] + pz * v[2]));
      }

  int nu = static_cast<int>(std::ceil(2.0 * halfU / su)) + 1;
  int nv = static_cast<int>(std::ceil(2.0 * halfV / sv)) + 1;
  std::printf("axis=%d size=%dx%d su=%.3f sv=%.3f half=%.1f,%.1f\n", axisCode, nu, nv, su, sv, halfU, halfV);

  double sum = 0, mx = -1e30;
  int nonzero = 0;
  const double x0 = -0.5 * (nu - 1) * su;
  const double y0 = -0.5 * (nv - 1) * sv;
  for (int j = 0; j < nv; ++j) {
    const double tv = y0 + j * sv;
    for (int i = 0; i < nu; ++i) {
      const double tu = x0 + i * su;
      const float val = dose->sampleWorld(cx + tu * u[0] + tv * v[0],
                                          cy + tu * u[1] + tv * v[1],
                                          cz + tu * u[2] + tv * v[2]);
      sum += val;
      mx = std::max(mx, (double)val);
      if (val > dose->minValue() + 1e-6f)
        ++nonzero;
    }
  }
  std::printf("  mean=%.4f max=%.4f nonzero=%d/%d\n", sum / (nu * nv), mx, nonzero, nu * nv);
}

int main(int argc, char** argv)
{
  const char* path = (argc > 1) ? argv[1]
    : "D:/testNewPhy/halcyon/examples/che1/0000137408/internal/doseOldAve.mhd";
  QString err;
  auto dose = DoseImage::loadFromFile(QString::fromLocal8Bit(path), &err);
  if (!dose) {
    std::printf("LOAD FAIL: %s\n", err.toLocal8Bit().constData());
    return 1;
  }
  std::printf("loaded %s min=%.4f max=%.4f dim=%zu %zu %zu\n",
              dose->name().toLocal8Bit().constData(), dose->minValue(), dose->maxValue(),
              dose->region().GetSize()[0], dose->region().GetSize()[1], dose->region().GetSize()[2]);
  sampleSlice(dose, 0);
  sampleSlice(dose, 1);
  sampleSlice(dose, 2);
  std::printf("OK\n");
  return 0;
}