multifok_scene_t *mfmi_make_scene
  ( interval_t dom[],
    char *sceneType,
    double pixelSize,
    double zDep,
    bool_t verbose
  );
  /* Generates a test scene with some number of objects, roughly
    spanning the box {dom[0]×dom[1]×dom[2]} in scene coordinates.
    
    If {sceneType} is "R", the scene have only a floor object
    of type {ot_RAMP}. See {mfmi_fill_scene_empty}.
    
    If {sceneType} is "D", "B", "C", or "P", the scene will have a flat floor (an
    object of type {ot_FLAT}) at {Z=0.0} plus a single foreground object
    of type {ot_DISK}, {ot_BALL}, {ot_CONE}, or {ot_PYRA}, respectively. 
    See {mfmi_fill_scene_singu}.
    
    If {sceneType} is "F" or "T", the scene will have a flat floor at
    {Z=0.0}, and there will be some number of foreground objects at
    random {XY} coordinates. The parameters {pixelSize} and {zDep} are
    relevant only for this kind of scene. See {mfmi_fill_scene_multi}.
    The {overlap} flag will be true if {sceneType} is "T", and false if
    {sceneType} is "F". */

void mfmi_fill_scene_empty
  ( multifok_scene_t *scene,
    multifok_scene_object_type_t type,
    bool_t verbose
  );
  /* The {scene} must have no objects.  Adds a floor object of the
    specified type ({ot_FLAT} or {ot_RAMP}). */

void mfmi_fill_scene_singu
  ( multifok_scene_t *scene,
    multifok_scene_object_type_t type,
    bool_t verbose
  );
  /* The {scene} must have no objects. Adds a floor object of the type
    {ot_FLAT} at {Z=0} and a single foreground object of the specified type
    (which must not be {ot_FLAT} or {ot_RAMP}). The object will be
    centered in the scene's domain {scene.dom}. */
  
void mfmi_fill_scene_multi
  ( multifok_scene_t *scene,
    double pixelSize,
    bool_t overlap,
    double zDep,
    bool_t verbose
  );
  /* The {scene} must have no objects.  Adds a floor object of the
    type {ot_FLAT} and zero or more foreground objects of random
    types, sizes, and positions.
    
    The {Z} spans of all these foreground objects will be strictly
    contained in the {Z} range {scene.dom[2]}. 
    
    If the {overlap} flag is true, the {XY} projections of the
    foreground objects will be disjoint and contained in the {XY}
    projection of {scene.dom}, even when out of focus (assuming that the
    depth of focus is {zDep}). Otherwise the {XY} projections of the
    foreground objects will probably overlap and extend outside the {XY}
    projection of {scene.dom}; only their centers will be in that
    rectangle.
    
    The {pixelSize} parameter is the size of one pixel in SCUs. 
    It affects the min size of objects. */

multifok_stack_t *mfmi_make_and_write_stack
  ( char *stackFolder,
    int32_t NX,
    int32_t NY, 
    multifok_scene_t *scene,
    char *patternFile,
    r3_t *lightVec,
    double ambient,
    double zFoc_min,
    double zFoc_max,
    double zFoc_step,
    double zDep,
    uint32_t HS,
    uint32_t KR,
    uint32_t NQ,
    i2_t iDeb[]
  );
  /* Makes a stack of {NI} frames {fr[0..NI-1]} with simulated views of the {scene},
    where {NI} is {o->NI + 1}.  
    
    The scene is texturized with a monochrome image to be read
    from the given  {patternFile}.
    
    The first {NI-1} frames {fr[0..NI-2]} will be blurred, with the {Z}
    coordinate {fr[i].zFoc} of the in-focus plane varying from
    {zFoc_min} by {zFoc_step}, up to but not exceeding {zFoc_max}. The
    depth of focus {fr[i].zDep} wil be {zDep}. For the meaning of {HS}
    and {KR}, see {multifok_scene_make_frame}.

    The last frame {fr[NI-1]} will be sharp, with {fr[NI-1].zDep = +INF}
    and {fr[NI-1].zFoc} arbitrary A single vertical ray will be used for
    each image sampling point.
    
    Also writes debuging information about the rays used to compute the
    pixels {iDeb[0..NQ-1]}, to files called
    "{stackFolder}/{frameFolder}/pixel-rays-{XXXX}-{YYYY}.txt", where
    {XXXX,YYYY} are the pixel indices formatted as "%04d". Here
    {stackFolder} is as explained in {PROG_INFO}. */

multifok_stack_t *mfmi_make_and_write_stack_from_pattern_function
  ( char *stackFolder,
    int32_t NX,
    int32_t NY,
    multifok_scene_t *scene, 
    multifok_pattern_t *pattern,
    r3_t *lightVec,
    double ambient,
    double zFoc_min,
    double zFoc_max,
    double zFoc_step,
    double zDep,
    uint32_t HS,
    uint32_t KR,
    uint32_t NQ,
    i2_t iDeb[] 
  );
  /* Same as {mfmi_make_stack}, execpt that the scene is colorized
    with the given pattern function. */

multifok_frame_t *mfmi_make_and_write_frame
  ( char *frameFolder,
    int32_t NX, 
    int32_t NY, 
    multifok_scene_t *scene,
    multifok_scene_tree_t *tree,
    multifok_pattern_t *pattern,
    r3_t *lightVec,
    double ambient,
    double zFoc,
    double zDep, 
    uint32_t HS,
    uint32_t KR, 
    bool_t verbose,
    uint32_t NQ,
    i2_t iDeb[]
  );
  /* Creates test images with size {NX,NY} showing the color, nominal
    blurring indicator, scene {Z} average, and scene {Z} deviation at each pixel
    for the given {scene}, as described in {multifok_scene_images_make}.
    
    As a special case, if {zDep} is infinite then the scene view {sVal}
    will be sharp everywhere, and {shrp} will be 1.0.
    
    Writes the images to "{frameFolder}/{tag}.png" were {tag} is "sVal",
    "hAvg", "hDev", "sNrm", and "shrp". Also writes pixel ray data files
    "{frameFolder}/pixel-rays-{XXXX}-{YYYY}.txt" */

void mfmi_write_pixel_profiles
  ( multifok_stack_t *stack,
    uint32_t NQ,
    i2_t iDeb[],
    multifok_scene_t *scene,
    char *stackFolder
  );
  /* Given a list of {NQ} of pixels {iDeb[0..NQ-1]} in the 
    image domain, writes a file "{stackFolder}/pixel-data-{XXXX}-{YYYY}.txt"
    for each pixel, where {XXXX} an {YYYY} are the pixel column and row indices.
    
    Each file will have up to {NI = stack.NI} lines with format

      "{ki} {zFoc} {zDep}  {hAvg} {hDev} {shrp} {sNrm.x} {sNrm.y} {sNrm.z} {sVal[0]} ...  {sVal[NC-1]}"

    where {ki} is a frame number in {0..NI-1} and {hAvg}, and 
    {hDev}, {shrp}, {sNrm}, and {sVal[0..NC-1]} are the values of
    that pixel in the images {hAvg} {hDev}, {shrp}, and {sNrm}, and 
    {sVal} of {stack.frame[ki]}. Note that {zDep} may be {+INF}.  */
    
void mfmi_select_debug_pixels
  ( int32_t NX,
    int32_t NY,
    multifok_scene_t *scene,
    uint32_t NQ_max,
    uint32_t *NQ_P,
    i2_t **iDeb_P
  );
  /* Selects a certain number {NQ <= NQ_max} of pixel index pairs {iDeb[0..NQ-1]}.
    Returns {NQ} and {iDeb} in {*NQ_P} and {ideb_P}. 
    
    The pixels will be preferrably selected over the scene's objects, if there 
    are any. */

float_image_t *mfmi_read_pattern_image(char *fname);
  /* Reads an image from file {fname}, to be used as pattern to paint the scene.
    If the image is in color, converts it to grayscale. */
    
FILE *mfmi_open_pixel_plot_data_file(char *stackFolder);
  /* Opens the file "{stackFolder}/pixplot.txt" for writing. */

          
r3_t morn_expected_normal(multifok_scene_object_t *obj, r3_t *pHit, r3_t *nHit)
  { 
    r3_t cObj, rObj;
    for (int32_t j = 0; j < 3; j++)
      { interval_mid_rad(&(obj->bbox[j]), &(cObj.c[j]), &(rObj.c[j])); }
    
    r3_t nExp;
    switch (obj->type)
      { case ot_DISK:
        case ot_FLAT:
          { /* Normal should be straight up: */
            nExp = (r3_t){{ 0.0, 0.0, 1.0 }};
          }
          break;

        case ot_RAMP:
          { /* Normal should point up or left but unclear how much: */
            /* Just check that {nHit} has the right shape: */
            nExp = (r3_t){{ -fabs(nHit->c[0]), 0.0, +fabs(nHit->c[2]) }};
            (void)r3_dir(&nExp, &nExp);
          }
          break;

        case ot_BALL:
          { /* Normal should be dir of {pHit-cObj}: */
            r3_t uHit; r3_sub(pHit, &cObj, &uHit);
            (void)r3_dir(&uHit, &nExp);
          }
          break;

        case ot_CONE:
          { /* Normal should be along {pHit-cObj} in {X,Y}: */
            r3_t aObj = (r3_t){{ cObj.c[0], cObj.c[1], cObj.c[2] + rObj.c[2] }}; /* Apex. */
            r3_t uHit; r3_sub(pHit, &aObj, &uHit);
            double r = hypot(uHit.c[0], uHit.c[1]);
            assert(uHit.c[2] < 0);
            double h = -uHit.c[2];
            nExp = (r3_t){{ uHit.c[0]*h, 2*uHit.c[1]*h, r }};
            (void)r3_dir(&nExp, &nExp);
          }
          break;

        case ot_PYRA:
          { /* Normal should be along {±X,Z} or {±Y,Z}: */
            r3_t aObj = (r3_t){{ cObj.c[0], cObj.c[1], cObj.c[2] + rObj.c[2] }}; /* Apex. */
            r3_t uHit; r3_sub(pHit, &aObj, &uHit);
            double h = HI(obj->bbox[2]) - LO(obj->bbox[2]);
            if (fabs(uHit.c[0]) > fabs(uHit.c[1]))
              { nExp = (r3_t){{ (uHit.c[0] > 0 ? +1 : -1)*rObj.c[0], 0.0, h }}; }
            else
              { nExp = (r3_t){{ 0.0, (uHit.c[1] > 0 ? +1 : -1)*rObj.c[1], h }}; }
            (void)r3_dir(&nExp, &nExp);
          }
          break;

        default:
          assert(FALSE);
      }
    return nExp;
  }