matrix_conv.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/matrix/matrix.h"
#include "jwsc/matrix/matrix_math.h"
#include "jwsc/matrix/matrix_fft.h"
#include "jwsc/matrix/matrix_conv.h"


#define  PI  3.14159265358979323846


Error* convolve_matrix_f
(
    Matrix_f**      m_out, 
    const Matrix_f* m_in, 
    const Matrix_f* h
)
{
    int64_t num_in_rows, num_in_cols; /* 64-bit since it's signed */
    int64_t num_h_rows, num_h_cols;
    int64_t in_row, in_col;
    int64_t h_row, h_col;
    int64_t h_row_center;
    int64_t h_col_center;
    int64_t i, j;

    float     sum;
    Matrix_f* m;

    num_in_rows = m_in->num_rows;
    num_in_cols = m_in->num_cols;
    num_h_rows  = h->num_rows;
    num_h_cols  = h->num_cols;

    /* Test if the filter is larger than m_in. */
    if (num_in_rows < num_h_rows || num_in_cols < num_h_cols)
    {
        if (*m_out != m_in)
        {
            free_matrix_f(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Convolution filter larger than input matrix");
    }

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_f(&m, num_in_rows, num_in_cols);

    h_row_center = num_h_rows / 2;
    h_col_center = num_h_cols / 2;

    for (in_row = 0; in_row < num_in_rows; in_row++)
    {
        for (in_col = 0; in_col < num_in_cols; in_col++)
        {
            sum = 0;
            for (h_row = 0; h_row < num_h_rows; h_row++)
            {
                for (h_col = 0; h_col < num_h_cols; h_col++)
                {
                    i = in_row + (h_row - h_row_center);
                    if (i < 0)
                    {
                        i = -i;
                    }
                    else if (i >= num_in_rows)
                    {
                        i = num_in_rows - (i - (num_in_rows - 1));
                    }

                    j = in_col + (h_col - h_col_center);
                    if (j < 0)
                    {
                        j = -j;
                    }
                    else if (j >= num_in_cols)
                    {
                        j = num_in_cols - (j - (num_in_cols - 1));
                    }

                    sum += m_in->elts[ i ][ j ] * h->elts[ h_row ][ h_col ];
                }
            }
            m->elts[ in_row ][ in_col ] = sum;
        }
    }

    if (*m_out == m_in)
    {
        copy_matrix_f(m_out, m);
        free_matrix_f(m);
    }
    else
    {
        *m_out = m;
    }

    return NULL;
}


Error* convolve_matrix_d
(
    Matrix_d**      m_out, 
    const Matrix_d* m_in, 
    const Matrix_d* h
)
{
    int64_t num_in_rows, num_in_cols; /* 64-bit since it's signed */
    int64_t num_h_rows, num_h_cols;
    int64_t in_row, in_col;
    int64_t h_row, h_col;
    int64_t h_row_center;
    int64_t h_col_center;
    int64_t i, j;

    double    sum;
    Matrix_d* m;

    num_in_rows = m_in->num_rows;
    num_in_cols = m_in->num_cols;
    num_h_rows  = h->num_rows;
    num_h_cols  = h->num_cols;

    /* Test if the filter is larger than m_in. */
    if (num_in_rows < num_h_rows || num_in_cols < num_h_cols)
    {
        if (*m_out != m_in)
        {
            free_matrix_d(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Convolution filter larger than input matrix");
    }

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_d(&m, num_in_rows, num_in_cols);

    h_row_center = num_h_rows / 2;
    h_col_center = num_h_cols / 2;

    for (in_row = 0; in_row < num_in_rows; in_row++)
    {
        for (in_col = 0; in_col < num_in_cols; in_col++)
        {
            sum = 0;
            for (h_row = 0; h_row < num_h_rows; h_row++)
            {
                for (h_col = 0; h_col < num_h_cols; h_col++)
                {
                    i = in_row + (h_row - h_row_center);
                    if (i < 0)
                    {
                        i = -i;
                    }
                    else if (i >= num_in_rows)
                    {
                        i = num_in_rows - (i - (num_in_rows - 1));
                    }

                    j = in_col + (h_col - h_col_center);
                    if (j < 0)
                    {
                        j = -j;
                    }
                    else if (j >= num_in_cols)
                    {
                        j = num_in_cols - (j - (num_in_cols - 1));
                    }

                    sum += m_in->elts[ i ][ j ] * h->elts[ h_row ][ h_col ];
                }
            }
            m->elts[ in_row ][ in_col ] = sum;
        }
    }

    if (*m_out == m_in)
    {
        copy_matrix_d(m_out, m);
        free_matrix_d(m);
    }
    else
    {
        *m_out = m;
    }

    return NULL;
}

#ifdef JWSC_HAVE_FFTW3F
Error* convolve_matrix_fft_f
(
    Matrix_f**      m_out, 
    const Matrix_f* m_in, 
    const Matrix_f* h
)
{
    uint32_t   num_in_rows, num_in_cols;
    uint32_t   num_h_rows, num_h_cols;
    Matrix_cf* dft_m_in;
    Matrix_cf* dft_h;

    num_in_rows = m_in->num_rows;
    num_in_cols = m_in->num_cols;
    num_h_rows  = h->num_rows;
    num_h_cols  = h->num_cols;

    /* Test if the filter is the same size as m_in. */
    if (num_in_rows != num_h_rows || num_in_cols != num_h_cols)
    {
        if (*m_out != m_in)
        {
            free_matrix_f(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Incompatible matrix sizes for FFT convolution");
    }

    dft_m_in = NULL;
    dft_h     = NULL;

    real_to_complex_matrix_dft_f(&dft_m_in, m_in);
    real_to_complex_matrix_dft_f(&dft_h, h);
    assert(multiply_matrices_elt_wise_cf(&dft_m_in, dft_m_in, dft_h) == NULL);
    complex_to_real_matrix_dft_f(m_out, dft_m_in, (num_in_cols % 2 == 0));
    scale_matrix_dft_f(m_out, *m_out);

    free_matrix_cf(dft_m_in);
    free_matrix_cf(dft_h);

    return NULL;
}
#endif

#ifdef JWSC_HAVE_FFTW3F
Error* convolve_matrix_fft_d
(
    Matrix_d**      m_out, 
    const Matrix_d* m_in, 
    const Matrix_d* h
)
{
    uint32_t   num_in_rows, num_in_cols;
    uint32_t   num_h_rows, num_h_cols;
    Matrix_cd* dft_m_in;
    Matrix_cd* dft_h;

    num_in_rows = m_in->num_rows;
    num_in_cols = m_in->num_cols;
    num_h_rows  = h->num_rows;
    num_h_cols  = h->num_cols;

    /* Test if the filter is the same size as m_in. */
    if (num_in_rows != num_h_rows || num_in_cols != num_h_cols)
    {
        if (*m_out != m_in)
        {
            free_matrix_d(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Incompatible matrix sizes for FFT convolution");
    }

    dft_m_in = NULL;
    dft_h     = NULL;

    real_to_complex_matrix_dft_d(&dft_m_in, m_in);
    real_to_complex_matrix_dft_d(&dft_h, h);
    assert(multiply_matrices_elt_wise_cd(&dft_m_in, dft_m_in, dft_h) == NULL);
    complex_to_real_matrix_dft_d(m_out, dft_m_in, (num_in_cols % 2 == 0));
    scale_matrix_dft_d(m_out, *m_out);

    free_matrix_cd(dft_m_in);
    free_matrix_cd(dft_h);

    return NULL;
}
#endif