psf_model.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 <cassert>
#include <cmath>
#include <iostream>
#include <fstream>
#include <sstream>

#include <inttypes.h>

#include <jwsc/config.h>
#include <jwsc/base/error.h>
#include <jwsc/base/limits.h>
#include <jwsc/prob/pdf.h>
#include <jwsc/vector/vector.h>
#include <jwsc/vector/vector_math.h>
#include <jwsc/matrix/matrix.h>
#include <jwsc/matrix/matrix_math.h>
#include <jwsc/matblock/matblock.h>
#include <jwsc/matblock/matblock_math.h>
#include <jwsc/matblock/matblock_fft.h>
#include <jwsc/imgblock/imgblock.h>
#include <jwsc/imgblock/imgblock_io.h>

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

#include "density.h"
#include "printable.h"
#include "imaging_model.h"
#include "psf_model.h"


#define  PADDING_CONV       8
#define  PADDING_FFT        16

#define  ALPHA_MU           0.88f
#define  ALPHA_SIGMA        0.1f
#define  ALPHA_MIN          1.0e-3f
#define  ALPHA_MAX          0.99f

#define  BETA_MU            0.5f
#define  BETA_SIGMA         0.1f
#define  BETA_MIN           0
#define  BETA_MAX           10.0f

#define  GAMMA_MU           0.8f
#define  GAMMA_SIGMA        0.1f
#define  GAMMA_MIN          1.0e-3f
#define  GAMMA_MAX          1.0f


using std::isinf;
using std::isnan;
using std::sqrt;
using std::log;
using std::ceil;
using std::ostream;
using std::istream;
using std::ofstream;
using std::ifstream;
using std::ostringstream;
using std::istringstream;
using namespace jwsc;
using jwscxx::base::IO_error;
using jwscxx::base::Arg_error;


/* ----------------------------  PSF_padding  ------------------------------- */

PSF_padding::PSF_padding()
{
    conv = PADDING_CONV;
    fft  = PADDING_FFT;
}

/* -------------------------------------------------------------------------- */










/* --------------------------------  PSF  ----------------------------------- */

PSF::PSF(const PSF_padding* padding)
{
    this->padding = padding;
    this->h       = 0;
}


PSF::~PSF()
{
    free_matblock_f(h);
}


PSF& PSF::operator= (const PSF& psf)
{
    padding = psf.padding;

    copy_matblock_f(&h, psf.h);

    return *this;
}


void PSF::write_h(const char* fname_fmt) const throw (IO_error)
{
    if (h == 0)
        return;

    uint32_t num_mats = h->num_mats - 2*(padding->fft + padding->conv);
    uint32_t num_rows = h->num_rows - 2*(padding->fft + padding->conv);
    uint32_t num_cols = h->num_cols - 2*(padding->fft + padding->conv);

    Matblock_f* h_unshifted = 0;
    unshift_matblock_f(&h_unshifted, h);

    uint32_t start_mat = padding->fft + padding->conv;
    uint32_t start_row = padding->fft + padding->conv;
    uint32_t start_col = padding->fft + padding->conv;

    Matblock_f* h_no_pad = 0;
    copy_matblock_block_f(&h_no_pad, h_unshifted, start_mat, start_row, 
            start_col, num_mats, num_rows, num_cols);
    free_matblock_f(h_unshifted);

    float max = h_no_pad->elts[0][0][0];
    for (uint32_t mat = 0; mat < num_mats; mat++)
    {
        for (uint32_t row = 0; row < num_rows; row++)
        {
            for (uint32_t col = 0; col < num_cols; col++)
            {
                if (max < h_no_pad->elts[ mat ][ row ][ col ])
                {
                    max = h_no_pad->elts[ mat ][ row ][ col ];
                }
            }
        }
    }

    if (max > 1.0e-8)
    {
        multiply_matblock_by_scalar_f(&h_no_pad, h_no_pad, 1.0f / max);
    }

    Imgblock_f* imgs = 0;
    create_imgblock_from_matblock_f(&imgs, h_no_pad);
    free_matblock_f(h_no_pad);

    Error* e;
    if ((e = write_imgblock_numbered_f(imgs, fname_fmt)) != 0)
    {
        IO_error ex(e->msg);
        free_error(e);
        throw ex;
    }
    free_imgblock_f(imgs);
}

/* -------------------------------------------------------------------------- */










/* -----------------------  PSF_model_density  ------------------------------ */

PSF_model_density::PSF_model_density() throw (Arg_error)
{
    float mu, sigma;
    float min, max;

    mu    = ALPHA_MU;
    sigma = ALPHA_SIGMA;
    min   = ALPHA_MIN;
    max   = ALPHA_MAX;
    set_alpha(mu, sigma, min, max);

    mu    = BETA_MU;
    sigma = BETA_SIGMA;
    min   = BETA_MIN;
    max   = BETA_MAX;
    set_beta(mu, sigma, min, max);

    mu    = GAMMA_MU;
    sigma = GAMMA_SIGMA;
    min   = GAMMA_MIN;
    max   = GAMMA_MAX;
    set_gamma(mu, sigma, min, max);
}


PSF_model_density* PSF_model_density::clone() const
{
    float mu, sigma;
    float min, max;

    PSF_model_density* d = new PSF_model_density();

    mu    = alpha.get_mu();
    sigma = alpha.get_sigma();
    min   = alpha.get_min();
    max   = alpha.get_max();
    d->set_alpha(mu, sigma, min, max);

    mu    = beta.get_mu();
    sigma = beta.get_sigma();
    min   = beta.get_min();
    max   = beta.get_max();
    d->set_beta(mu, sigma, min, max);

    mu    = gamma.get_mu();
    sigma = gamma.get_sigma();
    min   = gamma.get_min();
    max   = gamma.get_max();
    d->set_gamma(mu, sigma, min, max);

    return d;
}


void PSF_model_density::set_alpha(float mu, float sigma, float min, float max) 
    throw (Arg_error)
{
    if (min <= 0)
    {
        throw Arg_error("Minimum must be > 0");
    }
    if (max >= 1)
    {
        throw Arg_error("Maximum must be < 1");
    }
    alpha.set(mu, sigma, min, max);
}


void PSF_model_density::set_beta(float mu, float sigma, float min, float max) 
    throw (Arg_error)
{
    if (min < 0)
    {
        throw Arg_error("Minimum must be >= 0");
    }
    beta.set(mu, sigma, min, max);
}


void PSF_model_density::set_gamma(float mu, float sigma, float min, float max) 
    throw (Arg_error)
{
    if (min <= 0)
    {
        throw Arg_error("Minimum must be > 0");
    }
    if (max > 1)
    {
        throw Arg_error("Maximum must be <= 1");
    }
    gamma.set(mu, sigma, min, max);
}


float PSF_model_density::sample_alpha() const
{
    return (float)sample_gaussian_pdf_d(alpha.get_mu(), alpha.get_sigma(), 
            alpha.get_min(), alpha.get_max());
}


float PSF_model_density::sample_beta() const
{
    return (float)sample_gaussian_pdf_d(beta.get_mu(), beta.get_sigma(), 
            beta.get_min(), beta.get_max());
}


float PSF_model_density::sample_gamma() const
{
    return (float)sample_gaussian_pdf_d(gamma.get_mu(), gamma.get_sigma(), 
            gamma.get_min(), gamma.get_max());
}


PSF_model* PSF_model_density::sample
(
    const PSF_padding*   padding,
    const Imaging_model* imaging
)
const
{
    return new PSF_model(sample_alpha(), sample_beta(), sample_gamma(), 
            padding, this, imaging);
}

/* -------------------------------------------------------------------------- */










/* --------------------------  PSF_model  ----------------------------------- */

PSF_model::PSF_model
(
    float                    alpha,
    float                    beta,
    float                    gamma,
    const PSF_padding*       padding,
    const PSF_model_density* density,
    const Imaging_model*     imaging
) 
throw (Arg_error) : PSF(padding) 
{
    this->alpha   = alpha;
    this->beta    = beta;
    this->gamma   = gamma;
    this->density = density;

    check_alpha();
    check_beta();
    check_gamma();

    update_log_prob();

    update_h(imaging);
}


PSF_model::PSF_model(const PSF_model& model) : PSF(model.padding)
{
    alpha    = model.alpha;
    beta     = model.beta;
    gamma    = model.gamma;
    density  = model.density;
    log_prob = model.log_prob;

    if (model.h)
    {
        copy_matblock_f(&h, model.h);
    }
}


PSF_model::PSF_model
(
    const char*              fname, 
    const PSF_padding*       padding,
    const PSF_model_density* density,
    const Imaging_model*     imaging
)
throw (IO_error, Arg_error) : PSF(padding)
{
    ostringstream ost;
    istringstream ist;
    ifstream      in;

    this->density = density;

    in.open(fname);
    if (in.fail())
    {
        ost << fname << ": Could not open file";
        throw IO_error(ost.str());
    }

    const char* member_value;

    if (!(member_value = get_member_value("Alpha", in)))
    {
        ost << fname << ": Missing alpha";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> alpha;
    if (ist.fail())
    {
        ost << fname << ": Invalid alpha";
        throw Arg_error(ost.str());
    }
    ist.clear(std::ios_base::goodbit);

    if (!(member_value = get_member_value("Beta", in)))
    {
        ost << fname << ": Missing beta";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> beta;
    if (ist.fail())
    {
        ost << fname << ": Invalid beta";
        throw Arg_error(ost.str());
    }
    ist.clear(std::ios_base::goodbit);

    if (!(member_value = get_member_value("Gamma", in)))
    {
        ost << fname << ": Missing gamma";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> gamma;
    if (ist.fail())
    {
        ost << fname << ": Invalid gamma";
        throw Arg_error(ost.str());
    }
    ist.clear(std::ios_base::goodbit);

    if (!(member_value = get_member_value("Log Prob", in)))
    {
        ost << fname << ": Missing log probability";
        throw Arg_error(ost.str());
    }
    ist.str(member_value);
    ist >> log_prob;
    if (ist.fail())
    {
        ost << fname << ": Invalid log probability";
        throw Arg_error(ost.str());
    }
    ist.clear(std::ios_base::goodbit);

    in.close();
    if (in.fail())
    {
        ost << fname << ": Could not close file";
        throw IO_error(ost.str());
    }

    check_alpha();
    check_beta();
    check_gamma();

    update_log_prob();

    update_h(imaging);
}


PSF_model* PSF_model::clone()
{
    return new PSF_model(*this);
}


PSF_model& PSF_model::operator= (const PSF_model& psf)
{
    PSF::operator=(psf);

    alpha    = psf.alpha;
    beta     = psf.beta;
    gamma    = psf.gamma;
    log_prob = psf.log_prob;
    density  = psf.density;

    return *this;
}


void PSF_model::update_h(const Imaging_model* imaging)
{
    float x_scale = imaging->get_scale()->get_x();
    float y_scale = imaging->get_scale()->get_y();
    float z_scale = imaging->get_scale()->get_z();

    const Imaging_window* window = imaging->get_window();

    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    num_mats = window->get_num_imgs() + 2*(padding->conv + padding->fft);
    num_rows = window->get_num_rows() + 2*(padding->conv + padding->fft);
    num_cols = window->get_num_cols() + 2*(padding->conv + padding->fft);

    create_zero_matblock_f(&h, num_mats, num_rows, num_cols);

    /* Calculate the weights for (1/2) the Gaussians. */
    Vector_f* gauss_weight = 0;
    create_vector_f(&gauss_weight, static_cast<uint32_t>(ceil(num_mats / 2.0)));
    double gauss_weight_sum = 0.0;
    for (mat = 0; mat < (num_mats / 2); mat++)
    {
        gauss_weight->elts[ mat ] = pow(alpha, ceil(num_mats / 2.0) - mat);
        gauss_weight_sum += 2*gauss_weight->elts[ mat ];
    }
    if (num_mats % 2 == 1)
    {
        gauss_weight->elts[ mat ] = alpha;
        gauss_weight_sum += gauss_weight->elts[ mat ];
    }

    multiply_vector_by_scalar_f(&gauss_weight, gauss_weight, 
            1.0 / gauss_weight_sum);

    Matrix_f* m = 0;
    create_zero_matrix_f(&m, num_rows, num_cols);

    double sigma;
    double x_gauss, y_gauss;
    double gauss;
    double gauss_sum;
    double x, y, z;

    uint32_t start_mat = (uint32_t)ceil(num_mats / 2.0) - padding->fft;
    uint32_t start_row = (uint32_t)ceil(num_rows / 2.0) - padding->fft;
    uint32_t start_col = (uint32_t)ceil(num_cols / 2.0) - padding->fft;

    /* Compute the Gaussians in (1/8) the matrices of h. */
    for (mat = start_mat; mat < ceil(num_mats / 2.0); mat++)
    {
        z         = (ceil(num_mats / 2.0) - mat - 1)*z_scale;
        sigma     = sqrt(beta*fabs(z) + gamma*gamma);
        gauss_sum = 0.0;

        for (row = start_row; row < ceil(num_rows / 2.0); row++)
        {
            y       = (ceil(num_rows / 2.0) - row - 1)*y_scale;
#if defined JWSC_HAVE_SLATEC || defined JWSC_HAVE_ERF
            y_gauss = integrate_gaussian_pdf_d(0, sigma, y, y+y_scale);
#else
            y_gauss = gaussian_pdf_d(0, sigma, y+(y_scale*0.5));
#endif

            for (col = start_col; col < ceil(num_cols / 2.0); col++)
            {
                x       = (ceil(num_cols / 2.0) - col - 1)*x_scale;
#if defined JWSC_HAVE_SLATEC || defined JWSC_HAVE_ERF
                x_gauss = integrate_gaussian_pdf_d(0, sigma, x, x+x_scale);
#else
                x_gauss = gaussian_pdf_d(0, sigma, x+(x_scale*0.5));
#endif
                gauss   = x_gauss*y_gauss;

                uint32_t r = row; uint32_t c = col;
                m->elts[ r ][ c ] = gauss;

                r = row; c = num_cols-col-1;
                m->elts[ r ][ c ] = gauss;

                r = num_rows-row-1; c = col; 
                m->elts[ r ][ c ] = gauss;

                r = num_rows-row-1; c = num_cols-col-1; 
                m->elts[ r ][ c ] = gauss;

                if ((row == num_rows-row-1) && (col == num_cols-col-1))
                {
                    gauss_sum += gauss;
                }
                else
                {
                    if (row != num_rows-row-1) 
                        gauss_sum += 2*gauss;
                    if (col != num_cols-col-1) 
                        gauss_sum += 2*gauss;
                }
            }
        }

        multiply_matrix_by_scalar_f(&m, m, gauss_weight->elts[mat] / gauss_sum);

        copy_matrix_into_matblock_f(&h, h, m, MATBLOCK_MAT_MATRIX, mat);

        if (mat != num_mats-mat-1)
        {
            copy_matrix_into_matblock_f(&h, h, m, MATBLOCK_MAT_MATRIX, 
                    num_mats-mat-1);
        }
    }

    float sum = 0;

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                sum += h->elts[ mat ][ row ][ col ];
            }
        }
    }

    assert(sum >= 1.0e-12);
    multiply_matblock_by_scalar_f(&h, h, 1.0f/sum);

    shift_matblock_f(&h, h);

    free_vector_f(gauss_weight);
    free_matrix_f(m);
}


void PSF_model::set
(
    float                alpha, 
    float                beta, 
    float                gamma, 
    const Imaging_model* imaging
)
throw (Arg_error)
{
    float old_alpha = this->alpha;
    float old_beta  = this->beta;
    float old_gamma = this->gamma;

    try
    {
        this->alpha = alpha;
        this->beta  = beta;
        this->gamma = gamma;

        check_alpha();
        check_beta();
        check_gamma();

        update_log_prob();
        update_h(imaging);
    }
    catch (Arg_error e)
    {
        set(old_alpha, old_beta, old_gamma, imaging);
        throw e;
    }
}


void PSF_model::set_alpha(float alpha, const Imaging_model* imaging) 
    throw (Arg_error)
{
    float old_alpha = this->alpha;

    try
    {
        this->alpha = alpha;
        check_alpha();

        update_log_prob();
        update_h(imaging);
    }
    catch (Arg_error e)
    {
        set_alpha(old_alpha, imaging);
        throw e;
    }
}


void PSF_model::set_beta(float beta, const Imaging_model* imaging) 
    throw (Arg_error)
{
    float old_beta = this->beta;

    try
    {
        this->beta = beta;
        check_beta();

        update_log_prob();
        update_h(imaging);
    }
    catch (Arg_error e)
    {
        set_beta(old_beta, imaging);
        throw e;
    }
}


void PSF_model::set_gamma(float gamma, const Imaging_model* imaging) 
    throw (Arg_error)
{
    float old_gamma = this->gamma;

    try
    {
        this->gamma = gamma;
        check_gamma();

        update_log_prob();
        update_h(imaging);
    }
    catch (Arg_error e)
    {
        set_gamma(old_gamma, imaging);
        throw e;
    }
}


void PSF_model::print(const char* fname) const throw (IO_error)
{
    ofstream out;

    out.open(fname);
    if (out.fail())
    {
        ostringstream ost;
        ost << fname << ": Could not open file";
        throw IO_error(ost.str());
    }

    print(out);

    out.close();
    if (out.fail())
    {
        ostringstream ost;
        ost << fname << ": Could not close file";
        throw IO_error(ost.str());
    }
}


void PSF_model::print(std::ostream& out) const
{
    out << "    Alpha: " << alpha << '\n'
        << "     Beta: " << beta << '\n'
        << "    Gamma: " << gamma << '\n'
        << " Log Prob: " << log_prob << '\n';
}


void PSF_model::update_log_prob()
{
    float  mu, sigma;
    float  min, max;
    double delta;
    double p_1, p_2;

    log_prob = 0;

    mu        = density->get_alpha().get_mu();
    sigma     = density->get_alpha().get_sigma();
    min       = density->get_alpha().get_min();
    max       = density->get_alpha().get_max();
    delta     = (max - min) / 1000.0;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, alpha, alpha+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    mu        = density->get_beta().get_mu();
    sigma     = density->get_beta().get_sigma();
    min       = density->get_beta().get_min();
    max       = density->get_beta().get_max();
    delta     = (max - min) / 1000.0;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, beta, beta+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

    mu        = density->get_gamma().get_mu();
    sigma     = density->get_gamma().get_sigma();
    min       = density->get_gamma().get_min();
    max       = density->get_gamma().get_max();
    delta     = (max - min) / 1000.0;
    p_1       = integrate_gaussian_pdf_d(mu, sigma, gamma, gamma+delta); 
    p_2       = integrate_gaussian_pdf_d(mu, sigma, min, max);
    log_prob += log((p_1 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_1);
    log_prob -= log((p_2 < JWSC_MIN_LOG_ARG) ? JWSC_MIN_LOG_ARG : p_2);

#if defined ALTERNARIA_HAVE_ISINF && defined ALTERNARIA_HAVE_ISNAN
    assert(!isinf(log_prob) && !isnan(log_prob));
#endif
}


void PSF_model::check_alpha() const throw (Arg_error)
{
    if (density->get_alpha().get_min() > alpha || 
        density->get_alpha().get_max() < alpha)
    {
        ostringstream ost;
        ost << "Alpha " << alpha << " outside density range ["
            << density->get_alpha().get_min() << ','
            << density->get_alpha().get_max() << ']';
        throw Arg_error(ost.str());
    }
}


void PSF_model::check_beta() const throw (Arg_error)
{
    if (density->get_beta().get_min() > beta || 
        density->get_beta().get_max() < beta)
    {
        ostringstream ost;
        ost << "Beta " << beta << " outside density range ["
            << density->get_beta().get_min() << ','
            << density->get_beta().get_max() << ']';
        throw Arg_error(ost.str());
    }
}


void PSF_model::check_gamma() const throw (Arg_error)
{
    if (density->get_gamma().get_min() > gamma || 
        density->get_gamma().get_max() < gamma)
    {
        ostringstream ost;
        ost << "Gamma " << gamma << " outside density range ["
            << density->get_gamma().get_min() << ','
            << density->get_gamma().get_max() << ']';
        throw Arg_error(ost.str());
    }
}

/* -------------------------------------------------------------------------- */












/* -----------------------  PSF_delta_model  -------------------------------- */

PSF_delta_model::PSF_delta_model
(
    const PSF_padding*   padding,
    const Imaging_model* imaging
)
: PSF(padding)
{
    update_h(imaging);
}


void PSF_delta_model::update_h(const Imaging_model* imaging)
{
    const Imaging_window* window = imaging->get_window();

    uint32_t num_mats, num_rows, num_cols;
    num_mats = window->get_num_imgs() + 2*(padding->conv + padding->fft);
    num_rows = window->get_num_rows() + 2*(padding->conv + padding->fft);
    num_cols = window->get_num_cols() + 2*(padding->conv + padding->fft);

    create_zero_matblock_f(&h, num_mats, num_rows, num_cols);

    double h_val = 1.0;

    if (num_mats % 2 == 0) h_val *= 0.5;
    if (num_rows % 2 == 0) h_val *= 0.5;
    if (num_cols % 2 == 0) h_val *= 0.5;

    for (uint32_t mat = static_cast<uint32_t>(ceil(num_mats / 2.0) - 1.0); 
            mat <= (num_mats / 2); mat++)
    {
        for (uint32_t row = static_cast<uint32_t>(ceil(num_rows / 2.0) - 1.0); 
                row <= (num_rows / 2); row++)
        {
            for (uint32_t col = static_cast<uint32_t>(ceil(num_cols / 2.0) - 1.0); 
                    col <= (num_cols / 2); col++)
            {
                h->elts[ mat ][ row ][ col ] = h_val;
            }
        }
    }

    /* debug. */
    double delta_sum = 0;
    for (uint32_t mat = 0; mat < num_mats; mat++)
    {
        for (uint32_t row = 0; row < num_rows; row++)
        {
            for (uint32_t col = 0; col < num_cols; col++)
            {
                delta_sum += h->elts[ mat ][ row ][ col ];
            }
        }
    }
    assert(fabs(1.0 - delta_sum) < 1.0e-16);

    shift_matblock_f(&h, h);
}

/* -------------------------------------------------------------------------- */












/* ------------------------  PSF_gauss_model  ------------------------------- */

PSF_gauss_model::PSF_gauss_model
(
    float                sigma,
    const PSF_padding*   padding,
    const Imaging_model* imaging
)
throw (Arg_error) : PSF(padding)
{
    this->sigma = sigma;
    if (sigma < 1.0e-8)
    {
        throw Arg_error("Sigma is too small");
    }

    update_h(imaging);
}


PSF_gauss_model& PSF_gauss_model::operator= (const PSF_gauss_model& psf)
{
    PSF::operator=(psf);

    sigma = psf.sigma;

    return *this;
}


void PSF_gauss_model::update_h(const Imaging_model* imaging)
{
    float x_scale = imaging->get_scale()->get_x();
    float y_scale = imaging->get_scale()->get_y();
    float z_scale = imaging->get_scale()->get_z();

    const Imaging_window* window = imaging->get_window();

    uint32_t num_mats, num_rows, num_cols;
    num_mats = window->get_num_imgs() + 2*(padding->conv + padding->fft);
    num_rows = window->get_num_rows() + 2*(padding->conv + padding->fft);
    num_cols = window->get_num_cols() + 2*(padding->conv + padding->fft);

    create_zero_matblock_f(&h, num_mats, num_rows, num_cols);

    Matrix_f* m = 0;
    create_zero_matrix_f(&m, num_rows, num_cols);

    double x, y, z;
    double x_gauss, y_gauss, z_gauss;
    double gauss;

    uint32_t start_mat = (uint32_t)ceil(num_mats / 2.0) - padding->fft;
    uint32_t start_row = (uint32_t)ceil(num_rows / 2.0) - padding->fft;
    uint32_t start_col = (uint32_t)ceil(num_cols / 2.0) - padding->fft;

    /* Compute the Gaussian in (1/8) the matrices of h. */
    for (uint32_t mat = start_mat; mat < ceil(num_mats / 2.0); mat++)
    {
        z         = (ceil(num_mats / 2.0) - mat - 1)*z_scale;
#if defined JWSC_HAVE_SLATEC || defined JWSC_HAVE_ERF
        z_gauss = integrate_gaussian_pdf_d(0, sigma, z, z+z_scale);
#else
        z_gauss = gaussian_pdf_d(0, sigma, z+(z_scale*0.5));
#endif

        for (uint32_t row = start_row; row < ceil(num_rows / 2.0); row++)
        {
            y       = (ceil(num_rows / 2.0) - row - 1)*y_scale;
#if defined JWSC_HAVE_SLATEC || defined JWSC_HAVE_ERF
            y_gauss = integrate_gaussian_pdf_d(0, sigma, y, y+y_scale);
#else
            y_gauss = gaussian_pdf_d(0, sigma, y+(y_scale*0.5));
#endif

            for (uint32_t col = start_col; col < ceil(num_cols / 2.0); col++)
            {
                x       = (ceil(num_cols / 2.0) - col - 1)*x_scale;
#if defined JWSC_HAVE_SLATEC || defined JWSC_HAVE_ERF
                x_gauss = integrate_gaussian_pdf_d(0, sigma, x, x+x_scale);
#else
                x_gauss = gaussian_pdf_d(0, sigma, x+(x_scale*0.5));
#endif
                gauss   = x_gauss*y_gauss*z_gauss;

                uint32_t r = row; uint32_t c = col;
                m->elts[ r ][ c ] = gauss;

                r = row; c = num_cols-col-1;
                m->elts[ r ][ c ] = gauss;

                r = num_rows-row-1; c = col; 
                m->elts[ r ][ c ] = gauss;

                r = num_rows-row-1; c = num_cols-col-1; 
                m->elts[ r ][ c ] = gauss;
            }
        }

        copy_matrix_into_matblock_f(&h, h, m, MATBLOCK_MAT_MATRIX, mat);

        if (mat != num_mats-mat-1)
        {
            copy_matrix_into_matblock_f(&h, h, m, MATBLOCK_MAT_MATRIX, 
                    num_mats-mat-1);
        }
    }

    float sum = 0;

    for (uint32_t mat = 0; mat < num_mats; mat++)
    {
        for (uint32_t row = 0; row < num_rows; row++)
        {
            for (uint32_t col = 0; col < num_cols; col++)
            {
                sum += h->elts[ mat ][ row ][ col ];
            }
        }
    }

    assert(sum >= 1.0e-12);
    multiply_matblock_by_scalar_f(&h, h, 1.0f/sum);

    shift_matblock_f(&h, h);
}

/* -------------------------------------------------------------------------- */