// Moshade - Compute average cross section for 3D shapes
// Copyright (C) 2018-2026 Eric Bechet - bechet@cadxfem.org
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.

#include "moshadeconfig.h"
#ifdef MOSHADE_GRAPHICAL_INTERFACE
#include "ndisplay.h"
#endif //MOSHADE_GRAPHICAL_INTERFACE
#include "stl_io.h"
#include "geom_op.h"
#include "mesh_op.h"

#include <iostream>
#include <iomanip>
#include <sstream>
#include <cstdlib>
#include <cstring>
#ifdef OPENMP_FOUND
#include <omp.h>
#endif // OPENMP_FOUND

struct cmdline_info
{
  bool verbose; // print a lot of info
  bool quiet; // prints only the results
  bool naive; // choose naive integration quadrature
  bool override; // if .gm
  double szsph; // spheres size
  double szcyl; // cyl sizes
  int nele; // discretization of sph and cyl
  int N; // number of samples per half meridian
  int nprocs; // number of threads in parallel mode
  std::string file; // filename of geometry (.stl or descriptive .gm)
};

bool parse_cmdline(int argc, char *argv[], cmdline_info &data);
void help(void);
bool load_geometry(const cmdline_info &data,tri_mesh &mesh);
void display_geo(tri_mesh &mesh,bool compute_shades);
int compute_average_CS(tri_mesh& mesh, int N, bool naive, double& av_CS, bool verbose, int nprocs, bool quiet);

int main(int argc, char *argv[])
{
  // load sample data
  cmdline_info cmdline;
  bool ok=parse_cmdline(argc,argv,cmdline);
  if (ok)
  {
    if (!cmdline.quiet) std::cout << "moshade - compute average cross sections of complex geometrical shapes"<< std::endl;
    tri_mesh mesh;
    bool loadok=load_geometry(cmdline,mesh);
    if (loadok)
    {
      mesh.build_topology(true,cmdline.quiet); // to get normals and info
//      mesh.clear_topology();
      if (cmdline.N<1) // special case : display shape & perform some projections for testing
      {
        display_geo(mesh,cmdline.N<0);
      }
      else
      {
        if (!cmdline.quiet)
        if (cmdline.naive)
          std::cout << "Using a naive integration quadrature..."<< std::endl;
        else
          std::cout << "Using an optimized integration quadrature..."<< std::endl;
        std::cout.setf(std::ios::scientific);
        std::cout.precision(6);
        double avCS,tot_area=mesh.area(),volume=mesh.volume();
        int cnt=compute_average_CS(mesh,cmdline.N,cmdline.naive,avCS,cmdline.verbose,cmdline.nprocs,cmdline.quiet);
        std::cout.precision(6);
        if (!cmdline.quiet) std::cout << "  tot_area       volume       avg_CS       ratio    #samples "<< std::endl;
        std::cout << tot_area << " " << volume << " " << avCS << " " << tot_area/avCS << " " <<  cnt << std::endl;
      }
    }
    else std::cout << "Empty mesh !" << std::endl;
  }
  else
    help();
}

bool parse_cmdline(int argc, char *argv[], cmdline_info &data)
{
  data.naive=false;
  data.override=false;
  data.verbose=false;
  data.quiet=false;
  data.nprocs=0;
  if (argc>=3)
  {
    data.file=std::string(argv[1]);
    data.N=atoi(argv[2]);
    for (int i=3;i<argc;++i)
    {
      if (!mo_stricmp(argv[i],"naive"))
        data.naive=true;
      if (!mo_stricmp(argv[i],"override"))
      {
        data.override=true;
        data.szsph=atof(argv[++i]);
        data.szcyl=atof(argv[++i]);
        data.nele=atoi(argv[++i]);
      }
      if (!mo_stricmp(argv[i],"verbose"))
        data.verbose=true;
      if (!mo_stricmp(argv[i],"quiet"))
        data.quiet=true;
      if (!mo_stricmp(argv[i],"nproc"))
      {
        data.nprocs=atof(argv[++i]);
      }
    }
    return true;
  }
  else return false; // help !
}

bool load_geometry(const cmdline_info &data,tri_mesh &mesh)
{
  try
  {
    mesh.load_stl(data.file);
    return (mesh.size());
  }
  catch(int i)
  {
    try
    {
      mesh.load_gm(data.file,data.override,data.szsph,data.szcyl,data.nele);
      return (mesh.size());
    }
    catch (int i)
    {
      if (!data.quiet)
      {
        std::cerr << "File not read, failing..." << std::endl;
        std::cerr << "Testing with sample volume..." << std::endl;
      }
//      mesh.add_cylinder(npoint3(0,0,0),npoint3(10,10,10),1.0,true,true,10);
//      mesh.add_cylinder(npoint3(0,0,0),npoint3(0,10*sqrt(3)/sqrt(2),10*sqrt(3)/sqrt(2)),1.0,true,true,10);
//      mesh.add_cylinder(npoint3(0,0,0),npoint3(0,0,10*sqrt(3)),1.0,true,true,10);
      mesh.add_helix(1.0,0.3,100,5);
//      mesh.add_sphere(npoint3(0,0,0),1.0,4);
//     return(1);
      return (mesh.size());
    }
    return true;
  }
  return true;
}

void help(void)
{
  std::cerr << "Moshade - compute average cross sections of complex geometrical shapes." << std::endl;
  std::cerr << "This is free software (c) 2018-2026 Eric Bechet - eric@bechet.ca" << std::endl;
  std::cerr << "See the LICENCE file for license information." << std::endl;
  std::cerr << "Usage : moshade file value [naive] [nproc val] [verbose]" << std::endl;
  std::cerr << "where \"file\" is either a .stl ASCII STL file (stereolitography), or a composite .gm file" << std::endl;
  std::cerr << " \"value\" is an integer telling how many samples are used to compute " << std::endl;
  std::cerr << "the average cross section for each variable," << std::endl;
  std::cerr << "and \"naive\" tells to use a naive n*n integration quadrature "  << std::endl;
  std::cerr << "instead of an optimized one." << std::endl;
  std::cerr << "if \"value\" is zero, the part is simply displayed and no computations are done." << std::endl << std::endl;
  std::cerr << "\"nproc val\" sets the max number of threads for the computations (0=use all procs)" << std::endl << std::endl;
  std::cerr << "\"verbose\" will print all values of the orientations and cross sections " << std::endl;
  std::cerr << "on the standard output with the following structure :" << std::endl;
  std::cerr << "<topological info on the input file>" << std::endl;
  std::cerr << "theta1 phi1 #tris cross-section1 projected-area1 ratio1 %" << std::endl;
  std::cerr << "theta2 phi2 #tris cross-section2 projected-area2 ratio2 %" << std::endl;
  std::cerr << "..." << std::endl;
  std::cerr << "  min_CS        phimin       thmin"<< std::endl;
  std::cerr << "  max_CS        phimax       thmax"<< std::endl;
  std::cerr << " At the end, the final values are printed"<< std::endl;
  std::cerr << "total_area volume  average_cross_section   ratio   number_of_samples" << std::endl;
  std::cerr << std::endl;
}

void display_geo(tri_mesh &mesh,bool compute_shades)
{
#ifdef MOSHADE_GRAPHICAL_INTERFACE
  data_container data;
  ndisplay disp;
  properties p=data.getproptriangles();
  p.c.A=50;
  p.edgeon=true;
  p.edgecolor=p.c;
  data.setproptriangles(p);
  p=data.getproplines();
  p.c=color(255,255,255);
  data.setproplines(p);
  p=data.getproppoints();
  p.c=color(255,255,255);
  p.pointsize=3;
  data.setproppoints(p);
  mesh.display(data);
#endif //MOSHADE_GRAPHICAL_INTERFACE
  double tot_area=mesh.area();
  double volume=mesh.volume();
  if (compute_shades)
  {
#ifdef MOSHADE_GRAPHICAL_INTERFACE
    p.edgeon=false;
    data.setproptriangles(p);
#endif //MOSHADE_GRAPHICAL_INTERFACE
    double mattr[4][4];
    npoint3 bbmin,bbmax;
    mesh.bb(bbmin,bbmax);
    tri_mesh meshsh;
    tri_mesh meshtr;
    double area=0,prarea,sharea;
    int nbele;
    fill_tr_mat(2,bbmin[2]-(bbmax[2]-bbmin[2])*0.1,mattr);
    nbele=shade(0,0,mesh,sharea,prarea,&meshsh);
    area+=sharea*(2*Pi/3);
    std::cout << 0.0 << " " <<  0.0 << " " << nbele << " " << sharea << " " << prarea << " " << sharea/prarea << std::endl;
    meshtr.transform(meshsh,mattr); // translate projection to the side of the bb
#ifdef MOSHADE_GRAPHICAL_INTERFACE
    meshtr.display(data);
#endif //MOSHADE_GRAPHICAL_INTERFACE

    fill_tr_mat(0,bbmin[0]-(bbmax[0]-bbmin[0])*0.1,mattr);
    nbele=shade(Pi/2,0,mesh,sharea,prarea,&meshsh);
    area+=sharea*(2*Pi/3);
    std::cout << Pi/2 << " " <<  0.0 << " " << nbele << " " << sharea << " " << prarea << " " << sharea/prarea << std::endl;
    meshtr.transform(meshsh,mattr); // translate projection to the side of the bb
#ifdef MOSHADE_GRAPHICAL_INTERFACE
    meshtr.display(data);
#endif //MOSHADE_GRAPHICAL_INTERFACE

    fill_tr_mat(1,bbmin[1]-(bbmax[1]-bbmin[1])*0.1,mattr);
    nbele=shade(0,Pi/2,mesh,sharea,prarea,&meshsh);
    area+=sharea*(2*Pi/3);
    std::cout << 0.0 << " " <<  Pi/2 << " " << nbele << " " << sharea << " " << prarea << " " << sharea/prarea << std::endl;
    meshtr.transform(meshsh,mattr); // translate projection to the side of the bb
#ifdef MOSHADE_GRAPHICAL_INTERFACE
    meshtr.display(data);
#endif //MOSHADE_GRAPHICAL_INTERFACE
    std::cout << "  tot_area       volume       avg_CS       ratio    #samples "<< std::endl;
    std::cout << tot_area << " " << volume << " " << area*(0.5/Pi) << " " << tot_area/(area*(0.5/Pi)) << " " <<  "3" << std::endl;
  }
  else
  {
    std::cout << "  tot_area       volume"<< std::endl;
    std::cout << tot_area << " " << volume <<  std::endl;
  }
#ifdef MOSHADE_GRAPHICAL_INTERFACE
  disp.init_data(data);
  disp.display();
#endif //MOSHADE_GRAPHICAL_INTERFACE

}


int compute_average_CS(tri_mesh &mesh,int N,bool naive,double &av_CS,bool verbose, int nprocs, bool quiet)
{
  double area=0;
  double minCS=1e30;
  double maxCS=-1e30;
  double phmin,phmax;
  double thmin,thmax;
  int cnt=0;
  if (verbose)
    std::cout << "    theta          phi       #tri   shaded_area   proj_area     ratio     %done"<< std::endl;
  double deb=0;
  double limit=0;
  double step=0;
  if (naive)
  {
    deb=-Pi/2+Pi/(4*N);
    limit=Pi/2;
    step=Pi/(2*N);
  }
  else
  {
    deb=-1+1./(2.*N);
    limit=1.;
    step=1./N;
  }
  int effprocs=1;
  #ifdef OPENMP_FOUND
  if (nprocs>0)
    omp_set_num_threads(nprocs);
  effprocs=omp_get_num_threads();
  #endif // OPENMP_FOUND
#pragma omp parallel for reduction(+:area)
  for (int ii=0;ii<2*N;++ii)
  {
    double x=deb+ii*step;
    double th;
    if (naive) th=x;else th=asin(x);
    for (double phi=-Pi/2+Pi/(4*N);phi<Pi/2;phi+=Pi/(2*N))
//    for (double phi=-Pi+Pi/(4*N);phi<Pi;phi+=Pi/(2*N))
    {
      double sharea;
      double prarea;
      int nbele=shade(th,phi,mesh,sharea,prarea);
      double wh;
      if (naive)
        wh=cos(th)*(Pi/(2*N))*(Pi/(2*N));
      else
        wh=(Pi/(4*N))*(2./N);

//      wh/=2;
      area+=sharea*wh;
      #pragma omp critical
      {
        if (minCS>sharea)
        {
          minCS=sharea;
          thmin=th;
          phmin=phi;
        }
        if (maxCS<sharea)
        {
          maxCS=sharea;
          thmax=th;
          phmax=phi;
        }
        cnt++;
        if (verbose)
          std::cout << std::setw(13) << th << " " << std::setw(13) << phi << " "<< std::setw(6) << nbele << " " << std::setw(12)<< sharea << " "<< std::setw(12) << prarea << " " << std::setw(12)<< sharea/prarea << " " << std::setw(3)<< (int)(cnt/(4.0*N*N)*100) << "%" << std::endl;
        else if (!quiet)
          std::cerr << "Done " << (int)(cnt/(4.0*N*N)*100) << " %\r" <<  std::flush;
      }
    }
  }

  double mat1[4][4];
  double mat2[4][4];
  double matf[4][4];
  fill_rot_mat(0,phmin,mat1); // first rotate around x
  fill_rot_mat(1,thmin,mat2); // then rotate around y
  mat_prod(mat1,mat2,matf);
  if (verbose)
  {
    std::cout << "  min_CS        phimin       thmin           x            y           z"<< std::endl;
    std::cout << minCS << " " << phmin << " " << thmin  << " " << matf[0][2] << " " << matf[1][2] << " " << matf[2][2]   << std::endl;
  }
  fill_rot_mat(0,phmax,mat1); // first rotate around x
  fill_rot_mat(1,thmax,mat2); // then rotate around y
  mat_prod(mat1,mat2,matf);
  if (verbose)
  {
    double tot_area=mesh.area();
    std::cout << "  max_CS        phimax       thmax           x            y           z"<< std::endl;
    std::cout << maxCS << " " << phmax << " " << thmax  << " " << matf[0][2] << " " << matf[1][2] << " " << matf[2][2]   << std::endl;
  }
  av_CS=area*(0.5/Pi);
  return cnt;
}