alternaria.cpp

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <config.h>

#include <iostream>
#include <sstream>
#include <cstdlib>
#include <cstdio>
#include <cmath>

#include <inttypes.h>

#include <jwsc/config.h>
#include <jwsc/math/constants.h>
#include <jwsc/matblock/matblock.h>
#include <jwsc/imgblock/imgblock.h>
#include <jwsc/imgblock/imgblock_io.h>
#include <jwsc/imgblock/imgblock_util.h>

#ifdef JWSC_HAVE_FFTW3_H
#include <fftw3.h>
#endif

#include <jwsc++/base/exception.h>
#include <jwsc++/base/option.h>

#include "imaging_model.h"
#include "psf_model.h"
#include "structure.h"
#include "spore.h"
#include "hypha.h"
#include "alternaria_model.h"
#include "sampler.h"
#include "imaging_moves.h"
#include "psf_moves.h"
#include "alternaria_moves.h"
#include "dd_spore.h"


/* Program defaults. */
#define  SEED              2537984
#define  MOVE_INFO_FNAME   "move.info"
#define  BEST_INFO_FNAME   "best.info"
#define  ITERATIONS        10
#define  RES_IMG           1.0f
#define  RES_ROW           1.0f
#define  RES_COL           1.0f
#define  DATA_IMGS_NAMED   false
#define  IMG_IN_FNAME      0
#define  IMG_OUT_FNAME     "imaging.model"
#define  PSF_IN_FNAME      0
#define  PSF_OUT_FNAME     "psf.model"
#define  ALT_IN_FNAME      0
#define  ALT_OUT_FNAME     "alternaria.model"
#define  ALT_PRO_FMT       0
#define  ALT_SCENE_FMT     0
#define  PSF_H_FMT         0
#define  BASE_LEVEL        0
#define  DD_SPORES_FNAME   0


using std::cout;
using std::cerr;
using std::ostringstream;
using std::sscanf;
using namespace jwsc;
using jwscxx::base::Options;
using jwscxx::base::Exception;
using jwscxx::base::Arg_error;
using jwscxx::base::IO_error;


Options* opts = 0;

uint32_t seed = SEED;

const char* move_info_fname = MOVE_INFO_FNAME;

const char* best_info_fname = BEST_INFO_FNAME;

uint32_t iterations = ITERATIONS;

float res_img = RES_IMG;

float res_row = RES_ROW;

float res_col = RES_COL;

bool data_imgs_named = DATA_IMGS_NAMED;

const char* img_in_fname = IMG_IN_FNAME;

const char* img_out_fname = IMG_OUT_FNAME;

const char* psf_in_fname = PSF_IN_FNAME;

const char* psf_out_fname = PSF_OUT_FNAME;

const char* alt_in_fname = ALT_IN_FNAME;

const char* alt_out_fname = ALT_OUT_FNAME;

const char* alt_pro_fmt = ALT_PRO_FMT;

const char* alt_scene_fmt = ALT_SCENE_FMT;

const char* psf_h_fmt = PSF_H_FMT;

uint32_t base_level = BASE_LEVEL;

const char* dd_spores_fname = DD_SPORES_FNAME;

Imaging_window* img_window = 0;

Imaging_scale* img_scale = 0;

PSF_padding* padding = 0;

Imaging_model_density* imaging_d = 0;

PSF_model_density* psf_d = 0;

Apical_hypha_density* apical_hypha_d = 0;

Lateral_hypha_density* lateral_hypha_1_d = 0;

Lateral_hypha_density* lateral_hypha_n_d = 0;

Spore_density* spore_d = 0;

Alternaria_model_density* alt_d = 0;

bool img_window_set = false;


void read_data
(
    Matblock_f** data_out, 
    int          argc, 
    const char** argv, 
    int          argi
)
throw (Arg_error)
{
    uint32_t    num_imgs;
    Imgblock_f* imgs = NULL;
    Error*      err;

    /* Read the input Alternaria images. */
    if (data_imgs_named)
    {
        num_imgs = argc - argi;

        if (num_imgs < 1)
        {
            throw Arg_error("No input images");
        }

        if ((err = read_imgblock_named_f(&imgs, &(argv[ argi ]), num_imgs)) 
                != 0)
        {
            Arg_error ex = Arg_error(err->msg);
            delete err;
            throw ex;
        }
    }
    else
    {
        if ((argc - argi) < 2)
        {
            throw Arg_error("No input image format and count");
        }
        
        const char* input_fname = argv[ argi ];

        if (sscanf(argv[ argi+1 ], "%u", &num_imgs) != 1 || num_imgs <= 0)
        {
            throw Arg_error("No input image count");
        }       

        if ((err = read_imgblock_numbered_f(&imgs, input_fname, num_imgs)) != 0)
        {
            Arg_error ex = Arg_error(err->msg);
            delete err;
            throw ex;
        }

    }

    assert(down_sample_imgblock_f(&imgs, imgs, res_img, res_row, res_col) == 0);
    create_matblock_from_imgblock_f(data_out, imgs, 0.3f, 0.59f, 0.11f);
    free_imgblock_f(imgs);
}

void init_density_centroids() throw (Arg_error)
{
    float x_min = img_scale->get_x() * img_window->get_col();
    float x_max = x_min + img_scale->get_x()*(img_window->get_num_cols()-1);
    float y_min = img_scale->get_y() * img_window->get_row();
    float y_max = y_min + img_scale->get_y()*(img_window->get_num_rows()-1);
    float z_min = img_scale->get_z() * img_window->get_img();
    float z_max = z_min + img_scale->get_z()*(img_window->get_num_imgs()-1);

    spore_d->set_centroid_x(x_min, x_max);
    spore_d->set_centroid_y(y_min, y_max);
    spore_d->set_centroid_z(z_min, z_max);

    apical_hypha_d->set_centroid_x(x_min, x_max);
    apical_hypha_d->set_centroid_y(y_min, y_max);
    apical_hypha_d->set_centroid_z(z_min, z_max);

    lateral_hypha_1_d->set_centroid_x(x_min, x_max);
    lateral_hypha_1_d->set_centroid_y(y_min, y_max);
    lateral_hypha_1_d->set_centroid_z(z_min, z_max);

    lateral_hypha_n_d->set_centroid_x(x_min, x_max);
    lateral_hypha_n_d->set_centroid_y(y_min, y_max);
    lateral_hypha_n_d->set_centroid_z(z_min, z_max);
}

void init_sampler_moves(Sampler* sampler) throw (Arg_error)
{
    sampler->add_move(new Imaging_bg_move("Bg"));
    sampler->add_move(new PSF_move("PSF"));
    sampler->add_move(new Apical_hypha_birth_death_move("ApiBirth", 
                "ApiDeath"));
    sampler->add_move(new Spore_birth_death_move("SprBirth", "SprDeath"));
    sampler->add_move(new Spore_split_merge_spore_move("SprSplit", "SprMerge"));
    sampler->add_move(new Spore_transition_apical_move("ApiToSpr", "SprToApi"));
    sampler->add_move(new Spore_split_merge_2_apical_move("Spr2Api", "Api2Spr"));
    sampler->add_move(new Lateral_hypha_birth_death_move("LatBirth", 
                "LatDeath"));
    sampler->add_move(new Apical_hypha_split_merge_apical_move("ApiSplit", 
                "ApiMerge"));
    sampler->add_move(new Apical_hypha_split_merge_lateral_move("LatSplit", 
                "LatMerge"));
    sampler->add_move(new Structure_resize_move("Resize"));
    sampler->add_move(new Structure_rotate_move("Rotate"));
    sampler->add_move(new Structure_opacity_move("Opacity"));
    sampler->add_move(new Apical_structure_position_move("ApiPosn"));
    sampler->add_move(new Lateral_structure_lat_dist_move("LatDist"));

    sampler->normalize_move_probs();
}

void print_usage(void)
{
    cerr << "usage: alternaria OPTIONS [img-%%d num_imgs | img-1 ... img-n]\n";
    cerr << "where OPTIONS is one or more of the following:\n";
    opts->print(cerr, 33);
}

void process_help_opt(void)
{   
    print_usage();
    exit(EXIT_SUCCESS);
}

void process_version_opt(void)
{
    cout << ALTERNARIA_PACKAGE_STRING << "\n";
    exit(EXIT_SUCCESS);
}

void process_options_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'options' requires an argument");
    }
    try
    {
        opts->process(arg);
    }
    catch (IO_error e)
    {
        throw Arg_error(e.get_msg());
    }
}

void process_move_info_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'move-info' requires an argument");
    }
    move_info_fname = arg;
}

void process_best_info_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'best-info' requires an argument");
    }
    best_info_fname = arg;
}

#ifdef JWSC_HAVE_FFTW3F_THREADS
void process_threads_opt(const char* arg) throw (Arg_error)
{
    uint32_t threads;

    if (!arg)
    {
        throw Arg_error("Option 'threads' requires an argument");
    }
    if (sscanf(arg, "%u", &threads) != 1 || threads < 1)
    {
        throw Arg_error("Option 'threads' must be an integer > 0");
    }
    if (fftwf_init_threads() == 0)
    {
        throw Arg_error("Could not initialize FFTW threads");
    }
    fftwf_plan_with_nthreads(threads);
}
#endif

void process_iterations_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'iterations' requires an argument");
    }
    if (sscanf(arg, "%u", &iterations) != 1)
    {
        throw Arg_error("Option 'iterations' must be a non-negative integer");
    }
}

void process_resolution_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'resolution' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f", &res_col, &res_row, &res_img) != 3)
    {
        throw Arg_error("Option 'resolution' has format col,row,img");
    }
    if (res_col <= 0 || res_col > 1 || res_row <= 0 || res_row > 1
            || res_img <= 0 || res_img > 1)
    {
        throw Arg_error("Option 'resolution' must be in (0,1]");
    }
}

void process_seed_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'seed' requires an argument");
    }
    if (sscanf(arg, "%u", &seed) != 1)
    {
        throw Arg_error("Option 'seed' must be a non-negative integer");
    }
}

void process_img_in_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'img-in' requires an argument");
    }
    img_in_fname = arg;
}

void process_img_out_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'img-out' requires an argument");
    }
    img_out_fname = arg;
}

void process_psf_in_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'psf-in' requires an argument");
    }
    psf_in_fname = arg;
}

void process_psf_out_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'psf-out' requires an argument");
    }
    psf_out_fname = arg;
}

void process_alt_in_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'alt-in' requires an argument");
    }
    alt_in_fname = arg;
}

void process_alt_out_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'alt-out' requires an argument");
    }
    alt_out_fname = arg;
}

void process_alt_pro_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'alt-pro' requires an argument");
    }
    alt_pro_fmt = arg;
}

void process_dd_spores_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'dd-spores' requires an argument");
    }
    dd_spores_fname = arg;
}

void process_scene_fmt_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'scene-fmt' requires an argument");
    }
    alt_scene_fmt = arg;
}

void process_scene_scale_opt(const char* arg) throw (Arg_error)
{
    float x, y, z;

    if (!arg)
    {
        throw Arg_error("Option 'scene-scale' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f", &x, &y, &z) != 3)
    {
        throw Arg_error("Option 'scene-scale' has format x,y,z");
    }
    try
    {
        img_scale->set(x, y, z);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'scene-scale': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_scene_window_opt(const char* arg) throw (Arg_error)
{
    uint32_t col, row, img, ncols, nrows, nimgs;

    if (!arg)
    {
        throw Arg_error("Option 'scene-window' requires an argument");
    }
    if (sscanf(arg, "%u,%u,%u,%u,%u,%u", &col, &row, &img, &ncols, &nrows, &nimgs) != 6)
    {
        throw Arg_error("Option 'scene-window' has format col,row,img,ncols,nrows,nimgs");
    }
    try
    {
        img_window->set(col, row, img, ncols, nrows, nimgs);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'scene-window': " << e.get_msg();
        throw Arg_error(ost.str());
    }

    img_window_set = true;
}

void process_psf_fmt_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'psf-fmt' requires an argument");
    }
    psf_h_fmt = arg;
}

void process_apical_hypha_length_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'apical-hypha-length' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'apical-hypha-length' has format mu,sigma,min,max");
    }
    try
    {
        apical_hypha_d->set_length(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'apical-hypha-length': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_lateral_hypha_length_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'lateral-hypha-length' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'lateral-hypha-length' has format mu,sigma,min,max");
    }
    try
    {
        lateral_hypha_1_d->set_length(mu, sigma, min, max);
        lateral_hypha_n_d->set_length(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'lateral-hypha-length': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_apical_hypha_width_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'apical-hypha-width' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'apical-hypha-width' has format mu,sigma,min,max");
    }
    try
    {
        apical_hypha_d->set_width(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'apical-hypha-width': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_apical_hypha_dwidth_opt(const char* arg) throw (Arg_error)
{
    float sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'apical-hypha-dwidth' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f", &sigma, &min, &max) != 3)
    {
        throw Arg_error("Option 'apical-hypha-dwidth' has format sigma,min,max");
    }
    try
    {
        apical_hypha_d->set_dwidth(sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'apical-hypha-dwidth': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_lateral_hypha_width_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'lateral-hypha-width' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'lateral-hypha-width' has format mu,sigma,min,max");
    }
    try
    {
        lateral_hypha_1_d->set_width(mu, sigma, min, max);
        lateral_hypha_n_d->set_width(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'lateral-hypha-width': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_apical_hypha_theta_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'apical-hypha-theta' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'apical-hypha-theta' has format mu,sigma,min,max");
    }
    try
    {
        apical_hypha_d->set_theta(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'apical-hypha-theta': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_lateral_hypha_1_theta_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'lateral-hypha-1-theta' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'lateral-hypha-1-theta' has format mu,sigma,min,max");
    }
    try
    {
        lateral_hypha_1_d->set_theta(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'lateral-hypha-1-theta': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_lateral_hypha_n_theta_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'lateral-hypha-n-theta' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'lateral-hypha-n-theta' has format mu,sigma,min,max");
    }
    try
    {
        lateral_hypha_n_d->set_theta(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'lateral-hypha-n-theta': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_hypha_psi_opt(const char* arg) throw (Arg_error)
{
    float min, max;

    if (!arg)
    {
        throw Arg_error("Option 'hypha-psi' requires an argument");
    }
    if (sscanf(arg, "%f,%f", &min, &max) != 2)
    {
        throw Arg_error("Option 'hypha-psi' has format min,max");
    }
    try
    {
        apical_hypha_d->set_psi(min, max);
        lateral_hypha_1_d->set_psi(min, max);
        lateral_hypha_n_d->set_psi(min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'hypha-psi': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_hypha_opacity_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'hypha-opacity' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'hypha-opacity' has format mu,sigma,min,max");
    }
    try
    {
        apical_hypha_d->set_opacity(mu, sigma, min, max);
        lateral_hypha_1_d->set_opacity(mu, sigma, min, max);
        lateral_hypha_n_d->set_opacity(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'hypha-opacity': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_lateral_hypha_dist_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'lateral-hypha-dist' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'lateral-hypha-dist' has format mu,sigma,min,max");
    }
    try
    {
        lateral_hypha_1_d->set_lat_dist(mu, sigma, min, max);
        lateral_hypha_n_d->set_lat_dist(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'lateral-hypha-dist': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_spore_length_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'spore-length' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'spore-length' has format mu,sigma,min,max");
    }
    try
    {
        spore_d->set_length(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'spore-length': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_spore_width_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'spore-width' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'spore-width' has format mu,sigma,min,max");
    }
    try
    {
        spore_d->set_width(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'spore-width': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_spore_theta_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'spore-theta' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'spore-theta' has format mu,sigma,min,max");
    }
    try
    {
        spore_d->set_theta(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'spore-theta': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_spore_psi_opt(const char* arg) throw (Arg_error)
{
    float min, max;

    if (!arg)
    {
        throw Arg_error("Option 'spore-psi' requires an argument");
    }
    if (sscanf(arg, "%f,%f", &min, &max) != 2)
    {
        throw Arg_error("Option 'spore-psi' has format min,max");
    }
    try
    {
        spore_d->set_psi(min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'spore-psi': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_spore_opacity_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'spore-opacity' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'spore-opacity' has format mu,sigma,min,max");
    }
    try
    {
        spore_d->set_opacity(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'spore-opacity': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_psf_alpha_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'psf-alpha' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'psf-alpha' has format mu,sigma,min,max");
    }
    try
    {
        psf_d->set_alpha(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'psf-alpha': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_psf_beta_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'psf-beta' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'psf-beta' has format mu,sigma,min,max");
    }
    try
    {
        psf_d->set_beta(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'psf-beta': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_psf_gamma_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'psf-gamma' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'psf-gamma' has format mu,sigma,min,max");
    }
    try
    {
        psf_d->set_gamma(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'psf-gamma': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_img_bg_opt(const char* arg) throw (Arg_error)
{
    float mu, sigma, min, max;

    if (!arg)
    {
        throw Arg_error("Option 'img-bg' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &mu, &sigma, &min, &max) != 4)
    {
        throw Arg_error("Option 'img-bg' has format mu,sigma,min,max");
    }
    try
    {
        imaging_d->set_bg(mu, sigma, min, max);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'img-bg': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_levels_p_opt(const char* arg) throw (Arg_error)
{
    float p;

    if (!arg)
    {
        throw Arg_error("Option 'levels-p' requires an argument");
    }
    if (sscanf(arg, "%f", &p) != 1)
    {
        throw Arg_error("Option 'levels-p' must be in a number in (0,1)");
    }
    try
    {
        alt_d->set_num_levels_p(p);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'levels-p': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_structs_p_opt(const char* arg) throw (Arg_error)
{
    float spore_p, apical_hypha_p, lateral_hypha_1_p, lateral_hypha_n_p;

    if (!arg)
    {
        throw Arg_error("Option 'structs-p' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f,%f", &spore_p, &apical_hypha_p, &lateral_hypha_1_p, &lateral_hypha_n_p) != 4)
    {
        throw Arg_error("Option 'structs-p' has format spore,apical-hypha,lateral-hypha-1, lateral-hypha-n");
    }
    try
    {
        alt_d->set_structure_p(spore_p, apical_hypha_p, lateral_hypha_1_p,
                lateral_hypha_n_p);
    }
    catch (Arg_error e)
    {
        ostringstream ost;
        ost << "Option 'structs-p': " << e.get_msg();
        throw Arg_error(ost.str());
    }
}

void process_base_level_opt(const char* arg) throw (Arg_error)
{
    if (!arg)
    {
        throw Arg_error("Option 'base-level' requires an argument");
    }
    if (sscanf(arg, "%u", &base_level) != 1)
    {
        throw Arg_error("Option 'base-level' must be >= 0");
    }
}

void init_options()
{
    using jwscxx::base::Option_no_arg;
    using jwscxx::base::Option_with_arg;

    const char* l_name;
    const char* desc;

    opts = new Options();

    l_name = "help";
    desc   = "Program usage";
    opts->add(new Option_no_arg('h', l_name, desc, process_help_opt));

    l_name = "version";
    desc   = "Program version.";
    opts->add(new Option_no_arg('v', l_name, desc, process_version_opt));

    l_name = "options";
    desc   = "Process a file of program options.";
    opts->add(new Option_with_arg(0, l_name, desc, process_options_opt));

    l_name = "move-info";
    desc   = "File to print sampler move information to. Stdout by default.";
    opts->add(new Option_with_arg(0, l_name, desc, process_move_info_opt));

    l_name = "best-info";
    desc   = "File to print sampler best move information to. Stdout by default.";
    opts->add(new Option_with_arg(0, l_name, desc, process_best_info_opt));

#ifdef JWSC_HAVE_FFTW3F_THREADS
    l_name = "threads";
    desc   = "Number of threads to use.";
    opts->add(new Option_with_arg(0, l_name, desc, process_threads_opt));
#endif

    l_name = "iterations";
    desc   = "Number of iterations to run the sampler.";
    opts->add(new Option_with_arg(0, l_name, desc, process_iterations_opt));

    l_name = "resolution";
    desc   = "Fraction of the data to down sample to. Has format x,y,z. Must be in (0,1].";
    opts->add(new Option_with_arg(0, l_name, desc, process_resolution_opt));

    l_name = "seed";
    desc   = "Random seed to use in srand(). Must be a non-negative integer.";
    opts->add(new Option_with_arg(0, l_name, desc, process_seed_opt));

    l_name = "img-in";
    desc   = "Imaging input file to read from.";
    opts->add(new Option_with_arg(0, l_name, desc, process_img_in_opt));

    l_name = "img-out";
    desc   = "Imaging output file to write to.";
    opts->add(new Option_with_arg(0, l_name, desc, process_img_out_opt));

    l_name = "psf-in";
    desc   = "PSF input file to read from.";
    opts->add(new Option_with_arg(0, l_name, desc, process_psf_in_opt));

    l_name = "psf-out";
    desc   = "PSF output file to write to.";
    opts->add(new Option_with_arg(0, l_name, desc, process_psf_out_opt));

    l_name = "alt-in";
    desc   = "Alternaria input file to read from.";
    opts->add(new Option_with_arg(0, l_name, desc, process_alt_in_opt));

    l_name = "alt-out";
    desc   = "Alternaria output file to write to.";
    opts->add(new Option_with_arg(0, l_name, desc, process_alt_out_opt));

    l_name = "alt-pro";
    desc   = "Printf-fomatted file name for Alternaria propsals.";
    opts->add(new Option_with_arg(0, l_name, desc, process_alt_pro_opt));

    l_name = "dd-spores";
    desc   = "Data-driven spores file name.";
    opts->add(new Option_with_arg(0, l_name, desc, process_dd_spores_opt));

    l_name = "scene-fmt";
    desc   = "Printf-formatted file name for the optically sectioned scene. Must contain one integer conversion";
    opts->add(new Option_with_arg(0, l_name, desc, process_scene_fmt_opt));

    l_name = "scene-scale";
    desc   = "Number of alternaria world units per image pixel. Has format x,y,z.";
    opts->add(new Option_with_arg(0, l_name, desc, process_scene_scale_opt));

    l_name = "scene-window";
    desc   = "Window over the data to do the analysis on. Has format col,row,img,ncols,nrows,nimgs.";
    opts->add(new Option_with_arg(0, l_name, desc, process_scene_window_opt));

    l_name = "psf-fmt";
    desc   = "Printf-formatted file name for the optically sectioned PSF. Must contain one integer conversion.";
    opts->add(new Option_with_arg(0, l_name, desc, process_psf_fmt_opt));

    l_name = "apical-hypha-length";
    desc   = "Length of an apical hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_apical_hypha_length_opt));

    l_name = "lateral-hypha-length";
    desc   = "Length of a lateral hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_lateral_hypha_length_opt));

    l_name = "apical-hypha-width";
    desc   = "Width of an apical hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_apical_hypha_width_opt));

    l_name = "apical-hypha-dwidth";
    desc   = "Width difference between adjacent apical hypha. Has argument format sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_apical_hypha_dwidth_opt));

    l_name = "lateral-hypha-width";
    desc   = "Width of a lateral hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_lateral_hypha_width_opt));

    l_name = "apical-hypha-theta";
    desc   = "Theta Euler angle of an apical hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_apical_hypha_theta_opt));

    l_name = "lateral-hypha-1-theta";
    desc   = "Theta Euler angle of a lateral hypha at level 1. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_lateral_hypha_1_theta_opt));

    l_name = "lateral-hypha-n-theta";
    desc   = "Theta Euler angle of a lateral hypha at level > 1. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_lateral_hypha_n_theta_opt));

    l_name = "hypha-psi";
    desc   = "Psi Euler angle of a hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_hypha_psi_opt));

    l_name = "hypha-opacity";
    desc   = "Opacity of a hypha. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_hypha_opacity_opt));

    l_name = "lateral-hypha-dist";
    desc   = "Distance of a lateral hypha from its parent. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, 
                process_lateral_hypha_dist_opt));

    l_name = "spore-width";
    desc   = "Width of a spore. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_spore_width_opt));

    l_name = "spore-length";
    desc   = "Length of a spore. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_spore_length_opt));

    l_name = "spore-theta";
    desc   = "Theta Euler angle of a spore. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_spore_theta_opt));

    l_name = "spore-psi";
    desc   = "Psi Euler angle of a spore. Has argument format min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_spore_psi_opt));

    l_name = "spore-opacity";
    desc   = "Opacity of a spore. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_spore_opacity_opt));

    l_name = "psf-alpha";
    desc   = "PSF alpha parameter. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_psf_alpha_opt));

    l_name = "psf-beta";
    desc   = "PSF beta parameter. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_psf_beta_opt));

    l_name = "psf-gamma";
    desc   = "PSF gamma parameter. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_psf_gamma_opt));

    l_name = "img-bg";
    desc   = "Image background parameter. Has argument format mu,sigma,min,max.";
    opts->add(new Option_with_arg(0, l_name, desc, process_img_bg_opt));

    l_name = "levels-p";
    desc   = "Geometric distribution parameter over growth levels. Must be in (0,1)";
    opts->add(new Option_with_arg(0, l_name, desc, process_levels_p_opt));

    l_name = "structs-p";
    desc   = "Distribution over structure type. Has argument format spore,apical-hypha,lateral-hypha. All must sum to 1.";
    opts->add(new Option_with_arg(0, l_name, desc, process_structs_p_opt));

    l_name = "base-level";
    desc   = "Base level for the root Structure. Must be >= 0.";
    opts->add(new Option_with_arg(0, l_name, desc, process_base_level_opt));
}

void* operator new (size_t s) throw (std::bad_alloc)
{
    return malloc(s);
}

void operator delete (void* ptr) throw ()
{
    if (ptr)
        free(ptr);
}




int main(int argc, char** argv)
{
    Sampler*          sampler         = 0;
    Alternaria_model* alternaria      = 0;
    Alternaria_model* best_alternaria = 0;
    PSF_model*        psf             = 0;
    PSF_model*        best_psf        = 0;
    Imaging_model*    imaging         = 0;
    Imaging_model*    best_imaging    = 0;
    Matblock_f*       data            = 0;
    Sampler_move_parameters* params   = 0;
    DD_spore_set*     dd_spores       = 0;

    init_options();

    try
    {
        // Read data-driven.
        if (dd_spores_fname)
        {
            dd_spores = new DD_spore_set(dd_spores_fname);
        }

        sampler           = new Sampler();
        img_window        = new Imaging_window();
        img_scale         = new Imaging_scale();
        padding           = new PSF_padding();
        imaging_d         = new Imaging_model_density();
        psf_d             = new PSF_model_density();
        spore_d           = new Spore_density(dd_spores);
        apical_hypha_d    = new Apical_hypha_density();
        lateral_hypha_1_d = new Lateral_hypha_density();
        lateral_hypha_n_d = new Lateral_hypha_density();
        alt_d             = new Alternaria_model_density();

        int argi = argc - opts->process(argc, argv) + 1; 

        std::srand(seed);
        std::rand();
        std::rand();

        read_data(&data, argc, (const char**)argv, argi);

        // Update the image scale by the resolution.
        float x_scale = img_scale->get_x() / res_col;
        float y_scale = img_scale->get_y() / res_row;
        float z_scale = img_scale->get_z() / res_img;
        img_scale->set(x_scale, y_scale, z_scale);

        // Set the imaging window to the size of the data if not already set by 
        // an option.
        if (!img_window_set)
        {
            img_window->set(0, 0, 0, data->num_cols, data->num_rows, 
                    data->num_mats);
        }
        else
        {
            uint32_t img   = (uint32_t)ceil(res_img*img_window->get_img());
            uint32_t row   = (uint32_t)ceil(res_row*img_window->get_row());
            uint32_t col   = (uint32_t)ceil(res_col*img_window->get_col());
            uint32_t nimgs = (uint32_t)floor(res_img*
                                             img_window->get_num_imgs());
            uint32_t nrows = (uint32_t)floor(res_row*
                                             img_window->get_num_rows());
            uint32_t ncols = (uint32_t)floor(res_col*
                                             img_window->get_num_cols());
            img_window->set(col, row, img, ncols, nrows, nimgs);
        }

        // TEMPORARY
        Matblock_f* windowed_data = 0;
        copy_matblock_block_f(&windowed_data, data, 
                img_window->get_img(),
                img_window->get_row(),
                img_window->get_col(),
                img_window->get_num_imgs(),
                img_window->get_num_rows(),
                img_window->get_num_cols());
        Imgblock_f* imgs = 0;
        create_imgblock_from_matblock_f(&imgs, windowed_data);
        free_matblock_f(windowed_data);
        write_imgblock_numbered_f(imgs, "/tmp/data-%04u.tiff");
        free_imgblock_f(imgs);

        init_sampler_moves(sampler);
        init_density_centroids();

        // Read input models.
        if (img_in_fname)
        {
            imaging = new Imaging_model(img_in_fname, img_window, img_scale, 
                    imaging_d);
        }
        else
        {
            imaging = imaging_d->sample(img_window, img_scale);
        }
        if (psf_in_fname)
        {
            psf = new PSF_model(psf_in_fname, padding, psf_d, imaging);
        }
        else
        {
            psf = psf_d->sample(padding, imaging);
        }
        if (alt_in_fname)
        {
            alternaria = new Alternaria_model(alt_in_fname, alt_d, spore_d, 
                    apical_hypha_d, lateral_hypha_1_d, lateral_hypha_n_d, psf, 
                    imaging);
        }
        else
        {
            alternaria = alt_d->sample(0, 0, 0, 0, JWSC_PI_2, 0, base_level, 
                    spore_d, apical_hypha_d, lateral_hypha_1_d, 
                    lateral_hypha_n_d, psf, imaging);
        }

        params = new Sampler_move_parameters(alternaria, psf, imaging, data);

        // Run the sampler.
        sampler->run(&best_alternaria, &best_psf, &best_imaging, params, 
                iterations, move_info_fname, best_info_fname, alt_pro_fmt);


        // Write output models.
        if (img_out_fname)
        {
            best_imaging->print(img_out_fname);
        }
        if (psf_out_fname)
        {
            best_psf->print(psf_out_fname);
        }
        if (alt_out_fname)
        {
            best_alternaria->print(alt_out_fname);
        }
        if (psf_h_fmt)
        {
            best_psf->write_h(psf_h_fmt);
        }
        if (alt_scene_fmt)
        {
            best_alternaria->write_scene(alt_scene_fmt, best_psf);
        }
    }
    catch (Exception e)
    {
        cerr << "alternaria: ";
        e.print();
    }

    delete best_alternaria;
    delete best_psf;
    delete best_imaging;
    delete params;
    delete sampler;
    delete img_window;
    delete img_scale;
    delete padding;
    delete imaging_d;
    delete psf_d;
    delete spore_d;
    delete apical_hypha_d;
    delete lateral_hypha_1_d;
    delete lateral_hypha_n_d;
    delete alt_d;
    delete dd_spores;
    delete opts;
    free_matblock_f(data);

    return EXIT_SUCCESS;
}