#define PROG_NAME "balls_to_tomogram"
#define PROG_DESC "Create a tomogram of an arbitrary union of balls"
#define PROG_VERS "1.0"

#define balls_to_tomogram_C_COPYRIGHT \
  "Copyright © 2021 by the State University of Campinas (UNICAMP)"

/* Last edited on 2023-11-26 05:57:49 by stolfi */

#define PROG_HELP \
  "  " PROG_NAME " \\\n" \
  "    -voxelSize {VSIZE} \\\n" \
  "    [ -clip {XMIN} {XMAX} {YMIN} {YMAX} {ZMIN} {ZMAX} ] \\\n" \
  "    -smooth {SMRAD} {NPASSES} \\\n" \
  "    [ -modify { clip | drill | splat } {SHAPE} {CTRX} {CTRY} {CTRZ} {SZX} {SZY} {SZZ} ].. \\\n" \
  "    [ -blur {BLRAD} ] \\\n" \
  "    -outPrefix {OPREF} \\\n" \
  "    " argparser_help_info_HELP " \\\n" \
  "    < {INFILE}"

#define PROG_INFO \
  "NAME\n" \
  "  " PROG_NAME " - " PROG_DESC "\n" \
  "\n" \
  "SYNOPSIS\n" \
  PROG_HELP "\n" \
  "\n" \
  "DESCRIPTION\n" \
  "  The program writes a 3D binary tomogram (3D array of boolean voxels)" \
  " containing a model of union-of-balls model that is read from an" \
  " input file. This shape can be optionally smoothed by" \
  " erosion/dilation.  Optionally, the program writes also an" \
  " antialiased (grayscale) version of that binary tomogram.\n" \
  "\n" \
  "INPUT FILE FORMAT\n" \
  "  Each line of the input describes a ball.  It should have the" \
  " format \"{Si} {Xi} {Yi} {Ri}\" where {Si} is an arc length" \
  " value (irrelevant), {(Xi,Yi)} is the center of the ball, and" \
  " {Ri} its radius, all in millimeters. Blank lines and '#' comments" \
  " are ignored. \n" \
  "\n" \
  "OUTPUT FILE FORMAT\n" \
  "  The output file format is as follows:\n" \
  "\n" \
  "  " ppv_array_read_FORMAT_INFO "\n" \
  "\n" \
  "  The array will have three indices ({d = 3}), which are, in" \
  " order, {Z}, {Y}, and {X}.  This means that all voxels of each" \
  " horizontal plane occur in consecutive position, and ditto for" \
  " the voxels in any row parallel to the {X} axis.\n" \
  "\n" \
  "  The binary tomogram has {maxsmp=1} (1 bit per sample)\n" \
  " while the atialiased version, if any, has {maxsmp=255}\n" \
  " (8 bits per voxel).\n" \
  "\n" \
  "OPTIONS\n" \
  "  -voxelSize {VSIZE}\n" \
  "    This mandatory argument defines the size of a voxel in millimeters.\n" \
  "\n" \
  "  -clip {XMIN} {XMAX} {YMIN} {YMAX} {ZMIN} {ZMAX}\n" \
  "    This optional parameter specifies a clipping box for" \
  " the union of balls.  Balls or parts of balls that lie" \
  " outside the indicated coordinate ranges (mm) will be omitted from the model.\n" \
  "\n" \
  "  -smooth {SMRAD} {NPASSES}\n" \
  "    This mandatory argument specifies the radius (in millimeters) of the" \
  " erosion/dilation element that is used to smooth out corners" \
  " and grooves.  The smoothing will consist of the sequence {D E E D}, where" \
  " each {D} is {NPASSES} dilation steps by that element, and" \
  " each {E} is {NPASSES} ersion steps by the same. If either" \
  " parameter is zero, no smoothing is applied.\n" \
  "\n" \
  "    If parts of the object are clipped off, smoothing will usually be incorrect near" \
  " the clipped-off parts.  If clipping is used with smoothing, it is probably wise to expand" \
  " the clipping box beyond the desired area by {SMRAD*NPASSES} mm" \
  " in all directions.\n" \
  "\n" \
  "  -modify {OP} {SHAPE} {CTRX} {CTRY} {CTRZ} {SZX} {SZY} {SZZ}\n" \
  "    This optional argument specifies an additional modifier object to be combined" \
  " with the union of balls, AFTER any smoothing the shape.\n" \
  "\n" \
  "    If the {OP} is \"clip\", the modifier object is intersected with the" \
  " current thing.  If {OP} is \"drill\" the modifier is subtracted" \
  " from it.  If {OP} is \"-splat\", the modifier is united with it.\n" \
  "\n" \
  "    The {SHAPE} currently can be any of 'box', 'ball', 'cylX', 'cylY'," \
  " or 'cylZ'.  In any case, the geometric center of the object will be" \
  " at {(CTRX},{CTRY},CTRZ}).  For 'box', the parameters {SZX}, {SZY}," \
  " and {SZZ} are the dimensions of the box along the three coordinate" \
  " axes {X},{Y}, and {Z} (array axes 2, 1, and 0, respectively).  For 'ball', the" \
  " three dimensions must be equal t the radius.  For 'cylZ', a cylinder" \
  " with axis parallel to the {Z} axis, {SZX} and {SZY} must be equal" \
  " and specify the radius (not diameter) of the cylinder, while {SZZ} specifies" \
  " the length. The options 'cylX' and 'cylY' are analogous.  All these" \
  " coordinates and dimensions are in millimeters.\n" \
  "\n" \
  "  [ -blur {BLRAD} ]\n" \
  "    This optional argument defines the blurring radius to use" \
  " before writing out the tomogram.  The {BLUR-RADIUS} should be a non-negative" \
  " value in mm.  If omitted or 0, the output tomogram" \
  " will be the binary one.  If provided and greater" \
  " than 0, the output image size will be reduced by the" \
  " factor {1/(R+1)} along each axis, and each sample will" \
  " be replaced by a Hahn-weighted average of the binary voxels in a sindow with" \
  " size {2*R + 1} along each axis; where {R} is {BLRAD/VSIZE} rounded to" \
  " an integer.  This will smooth out corners" \
  " and erasse fine detail, but should result" \
  " in better STL meshes.\n" \
  "\n" \
  "  -outPrefix {OPREF}\n" \
  "    This mandatory option specifies the prefix for output file" \
  " names.  The binary tomogram output will be" \
  " named \"{OPREF}-b.tom\" and the antialiased one \"{OPREF}-g.tom\".\n" \
  "\n" \
  "DOCUMENTATION OPTIONS\n" \
  argparser_help_info_HELP_INFO "\n" \
  "\n" \
  "SEE ALSO\n" \
  "  tomo_to_stl(1), salamic(1).\n" \
  "\n" \
  "AUTHOR\n" \
  "  Created 2021-06-05 by Jorge Stolfi, IC-UNICAMP.\n" \
  "\n" \
  "MODIFICATION HISTORY\n" \
  " 2021-06-05 J. Stolfi: Created, by adapting {make_tesla_valve}.\n" \
  " 2021-06-08 J. Stolfi: Reformed to a generic union-of-balls tomogram maker.\n" \
  " 2021-06-11 J. Stolfi: Changed smoothing to be {D^n E^n E^n D^n}.\n" \
  " 2021-06-14 J. Stolfi: Added \"modify\" option.\n" \
  " 2021-06-15 J. Stolfi: Added \"-clip\" option, revised the clipping logic.\n" \
  " 2021-07-04 J. Stolfi: Added \"-blur\" option.\n" \
  " 2021-07-15 J. Stolfi: Added \"-outPrefix\" option, writes both tomograms.\n" \
  "\n" \
  "WARRANTY\n" \
  argparser_help_info_NO_WARRANTY "\n" \
  "\n" \
  "RIGHTS\n" \
  "  " balls_to_tomogram_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 <fget.h>
#include <argparser.h>
#include <r3.h>
#include <r4.h>
#include <r3x3.h>
#include <r3_motion.h>
#include <box.h>
#include <interval.h>
#include <wt_table.h>
#include <wt_table_hann.h>
#include <ppv_array.h>
#include <ppv_array_read.h>
#include <ppv_array_write.h>
#include <ppv_brush.h>
#include <ppv_array_blur.h>

#include <voxb_obj.h>
#include <voxb_splat.h>
#include <voxb_erolate.h>

/* COMMAND-LINE OPTIONS */

typedef struct btt_modifier_t
  { voxb_op_t op;   /* The operation to use for splatting. */
    char* shape;    /* The {SHAPE} argument of "-clip", or {NULL} if not given. */
    r3_t ctr;       /* The {CTRX} {CTRY} {CTRZ} arguments of "-clip" (mm). */
    r3_t sz;        /* The {SZX} {SZY} {SZZ} arguments of "-clip" (mm). */
  } btt_modifier_t;
  /* Parmeters from a "-clip", "-drill", or "-splat" option. */
  
vec_typedef(btt_modifier_vec_t,btt_modifier_vec,btt_modifier_t);

typedef struct btt_options_t
  { double voxelSize;         /* Size of a voxel in millimeters. */
    r3_t pmin;                /* Low point of the clipping box. */
    r3_t pmax;                /* High point of the clipping box. */
    double smooth_radius;     /* Radius of smothing element for erosion/dilation (mm). */
    int32_t smooth_passes;    /* Number of iterations of each erosion/dilation. */
    btt_modifier_vec_t mods;  /* Array of modifiers to apply. */
    double blur;              /* Blurring radius (mm). */
    char *outPrefix;          /* Output filename prefix. */
  } btt_options_t;

#define btt_voxelSize_MIN (0.10)
  /* Min voxel size (mm). */

#define btt_voxelSize_MAX (2.00)
  /* Max voxel  size (mm). */

#define btt_object_size_MIN (3.0)
  /* Min object size on any axis (mm). */

#define btt_object_size_MAX (500.0)
  /* Max object size on any axis (mm). */

#define btt_array_margin (3)
  /* Extra voxels to leave all around object in tomogram. */

#define btt_array_size_MIN (10+2*(btt_array_margin))
  /* Min tomogram size along any axis. */

#define btt_array_size_MAX (10000)
  /* Max tomogram size along any axis. */

#define btt_array_tot_voxels_MAX (1024*1024*(int64_t)1024)
  /* Max tomogram size along any axis. */

#define btt_smooth_radius_MAX (7.0)
  /* Max smoohing radius (vx). */

#define btt_smooth_passes_MAX (10)
  /* Max number of smoothing passes. */

#define btt_blur_MAX (10)
  /* Max blur window radius (vx). */

/* INTERNAL PROTOTYPES */

btt_options_t *btt_parse_options(int32_t argc, char **argv);
  /* Parses the command line arguments and packs them as an {btt_options_t}. */

r4_vec_t btt_read_union_of_balls(FILE *rd, int32_t nskip);
  /* Reads from file {rd} a set of balls. Each ball is specifed by an
    {r4_t} value where the first three coordinates are its center and
    the last one is the radius. In the file, each line should have at
    least {nskip+4} numbers; the first {nskip} ones are ignored. Blank
    lines, and comments are ignored. */

void btt_get_bounding_box(r4_vec_t ball, r3_t *pmin, r3_t *pmax);
  /* On entry, {*pmin} and {*pmax} should be the low and high corners 
    of the clipping box. On exit, they will be the corners of the
    intersection of that box with the bounding box of the union of all
    balls in {ball}, including their radii. All coordinates
    sould be in mm from some arbitrary origin. */

ppv_array_t *btt_allocate_tomogram(r3_t *szobj, ppv_size_t mrg, double vsz);
  /* Allocates a tomogram with enough space for a thing with dimensions {szobj} (mm)
    plus a safety margin of {mrg} voxels all around.  Assumes the
    voxel size is {vsz} (mm).  */

ppv_array_t *btt_crop_tomogram(ppv_array_t *A, ppv_size_t mrg);
  /* Returns an array descriptor that shares the same samples as {A} 
    minus a layer of {mrg} voxels all aroud. */

void btt_splat_balls
  ( ppv_array_t *A, 
    double vsz, 
    r4_vec_t ball, 
    r3_t *shift
  );
  /* Splat into {A} the union of all the balls in {ball}. The balls are implicitly
    shifted by the vector {shift}, rescaled from mm to vx, and shifted again so that the origin
    is at the center of {A}'s domain. */

void btt_smooth_object(ppv_array_t *A, double vsz, double rad, int32_t passes);
  /* Applies smoothing to the object in the tomogram {A} by performing
    morphological closing and opening with a digital ball of radius
    {rad} (mm). Assumes the voxel size is {vsz} (mm). */

void btt_modify_object
  ( ppv_array_t *A, 
    double vsz, 
    btt_modifier_t *m, 
    r3_t *shift
  );
  /* Applies to the object currently in {A} the modifier described by
    {m}, implicitly shifted by {shift}, and converted from mm to vx. */

void btt_write_tomogram(char *outFile, ppv_array_t *A); 
  /* Writes the tomogram to a file named "{outFile}", in the format
    described by {ppv_array_read_FORMAT_INFO}.  */

void btt_splat_ball(ppv_array_t *A, r3_t *ctr, double rad);
  /* Splats into the array {A} a ball with center {ctr} and radius {rad}.
    All coordinates should be in {vx} relative to the low corner of {A}'s domain. */

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

/* IMPLEMENTATIONS */

bool_t voxb_splat_debug = FALSE;

int32_t main(int32_t argc, char** argv)
  {
    voxb_splat_debug = FALSE;
    
    /* ppv_dim_t d = 3; */ /* For now. */
    btt_options_t *o = btt_parse_options(argc, argv);
    
    double vsz = o->voxelSize;
    
    /* Read the union-of-balls model: */
    FILE *rd = stdin;
    int32_t nskip = 1; /* Skip the arc length. */
    r4_vec_t ball = btt_read_union_of_balls(rd, nskip);
    fclose(rd);

    /* Get the bounding box: */
    r3_t pmin = o->pmin;
    r3_t pmax = o->pmax;
    btt_get_bounding_box(ball, &pmin, &pmax);
    
    /* Compute the alloed object size range after clipping: */
    double szobj_arr_MAX = (btt_array_size_MAX - 2*btt_array_margin)*vsz; /* (mm) */
    double szobj_MAX = fmin(szobj_arr_MAX, btt_object_size_MAX); /* (mm) */
    double szobj_arr_MIN = btt_array_size_MIN*vsz; /* (mm) */
    double szobj_MIN = fmax(szobj_arr_MIN, btt_object_size_MIN); /* (mm) */

    /* Compute the object size {szobj} and the {shift} vector: */
    r3_t szobj, shift;
    for (int32_t i = 0; i < 3; i++) 
      { szobj.c[i] = pmax.c[i] - pmin.c[i];
        shift.c[i] = - (pmin.c[i] + pmax.c[i])/2;
        demand(szobj.c[i] >= szobj_MIN, "clipped object too small");
        demand(szobj.c[i] <= szobj_MAX, "clipped object too big");
      }
    r3_gen_print(stderr, &szobj, "%8.3f", "object size = ( ", " ", " )\n");
    r3_gen_print(stderr, &shift, "%8.3f", "shift = ( ", " ", " )\n");
    
    /* Allocate the tomogram with a safety margin and crop it out: */
    ppv_size_t mrg = btt_array_margin;
    ppv_array_t *A = btt_allocate_tomogram(&szobj, mrg, vsz);
    ppv_array_t *Acrop = btt_crop_tomogram(A, mrg);

    /* Splat the balls: */
    btt_splat_balls(Acrop, vsz, ball, &shift);
    
    /* Apply smoothing: */
    if (o->smooth_radius > 0)
       { btt_smooth_object(Acrop, vsz, o->smooth_radius, o->smooth_passes); }
    
    /* Modify object: */
    for (int32_t km = 0; km < o->mods.ne;  km++)
       { btt_modifier_t m = o->mods.e[km];
         btt_modify_object(Acrop, vsz, &m, &shift);
       }

    /* Write it out: */
    char *fname = NULL;
    asprintf(&fname, "%s-b.tom", o->outPrefix);
    btt_write_tomogram(fname, A);
    free(fname);
       
    if (o->blur > 0.0)
      { /* Apply blurring: */
        int32_t radius = (int32_t)ceil(o->blur/vsz);
        int32_t szw = 2*radius +1;
        double wt[szw]; 
        double flat = 0.0; /* No flat top. */
        int32_t stride;
        wt_table_hann_fill(szw, flat, wt, &stride);
        assert(stride == radius + 1);
        fprintf(stderr, "blurring with radius = %d stride = %d (vx)\n", radius, stride);
        ppv_sample_t mB = 255; /* Max sample value of blurred array. */
        ppv_array_t *B = ppv_array_blur(A, NULL, mB, radius, wt, stride, NULL);
        free(A->el); free(A);
        
        /* Write it out: */
        fname = NULL;
        asprintf(&fname, "%s-g.tom", o->outPrefix);
        btt_write_tomogram(fname, B);
      }

    return 0;
  }
  
r4_vec_t btt_read_union_of_balls(FILE *rd, int32_t nskip)
  { r4_vec_t ball = r4_vec_new(1000);
    int32_t nb = 0;   /* Number of balls found. */
    int32_t nlin = 0; /* Number of lines read. */
    while (TRUE)
      {
        fget_skip_spaces(rd);
        if (fget_test_char(rd, EOF)) { break; }
        nlin++;
        if (fget_test_char(rd, '#'))
          { /* Skip to eol, including it: */
            fget_skip_to_eol(rd);
          }
        else if (fget_test_char(rd, '\n'))
          { /* Nothing to do */
          }
        else
          { /* Assume that it is a data line: */
            for (int32_t i = 0; i < nskip; i++) { (void)fget_double(rd); }
            r4_t pr;
            for (int32_t i = 0; i < 4; i++) { pr.c[i] = fget_double(rd); }
            fget_skip_to_eol(rd);
            r4_vec_expand(&ball, nb);
            ball.e[nb] = pr;
            nb++;
          }
      }
    r4_vec_trim(&ball, nb);
    fprintf(stderr, "read %d balls\n", nb);
    return ball;
  }
  
void btt_get_bounding_box(r4_vec_t ball, r3_t *pmin, r3_t *pmax)
  { int32_t nb = ball.ne;
    
    /* Find bounding box {bmin,bmax}: */
    r3_t bmin = (r3_t){{ +INF, +INF, +INF }};
    r3_t bmax = (r3_t){{ -INF, -INF, -INF }};
    for (int32_t kb = 0; kb < nb; kb++)
      { r4_t *ballk = &(ball.e[kb]);
        double rk = ballk->c[3];
        for (int32_t i = 0; i < 3; i++)
          { double xki = ballk->c[i];
            double xki_min = xki - rk; 
            bmin.c[i] = fmin(xki_min, bmin.c[i]);
            double xki_max = xki + rk; 
            bmax.c[i] = fmax(xki_max, bmax.c[i]);
          }
      }
    /* Interset with client's clip box {pmin,pmax}: */
    for (int32_t i = 0; i < 3; i++)
      { pmin->c[i] = fmax(pmin->c[i], bmin.c[i]);
        pmax->c[i] = fmin(pmax->c[i], bmax.c[i]);
        demand(pmin->c[i] < pmax->c[i], "clip box is dsjoint from balls bounding box");
      }
  }

ppv_array_t *btt_allocate_tomogram(r3_t *szobj, ppv_size_t mrg, double vsz)
  {
    ppv_dim_t d = 3;
    ppv_sample_t maxsmp = 1;
    ppv_size_t sz[d];
    int64_t nvt = 1; /* Total number of volxels. */
    fprintf(stderr, "array size =");
    for (ppv_axis_t ax = 0; ax < d; ax++) 
      { /* Get the object size (mm) along tomogram axis {ax}, with space for smoothing: */
        ppv_size_t szk = (ppv_size_t)(ceil(szobj->c[2-ax]/vsz) + 2*(double)mrg);
        fprintf(stderr, " " ppv_size_t_FMT, szk);
        affirm(szk >= btt_array_size_MIN, "voxel array dimension is too small");
        affirm(szk <= btt_array_size_MAX, "voxel array dimension is too big");
        sz[ax] = szk;
        nvt = nvt*szk;
      }
    fprintf(stderr, "\n");
    fprintf(stderr, "total %ld voxels", nvt);
    demand(nvt <= btt_array_tot_voxels_MAX, "voxel array is too big");

    ppv_array_t *A = ppv_array_new(d, sz, maxsmp);

    ppv_nbits_t bps = A->bps;
    int64_t nMb = (int64_t)ceil((double)(nvt*bps)/1000000.0);
    int64_t nMB = (int64_t)ceil((double)(nvt*bps)/8000000.0);
    fprintf(stderr, " (%ld Mb, %ld MB)\n", nMb, nMB);
    return A;
  }

ppv_array_t *btt_crop_tomogram(ppv_array_t *A, ppv_size_t mrg)
  {
    ppv_dim_t d = A->d;
    ppv_array_t *Acrop = ppv_array_clone(A);
    for (ppv_axis_t ax = 0; ax < d; ax++)
      { ppv_size_t keep = A->size[ax] - 2*mrg;
        ppv_crop(Acrop, ax, mrg, keep);
      }
    return Acrop;
  }

void btt_splat_balls
  ( ppv_array_t *A, 
    double vsz, 
    r4_vec_t ball, 
    r3_t *shift
  )
  { 
    /* Corners and centers of {A}'s domain, in vx but {XYZ} order: */
    double NX = (double)A->size[2]; 
    double NY = (double)A->size[1]; 
    double NZ = (double)A->size[0]; 
    r3_t loA = (r3_t){{ 0, 0, 0 }};
    r3_t hiA = (r3_t){{ NX, NY, NZ }};
    r3_t ctrA; r3_mix(0.5, &loA, 0.5, &hiA, &ctrA);

    /* Splat the balls, shifted and scaled: */
    int32_t nb = ball.ne;
    int32_t per = 1;  while (100*per < nb) { per = 10*per; }
    fprintf(stderr, "splatting %d balls ('.' = %d balls)\n", nb, per);
    for (int32_t ib = 0; ib < nb; ib++)
      { if ((ib % per) == 0) { fputs(".", stderr); }
        r4_t pr = ball.e[ib];
        /* Scale ball center coordinates to voxels and shift to center: */
        double xb = (pr.c[0] + shift->c[0])/vsz + ctrA.c[0];
        double yb = (pr.c[1] + shift->c[1])/vsz + ctrA.c[1];
        double zb = (pr.c[2] + shift->c[2])/vsz + ctrA.c[2];
        r3_t ctrb = (r3_t){{ xb, yb, zb }};
        /* Scale the ball radius: */
        double rb = pr.c[3]/vsz;
        /* Check containment: */
        bool_t ok = TRUE; /* Does this box enter the clip box? */
        for (int32_t i = 0; i < 3; i++)
          { if (ctrb.c[i] + rb <= loA.c[i]) { ok = FALSE; }
            if (ctrb.c[i] - rb >= hiA.c[i]) { ok = FALSE; }
          }
        if (ok) { btt_splat_ball(A, &ctrb, rb); }
    } 
    fputs("\n", stderr);
  }
  
void btt_smooth_object(ppv_array_t *A, double vsz, double rad, int32_t passes)
  {
    fprintf(stderr, "smoothing with rad = %.3f passes = %d\n", rad, passes);

    assert(rad > 0);
    double smv = rad/vsz;          /* Erosion/dilation radius of each pass (vx). */

    fprintf(stderr, "smoothing with rad = %.3f mm (%.4f vx) passes = %d\n", rad, smv, passes);

    /* Smooth the shape by erosion/dilation. */
    ppv_brush_t *b = ppv_brush_make_ball(A->d, smv);

    for (int32_t k = 0; k < passes; k++) { voxb_erolate_with_brush(A, b, FALSE); }  /* Dilate. */
    for (int32_t k = 0; k < 2*passes; k++) { voxb_erolate_with_brush(A, b, TRUE); } /* Erode. */
    for (int32_t k = 0; k < passes; k++) { voxb_erolate_with_brush(A, b, FALSE); }  /* Dilate. */

    ppv_brush_free(b);
  }

void btt_modify_object
  ( ppv_array_t *A, 
    double vsz, 
    btt_modifier_t *m, 
    r3_t *shift
  )
  {
    fprintf(stderr, "applying modifier op = %d shape = %s ...\n", m->op, m->shape);
    r3_pred_t *modpred = NULL;  /* Predicate that defines the modifier object. */
    
    /* Corners and centers of {A}'s domain, in vx but {XYZ} order: */
    r3_t loA = (r3_t){{ 0, 0, 0 }};
    r3_t hiA = (r3_t){{ (double)A->size[2], (double)A->size[1], (double)A->size[0] }};
    r3_t ctrA; r3_mix(0.5, &loA, 0.5, &hiA, &ctrA);
    
    /* Convert modifier's center to voxel units {ctrC}, shifted, in {XYZ} order: */
    r3_t ctrC; 
    r3_add(&(m->ctr), shift, &ctrC);
    r3_scale(1/vsz, &ctrC, &ctrC);
    r3_add(&ctrC, &ctrA, &ctrC);
    r3_gen_print(stderr, &ctrC, "%8.3f", "modifier center = ( ", " ", " ) (vx)\n");
   
    /* Convert modifier {sz} to voxel units {szC}, in {XYZ} order: */
    r3_t szC; r3_scale(1/vsz, (&m->sz), &szC);
    r3_gen_print(stderr, &szC, "%8.3f", "modifier size = ( ", " ", " ) (vx)\n");
    
    /* Parameters used by the {modpred_XXX} procs: */
    ppv_axis_t clax;  /* Array axis parallel to cylinder axis. */
    r3_t hbox;        /* Half-size of box (vx, {XYZ} order). */
    double rad;       /* Radius for ball or cylinder (vx). */
    double hlen;      /* Half-length of cylinder (vx). */
    double maxR;      /* Max extent of modifier from its center. */
    
    auto bool_t modpred_box(r3_t *p);
    auto bool_t modpred_ball(r3_t *p);
    auto bool_t modpred_cyl(r3_t *p);
      /* Predicates that define the modifier objects. */
    r3_motion_state_t S;
    S.p = ctrC;
    r3x3_ident(&(S.M));

    if (strcmp(m->shape, "box") == 0)
      { modpred = &modpred_box; 
        r3_scale(0.5, &szC, &hbox); 
        maxR = r3_norm(&hbox);
      }
    else if (strcmp(m->shape, "ball") == 0)
      { modpred = &modpred_ball; 
        rad = szC.c[0]; 
        maxR = rad;
      }
    else if (strcmp(m->shape, "cylX") == 0)
      { modpred = &modpred_cyl; 
        clax = 0; 
        hlen = 0.5*szC.c[0]; 
        rad = szC.c[1]; 
        maxR = hypot(hlen,rad);
      }
    else if (strcmp(m->shape, "cylY") == 0)
      { modpred = &modpred_cyl; 
        clax = 1; 
        hlen = 0.5*szC.c[1]; 
        rad = szC.c[2]; 
        maxR = hypot(hlen,rad);
      }
    else if (strcmp(m->shape, "cylZ") == 0)
      { modpred = &modpred_cyl; 
        clax = 2; 
        hlen = 0.5*szC.c[2]; 
        rad = szC.c[0]; 
        maxR = hypot(hlen,rad);
      }
    else
      { assert(FALSE); }
      
    for (int32_t i = 0; i < 3; i++)
       { if ((ctrC.c[i] + maxR <= loA.c[i]) ||  (ctrC.c[i] - maxR >= hiA.c[i]))
           { fprintf(stderr, "modifier entirely outside domain, skipped\n");
             return;
           }
       }

    voxb_splat_object(A, modpred, &S, maxR+1.0, m->op, (m->op != voxb_op_AND));
    
    return;
    
    bool_t modpred_box(r3_t *p)
      { for (int32_t i = 0; i < 3; i++) 
          { if (fabs(p->c[i]) > hbox.c[i]) { return FALSE; } }
        return TRUE;
      }
      
    bool_t modpred_ball(r3_t *p)
      { return r3_norm_sqr(p) <= rad*rad; }
      
    bool_t modpred_cyl(r3_t *p)
      { double d2 = 0; /* Distance from axis, squared. */
        for (int32_t i = 0; i < 3; i++) 
          { double di = fabs(p->c[i]);
            if (i == clax)
              { if (di > hlen) { return FALSE; } }
            else
              { if (di > rad) { return FALSE; }
                d2 += di*di;
              }
          }
        return d2 <= rad*rad;
      }
  }
    
void btt_splat_ball(ppv_array_t *A, r3_t *ctr, double rad)
  { 
    auto bool_t ball(r3_t *p);
      /* Indicator function for the ball. */
       
    /* Position and orientation of ball: */
    r3_motion_state_t S;
    S.p = (*ctr);
    r3x3_ident(&(S.M));

    voxb_splat_debug = TRUE;
    voxb_splat_object(A, ball, &S, rad, voxb_op_OR, FALSE);
    voxb_splat_debug = FALSE;

    return;
      
    bool_t ball(r3_t *p)
      { return voxb_obj_ball(p, rad); }
  }
  
void btt_write_tomogram(char *outFile, ppv_array_t *A)
  {
    FILE *wr = open_write(outFile, TRUE);

    bool_t plain = FALSE;
    ppv_array_write_file(wr, A, plain);
    
    fclose(wr); 
  }

btt_options_t *btt_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: */
    btt_options_t *o = (btt_options_t *)malloc(sizeof(btt_options_t)); 
    
    /* Parse keyword parameters: */
    
    argparser_get_keyword(pp, "-voxelSize");
    o->voxelSize = argparser_get_next_double(pp, btt_voxelSize_MIN, btt_voxelSize_MAX);

    if (argparser_keyword_present(pp, "-clip"))
      { for (int32_t i = 0; i < 3; i++)
          { double coord_MAX = btt_object_size_MAX/2;
            o->pmin.c[i] = argparser_get_next_double(pp, -coord_MAX, +coord_MAX);
            o->pmax.c[i] = argparser_get_next_double(pp, o->pmin.c[i]+1.0, +coord_MAX);
          }
      }
    else
      { o->pmin = (r3_t){{ -INF, -INF, -INF }};
        o->pmax = (r3_t){{ +INF, +INF, +INF }};
      }
    
    argparser_get_keyword(pp, "-smooth");
    o->smooth_radius = argparser_get_next_double(pp, 0, btt_smooth_radius_MAX);
    o->smooth_passes = (int32_t)argparser_get_next_int(pp, 0, btt_smooth_passes_MAX);
    
    int32_t nmods = 0; /* Number of modifiers specified. */
    o->mods = btt_modifier_vec_new(10);
    while (argparser_keyword_present(pp, "-modify")) 
      { btt_modifier_t m;
        char *op = argparser_get_next_non_keyword(pp);
        if (strcmp(op, "clip") == 0)
          { m.op = voxb_op_AND; }
        else if (strcmp(op, "drill") == 0)
          { m.op = voxb_op_SUB; }
        else if (strcmp(op, "splat") == 0)
          { m.op = voxb_op_OR; }
        else
          { argparser_error(pp, "invalid \"-modify\" operation"); }
        fprintf(stderr, "op = %s %d\n", op, m.op);
        m.shape = argparser_get_next_non_keyword(pp);
        for (int32_t i = 0; i < 3; i++) 
          { m.ctr.c[i] = argparser_get_next_double(pp, -1000.0, 1000.0); }
        for (int32_t i = 0; i < 3; i++) 
          { m.sz.c[i] = argparser_get_next_double(pp, 1.0, 1000.0); }
        /* Consistency of sizes: */
        if (strcmp(m.shape, "cylZ") == 0) 
          { demand(m.sz.c[0] == m.sz.c[1], "bad 'cylZ' clip size"); }
        else if (strcmp(m.shape, "cylX") == 0) 
          { demand(m.sz.c[1] == m.sz.c[2], "bad 'cylX' clip size"); }
        else if (strcmp(m.shape, "cylY") == 0) 
          { demand(m.sz.c[0] == m.sz.c[2], "bad 'cylY' clip size"); }
        else  if (strcmp(m.shape, "ball") == 0) 
          { demand(m.sz.c[0] == m.sz.c[1], "bad 'ball' clip size");
            demand(m.sz.c[1] == m.sz.c[2], "bad 'ball' clip size");
          }
        else
          { argparser_error(pp, "invalid \"-modify\" shape"); }
         btt_modifier_vec_expand(&(o->mods), nmods);
         o->mods.e[nmods] = m;
         nmods++;
      }
    btt_modifier_vec_trim(&(o->mods), nmods);

    if (argparser_keyword_present(pp, "-blur"))
      { o->blur = argparser_get_next_double(pp, 0, btt_blur_MAX*o->voxelSize); }
    else
      { o->blur = 0.0; }
    
    argparser_get_keyword(pp, "-outPrefix");
    o->outPrefix = argparser_get_next_non_keyword(pp);

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

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

vec_typeimpl(btt_modifier_vec_t,btt_modifier_vec,btt_modifier_t);