// Gnurbs - A curve and surface library
// Copyright (C) 2008-2026 Eric Bechet
//
// See the LICENSE file for contributions and license information.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//

#include <iostream>
#include <fstream>
#include <sstream>
#include "nbspline.h"
#include "ndisplay.h"
#include "nutil.h"
#include "point_set.h"

int main(void)
{
  PointSet pts,dest,dest2;
//  pts.reserve(1000000);
//  dest.reserve(1000000);
//  dest2.reserve(1000000);
  data_container data;
  double X=0,Y=0,Z=0;
  properties p;
  p=data.getproppoints();
  p.c=color(255,255,255,255);
  p.pointsize=2.0;
  data.setproppoints(p);
  p=data.getproplines();
  p.c=color(255,255,255);
  p.thickness=2.0;
  data.setproplines(p);
  std::cout.precision(15);
  std::cout << std::showpos;
  std::cout << std::scientific;
//  read_data_bin_gz("Pts_Aliantus_jpgTif_20121113_Net_GC3_FiltreCC.bin.gz",pts);
//  pts.read_bin_gz("01_2-nettoye.xyz.bin.gz");
//  read_data_bin_gz("02_3-nettoye.xyz.bin.gz",pts);
//pts.insert(scanpt());

//  pts.read_bin_gz("Arbe1_h.bin.gz");
//  pts.read_bin("Arbe1_h.bin");
//  pts.read_gz("Arbe1_h.txt.gz");
//  pts.read_bin("arbre2_hExport.bin");
  pts.read_bin("arbre3_hexport.bin");
//  pts.read_gz("arbre3_hexport.txt.gz");
//  pts.translate(true,X,Y,Z);

  npoint3 min,max;
  pts.bbox(min,max);
  std::cout << "Bounding box : " << std::endl;
  min.print(std::cout);
  max.print(std::cout);
  const int num_slice=20;
  nbspline courbe_MC[num_slice]; // la courbe sous forme b-spline
/*  double Zplan=152.25;

  double low=Zplan-Z;
  double up=Zplan-Z+2.75;
  double deltaz=up-low;
  double DeltaZplan=0.1*deltaz/((num_slice==1)?1:(num_slice));
*/

  double Zplan=-2750;

  double low=Zplan;
  double up=Zplan+8000;
  double deltaz=up-low;
  double DeltaZplan=0.019*deltaz/((num_slice==1)?1:(num_slice));
  p=data.getproppoints();
  p.pointsize=1;
  data.setproppoints(p);
  for (int i=0;i<num_slice;++i)
  {
    std::cout << "slice " << i <<  " : ";
    std::cout << low+deltaz*i/((num_slice==1)?1:(num_slice)) << " -> " << low+deltaz*i/((num_slice==1)?1:(num_slice))+DeltaZplan << std::endl;
    dest.resize(0);
//    dest2.resize(0);
//    std::cout << "filter" << std::endl;
    pts.filter(dest,low+deltaz*i/((num_slice==1)?1:(num_slice)),DeltaZplan);
//    double XX=0;
//    for (int i=0;i<dest.size();++i) XX+=dest[i].x();
//    XX/=dest.size();
//    dest.filterX(dest2,XX);
//    std::cout << "compute parameters" << std::endl;
    dest.compute_parameters();
    std::cout << "least_squares_bspline" << std::endl;
    int degree=3;
    int nbCP = (dest.size()/1000>(3*degree))?dest.size()/1000:(3*degree);
    dest.least_squares_bspline(degree,nbCP,courbe_MC[i]); // appel a la fonction de d'approximation
//    std::cout << "display" << std::endl;
//    std::cout << "min_u=" << courbe_MC[i].min_u() << std::endl;
//    std::cout << "max_u=" << courbe_MC[i].max_u() << std::endl;
    npoint pt;
    npoint3 pt3;
/*    courbe_MC[i].P(0.5,pt);
    pt3=pt;
    data.add_point(pt3);*/
//    for (int cp=0;cp<courbe_MC[i].nb_CP();++cp) courbe_MC[i].CP(cp)[0]+=i*1500;
//    for (int ii=0;ii<dest.size();++ii) dest[ii][0] +=i*1500;
if (1)   { 
    p=data.getproppoints();
  p.pointsize=2;
  data.setproppoints(p);
courbe_MC[i].Display(data);
    p=data.getproppoints();
  p.pointsize=1;
  data.setproppoints(p);

    dest.Display(data,1);}
  }
  p=data.getproppoints();
  p.pointsize=1;
  p.c=color(0,255,0,63);
  data.setproppoints(p);
  pts.Display(data,500);
  ndisplay display;
  display.init_data(data);
  display.display();
  return 0;
}