surftopgm.cpp

Go to the documentation of this file.

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_euclidean_traits_xy_3.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Bbox_2.h>
#include <CGAL/intersection_3_1.h>

#include <string>
#include <iostream>
#include <istream>
#include <ostream>
#include <fstream>

#include "parsecl.hpp"

struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};

typedef K::Plane_3 Plane;
typedef K::Line_3 Line;
typedef CGAL::Bbox_2 BoundingBox;

typedef CGAL::Triangulation_euclidean_traits_xy_3<K>  Gt;
typedef CGAL::Delaunay_triangulation_2<Gt> Delaunay;
typedef Delaunay::Point Point;
typedef Delaunay::Point_iterator PointIterator;
typedef Delaunay::Vertex_handle VertexHandle;
typedef Delaunay::Face_handle FaceHandle;
typedef Delaunay::Vertex_circulator VertexCirculator;

int main(int argc, char** argv) {
  using namespace Arak::Util;
  // Parse the command line.
  CommandLine cl;
  CommandLine::Parameter<double> 
    gamma("-g", "--gamma", "the gamma correction value", 1.0 / 0.45);
  cl.add(gamma);
  CommandLine::Parameter<int> 
    widthSpec("-w", "--width", "the image width", -1);
  cl.add(widthSpec);
  CommandLine::Parameter<int> 
    heightSpec("-h", "--height", "the image height", -1);
  cl.add(heightSpec);
  CommandLine::Parameter<int> 
    depth("-d", "--depth", "the image depth (# gray values)", 256);
  cl.add(depth);
  CommandLine::Parameter<double> 
    max("-min", "--min-intensity", "the surface value mapped to white", 0.0);
  cl.add(max);
  CommandLine::Parameter<double> 
    min("-max", "--max-intensity", "the surface value mapped to black", 1.0);
  cl.add(min);
  CommandLine::Option
    linear("-l", "--linear-interpolation", 
     "use linear instead of nearest-neighbor interpolation");
  cl.add(linear);
  CommandLine::Parameter<int> 
    def("-n", "--nan", 
  "the value to use for pixels outside the convex hull", 255);
  cl.add(def);
  CommandLine::Parameter<std::string> 
    inputPath("-i", "--input", 
        "the input file of point color estimates", std::string());
  cl.add(inputPath);
  CommandLine::Parameter<std::string> 
    outputPath("-o", "--output", "the output PGM file", std::string());
  cl.add(outputPath);
  if (!cl.parse(argc, argv, std::cerr)) {
    cl.printUsage(std::cerr);
    exit(1);
  }

  // Grab the input and output streams.
  std::istream* in = &std::cin;
  if (inputPath.supplied()) {
    in = new std::ifstream(inputPath.value().data());
    if (!in->good()) {
      std::cerr << "Error: cannot open file " << inputPath.value()
    << " for reading." << std::endl;
      exit(1);
    }
  }
  std::ostream* out = &std::cout;
  if (outputPath.supplied()) {
    out = new std::ofstream(outputPath.value().data());
    if (!out->good()) {
      std::cerr << "Error: cannot open file " << outputPath.value()
    << " for writing." << std::endl;
      exit(1);
    }
  }

  // Compute the gamma corrected values.
  int* gcv = new int[depth.value()];
  gcv[0] = 0;
  for (int i = 1; i < depth.value(); i++) {
    double raw = double(i) / double(depth.value() - 1);
    double corr = pow(raw, gamma.value());
    gcv[i] = int(round(corr * double(depth.value() - 1)));
  }

  // Read the point data.
  std::istream_iterator<Point> begin(*in);
  std::istream_iterator<Point> end;
  Delaunay dt;
  dt.insert(begin, end);
  std::cerr << "Read " << dt.number_of_vertices() 
      << " points." << std::endl;
  
  // Compute the bounding box by traversing the convex hull of the points.
  VertexHandle infVertexHandle = dt.infinite_vertex();
  VertexCirculator circ = dt.incident_vertices(infVertexHandle);
  VertexCirculator first = dt.incident_vertices(infVertexHandle);
  Point p = circ->point();
  BoundingBox bbox(CGAL::to_double(p.x()), CGAL::to_double(p.y()), 
       CGAL::to_double(p.x()), CGAL::to_double(p.y()));
  while (++circ != first) {
    p = circ->point();
    bbox = bbox + BoundingBox(CGAL::to_double(p.x()), CGAL::to_double(p.y()), 
            CGAL::to_double(p.x()), CGAL::to_double(p.y()));
  }
  std::cerr << "Bounding box: " << bbox << std::endl;

  // Compute the height and width from the specs and the aspect ratio
  // of the bounding box.
  int height, width;
  if (widthSpec.supplied() && heightSpec.supplied()) {
    // Use the specified values, ignoring the aspect ratio.
    height = heightSpec.value();
    width = widthSpec.value();
  } else if (!widthSpec.supplied() && !heightSpec.supplied()) {
    // Infer a good set of values from the aspect ratio and sampling
    // density.
    int size = dt.number_of_vertices();
    double aspectRatio = 
      (bbox.ymax() - bbox.ymin()) / (bbox.xmax() - bbox.xmin());
    width = int(round(sqrt(double(size) / aspectRatio)));
    height = int(aspectRatio * double(width));
  } else if (!heightSpec.supplied()) {
    // Use the supplied width and scale the height to maintain the
    // aspect ratio.
    width = widthSpec.value();
    double aspectRatio = 
      (bbox.ymax() - bbox.ymin()) / (bbox.xmax() - bbox.xmin());
    height = int(aspectRatio * double(width));
  } else {
    // Use the supplied height and scale the width to maintain the
    // aspect ratio.
    height = heightSpec.value();
    double aspectRatio = 
      (bbox.ymax() - bbox.ymin()) / (bbox.xmax() - bbox.xmin());
    width = int(double(height) / aspectRatio);
  }

  // Compute the units per dot (UPD) and report the resolution.
  K::FT upd((bbox.xmax() - bbox.xmin()) / double(width));
  std::cerr << "Output resolution: " << width << "x" << height << std::endl;

  // Compute the interpolated values and output the file.
  *out << "P2" << std::endl;                       // magic number
  *out << "# " << outputPath.value() << std::endl; // filename
  *out << width << " " << height << std::endl;     // resolution
  *out << (depth.value() - 1) << std::endl;        // max-value
  K::FT y = bbox.ymax();
  int num_def = 0;
  FaceHandle fh = dt.infinite_face();
  for (int i = height - 1; i >= 0; i--) {
    K::FT x = bbox.xmin();
    for (int j = 0; j < width; j++) {
      Point q(x, y, K::FT(0));
      fh = dt.locate(q, fh);
      if ((fh == NULL) || dt.is_infinite(fh)) {
  // Point lies outside convex hull.
  *out << def.value() << " ";
  num_def++;
      } else if (linear.supplied()) {
  // Use linear interpolation.  First form the plane containing
  // the vertices of the face.
  Plane plane(fh->vertex(0)->point(),
        fh->vertex(1)->point(),
        fh->vertex(2)->point());
  // Now form the vertical line at the query point.
  Line line(q, Point(x, y, K::FT(1.0)));
  // Compute the intersection point in 3-space.
  Point interp;
  CGAL::Object result = intersection(plane, line);
  assert(CGAL::assign(interp, result));
  // The interpolated value is the z-coordinate.
  double value = CGAL::to_double(interp.z());
  // Compute the intensity associated with this value.
  double intensity = 
    (value - min.value()) / (max.value() - min.value());
  intensity = (intensity < 0.0) ? 0.0 : intensity;
  intensity = (intensity > 1.0) ? 1.0 : intensity;
  // Compute the pixel value and then use gamma correction.
  *out << gcv[int(round(intensity * double(depth.value() - 1)))] << " ";
      } else {
  // Use nearest neighbor interpolation.
  Point p0 = fh->vertex(0)->point();
  Point p1 = fh->vertex(1)->point();
  Point p2 = fh->vertex(2)->point();
  // The interpolated value is the z-coordinate of the nearest
  // vertex.
  double value;
  if (CGAL::compare_distance_to_point(q, p0, p1) ==
      CGAL::SMALLER) {
    if (CGAL::compare_distance_to_point(q, p0, p2) ==
        CGAL::SMALLER) 
      value = CGAL::to_double(p0.z());
    else
      value = CGAL::to_double(p2.z());
  } else {
    if (CGAL::compare_distance_to_point(q, p1, p2) ==
        CGAL::SMALLER) 
      value = CGAL::to_double(p1.z());
    else
      value = CGAL::to_double(p2.z());
  }
  // Compute the intensity associated with this value.
  double intensity = 
    (value - min.value()) / (max.value() - min.value());
  intensity = (intensity < 0.0) ? 0.0 : intensity;
  intensity = (intensity > 1.0) ? 1.0 : intensity;
  // Compute the pixel value and then use gamma correction.
  *out << gcv[int(round(intensity * double(depth.value() - 1)))] << " ";
      }
      x += upd;
    }
    y -= upd;
    *out << std::endl;
  }
  *out << std::endl;
  out->flush();  

  // Delete any allocated streams.
  if (inputPath.supplied()) delete in;
  if (outputPath.supplied()) delete out;
  
  // Report on the number of default pixel values.
  if (num_def > 0)
    std::cerr << "Warning: " << num_def 
        << " pixel values defaulted to " 
        << def.value() << std::endl;

  return 0;
}

---