#define PROG_NAME "raspol"
#define PROG_DESC "Converts an unstructured polygonal drawing to a raster"
#define PROG_VERS "1.0"

#define raspol_C_COPYRIGHT \
  "Copyright © 2016 by the State University of Campinas (UNICAMP) and" \
  " the Federal Technological University of Parana (UTFPR)"

/* Last edited on 2021-07-09 01:06:00 by jstolfi */

#define PROG_HELP \
  "  " PROG_NAME " \\\n" \
  "    -polyFile {POLY_FILE} \\\n" \
  "    -format { ascii | binary } \\\n" \
  "    -pixelSize {PIXEL_SIZE} \\\n" \
  "    -start {X_LO} {Y_LO} \\\n" \
  "    -stop {X_HI} {Y_HI} \\\n" \
  "    [ -neGuess {GUESSED_NUM_SEGS} ] \\\n" \
  "    " argparser_help_info_HELP " \\\n" \
  "    > {OUTFILE}"
  
#define PROG_INFO \
  "NAME\n" \
  "  " PROG_NAME " - " PROG_DESC "\n" \
  "\n" \
  "SYNOPSIS\n" \
  PROG_HELP "\n" \
  "\n" \
  "DESCRIPTION\n" \
  "  The program reads from {POLY_FILE} a list of line segments" \
  " that are supposed to be the outline of a polygonal figure, and" \
  " outputs a bilevel raster image where each pixel is 1 or 0 depending" \
  " on whether its center is inside or outside the figure, respectively.\n" \
  "\n" \
  "  The figure may have several disconnected components of arbitrarily" \
  " complicated shape, and each component may have zero or more" \
  " holes.  Components nested in the holes of other components, many" \
  " levels deep.  The segments may cross or partly overlap each" \
  " other.  The figure is rasterized only within the 'region of" \
  " interest', a user-specifed axiis-aligned rectangle of the plane.\n" \
  "\n" \
  "INPUT\n" \
  "  The input file must be in a" \
  " special 'STereolitography Poligon' (STP) format.  There are two variants, ASCII and binary.\n" \
  "\n" \
  "  " stpoly_STP_format_ascii_INFO "\n" \
  "\n" \
  "  " stpoly_STP_format_binary_INFO "\n" \
  "\n" \
  "OUTPUT\n" \
  "  The program writes to standard output a binary image file in" \
  " PNG format with 1 bit per pixel, that is assumed to represent" \
  " the axis-aligned rectangle of the plane with" \
  " corners  {P_LO = (X_LO,Y_LO)} and {P_HI = (X_HI,Y_HI)}.  More" \
  " precisely, each pixel is assumed to represent a square with" \
  " {PIXEL_SIZE} millimeters on each side.  The" \
  " image then represents a mesh of such squares, whose inferior" \
  " corner is at the point with coordinates {(X_LO,Y_LO)} in" \
  " millimeters, rounded down to an integer multiple of {PIXEL_SIZE}.  The" \
  " number of columns {NX} and the number of rows {NY} will be such" \
  " that the mesh covers completely the rectangle with corners {P_LO} and {P_HI}.\n" \
  "\n" \
  "  Each pixel of the resulting image will be set to 1 or 0 depending" \
  " on whether its center is inside or outside the figure," \
  " respectively.  By definition, a point {p} is inside the" \
  " figure if a ray {R(p)} starting at {p} and directed horizontally" \
  " towards the {-X} direction crosses transversally an odd number" \
  " of segments.\n" \
  "\n" \
  "  In order to avoid ambiguous cases, all segment endpoint" \
  " coordinates from the input file will be rounded" \
  " to even integer multiples of {EPS = PIXEL_SIZE/10} before further" \
  " processing.  Then the the coordinates of pixel centers will" \
  " be odd multiples of {EPS}.  This expedient ensures that the" \
  " ray {R(p)} of any pixel center {p} will not go through any" \
  " vertex of the figure, and will not be collinear with" \
  " any of the segments.\n" \
  "\n" \
  "OPTIONS\n" \
  "  -polyFile {POLY_FILE}\n" \
  "    This mandatory argument specifies the nameof the input STP file.\n" \
  "\n" \
  "  -format ascii\n" \
  "  -format binary\n" \
  "    This mandatory argument specifies the format of the input STP file.\n" \
  "\n" \
  "  -neGuess {GUESSED_NUM_SEGS} \n" \
  "    This optional argument is a hint for the number of" \
  " segments that are expected to be in the input file.  It" \
  " is used by the program to preallocate" \
  " internal tables, so that they don't have to be expanded dynamically while" \
  " the file is being read.  The program is more efficient if {GUESSED_NUM_SEGS} is the" \
  " actual number of segments, or somewhat higher.  If omitted, the" \
  " program assumes \"-neGuess " stringify(raspol_neGuess_DEFAULT) "\".\n" \
  "\n" \
  "  -pixelSize {PIXEL_SIZE}\n" \
  "    This mandatory argument specifies the size of pixels of the output image" \
  " (in millimeters).  It had better be a fraction of the desired" \
  " accuracy, so that the rasterization artifacts are too small to" \
  " matter.  Note however that the size of the output image has" \
  " approximately {X_HI/PIXEL_SIZE} pixel columns and {Y_HI/PIXEL_SIZE} rows.  If" \
  " the desired area is {SZ} millimiters wide, then {PIXEL_SIZE} had" \
  " better be larger than {SZ/65535}.\n" \
  "\n" \
  "  -start {X_LO} {Y_LO}\n" \
  "  -stop {X_HI} {Y_HI}\n" \
  "    These mandatory arguments specify the coordinates of the inferior" \
  " corner {(X_LO,Y_LO)} and the superior corner {(X_HI,Y_HI)} of" \
  " the region of interest, all in millimeters.\n" \
  "\n" \
  "DOCUMENTATION OPTIONS\n" \
  argparser_help_info_HELP_INFO "\n" \
  "\n" \
  "SEE ALSO\n" \
  "  ImageMagick(1).\n" \
  "\n" \
  "AUTHOR\n" \
  "  Created 2016-04-10 by Jorge Stolfi, IC-UNICAMP.\n" \
  "\n" \
  "MODIFICATION HISTORY\n" \
  " 2016-04-13 J. Stolfi: Created.\n" \
  "\n" \
  "WARRANTY\n" \
  argparser_help_info_NO_WARRANTY "\n" \
  "\n" \
  "RIGHTS\n" \
  "  " raspol_C_COPYRIGHT ".\n" \
  "\n" \
  argparser_help_info_STANDARD_RIGHTS

#define stringify(x) strngf(x)
#define strngf(x) #x
  /* Hack to include non-string defines in strings.  Don't ask why 2 levels. */ 

#define _GNU_SOURCE
#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <math.h>

#include <bool.h>
#include <affirm.h>
#include <jsfile.h>
#include <jsmath.h>
#include <cmp.h>
#include <argparser.h>

#include <r2.h>

#include <ppv_array.h>
#include <ppv_array_read.h>
#include <ppv_array_write.h>

#include <jspnm.h>
#include <uint16_image.h>
#include <uint16_image_write_png.h>

#include <stpoly.h>
#include <stpoly_STP.h>

#include <raspol_stats.h>
#include <raspol_slicer.h>

/* COMMAND-LINE OPTIONS */

typedef struct raspol_options_t
  { char *polyFile;  /* Name of the STP input file. */
    bool_t binary;   /* FALSE if ASCII STP format, TRUE if binary STP format. */
    double pixelSize;     /* Pixel size (mm). */
    r2_t start;      /* Inferior corner of region of interest (mm). */
    r2_t stop;       /* Superior corner of region of interest (mm).  */
    int neGuess;     /* Estimated number of segments in input file. */
  } raspol_options_t; 

#define raspol_image_size_MAX 65356
  /* Maximum image size along any axis (2^16), for sanity checks. */

/* INTERNAL PROTOTYPES */

typedef struct raspol_vert_t 
  { i2_t p;  /* Quantized coordinates of a vertex. */
  } raspol_vert_t;
  /* A quantized vertex record. */

typedef struct raspol_seg_t 
  { int32_t end[2]; /* Vertex indices. */ 
  } raspol_seg_t;
  /* A segment in a polygonal figure, represented by the indices of its two vertices. */

void raspol_fill_row(ppv_array_t *img, int32_t iy, int32_t na, int32_t ixa[]);
  /* Sets the pixels of row {iy} of the image {img} according to the 
    edge intersections described by {ixa[0..na-1]}.
    The elements of {ixa} must be sorted in increasing order.
    
    Specifically, let /the scanline/ be the horizontal line that goes
    through the centers of the pixels of row {iy}. Assumes that for each
    {k} in {0..na-1}, there is an edge of the figure that crosses the
    scanline strictly between the centers of the pixels {ix-1} and {ix},
    where {ix=ixa[k]}. Each pixel of the row is then 1 if and only if
    there is an odd number of elements of {ixa} that are less than or
    equal to {ix}. */

void raspol_write_image(FILE *wr, ppv_array_t *img);
  /* Writes the array {img} to {wr} as a bilevel PNG image. */

raspol_options_t *raspol_parse_options(int32_t argc, char **argv);

int32_t main(int32_t argc,char** argv);

/* IMPLEMENTATIONS */

int32_t main(int32_t argc, char** argv)
  {
    raspol_options_t *o = raspol_parse_options(argc, argv);
    
    /* Choose the pixel size in multiples of the quantization unit: */
    int32_t img_psz = 10;   /* Pixel size in quantization units. */
    assert((img_psz % 4) == 2); /* So that pixel centers are odd. */

    /* Choose the quantization unit {eps}: */
    double ulp = 1.0e-6;    /* Round eps to a multiple of {ulp} to reduce rounding when printing. */
    double eps = ulp*floor(o->pixelSize/img_psz/ulp + 0.5); /* Quantization unit of length. */

    double img_psz_mm = eps*img_psz;    /* A more accurate value for the pixel size. */
    fprintf(stderr, "specified pixel size (%.8f) rounded to %.8f mm\n", o->pixelSize, img_psz_mm);
    fprintf(stderr, "fundamental length unit = %.8f mm (1/%d of pixel size)\n", eps, img_psz);
    
    /* Quantize the corners of the region of interest, rounding to multiples of {img_psz}: */
    /* Note that decimal fractions are usually inexact {double} values. */
    /* Still, the user should be able to specify the exact start and stop. */

    i2_t img_start; /* Quantized inferior corner in {eps} units. */
    i2_t img_stop;  /* Quantized superior corner in {eps} units. */
    i2_t img_size;  /* Size of image in pixels. */
    for (int j = 0; j < 2; j++)
      { 
        /* Compute image bounds and size along axiz {j}: */
        img_start.c[j] = (int32_t)iroundfrac(o->start.c[j], eps, img_psz, 0, INT32_MAX); 
        img_stop.c[j] = (int32_t)iroundfrac(o->stop.c[j], eps, img_psz, 0, INT32_MAX); 
        img_size.c[j] = (img_stop.c[j] - img_start.c[j])/img_psz;
        
        /* Paranoia: */
        assert(img_start.c[j] % img_psz == 0);
        assert(img_stop.c[j] % img_psz == 0);
        assert(img_size.c[j] > 0);
        
        /* Save the quantized coords and recompute the float ones: */
        double img_size_mm = img_psz_mm*img_size.c[j]; 
        double img_start_mm = img_psz_mm*img_start.c[j];
        double img_stop_mm = img_psz_mm*img_stop.c[j];
        fprintf(stderr, "image spans %.6f mm (%d pixels) in %c", img_size_mm, img_size.c[j], "XY"[j]); 
        fprintf(stderr, " from %.6f to %.6f mm", img_start_mm, img_stop_mm);
        fprintf(stderr, " (pixel grid coords %d to %d)", img_start.c[j], img_stop.c[j]);
        fprintf(stderr, "\n");
        demand (img_size.c[j] <= raspol_image_size_MAX, "too many pixel rows or columns");
      }
    
    /* Allocate the image (1 bit per pixel, 32 pixels per word): */
    /* Index 0 is {X} (to the left), index 1 is {Y} (up). */
    ppv_dim_t d = 2;
    ppv_size_t sz[d];
    for (int ia = 0; ia < d; ia++) { sz[ia] = (ia < 2 ? (ppv_size_t)img_size.c[ia] : 1); }
    ppvsample_t maxsmp = 1;
    ppv_array_t *img = ppv_array_new(d, sz, maxsmp);
    
    /* Read the polygonal outline, rounding the vertices to even multiples of {eps}: */
    bool_t even = TRUE; /* Round all vertex coordinates to even multiples of {eps}. */
    stpoly_t poly = stpoly_read_STP(o->polyFile, o->binary, (float)eps, o->neGuess, even);
    uint32_t ne = stpoly_edge_count(poly);

    /* Process it: */
    int32_t *ixa = notnull(malloc(ne*sizeof(int32_t)), "no mem"); /* Work area for {process_slice}. */
    
    auto void process_slice(int32_t Y, uint32_t na, stpoly_edge_t a[]);
      /* Called by {raspol_slicer_slice} to process each slice
        by a slicing line at ordinate {Y}, that crosses the edges
        {a[0..na-1]}. */ 
    
    /* Choose the ordinates of the slicing lines: */
    int32_t nl = img_size.c[1]; /* One slice per image row. */
    int32_t startY = img_start.c[1] + img_psz/2; /* Start at pixel centers. */
    int32_t stepY = img_psz;    
    raspol_slicer_slice(poly, nl, startY, stepY, process_slice);
    
    /* Write the image out: */
    raspol_write_image(stdout, img);
    
    free(img->el);
    stpoly_free(poly);

    return 0;
    
    /* INTERNAL IMPLEMENTATIONS */

    void process_slice(int32_t Y, uint32_t na, stpoly_edge_t a[])
      { 
        int32_t NX = img_size.c[0];
        
        auto int cmp32(const void *a, const void *b);
          /* Compares the two {int32_t} pointed to by {a} and {b}. */
        
        /* Store the quantized {X}-coordinates of crossings in {ix[0..na-1]}: */
        assert(na <= ne);
        assert((img_psz % 2) == 0);
        int ia;
        for (ia = 0; ia < na; ia++)
          { /* Compute the {X} of the intersection of edge {a[ia]}: */
            double aX = stpoly_edge_line_intersection(a[ia], Y, eps);
            
            /* Quantize to even multiple of {eps}, since pixel centers are odd: */
            int32_t qX = (int32_t)iroundfrac(aX, eps, 2, 0, INT32_MAX);
            
            /* Convert {qX} to pixel displacement from boredr: */
            int32_t dX = qX - img_start.c[0] + img_psz/2;
            assert((dX % 2) == 1);
            int32_t ix = dX/img_psz;
            if (ix < 0) { ix = 0; }
            if (ix > NX) { ix = NX; }
            ixa[ia] = ix;
          }
        
        /* Sort the intersections by increasing {X}: */
        qsort(ixa, na, sizeof(int32_t), cmp32);
        
        /* Fill the array row: */
        int32_t iy = (Y - startY)/stepY; /* Row index. */
        raspol_fill_row(img, iy, na, ixa);
        
        return;
        
        int cmp32(const void *a, const void *b)
          { return cmp_int32((int32_t *)a, (int32_t *)b); }
        
      }
      
      
  }

void raspol_fill_row(ppv_array_t *img, int32_t iy, int32_t na, int32_t ixa[])
  { 
    int32_t NX = (int32_t)img->size[0];
    int32_t NY = (int32_t)img->size[1];
    assert((iy >= 0) && (iy < NY));
    
    int32_t d = img->d; /* Axes in a {ppv-array_t}. */
    assert(d == 2);
    ppv_index_t imgix[d];
    for (int k = 0; k < d; k++) { imgix[k] = 0; }
    int32_t ia = 0; /* Next intersection to consider is {ixa[ia]}. */
    int32_t ix = 0; /* Current pixel column. */
    ppv_sample_t smp = 0; /* Current pixel value. */
    while (ix <= NX)
      { /* Skip all intersections before pixel {ix}: */
        while ((ia < na) && (ixa[ia] <= ix)) { ia++; smp = 1 - smp; }
        if (ix >= NX) { break; }
        assert(ia < na);
        int32_t ix_next = ixa[ia]; /* Column of next interesting pixel. */
        if (smp > 0)
          { /* Fill pixels from {ix} to {ix_next-1} with 1: */
            imgix[0] = ix;
            imgix[1] = iy;
            ppv_pos_t vpos = ppv_sample_pos(img, imgix);
            while (ix < ix_next)
              { ppv_set_sample_at_pos(img, img->bps, img->bpw, vpos, smp);
                ix++;
                vpos += img->step[0];
              }
          }
        else
          { /* Skip pixels from {ix} to {ix_next-1}: */
            ix = ix_next;
          }
      }
    assert(ia == na);
    /* Check for closure: */
    if (smp > 0)
      { fprintf(stderr, "%s: row %d: polygon is not closed\n", __FUNCTION__, iy); }
  }

void raspol_write_image(FILE *wr, ppv_array_t *img)
  {
    bool_t plain = TRUE;
    ppv_array_write_file(wr, img, plain);
  }

#define raspol_neGuess_DEFAULT 100000
  /* Maximum guess for the number of faces in the mesh. */

#define raspol_pixelSize_MIN 0.00001
#define raspol_pixelSize_MAX 10.0
  /* Max and min pixel sizes, for safety. */

#define raspol_coord_MIN -1000.0
#define raspol_coord_MAX +1000.0
  /* Max and min coordinates of a point, for safety. */
  
#define raspol_neGuess_DEFAULT 100000

#define raspol_ne_MAX stpoly_ne_MAX
  /* Default guess for number of segments in outline. */

raspol_options_t *raspol_parse_options(int32_t argc, char **argv)
  {
    /* Initialize argument parser: */
    argparser_t *pp = argparser_new(stderr, argc, argv);
    argparser_set_help(pp, PROG_NAME " version " PROG_VERS ", usage:\n" PROG_HELP);
    argparser_set_info(pp, PROG_INFO);
    argparser_process_help_info_options(pp);

    /* Allocate the command line argument record: */
    raspol_options_t *o = (raspol_options_t *)malloc(sizeof(raspol_options_t));

    /* Parse keyword parameters: */

    argparser_get_keyword(pp, "-polyFile");
    o->polyFile = argparser_get_next(pp);
    
    argparser_get_keyword(pp, "-format");
    if (argparser_keyword_present_next(pp, "ascii"))
      { o->binary = FALSE; }
    else if (argparser_keyword_present_next(pp, "binary"))
      { o->binary = TRUE; }
    else
      { argparser_error(pp, "unrecognized file format"); }
      
    argparser_get_keyword(pp, "-pixelSize");
    o->pixelSize = argparser_get_next_double(pp, raspol_pixelSize_MIN, raspol_pixelSize_MAX);

    argparser_get_keyword(pp, "-start");
    o->start.c[0] = argparser_get_next_double(pp, raspol_coord_MIN, raspol_coord_MAX);
    o->start.c[1] = argparser_get_next_double(pp, raspol_coord_MIN, raspol_coord_MAX);

    argparser_get_keyword(pp, "-stop");
    o->stop.c[0] = argparser_get_next_double(pp, o->start.c[0] + o->pixelSize, raspol_coord_MAX);
    o->stop.c[1] = argparser_get_next_double(pp, o->start.c[1] + o->pixelSize, raspol_coord_MAX);

    if (argparser_keyword_present(pp, "-neGuess"))
      { o->neGuess = (int)argparser_get_next_int(pp, 1, raspol_ne_MAX); }
    else
      { o->neGuess = raspol_neGuess_DEFAULT; }

    /* Parse positional arguments: */
    argparser_skip_parsed(pp);

    /* Check for spurious arguments: */
    argparser_finish(pp);

    return o;
  }