matblock_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/matblock/matblock.h"
#include "jwsc/matblock/matblock_math.h"
#include "jwsc/matblock/matblock_fft.h"
#include "jwsc/matblock/matblock_conv.h"


#define  PI  3.14159265358979323846


Error* convolve_matblock_f
(
    Matblock_f**      m_out, 
    const Matblock_f* m_in, 
    const Matblock_f* h
)
{
    int64_t num_in_mats, num_in_rows, num_in_cols; /* 64-bit since signed */
    int64_t num_h_mats, num_h_rows, num_h_cols;
    int64_t in_mat, in_row, in_col;
    int64_t h_mat, h_row, h_col;
    int64_t h_mat_center;
    int64_t h_row_center;
    int64_t h_col_center;
    int64_t i, j, k;

    float       sum;
    Matblock_f* m;

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

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

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

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

    for (in_mat = 0; in_mat < num_in_mats; in_mat++)
    {
        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_mat = 0; h_mat < num_h_mats; h_mat++)
                {
                    for (h_row = 0; h_row < num_h_rows; h_row++)
                    {
                        for (h_col = 0; h_col < num_h_cols; h_col++)
                        {
                            i = in_mat + (h_mat - h_mat_center);
                            if (i < 0)
                            {
                                i = -i;
                            }
                            else if (i >= num_in_mats)
                            {
                                i = num_in_mats - (i - (num_in_mats - 1));
                            }

                            j = in_row + (h_row - h_row_center);
                            if (j < 0)
                            {
                                j = -j;
                            }
                            else if (j >= num_in_rows)
                            {
                                j = num_in_rows - (j - (num_in_rows - 1));
                            }

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

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

    if (*m_out == m_in)
    {
        copy_matblock_f(m_out, m);
        free_matblock_f(m);
    }
    else
    {
        *m_out = m;
    }

    return NULL;
}

Error* convolve_matblock_d
(
    Matblock_d**      m_out, 
    const Matblock_d* m_in, 
    const Matblock_d* h
)
{
    int64_t num_in_mats, num_in_rows, num_in_cols; /* 64-bit since signed */
    int64_t num_h_mats, num_h_rows, num_h_cols;
    int64_t in_mat, in_row, in_col;
    int64_t h_mat, h_row, h_col;
    int64_t h_mat_center;
    int64_t h_row_center;
    int64_t h_col_center;
    int64_t i, j, k;

    double      sum;
    Matblock_d* m;

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

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

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

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

    for (in_mat = 0; in_mat < num_in_mats; in_mat++)
    {
        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_mat = 0; h_mat < num_h_mats; h_mat++)
                {
                    for (h_row = 0; h_row < num_h_rows; h_row++)
                    {
                        for (h_col = 0; h_col < num_h_cols; h_col++)
                        {
                            i = in_mat + (h_mat - h_mat_center);
                            if (i < 0)
                            {
                                i = -i;
                            }
                            else if (i >= num_in_mats)
                            {
                                i = num_in_mats - (i - (num_in_mats - 1));
                            }

                            j = in_row + (h_row - h_row_center);
                            if (j < 0)
                            {
                                j = -j;
                            }
                            else if (j >= num_in_rows)
                            {
                                j = num_in_rows - (j - (num_in_rows - 1));
                            }

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

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

    if (*m_out == m_in)
    {
        copy_matblock_d(m_out, m);
        free_matblock_d(m);
    }
    else
    {
        *m_out = m;
    }

    return NULL;
}

#ifdef JWSC_HAVE_FFTW3F
Error* convolve_matblock_fft_f
(
    Matblock_f**      m_out, 
    const Matblock_f* m_in, 
    const Matblock_f* h
)
{
    Matblock_cf* s_1 = NULL;
    Matblock_cf* s_2 = NULL;
    Matblock_f*  s_3 = NULL;
    Matblock_cf* s_4 = NULL;

    convolve_matblock_fft_scratch_f(m_out, m_in, h, &s_1, &s_2, &s_3, &s_4);

    free_matblock_cf(s_1);
    free_matblock_cf(s_2);
    free_matblock_f(s_3);
    free_matblock_cf(s_4);

    return NULL;
}
#endif

#ifdef JWSC_HAVE_FFTW3F
Error* convolve_matblock_fft_d
(
    Matblock_d**      m_out, 
    const Matblock_d* m_in, 
    const Matblock_d* h
)
{
    Matblock_cd* s_1 = NULL;
    Matblock_cd* s_2 = NULL;
    Matblock_d*  s_3 = NULL;
    Matblock_cd* s_4 = NULL;

    convolve_matblock_fft_scratch_d(m_out, m_in, h, &s_1, &s_2, &s_3, &s_4);

    free_matblock_cd(s_1);
    free_matblock_cd(s_2);
    free_matblock_d(s_3);
    free_matblock_cd(s_4);

    return NULL;
}
#endif

#ifdef JWSC_HAVE_FFTW3F
Error* convolve_matblock_fft_scratch_f
(
    Matblock_f**      m_out, 
    const Matblock_f* m_in, 
    const Matblock_f* h,
    Matblock_cf**     s_1,
    Matblock_cf**     s_2,
    Matblock_f**      s_3,
    Matblock_cf**     s_4
)
{
    uint32_t     num_in_mats, num_in_rows, num_in_cols;
    uint32_t     num_h_mats, num_h_rows, num_h_cols;

    num_in_mats = m_in->num_mats;
    num_in_rows = m_in->num_rows;
    num_in_cols = m_in->num_cols;
    num_h_mats  = h->num_mats;
    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_mats != num_h_mats ||
        num_in_rows != num_h_rows || 
        num_in_cols != num_h_cols)
    {
        if (*m_out != m_in)
        {
            free_matblock_f(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Incompatible matblock sizes for FFT convolution");
    }

    real_to_complex_matblock_dft_scratch_f(s_1, m_in, s_3);
    real_to_complex_matblock_dft_scratch_f(s_2, h, s_3);
    assert(multiply_matblocks_elt_wise_cf(s_1, *s_1, *s_2) == NULL);
    complex_to_real_matblock_dft_scratch_f(m_out, *s_1, 
            (num_in_cols % 2 == 0), s_4);
    scale_matblock_dft_f(m_out, *m_out);

    return NULL;
}
#endif

#ifdef JWSC_HAVE_FFTW3F
Error* convolve_matblock_fft_scratch_d
(
    Matblock_d**      m_out, 
    const Matblock_d* m_in, 
    const Matblock_d* h,
    Matblock_cd**     s_1,
    Matblock_cd**     s_2,
    Matblock_d**      s_3,
    Matblock_cd**     s_4
)
{
    uint32_t     num_in_mats, num_in_rows, num_in_cols;
    uint32_t     num_h_mats, num_h_rows, num_h_cols;

    num_in_mats = m_in->num_mats;
    num_in_rows = m_in->num_rows;
    num_in_cols = m_in->num_cols;
    num_h_mats  = h->num_mats;
    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_mats != num_h_mats ||
        num_in_rows != num_h_rows || 
        num_in_cols != num_h_cols)
    {
        if (*m_out != m_in)
        {
            free_matblock_d(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Incompatible matblock size for FFT convolution");
    }

    real_to_complex_matblock_dft_scratch_d(s_1, m_in, s_3);
    real_to_complex_matblock_dft_scratch_d(s_2, h, s_3);
    assert(multiply_matblocks_elt_wise_cd(s_1, *s_1, *s_2) == NULL);
    complex_to_real_matblock_dft_scratch_d(m_out, *s_1, 
            (num_in_cols % 2 == 0), s_4);
    scale_matblock_dft_d(m_out, *m_out);

    return NULL;
}
#endif