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


Error* down_sample_image_f
(
    Image_f**      img_out, 
    const Image_f* img_in, 
    float          row_factor,
    float          col_factor
)
{
    int64_t num_in_rows, num_in_cols;
    int64_t num_out_rows, num_out_cols;
    int64_t num_h_rows, num_h_cols;
    int64_t in_row, in_col;
    int64_t out_row, out_col;
    int64_t h_row, h_col;
    int64_t h_row_center;
    int64_t h_col_center;
    int64_t i, j;

    float col_delta_sum, row_delta_sum;
    float r_sum, g_sum, b_sum;
    float row_delta;
    float col_delta;
    float row_sigma;
    float col_sigma;
    float weight;
    float weight_sum;

    Image_f*   img = NULL;
    Matrix_f*  h   = NULL;

    /* Test if the factors are in (0, 1]. */
    if (0 >= row_factor || row_factor > 1.0 ||
        0 >= col_factor || col_factor > 1.0)
    {
        if (*img_out != img_in)
        {
            free_image_f(*img_out);
            *img_out = NULL;
        }
        return JWSC_EARG("Down sampling factors must be in (0, 1]");
    }

    num_in_rows = img_in->num_rows;
    num_in_cols = img_in->num_cols;

    num_out_rows = ceil(row_factor*num_in_rows);
    num_out_cols = ceil(col_factor*num_in_cols);

    if (num_in_rows == num_out_rows && num_in_cols == num_out_cols &&
            *img_out != img_in)
    {
        copy_image_f(img_out, img_in);
        return  NULL;
    }

    row_delta = 1.0 / row_factor;
    col_delta = 1.0 / col_factor;

    row_sigma = 0.35 * row_delta;
    col_sigma = 0.35 * col_delta;

    assert(create_auto_2d_gaussian_filter_f(&h, row_sigma, col_sigma) == NULL);
    num_h_rows = h->num_rows;
    num_h_cols = h->num_cols;

    img = (*img_out == img_in) ? NULL : *img_out;
    create_image_f(&img, num_out_rows, num_out_cols);

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

    row_delta_sum = 0.5*row_delta;
    in_row = ceil(row_delta_sum) - 1;
    for (out_row = 0; out_row < num_out_rows; out_row++)
    {
        col_delta_sum = 0.5*col_delta;
        in_col = ceil(col_delta_sum) - 1;
        for (out_col = 0; out_col < num_out_cols; out_col++)
        {
            r_sum = 0;
            g_sum = 0;
            b_sum = 0;
            weight_sum = 0;

            i = -1;
            j = -1;

            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 = 0;
                    }
                    else if (i >= num_in_rows)
                    {
                        i = num_in_rows - 1;
                    }

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

                    weight = h->elts[ h_row ][ h_col ];

                    r_sum += img_in->pxls[ i ][ j ].r * weight;
                    g_sum += img_in->pxls[ i ][ j ].g * weight;
                    b_sum += img_in->pxls[ i ][ j ].b * weight;

                    weight_sum += weight;
                }
            }

            if (weight_sum > 0)
            {
                r_sum /= weight_sum;
                g_sum /= weight_sum;
                b_sum /= weight_sum;

                img->pxls[ out_row ][ out_col ].r = r_sum;
                img->pxls[ out_row ][ out_col ].g = g_sum;
                img->pxls[ out_row ][ out_col ].b = b_sum;
            }
            else
            {
                img->pxls[ out_row ][ out_col ].r = img_in->pxls[ i ][ j ].r;
                img->pxls[ out_row ][ out_col ].g = img_in->pxls[ i ][ j ].g;
                img->pxls[ out_row ][ out_col ].b = img_in->pxls[ i ][ j ].b;
            }

            col_delta_sum += col_delta;
            in_col = floor(col_delta_sum);
            if (in_col >= num_in_cols) 
                in_col = num_in_cols - 1;
        }

        row_delta_sum += row_delta;
        in_row = floor(row_delta_sum);
        if (in_row >= num_in_rows) 
            in_row = num_in_rows - 1;
    }

    if (*img_out == img_in)
    {
        copy_image_f(img_out, img);
        free_image_f(img);
    }
    else
    {
        *img_out = img;
    }

    free_matrix_f(h);

    return NULL;
}

void extend_image_f
(
    Image_f**      img_out, 
    const Image_f* img_in, 
    uint32_t       row_padding,
    uint32_t       col_padding
)
{
    uint32_t num_rows, num_cols;
    uint32_t num_prows, num_pcols;
    uint32_t prow, pcol;

    Image_f* img;

    num_rows = img_in->num_rows;
    num_cols = img_in->num_cols;

    num_prows = num_rows + 2*row_padding;
    num_pcols = num_cols + 2*col_padding;

    img = (*img_out == img_in) ? NULL : *img_out;
    create_image_f(&img, num_prows, num_pcols);

    for (prow = 0; prow < (num_prows - row_padding); prow++)
    {
        if (prow < row_padding)
        {
            for (pcol = 0; pcol < col_padding; pcol++)
            {
                img->pxls[ prow ][ pcol ] 
                    = img_in->pxls[ 0 ][ 0 ];

                img->pxls[ prow ][ num_pcols-pcol-1 ] 
                    = img_in->pxls[ 0 ][ num_cols-1 ];

                img->pxls[ num_prows-prow-1 ][ pcol ] 
                    = img_in->pxls[ num_rows-1 ][ 0 ];

                img->pxls[ num_prows-prow-1 ][ num_pcols-pcol-1 ] 
                    = img_in->pxls[ num_rows-1 ][ num_cols-1 ];
            }
            for (pcol = col_padding; pcol < (num_pcols-col_padding); pcol++)
            {
                img->pxls[ prow ][ pcol ] 
                    = img_in->pxls[ 0 ][ pcol-col_padding ];
                img->pxls[ num_prows-prow-1 ][ pcol ] 
                    = img_in->pxls[ num_rows-1 ][ pcol-col_padding ];
            }
        }
        else
        {
            for (pcol = 0; pcol < col_padding; pcol++)
            {
                img->pxls[ prow ][ pcol ] 
                    = img_in->pxls[ prow-row_padding ][ 0 ];
                img->pxls[ prow ][ num_pcols-pcol-1 ] 
                    = img_in->pxls[ prow-row_padding ][ num_cols-1 ];
            }
            for (pcol = col_padding; pcol < (num_pcols-col_padding); pcol++)
            {
                img->pxls[ prow ][ pcol ] 
                    = img_in->pxls[ prow-row_padding][ pcol-col_padding ];
            }
        }
    }

    *img_out = img;
}

Error* crop_image_f
(
    Image_f**      img_out, 
    const Image_f* img_in, 
    uint32_t       row,
    uint32_t       col,
    uint32_t       num_rows,
    uint32_t       num_cols
)
{
    uint32_t num_in_rows;
    uint32_t num_in_cols;
    uint32_t rrow, ccol;
    uint32_t r, c;

    Image_f* img = NULL;

    num_in_rows = img_in->num_rows;
    num_in_cols = img_in->num_cols;

    if (row == num_in_rows || col == num_in_cols ||
        row+num_rows > num_in_rows || col+num_cols > num_in_cols)
    {
        free_image_f(*img_out);
        *img_out = NULL;
        return JWSC_EARG("Invalid image region to crop");
    }

    img = (*img_out == img_in) ? NULL : *img_out;
    create_image_f(&img, num_rows, num_cols);

    r = 0;

    for (rrow = row; rrow < row+num_rows; rrow++)
    {
        c = 0;

        for (ccol = col; ccol < col+num_cols; ccol++)
        {
            img->pxls[ r ][ c++ ] = img_in->pxls[ rrow ][ ccol ];
        }

        r++;
    }

    if (*img_out == img_in)
    {
        copy_image_f(img_out, img);
    }
    else
    {
        *img_out = img;
    }

    return NULL;
}

void threshold_image_f
(
    Matrix_f**     mat_out, 
    const Image_f* img_in,
    float          t
)
{
    create_matrix_from_image_f(mat_out, img_in, 0.3, 0.59, 0.11);
    threshold_matrix_f(mat_out, *mat_out, t);
}

void threshold_matrix_f
(
    Matrix_f**      mat_out, 
    const Matrix_f* mat_in,
    float           t
)
{
    uint32_t row, num_rows;
    uint32_t col, num_cols;

    num_rows = mat_in->num_rows;
    num_cols = mat_in->num_cols;

    create_matrix_f(mat_out, num_rows, num_cols);

    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            if (mat_in->elts[ row ][ col ] >= t)
            {
                (*mat_out)->elts[ row ][ col ] = 1.0f;
            }
            else
            {
                (*mat_out)->elts[ row ][ col ] = 0.0f;
            }
        }
    }
}

void dynamically_threshold_image_f
(
    Matrix_f**     mat_out, 
    const Image_f* img_in
)
{
    create_matrix_from_image_f(mat_out, img_in, 0.3, 0.59, 0.11);
    dynamically_threshold_matrix_f(mat_out, *mat_out);
}

void dynamically_threshold_matrix_f
(
    Matrix_f**      mat_out, 
    const Matrix_f* mat_in
)
{
    static const float epsilon = 1.0e-12;

    uint32_t row, num_rows;
    uint32_t col, num_cols;
    uint32_t num_fg;
    uint32_t num_bg;
    float min, max;
    float x;
    float t;
    float t_prev;
    float mean_fg;
    float mean_bg;

    num_rows = mat_in->num_rows;
    num_cols = mat_in->num_cols;

    max = min = mat_in->elts[0][0];

    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            if (mat_in->elts[ row ][ col ] > max)
            {
                max = mat_in->elts[ row ][ col ];
            }
            else if (mat_in->elts[ row ][ col ] < min)
            {
                min = mat_in->elts[ row ][ col ];
            }
        }
    }
    assert((max - min) >= 0.0);

    t = 0.5f*(max - min) + min;
    t_prev = t - min;

    while (fabs(t - t_prev) > epsilon)
    {
        t_prev = t;

        num_fg = num_bg = 0;
        mean_fg = mean_bg = 0.0;

        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                x = mat_in->elts[ row ][ col ];

                if (x < t)
                {
                    num_bg++;
                    mean_bg += x;
                }
                else
                {
                    num_fg++;
                    mean_fg += x;
                }
            }
        }

        if (num_bg)
        {
            mean_bg *= 1.0f / (float)num_bg;
        }
        if (num_fg)
        {
            mean_fg *= 1.0f / (float)num_fg;
        }

        t = 0.5f * (mean_fg + mean_bg);
    }

    threshold_matrix_f(mat_out, mat_in, t);
}