2d.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/matrix/matrix.h"
#include "jwsc/matrix/matrix_math.h"
#include "jwsc/filter/2d.h"


#define  PI  3.14159265358979323846


void create_sobel_x_filter_f(Matrix_f** sobel_out)
{
    Matrix_f* sobel;

    create_zero_matrix_f(sobel_out, 3, 3);
    sobel = *sobel_out;

    sobel->elts[ 0 ][ 0 ] = -1;
    sobel->elts[ 0 ][ 2 ] =  1;
    sobel->elts[ 1 ][ 0 ] = -2;
    sobel->elts[ 1 ][ 2 ] =  2;
    sobel->elts[ 2 ][ 0 ] = -1;
    sobel->elts[ 2 ][ 2 ] =  1;
}

void create_sobel_x_filter_d(Matrix_d** sobel_out)
{
    Matrix_d* sobel;

    create_zero_matrix_d(sobel_out, 3, 3);
    sobel = *sobel_out;

    sobel->elts[ 0 ][ 0 ] = -1;
    sobel->elts[ 0 ][ 2 ] =  1;
    sobel->elts[ 1 ][ 0 ] = -2;
    sobel->elts[ 1 ][ 2 ] =  2;
    sobel->elts[ 2 ][ 0 ] = -1;
    sobel->elts[ 2 ][ 2 ] =  1;
}

void create_sobel_y_filter_f(Matrix_f** sobel_out)
{
    Matrix_f* sobel;

    create_zero_matrix_f(sobel_out, 3, 3);
    sobel = *sobel_out;

    sobel->elts[ 0 ][ 0 ] =  1;
    sobel->elts[ 2 ][ 0 ] = -1;
    sobel->elts[ 0 ][ 1 ] =  2;
    sobel->elts[ 2 ][ 1 ] = -2;
    sobel->elts[ 0 ][ 2 ] =  1;
    sobel->elts[ 2 ][ 2 ] = -1;
}

void create_sobel_y_filter_d(Matrix_d** sobel_out)
{
    Matrix_d* sobel;

    create_zero_matrix_d(sobel_out, 3, 3);
    sobel = *sobel_out;

    sobel->elts[ 0 ][ 0 ] =  1;
    sobel->elts[ 2 ][ 0 ] = -1;
    sobel->elts[ 0 ][ 1 ] =  2;
    sobel->elts[ 2 ][ 1 ] = -2;
    sobel->elts[ 0 ][ 2 ] =  1;
    sobel->elts[ 2 ][ 2 ] = -1;
}

Error* create_LoG_filter_f(Matrix_f** LoG_out, float sigma)
{
    uint32_t width, center;
    uint32_t r, c;

    float   x_2, y_2;
    float   sigma_sqr;

    Matrix_f* LoG;

    if (sigma < 0.01)
    {
        return JWSC_EARG("Sigma for LoG is too small");
    }

    if ((width = ceil(6.0 * sqrt(2.0) * sigma)) % 2 == 0) width++;
    center = width / 2;

    create_matrix_f(LoG_out, width, width);
    LoG = *LoG_out;

    sigma_sqr = sigma*sigma;

    for (r = 0; r < width; r++)
    {
        y_2 = (r - center)*(r - center);
        for (c = 0; c < width; c++)
        {
            x_2 = (c - center)*(c - center);

            LoG->elts[ r ][ c ] = (x_2 + y_2 - 2*sigma_sqr) *
                exp(-0.5*(x_2 + y_2)/sigma_sqr) / (sigma_sqr*sigma_sqr);
        }
    }

    return NULL;
}

Error* create_auto_2d_gaussian_filter_f
(
    Matrix_f** h_out, 
    float      row_sigma, 
    float      col_sigma
)
{
    uint32_t num_rows, num_cols;

    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_2d_gaussian_filter_f(h_out, row_sigma, col_sigma, 
            num_rows, num_cols);
}

Error* create_auto_2d_gaussian_filter_d
(
    Matrix_d** h_out, 
    double     row_sigma, 
    double     col_sigma
)
{
    uint32_t num_rows, num_cols;

    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_2d_gaussian_filter_d(h_out, row_sigma, col_sigma, 
            num_rows, num_cols);
}

Error* create_2d_gaussian_filter_f
(
    Matrix_f** h_out, 
    float      row_sigma,
    float      col_sigma,
    uint32_t   num_rows,
    uint32_t   num_cols
)
{
    uint32_t  row, col;
    uint32_t  r, c;
    float     x, y;
    float     x_gauss, y_gauss;
    float     gauss;
    Matrix_f* h;

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

    assert(num_rows > 0 && num_cols > 0);

    create_matrix_f(h_out, num_rows, num_cols);
    h = *h_out;

#if defined JWSC_HAVE_SLATEC || defined JWSC_HAVE_ERF
    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;

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

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

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

            r = num_rows-row-1; c = num_cols-col-1;
            h->elts[ r ][ c ] = gauss;
        }
    }
#else
    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;

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

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

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

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

    return NULL;
}

Error* create_2d_gaussian_filter_d
(
    Matrix_d** h_out, 
    double     row_sigma,
    double     col_sigma,
    uint32_t   num_rows,
    uint32_t   num_cols
)
{
    uint32_t  row, col;
    uint32_t  r, c;
    double    x, y;
    double    x_gauss, y_gauss;
    double    gauss;
    Matrix_d* h;

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

    assert(num_rows > 0 && num_cols > 0);

    create_matrix_d(h_out, num_rows, num_cols);
    h = *h_out;

#if defined JWSC_HAVE_SLATEC || defined JWSC_HAVE_ERF
    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;

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

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

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

            r = num_rows-row-1; c = num_cols-col-1;
            h->elts[ r ][ c ] = gauss;
        }
    }
#else
    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;

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

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

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

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

    return NULL;
}

Error* create_auto_2d_gaussian_dx_filter_f
(
    Matrix_f** h_out, 
    float      row_sigma,
    float      col_sigma
)
{
    uint32_t num_rows, num_cols;

    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_2d_gaussian_dx_filter_f(h_out, row_sigma, col_sigma, 
            num_rows, num_cols);
}

Error* create_auto_2d_gaussian_dx_filter_d
(
    Matrix_d** h_out, 
    double     row_sigma,
    double     col_sigma
)
{
    uint32_t num_rows, num_cols;

    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_2d_gaussian_dx_filter_d(h_out, row_sigma, col_sigma, 
            num_rows, num_cols);
}

Error* create_2d_gaussian_dx_filter_f
(
    Matrix_f** h_out, 
    float      row_sigma,
    float      col_sigma,
    uint32_t   num_rows,
    uint32_t   num_cols
)
{
    int64_t row, col;
    int64_t rows_center;
    int64_t cols_center;

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

    Matrix_f* h;

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

    create_matrix_f(h_out, num_rows, num_cols);
    h = *h_out;

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

    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            val       = (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[ row ][ col ] = gauss_val;
        }
    }

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

    return NULL;
}

Error* create_2d_gaussian_dx_filter_d
(
    Matrix_d** h_out, 
    double     row_sigma,
    double     col_sigma,
    uint32_t   num_rows,
    uint32_t   num_cols
)
{
    int64_t row, col;
    int64_t rows_center;
    int64_t cols_center;

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

    Matrix_d* h;

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

    create_matrix_d(h_out, num_rows, num_cols);
    h = *h_out;

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

    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            val       = (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[ row ][ col ] = gauss_val;
        }
    }

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


    return NULL;
}

Error* create_auto_2d_gaussian_dy_filter_f
(
    Matrix_f** h_out, 
    float      row_sigma,
    float      col_sigma
)
{
    uint32_t num_rows, num_cols;

    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_2d_gaussian_dy_filter_f(h_out, row_sigma, col_sigma, 
            num_rows, num_cols);
}

Error* create_auto_2d_gaussian_dy_filter_d
(
    Matrix_d** h_out, 
    double     row_sigma,
    double     col_sigma
)
{
    uint32_t num_rows, num_cols;

    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_2d_gaussian_dy_filter_d(h_out, row_sigma, col_sigma, 
            num_rows, num_cols);
}

Error* create_2d_gaussian_dy_filter_f
(
    Matrix_f** h_out, 
    float      row_sigma,
    float      col_sigma,
    uint32_t   num_rows,
    uint32_t   num_cols
)
{
    int64_t row, col;
    int64_t rows_center;
    int64_t cols_center;

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

    Matrix_f* h;

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

    create_matrix_f(h_out, num_rows, num_cols);
    h = *h_out;

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

    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            val       = (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[ row ][ col ] = gauss_val;
        }
    }

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

    return NULL;
}

Error* create_2d_gaussian_dy_filter_d
(
    Matrix_d** h_out, 
    double     row_sigma,
    double     col_sigma,
    uint32_t   num_rows,
    uint32_t   num_cols
)
{
    int64_t row, col;
    int64_t rows_center;
    int64_t cols_center;

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

    Matrix_d* h;

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

    create_matrix_d(h_out, num_rows, num_cols);
    h = *h_out;

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

    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            val       = (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[ row ][ col ] = gauss_val;
        }
    }

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

    return NULL;
}