3d.c

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 <jwsc/config.h>

#include <stdlib.h>
#include <inttypes.h>
#include <math.h>
#include <assert.h>

#include "jwsc/base/error.h"
#include "jwsc/prob/pdf.h"
#include "jwsc/matblock/matblock.h"
#include "jwsc/matblock/matblock_math.h"
#include "jwsc/filter/3d.h"


#define  PI  3.14159265358979323846


Error* create_auto_3d_gaussian_filter_f
(
    Matblock_f** h_out, 
    float        mat_sigma,
    float        row_sigma,
    float        col_sigma
)
{
    uint32_t num_mats, num_rows, num_cols;

    if ((num_mats = ceil(sqrt(2) * 6.0 * mat_sigma)) % 2 == 0) num_mats++;
    if ((num_rows = ceil(sqrt(2) * 6.0 * row_sigma)) % 2 == 0) num_rows++;
    if ((num_cols = ceil(sqrt(2) * 6.0 * col_sigma)) % 2 == 0) num_cols++;

    return create_3d_gaussian_filter_f(h_out, mat_sigma, row_sigma, col_sigma, 
            num_mats, num_rows, num_cols);
}

Error* create_auto_3d_gaussian_filter_d
(
    Matblock_d** h_out, 
    double       mat_sigma,
    double       row_sigma,
    double       col_sigma
)
{
    uint32_t num_mats, num_rows, num_cols;

    if ((num_mats = ceil(sqrt(2) * 6.0 * mat_sigma)) % 2 == 0) num_mats++;
    if ((num_rows = ceil(sqrt(2) * 6.0 * row_sigma)) % 2 == 0) num_rows++;
    if ((num_cols = ceil(sqrt(2) * 6.0 * col_sigma)) % 2 == 0) num_cols++;

    return create_3d_gaussian_filter_d(h_out, mat_sigma, row_sigma, col_sigma, 
            num_mats, num_rows, num_cols);
}

Error* create_3d_gaussian_filter_f
(
    Matblock_f** h_out, 
    float        mat_sigma,
    float        row_sigma,
    float        col_sigma,
    uint32_t     num_mats,
    uint32_t     num_rows,
    uint32_t     num_cols
)
{
    uint32_t    mat, row, col;
    uint32_t    m, r, c;
    float       x, y, z;
    float       x_gauss, y_gauss, z_gauss;
    float       gauss;
    Matblock_f* h;

    if (mat_sigma < 0.01 || row_sigma < 0.01 || col_sigma < 0.01)
    {
        return JWSC_EARG("Sigma for Gaussian is too small");
    }

    create_matblock_f(h_out, num_mats, num_rows, num_cols);
    h = *h_out;

#if defined JWSC_HAVE_SLATEC || defined JWSC_HAVE_ERF
    for (mat = 0; mat < ceil(0.5*num_mats); mat++)
    {
        z       = (ceil(0.5*num_mats) - mat - 1);
        z_gauss = integrate_gaussian_pdf_d(0, mat_sigma, z, z+1);
        for (row = 0; row < ceil(0.5*num_rows); row++)
        {
            y       = (ceil(0.5*num_rows) - row - 1);
            y_gauss = integrate_gaussian_pdf_d(0, row_sigma, y, y+1);
            for (col = 0; col < ceil(0.5*num_cols); col++)
            {
                x       = (ceil(0.5*num_cols) - col - 1);
                x_gauss = integrate_gaussian_pdf_d(0, col_sigma, x, x+1);
                gauss   = x_gauss*y_gauss*z_gauss;

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

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

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

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

                m = num_mats-mat-1; r = row; c = col;
                h->elts[ m ][ r ][ c ] = gauss;

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

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

                m = num_mats-mat-1; r = num_rows-row-1; c = num_cols-col-1;
                h->elts[ m ][ r ][ c ] = gauss;
            }
        }
    }
#else
    for (mat = 0; mat < ceil(0.5*num_mats); mat++)
    {
        z       = (ceil(0.5*num_mats) - mat - 1);
        z_gauss = gaussian_pdf_d(0, row_sigma, z+0.5);
        for (row = 0; row < ceil(0.5*num_rows); row++)
        {
            y       = (ceil(0.5*num_rows) - row - 1);
            y_gauss = gaussian_pdf_d(0, row_sigma, y+0.5);
            for (col = 0; col < ceil(0.5*num_cols); col++)
            {
                x       = (ceil(0.5*num_cols) - col - 1);
                x_gauss = gaussian_pdf_d(0, col_sigma, x+0.5);
                gauss   = x_gauss*y_gauss*z_gauss;

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

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

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

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

                m = num_mats-mat-1; r = row; c = col;
                h->elts[ m ][ r ][ c ] = gauss;

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

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

                m = num_mats-mat-1; r = num_rows-row-1; c = num_cols-col-1;
                h->elts[ m ][ r ][ c ] = gauss;
            }
        }
    }
#endif

    return NULL;
}

Error* create_3d_gaussian_filter_d
(
    Matblock_d** h_out, 
    double       mat_sigma,
    double       row_sigma,
    double       col_sigma,
    uint32_t     num_mats,
    uint32_t     num_rows,
    uint32_t     num_cols
)
{
    uint32_t    mat, row, col;
    uint32_t    m, r, c;
    double      x, y, z;
    double      x_gauss, y_gauss, z_gauss;
    double      gauss;
    Matblock_d* h;

    if (mat_sigma < 0.01 || row_sigma < 0.01 || col_sigma < 0.01)
    {
        return JWSC_EARG("Sigma for Gaussian is too small");
    }

    create_matblock_d(h_out, num_mats, num_rows, num_cols);
    h = *h_out;

#if defined JWSC_HAVE_SLATEC || defined JWSC_HAVE_ERF
    for (mat = 0; mat < ceil(0.5*num_mats); mat++)
    {
        z       = (ceil(0.5*num_mats) - mat - 1);
        z_gauss = integrate_gaussian_pdf_d(0, mat_sigma, z, z+1);
        for (row = 0; row < ceil(0.5*num_rows); row++)
        {
            y       = (ceil(0.5*num_rows) - row - 1);
            y_gauss = integrate_gaussian_pdf_d(0, row_sigma, y, y+1);
            for (col = 0; col < ceil(0.5*num_cols); col++)
            {
                x       = (ceil(0.5*num_cols) - col - 1);
                x_gauss = integrate_gaussian_pdf_d(0, col_sigma, x, x+1);
                gauss   = x_gauss*y_gauss*z_gauss;

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

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

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

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

                m = num_mats-mat-1; r = row; c = col;
                h->elts[ m ][ r ][ c ] = gauss;

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

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

                m = num_mats-mat-1; r = num_rows-row-1; c = num_cols-col-1;
                h->elts[ m ][ r ][ c ] = gauss;
            }
        }
    }
#else
    for (mat = 0; mat < ceil(0.5*num_mats); mat++)
    {
        z       = (ceil(0.5*num_mats) - mat - 1);
        z_gauss = gaussian_pdf_d(0, row_sigma, z+0.5);
        for (row = 0; row < ceil(0.5*num_rows); row++)
        {
            y       = (ceil(0.5*num_rows) - row - 1);
            y_gauss = gaussian_pdf_d(0, row_sigma, y+0.5);
            for (col = 0; col < ceil(0.5*num_cols); col++)
            {
                x       = (ceil(0.5*num_cols) - col - 1);
                x_gauss = gaussian_pdf_d(0, col_sigma, x+0.5);
                gauss   = x_gauss*y_gauss*z_gauss;

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

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

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

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

                m = num_mats-mat-1; r = row; c = col;
                h->elts[ m ][ r ][ c ] = gauss;

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

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

                m = num_mats-mat-1; r = num_rows-row-1; c = num_cols-col-1;
                h->elts[ m ][ r ][ c ] = gauss;
            }
        }
    }
#endif

    return NULL;
}

Error* create_auto_3d_gaussian_dx_filter_f
(
    Matblock_f** h_out, 
    float        mat_sigma,
    float        row_sigma,
    float        col_sigma
)
{
    uint32_t num_mats, num_rows, num_cols;

    if ((num_mats = ceil(sqrt(2) * 6.0 * mat_sigma)) % 2 == 0) num_mats++;
    if ((num_rows = ceil(sqrt(2) * 6.0 * row_sigma)) % 2 == 0) num_rows++;
    if ((num_cols = ceil(sqrt(2) * 6.0 * col_sigma)) % 2 == 0) num_cols++;

    return create_3d_gaussian_dx_filter_f(h_out, mat_sigma, row_sigma, 
            col_sigma, num_mats, num_rows, num_cols);
}

Error* create_auto_3d_gaussian_dx_filter_d
(
    Matblock_d** h_out, 
    double       mat_sigma,
    double       row_sigma,
    double       col_sigma
)
{
    uint32_t num_mats, num_rows, num_cols;

    if ((num_mats = ceil(sqrt(2) * 6.0 * mat_sigma)) % 2 == 0) num_mats++;
    if ((num_rows = ceil(sqrt(2) * 6.0 * row_sigma)) % 2 == 0) num_rows++;
    if ((num_cols = ceil(sqrt(2) * 6.0 * col_sigma)) % 2 == 0) num_cols++;

    return create_3d_gaussian_dx_filter_d(h_out, mat_sigma, row_sigma, 
            col_sigma, num_mats, num_rows, num_cols);
}

Error* create_3d_gaussian_dx_filter_f
(
    Matblock_f** h_out, 
    float        mat_sigma,
    float        row_sigma,
    float        col_sigma,
    uint32_t     num_mats,
    uint32_t     num_rows,
    uint32_t     num_cols
)
{
    int64_t mat, row, col;
    int64_t mats_center;
    int64_t rows_center;
    int64_t cols_center;

    float inv_mat_sigma_sqr;
    float inv_row_sigma_sqr;
    float inv_col_sigma_sqr;
    float gauss_val;
    float val;

    Matblock_f* h;

    if (mat_sigma < 0.01 || row_sigma < 0.01 || col_sigma < 0.01)
    {
        return JWSC_EARG("Sigma for Gaussian is too small");
    }

    create_matblock_f(h_out, num_mats, num_rows, num_cols);
    h = *h_out;

    mats_center       = num_mats / 2;
    rows_center       = num_rows / 2;
    cols_center       = num_cols / 2;
    inv_mat_sigma_sqr = 1.0 / (mat_sigma*mat_sigma);
    inv_row_sigma_sqr = 1.0 / (row_sigma*row_sigma);
    inv_col_sigma_sqr = 1.0 / (col_sigma*col_sigma);

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                val = (mat-mats_center)*(mat-mats_center)*inv_mat_sigma_sqr +
                      (row-rows_center)*(row-rows_center)*inv_row_sigma_sqr + 
                      (col-cols_center)*(col-cols_center)*inv_col_sigma_sqr;
                gauss_val = -1*(col - cols_center)*exp(-0.5 * val);
                h->elts[ mat ][ row ][ col ] = gauss_val;
            }
        }
    }

    val = (inv_mat_sigma_sqr*inv_row_sigma_sqr*inv_col_sigma_sqr) / (2 * PI);
    multiply_matblock_by_scalar_f(h_out, h, val);

    return NULL;
}

Error* create_3d_gaussian_dx_filter_d
(
    Matblock_d** h_out, 
    double       mat_sigma,
    double       row_sigma,
    double       col_sigma,
    uint32_t     num_mats,
    uint32_t     num_rows,
    uint32_t     num_cols
)
{
    int64_t mat, row, col;
    int64_t mats_center;
    int64_t rows_center;
    int64_t cols_center;

    double inv_mat_sigma_sqr;
    double inv_row_sigma_sqr;
    double inv_col_sigma_sqr;
    double gauss_val;
    double val;

    Matblock_d* h;

    if (mat_sigma < 0.01 || row_sigma < 0.01 || col_sigma < 0.01)
    {
        return JWSC_EARG("Sigma for Gaussian is too small");
    }

    create_matblock_d(h_out, num_mats, num_rows, num_cols);
    h = *h_out;

    mats_center       = num_mats / 2;
    rows_center       = num_rows / 2;
    cols_center       = num_cols / 2;
    inv_mat_sigma_sqr = 1.0 / (mat_sigma*mat_sigma);
    inv_row_sigma_sqr = 1.0 / (row_sigma*row_sigma);
    inv_col_sigma_sqr = 1.0 / (col_sigma*col_sigma);

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                val = (mat-mats_center)*(mat-mats_center)*inv_mat_sigma_sqr +
                      (row-rows_center)*(row-rows_center)*inv_row_sigma_sqr + 
                      (col-cols_center)*(col-cols_center)*inv_col_sigma_sqr;
                gauss_val = -1*(col - cols_center)*exp(-0.5 * val);
                h->elts[ mat ][ row ][ col ] = gauss_val;
            }
        }
    }

    val = (inv_mat_sigma_sqr*inv_row_sigma_sqr*inv_col_sigma_sqr) / (2 * PI);
    multiply_matblock_by_scalar_d(h_out, h, val);

    return NULL;
}

Error* create_auto_3d_gaussian_dy_filter_f
(
    Matblock_f** h_out, 
    float        mat_sigma,
    float        row_sigma,
    float        col_sigma
)
{
    uint32_t num_mats, num_rows, num_cols;

    if ((num_mats = ceil(sqrt(2) * 6.0 * mat_sigma)) % 2 == 0) num_mats++;
    if ((num_rows = ceil(sqrt(2) * 6.0 * row_sigma)) % 2 == 0) num_rows++;
    if ((num_cols = ceil(sqrt(2) * 6.0 * col_sigma)) % 2 == 0) num_cols++;

    return create_3d_gaussian_dy_filter_f(h_out, mat_sigma, row_sigma, 
            col_sigma, num_mats, num_rows, num_cols);
}

Error* create_auto_3d_gaussian_dy_filter_d
(
    Matblock_d** h_out, 
    double       mat_sigma,
    double       row_sigma,
    double       col_sigma
)
{
    uint32_t num_mats, num_rows, num_cols;

    if ((num_mats = ceil(sqrt(2) * 6.0 * mat_sigma)) % 2 == 0) num_mats++;
    if ((num_rows = ceil(sqrt(2) * 6.0 * row_sigma)) % 2 == 0) num_rows++;
    if ((num_cols = ceil(sqrt(2) * 6.0 * col_sigma)) % 2 == 0) num_cols++;

    return create_3d_gaussian_dy_filter_d(h_out, mat_sigma, row_sigma, 
            col_sigma, num_mats, num_rows, num_cols);
}

Error* create_3d_gaussian_dy_filter_f
(
    Matblock_f** h_out, 
    float        mat_sigma,
    float        row_sigma,
    float        col_sigma,
    uint32_t     num_mats,
    uint32_t     num_rows,
    uint32_t     num_cols
)
{
    int64_t mat, row, col;
    int64_t mats_center;
    int64_t rows_center;
    int64_t cols_center;

    float inv_mat_sigma_sqr;
    float inv_row_sigma_sqr;
    float inv_col_sigma_sqr;
    float gauss_val;
    float val;

    Matblock_f* h;

    if (mat_sigma < 0.01 || row_sigma < 0.01 || col_sigma < 0.01)
    {
        return JWSC_EARG("Sigma for Gaussian is too small");
    }

    create_matblock_f(h_out, num_mats, num_rows, num_cols);
    h = *h_out;

    mats_center       = num_mats / 2;
    rows_center       = num_rows / 2;
    cols_center       = num_cols / 2;
    inv_mat_sigma_sqr = 1.0 / (mat_sigma*mat_sigma);
    inv_row_sigma_sqr = 1.0 / (row_sigma*row_sigma);
    inv_col_sigma_sqr = 1.0 / (col_sigma*col_sigma);

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                val = (mat-mats_center)*(mat-mats_center)*inv_mat_sigma_sqr +
                      (row-rows_center)*(row-rows_center)*inv_row_sigma_sqr + 
                      (col-cols_center)*(col-cols_center)*inv_col_sigma_sqr;
                gauss_val = -1*(row - rows_center)*exp(-0.5 * val);
                h->elts[ mat ][ row ][ col ] = gauss_val;
            }
        }
    }

    val = (inv_mat_sigma_sqr*inv_row_sigma_sqr*inv_col_sigma_sqr) / (2 * PI);
    multiply_matblock_by_scalar_f(h_out, h, val);

    return NULL;
}

Error* create_3d_gaussian_dy_filter_d
(
    Matblock_d** h_out, 
    double       mat_sigma,
    double       row_sigma,
    double       col_sigma,
    uint32_t     num_mats,
    uint32_t     num_rows,
    uint32_t     num_cols
)
{
    int64_t mat, row, col;
    int64_t mats_center;
    int64_t rows_center;
    int64_t cols_center;

    double inv_mat_sigma_sqr;
    double inv_row_sigma_sqr;
    double inv_col_sigma_sqr;
    double gauss_val;
    double val;

    Matblock_d* h;

    if (mat_sigma < 0.01 || row_sigma < 0.01 || col_sigma < 0.01)
    {
        return JWSC_EARG("Sigma for Gaussian is too small");
    }

    create_matblock_d(h_out, num_mats, num_rows, num_cols);
    h = *h_out;

    mats_center       = num_mats / 2;
    rows_center       = num_rows / 2;
    cols_center       = num_cols / 2;
    inv_mat_sigma_sqr = 1.0 / (mat_sigma*mat_sigma);
    inv_row_sigma_sqr = 1.0 / (row_sigma*row_sigma);
    inv_col_sigma_sqr = 1.0 / (col_sigma*col_sigma);

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                val = (mat-mats_center)*(mat-mats_center)*inv_mat_sigma_sqr +
                      (row-rows_center)*(row-rows_center)*inv_row_sigma_sqr + 
                      (col-cols_center)*(col-cols_center)*inv_col_sigma_sqr;
                gauss_val = -1*(row - rows_center)*exp(-0.5 * val);
                h->elts[ mat ][ row ][ col ] = gauss_val;
            }
        }
    }

    val = (inv_mat_sigma_sqr*inv_row_sigma_sqr*inv_col_sigma_sqr) / (2 * PI);
    multiply_matblock_by_scalar_d(h_out, h, val);

    return NULL;
}

Error* create_auto_3d_gaussian_dz_filter_f
(
    Matblock_f** h_out, 
    float        mat_sigma,
    float        row_sigma,
    float        col_sigma
)
{
    uint32_t num_mats, num_rows, num_cols;

    if ((num_mats = ceil(sqrt(2) * 6.0 * mat_sigma)) % 2 == 0) num_mats++;
    if ((num_rows = ceil(sqrt(2) * 6.0 * row_sigma)) % 2 == 0) num_rows++;
    if ((num_cols = ceil(sqrt(2) * 6.0 * col_sigma)) % 2 == 0) num_cols++;

    return create_3d_gaussian_dz_filter_f(h_out, mat_sigma, row_sigma, 
            col_sigma, num_mats, num_rows, num_cols);
}

Error* create_auto_3d_gaussian_dz_filter_d
(
    Matblock_d** h_out, 
    double       mat_sigma,
    double       row_sigma,
    double       col_sigma
)
{
    uint32_t num_mats, num_rows, num_cols;

    if ((num_mats = ceil(sqrt(2) * 6.0 * mat_sigma)) % 2 == 0) num_mats++;
    if ((num_rows = ceil(sqrt(2) * 6.0 * row_sigma)) % 2 == 0) num_rows++;
    if ((num_cols = ceil(sqrt(2) * 6.0 * col_sigma)) % 2 == 0) num_cols++;

    return create_3d_gaussian_dz_filter_d(h_out, mat_sigma, row_sigma, 
            col_sigma, num_mats, num_rows, num_cols);
}

Error* create_3d_gaussian_dz_filter_f
(
    Matblock_f** h_out, 
    double       mat_sigma,
    double       row_sigma,
    double       col_sigma,
    uint32_t     num_mats,
    uint32_t     num_rows,
    uint32_t     num_cols
)
{
    int64_t mat, row, col;
    int64_t mats_center;
    int64_t rows_center;
    int64_t cols_center;

    float inv_mat_sigma_sqr;
    float inv_row_sigma_sqr;
    float inv_col_sigma_sqr;
    float gauss_val;
    float val;

    Matblock_f* h;

    if (mat_sigma < 0.01 || row_sigma < 0.01 || col_sigma < 0.01)
    {
        return JWSC_EARG("Sigma for Gaussian is too small");
    }

    create_matblock_f(h_out, num_mats, num_rows, num_cols);
    h = *h_out;

    mats_center       = num_mats / 2;
    rows_center       = num_rows / 2;
    cols_center       = num_cols / 2;
    inv_mat_sigma_sqr = 1.0 / (mat_sigma*mat_sigma);
    inv_row_sigma_sqr = 1.0 / (row_sigma*row_sigma);
    inv_col_sigma_sqr = 1.0 / (col_sigma*col_sigma);

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                val = (mat-mats_center)*(mat-mats_center)*inv_mat_sigma_sqr +
                      (row-rows_center)*(row-rows_center)*inv_row_sigma_sqr + 
                      (col-cols_center)*(col-cols_center)*inv_col_sigma_sqr;
                gauss_val = -1*(mat - mats_center)*exp(-0.5 * val);
                h->elts[ mat ][ row ][ col ] = gauss_val;
            }
        }
    }

    val = (inv_mat_sigma_sqr*inv_row_sigma_sqr*inv_col_sigma_sqr) / (2 * PI);
    multiply_matblock_by_scalar_f(h_out, h, val);

    return NULL;
}

Error* create_3d_gaussian_dz_filter_d
(
    Matblock_d** h_out, 
    double       mat_sigma,
    double       row_sigma,
    double       col_sigma,
    uint32_t     num_mats,
    uint32_t     num_rows,
    uint32_t     num_cols
)
{
    int64_t mat, row, col;
    int64_t mats_center;
    int64_t rows_center;
    int64_t cols_center;

    double inv_mat_sigma_sqr;
    double inv_row_sigma_sqr;
    double inv_col_sigma_sqr;
    double gauss_val;
    double val;

    Matblock_d* h;

    if (mat_sigma < 0.01 || row_sigma < 0.01 || col_sigma < 0.01)
    {
        return JWSC_EARG("Sigma for Gaussian is too small");
    }

    create_matblock_d(h_out, num_mats, num_rows, num_cols);
    h = *h_out;

    mats_center       = num_mats / 2;
    rows_center       = num_rows / 2;
    cols_center       = num_cols / 2;
    inv_mat_sigma_sqr = 1.0 / (mat_sigma*mat_sigma);
    inv_row_sigma_sqr = 1.0 / (row_sigma*row_sigma);
    inv_col_sigma_sqr = 1.0 / (col_sigma*col_sigma);

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                val = (mat-mats_center)*(mat-mats_center)*inv_mat_sigma_sqr +
                      (row-rows_center)*(row-rows_center)*inv_row_sigma_sqr + 
                      (col-cols_center)*(col-cols_center)*inv_col_sigma_sqr;
                gauss_val = -1*(mat - mats_center)*exp(-0.5 * val);
                h->elts[ mat ][ row ][ col ] = gauss_val;
            }
        }
    }

    val = (inv_mat_sigma_sqr*inv_row_sigma_sqr*inv_col_sigma_sqr) / (2 * PI);
    multiply_matblock_by_scalar_d(h_out, h, val);

    return NULL;
}