edge.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 <stdio.h>
#include <inttypes.h>
#include <assert.h>
#include <math.h>
#include <string.h>
#include <errno.h>

#include "jwsc/base/error.h"
#include "jwsc/math/constants.h"
#include "jwsc/vector/vector.h"
#include "jwsc/vector/vector.h"
#include "jwsc/matrix/matrix.h"
#include "jwsc/matrix/matrix_fft.h"
#include "jwsc/matrix/matrix_conv.h"
#include "jwsc/matrix/matrix_util.h"
#include "jwsc/filter/2d.h"
#include "jwsc/prob/pmf.h"
#include "jwsc/image/edge.h"

#ifdef JWSC_HAVE_DMALLOC
#include <dmalloc.h>
#endif


typedef struct Edge_list_f
{
    Edge_point_f* pt;

    struct Edge_list_f* next;
}
Edge_list_f;


static void free_edge_list_f(Edge_list_f* edges)
{
    Edge_list_f* next;

    next = edges->next;

    free(edges->pt);
    free(edges);

    if (next)
    {
        free_edge_list_f(next);
    }
}


static void copy_edge_list_into_set_f(Edge_set_f* set, Edge_list_f* list)
{
    uint32_t      num_pts;
    uint32_t      e;
    Edge_list_f*  iter;
    Edge_point_f* pt;
    Edge_point_f* pts;

    set->num_edges++;

    assert(set->num_pts = realloc(set->num_pts, 
                set->num_edges * sizeof(uint32_t)));

    num_pts = 0;
    iter = list;
    while (iter != NULL)
    {
        num_pts++;
        iter = iter->next;
    }
    set->num_pts[ set->num_edges - 1 ] = num_pts;
    set->total_num_pts += num_pts;

    if (set->num_edges == 1)
    {
        assert(set->pts = malloc(set->num_edges * sizeof(Edge_point_f*)));
        assert(*(set->pts) = malloc(set->total_num_pts * sizeof(Edge_point_f)));
    }
    else
    {
        assert(set->pts = realloc(set->pts, 
                    set->num_edges * sizeof(Edge_point_f*)));
        assert(*(set->pts) = realloc(*(set->pts), 
                    set->total_num_pts * sizeof(Edge_point_f)));

        pts = *(set->pts);
        for (e = 1; e < set->num_edges; e++)
        {
            pts += set->num_pts[ e-1 ];
            set->pts[ e ] = pts;
        }
    }

    pt   = set->pts[ set->num_edges - 1 ];
    iter = list;
    while (iter != NULL)
    {
        *pt = *(iter->pt);

        pt++;
        iter = iter->next;
    }
}


static Edge_list_f* merge_edge_lists_f
(
    Edge_list_f* left, 
    Edge_list_f* right
)
{
    Edge_list_f* iter;
    Edge_list_f* next;
    Edge_list_f* prev;

    assert(left != NULL || right != NULL);

    if (left == NULL)
    {
        return right;
    }
    else if (right == NULL)
    {
        return left;
    }

    prev = right;
    iter = left;
    while (iter != NULL)
    {
        next = iter->next;
        iter->next = prev;
        prev = iter;
        iter = next;
    }

    return prev;
}



typedef struct
{
    float dcol;

    float drow;

    float mag;

    uint8_t marked;
}
Gradient_f;


static Error* create_gradient_map_f
(
    Gradient_f***   map_out,
    const Matrix_f* m,
    float           sigma
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    float    dcol, drow;

    Matrix_f*    gauss  = NULL;
    Matrix_f*    m_dcol = NULL;
    Matrix_f*    m_drow = NULL;
    Gradient_f** map;
    Gradient_f*  map_elts;
    Error*       e;

    num_rows = m->num_rows;
    num_cols = m->num_cols;

    assert(*map_out == NULL);

#ifdef JWSC_HAVE_FFTW3F
    if ((e = create_2d_gaussian_dx_filter_f(&gauss, sigma, sigma,
                    num_rows, num_cols)) != NULL)
    {
        return e;
    }
    shift_matrix_f(&gauss, gauss);
    assert(convolve_matrix_fft_f(&m_dcol, m, gauss) == NULL);

    if ((e = create_2d_gaussian_dy_filter_f(&gauss, sigma, sigma,
                    num_rows, num_cols)) != NULL)
    {
        return e;
    }
    shift_matrix_f(&gauss, gauss);
    assert(convolve_matrix_fft_f(&m_drow, m, gauss) == NULL);
#else
    if ((e = create_auto_2d_gaussian_dx_filter_f(&gauss, sigma, sigma)) 
            != NULL)
    {
        return e;
    }
    assert(convolve_matrix_f(&m_dcol, m, gauss) == NULL);

    if ((e = create_auto_2d_gaussian_dy_filter_f(&gauss, sigma, sigma)) 
            != NULL)
    {
        return e;
    }
    assert(convolve_matrix_f(&m_drow, m, gauss) == NULL);
#endif

    free_matrix_f(gauss);

    assert(map_elts = malloc(num_rows*num_cols*sizeof(Gradient_f)));
    assert(*map_out = malloc(num_rows*sizeof(Gradient_f*)));
    map = *map_out;

    for (row = 0; row < num_rows; row++)
    {
        map[ row ] = &(map_elts[ row*num_cols ]);
    }

    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            dcol = m_dcol->elts[ row ][ col ];
            drow = m_drow->elts[ row ][ col ];

            map[ row ][ col ].dcol   = dcol;
            map[ row ][ col ].drow   = drow;
            map[ row ][ col ].mag    = sqrt(dcol*dcol + drow*drow);
            map[ row ][ col ].marked = 0;

            if (map[ row ][ col ].mag >= 1.0e-8)
            {
                map[ row ][ col ].dcol /= map[ row ][ col ].mag;
                map[ row ][ col ].drow /= map[ row ][ col ].mag;
            }
        }
    }

    free_matrix_f(m_dcol);
    free_matrix_f(m_drow);

    return NULL;
}


static Gradient_f* follow_gradient_f
(
    int32_t*     grad_row_in_out, 
    int32_t*     grad_col_in_out,
    float        dir_x,
    float        dir_y,
    Gradient_f** map,
    uint32_t     num_rows,
    uint32_t     num_cols,
    uint32_t     padding
)
{
    static const int neighbor_positions[ 8 ][ 2 ] = { {  1,  0 }, 
                                                {  1,  1 }, 
                                                {  0,  1 }, 
                                                { -1,  1 }, 
                                                { -1,  0 }, 
                                                { -1, -1 }, 
                                                {  0, -1 }, 
                                                {  1, -1 } };
    static const float neighbor_vectors[ 8 ][ 2 ] = { {  1, 0 }, 
                                                {  0.70710678f,  0.70710678f },
                                                {  0,            1 },
                                                { -0.70710678f,  0.70710678f },
                                                { -1,            0 },
                                                { -0.70710678f, -0.70710678f },
                                                {  0,           -1 },
                                                {  0.70710678f, -0.70710678f }};

    uint8_t     i;
    uint8_t     max_neighbor;
    int32_t     grad_col, grad_row;
    float       max_dp;
    float       dp;

    Gradient_f* grad;

    grad_col = *grad_col_in_out;
    grad_row = *grad_row_in_out;

    max_neighbor = 0;
    max_dp       = -1.0f;
    grad         = &map[ grad_row ][ grad_col ];

    for (i = 0; i < 8; i++)
    {
        dp = dir_x * neighbor_vectors[ i ][ 0 ] +
             dir_y * neighbor_vectors[ i ][ 1 ];

        /* Should have some stochastic tie-breaker here. */
        if (dp >= max_dp)
        {
            max_neighbor = i;
            max_dp       = dp;
        }
    }

    grad_row += neighbor_positions[ max_neighbor ][ 1 ];
    grad_col += neighbor_positions[ max_neighbor ][ 0 ];

    if (grad_row == padding-1) 
        grad_row = padding;
    else if (grad_row == num_rows-padding) 
        grad_row = num_rows-padding-1;
   
    if (grad_col == padding-1) 
        grad_col = padding;
    else if (grad_col == num_cols-padding) 
        grad_col = num_cols-padding-1;

    *grad_col_in_out = grad_col;
    *grad_row_in_out = grad_row;

    return &map[ grad_row ][ grad_col ];
}


static Gradient_f* get_maximal_along_gradient_f
(
    int32_t*     grad_row_in_out, 
    int32_t*     grad_col_in_out, 
    Gradient_f** map,
    uint32_t     num_rows,
    uint32_t     num_cols,
    uint32_t     padding
) 
{
    int32_t grad_row, grad_col;
    int32_t curr_grad_row, curr_grad_col;
    int32_t next_grad_row, next_grad_col;
    int32_t prev_grad_row, prev_grad_col;

    Gradient_f* curr_grad;
    Gradient_f* next_grad;
    Gradient_f* prev_grad;
    Gradient_f* max_grad;

    grad_row = *grad_row_in_out;
    grad_col = *grad_col_in_out;

    curr_grad_row = grad_row;
    curr_grad_col = grad_col;
    curr_grad     = &map[ curr_grad_row ][ curr_grad_col ];

    prev_grad_row = curr_grad_row;
    prev_grad_col = curr_grad_col;
    prev_grad     = follow_gradient_f(&prev_grad_row, &prev_grad_col,
                                      -(curr_grad->dcol), -(curr_grad->drow),
                                      map, num_rows, num_cols, padding);

    next_grad_row = curr_grad_row;
    next_grad_col = curr_grad_col;
    next_grad     = follow_gradient_f(&next_grad_row, &next_grad_col,
                                      curr_grad->dcol, curr_grad->drow,
                                      map, num_rows, num_cols, padding);

    max_grad = NULL;
    while (max_grad == NULL)
    {
        if (curr_grad->mag >= next_grad->mag && 
            curr_grad->mag >= prev_grad->mag)
        {
            max_grad = curr_grad;
        }
        else if (curr_grad->mag >= prev_grad->mag)
        {
            curr_grad_row = next_grad_row;
            curr_grad_col = next_grad_col;
            curr_grad     = next_grad;
        }
        else
        {
            curr_grad_row = prev_grad_row;
            curr_grad_col = prev_grad_col;
            curr_grad     = prev_grad;
        }

        prev_grad_row = curr_grad_row;
        prev_grad_col = curr_grad_col;
        prev_grad     = follow_gradient_f(&prev_grad_row, &prev_grad_col,
                -(curr_grad->dcol), -(curr_grad->drow), map, 
                num_rows, num_cols, padding);

        next_grad_row = curr_grad_row;
        next_grad_col = curr_grad_col;
        next_grad     = follow_gradient_f(&next_grad_row, &next_grad_col,
                                          curr_grad->dcol, curr_grad->drow, 
                                          map, num_rows, num_cols, padding);
    }

    *grad_row_in_out = curr_grad_row;
    *grad_col_in_out = curr_grad_col;

    return max_grad;
}


static Edge_list_f* follow_edge_right_f
(
    float        dir_x, 
    float        dir_y, 
    int32_t      grad_row,
    int32_t      grad_col, 
    Gradient_f** map,
    uint32_t     num_rows,
    uint32_t     num_cols,
    uint32_t     padding,
    float        thresh
)
{
    Gradient_f*  grad_curr = NULL;
    Gradient_f*  grad_next = NULL;
    Edge_list_f* node = NULL;

    grad_curr = &(map[ grad_row ][ grad_col ]);
    
    follow_gradient_f(&grad_row, &grad_col, dir_x, dir_y, map, 
            num_rows, num_cols, padding);

    grad_next = get_maximal_along_gradient_f(&grad_row, &grad_col, map, 
            num_rows, num_cols, padding);

    if (!grad_next->marked && grad_next->mag >= thresh &&
        grad_row >= padding && grad_row < num_rows-padding &&
        grad_col >= padding && grad_col < num_cols-padding)
    {
        grad_next->marked = 1;

        assert(node = malloc(sizeof(Edge_list_f)));

        assert(node->pt = malloc(sizeof(Edge_point_f)));
        node->pt->row  = grad_row - padding;
        node->pt->col  = grad_col - padding;
        node->pt->dcol = grad_next->dcol;
        node->pt->drow = grad_next->drow;
        node->pt->mag  = grad_next->mag;

        dir_x =  grad_next->drow;
        dir_y = -grad_next->dcol;

        node->next = follow_edge_right_f(dir_x, dir_y, grad_row, grad_col, map, 
                num_rows, num_cols, padding, thresh);
    }

    return node;
}


static Edge_list_f* follow_edge_left_f
(
    float        dir_x, 
    float        dir_y, 
    int32_t      grad_row,
    int32_t      grad_col, 
    Gradient_f** map,
    uint32_t     num_rows,
    uint32_t     num_cols,
    uint32_t     padding,
    float        thresh
)
{
    Gradient_f*  grad_curr = NULL;
    Gradient_f*  grad_next = NULL;
    Edge_list_f* node = NULL;
    
    grad_curr = &(map[ grad_row ][ grad_col ]);
    
    follow_gradient_f(&grad_row, &grad_col, dir_x, dir_y, map, 
            num_rows, num_cols, padding);

    grad_next = get_maximal_along_gradient_f(&grad_row, &grad_col, map, 
            num_rows, num_cols, padding);

    if (!grad_next->marked && grad_next->mag >= thresh &&
        grad_row >= padding && grad_row < num_rows-padding &&
        grad_col >= padding && grad_col < num_cols-padding)
    {
        grad_next->marked = 1;

        assert(node = malloc(sizeof(Edge_list_f)));

        assert(node->pt = malloc(sizeof(Edge_point_f)));
        node->pt->row  = grad_row - padding;
        node->pt->col  = grad_col - padding;
        node->pt->dcol = grad_next->dcol;
        node->pt->drow = grad_next->drow;
        node->pt->mag  = grad_next->mag;

        dir_x = -grad_next->drow;
        dir_y =  grad_next->dcol;

        node->next = follow_edge_left_f(dir_x, dir_y, grad_row, grad_col, map, 
                num_rows, num_cols, padding, thresh);
    }

    return node;
}


static Edge_set_f* get_edge_set_from_gradient_map_f
(
    Gradient_f**   map, 
    uint32_t       num_rows, 
    uint32_t       num_cols,
    uint32_t       padding,
    float          begin_thresh, 
    float          end_thresh
)
{
    int32_t  grad_row;
    int32_t  grad_col;
    uint32_t row, col;
    float    dir_x;
    float    dir_y;

    Gradient_f*   grad;
    Gradient_f*   max_grad;
    Edge_list_f*  node;
    Edge_list_f*  right;
    Edge_list_f*  left;
    Edge_list_f*  edge_list;
    Edge_set_f*   edge_set;

    assert(edge_set = malloc(sizeof(Edge_set_f)));
    edge_set->num_edges = 0;
    edge_set->total_num_pts = 0;
    edge_set->num_pts = NULL;
    edge_set->pts     = NULL;

    for (row = padding; row < num_rows-padding; row++)
    {
        for (col = padding; col < num_cols-padding; col++)
        {
            grad = &map[ row ][ col ];

            if (!grad->marked && grad->mag >= begin_thresh)
            {
                grad_row = row;
                grad_col = col;

                max_grad = get_maximal_along_gradient_f(&grad_row, &grad_col, 
                        map, num_rows, num_cols, padding);

                if (!max_grad->marked && 
                    grad_row >= padding && grad_row < num_rows-padding &&
                    grad_col >= padding && grad_col < num_cols-padding)
                {
                    max_grad->marked = 1;

                    assert(node = malloc(sizeof(Edge_list_f)));
                    node->next = NULL;

                    assert(node->pt = malloc(sizeof(Edge_point_f)));
                    node->pt->row   = grad_row - padding;
                    node->pt->col   = grad_col - padding;
                    node->pt->dcol  = grad->dcol;
                    node->pt->drow  = grad->drow;
                    node->pt->mag   = grad->mag;

                    dir_x =  max_grad->drow;
                    dir_y = -max_grad->dcol;
                    left  = follow_edge_right_f(dir_x, dir_y, grad_row,
                            grad_col, map, num_rows, num_cols, padding,
                            end_thresh);

                    dir_x = -max_grad->drow;
                    dir_y =  max_grad->dcol;
                    right = follow_edge_left_f(dir_x, dir_y, grad_row, grad_col,
                            map, num_rows, num_cols, padding, end_thresh);

                    edge_list = merge_edge_lists_f(node, right);
                    edge_list = merge_edge_lists_f(left, edge_list);
                    copy_edge_list_into_set_f(edge_set, edge_list);
                    free_edge_list_f(edge_list);
                }
            }
        }
    }

    return edge_set;
}




Error* detect_image_edge_set_f
(
    Edge_set_f**   edges_out,
    const Image_f* img,
    float          sigma, 
    float          begin_thresh, 
    float          end_thresh,
    uint32_t       padding
)
{
    Matrix_f* m = NULL;
    Error*    e;

    create_matrix_from_image_f(&m, img, 0.3, 0.59, 0.11);
    if ((e = detect_matrix_edge_set_f(edges_out, m, sigma, begin_thresh,
                    end_thresh, padding)) != NULL)
    {
        return e;
    }

    free_matrix_f(m);
    
    return NULL;
}

Error* DEPRECATED_detect_image_edge_set_f
(
    Edge_set_f**   r_edges_out,
    Edge_set_f**   g_edges_out,
    Edge_set_f**   b_edges_out,
    const Image_f* img,
    float          sigma, 
    float          begin_thresh, 
    float          end_thresh,
    uint32_t       padding
)
{
    Matrix_f* m = NULL;
    Error*    e;

    create_matrix_from_image_f(&m, img, 1.0, 0.0, 0.0);
    if ((e = detect_matrix_edge_set_f(r_edges_out, m, sigma, begin_thresh,
                    end_thresh, padding)) != NULL)
    {
        return e;
    }

    create_matrix_from_image_f(&m, img, 0.0, 0.0, 0.0);
    if ((e = detect_matrix_edge_set_f(g_edges_out, m, sigma, begin_thresh,
                    end_thresh, padding)) != NULL)
    {
        return e;
    }

    create_matrix_from_image_f(&m, img, 0.0, 0.0, 1.0);
    if ((e = detect_matrix_edge_set_f(b_edges_out, m, sigma, begin_thresh,
                    end_thresh, padding)) != NULL)
    {
        return e;
    }

    free_matrix_f(m);
    
    return NULL;
}

Error* detect_matrix_edge_set_f
(
    Edge_set_f**    edges_out,
    const Matrix_f* m,
    float           sigma, 
    float           begin_thresh, 
    float           end_thresh,
    uint32_t        padding
)
{
    Matrix_f*    padded_m = NULL;
    Gradient_f** map = NULL;
    Error*       e;

    if (*edges_out != NULL)
    {
        free_edge_set_f(*edges_out);
        *edges_out = NULL;
    }

    extend_matrix_f(&padded_m, m, padding, padding);

    if ((e = create_gradient_map_f(&map, padded_m, sigma)) != NULL)
    {
        free_matrix_f(padded_m);
        return e;
    }

    *edges_out = get_edge_set_from_gradient_map_f(map, padded_m->num_rows,
            padded_m->num_cols, padding, begin_thresh, end_thresh);

    (*edges_out)->num_rows = m->num_rows;
    (*edges_out)->num_cols = m->num_cols;

    free_matrix_f(padded_m);
    free(*map);
    free(map);

    return NULL;
}

void free_edge_set_f(Edge_set_f* edges)
{
    if (edges == NULL)
        return;

    if (edges->total_num_pts > 0)
    {
        free(*(edges->pts));
    }
    free(edges->pts);
    free(edges->num_pts);
    free(edges);
}

void remove_short_edges_f(Edge_set_f* edges, uint32_t min_len)
{
    uint32_t  i, j, k;
    uint32_t  num_pts_to_remove;
    uint32_t  num_edges_to_remove;
    uint32_t* old_num_pts;

    Edge_point_f*  pt;
    Edge_point_f** old_pts;

    num_pts_to_remove = 0;
    num_edges_to_remove = 0;

    for (i = 0; i < edges->num_edges; i++)
    {
        if (edges->num_pts[ i ] < min_len)
        {
            num_pts_to_remove += edges->num_pts[ i ];
            num_edges_to_remove++;
        }
    }

    if (num_edges_to_remove == 0)
    {
        assert(num_pts_to_remove == 0);
        return;
    }

    old_pts = edges->pts;
    old_num_pts = edges->num_pts;

    edges->num_edges -= num_edges_to_remove;
    edges->total_num_pts -= num_pts_to_remove;

    edges->pts = malloc(edges->num_edges * sizeof(Edge_point_f*));
    assert(edges->pts);

    *(edges->pts) = malloc(edges->total_num_pts * sizeof(Edge_point_f));
    assert(*(edges->pts));

    edges->num_pts = malloc(edges->num_edges * sizeof(uint32_t));
    assert(edges->num_pts);

    i = 0;
    pt = *(edges->pts);

    for (j = 0; j < edges->num_edges + num_edges_to_remove; j++)
    {
        if (old_num_pts[ j ] >= min_len)
        {
            edges->num_pts[ i ] = old_num_pts[ j ];

            edges->pts[ i ] = pt;

            for (k = 0; k < edges->num_pts[ i ]; k++)
            {
                *pt = old_pts[ j ][ k ];
                pt++;
            }

            i++;
        }
    }

    free(*old_pts);
    free(old_pts);
    free(old_num_pts);
}

static void recursively_break_edges_at_corners_f
(
    Edge_set_f* edges, 
    uint32_t    i, 
    float       thresh, 
    uint32_t    num_avg
)
{
    uint32_t j, k;
    uint32_t min_j;

    float  dcol_1, drow_1;
    float  dcol_2, drow_2;
    float  dp;
    float  min_dp;

    assert(num_avg > 0);

    if (edges->num_pts[ i ] < 2*num_avg)
        return;

    min_dp = thresh;

    for (j = num_avg-1; j < edges->num_pts[ i ] - num_avg; j++)
    {
        dcol_1 = drow_1 = 0;
        dcol_2 = drow_2 = 0;

        for (k = 0; k < num_avg; k++)
        {
            dcol_1 += edges->pts[ i ][ j-k ].dcol;
            drow_1 += edges->pts[ i ][ j-k ].drow;
            dcol_2 += edges->pts[ i ][ j+k+1 ].dcol;
            drow_2 += edges->pts[ i ][ j+k+1 ].drow;
        }

        dcol_1 /= num_avg;
        drow_1 /= num_avg;
        dcol_2 /= num_avg;
        drow_2 /= num_avg;

        dp = fabs(dcol_1*dcol_2 + drow_1*drow_2);
        if (dp < min_dp)
        {
            min_j = j;
            min_dp = dp;
        }
    }

    if (min_dp < thresh)
    {
        edges->num_edges++;

        edges->num_pts = realloc(edges->num_pts, 
                edges->num_edges * sizeof(uint32_t));
        assert(edges->num_pts);

        edges->pts = realloc(edges->pts, 
                edges->num_edges * sizeof(Edge_point_f*));
        assert(edges->pts);

        for (j = edges->num_edges-1; j > i; j--)
        {
            edges->num_pts[ j ] = edges->num_pts[ j-1 ];
            edges->pts[ j ]     = edges->pts[ j-1 ];
        }

        edges->num_pts[ i ]    = min_j+1;
        edges->num_pts[ i+1 ] -= min_j+1;
        edges->pts[ i+1 ]     += min_j+1;

        recursively_break_edges_at_corners_f(edges, i+1, thresh, num_avg);
        recursively_break_edges_at_corners_f(edges, i, thresh, num_avg);
    }
}

void break_edges_at_corners_f(Edge_set_f* edges, float thresh, uint32_t num_avg)
{
    uint32_t i;

    if (num_avg == 0)
        return;

    for (i = 0; i < edges->num_edges; i++)
    {
        recursively_break_edges_at_corners_f(edges, i, thresh, num_avg);
    }
}

void sample_edge_set_f
(
    Edge_set_f**      edges_out, 
    const Edge_set_f* edges_in, 
    float             p
)
{
    uint32_t      e, pt;
    uint32_t*     num_pts;
    Edge_point_f* pts;
    Edge_set_f*   edges;

    assert(edges_in->total_num_pts > 0);

    if (*edges_out == NULL || *edges_out == edges_in)
    {
        assert(edges = malloc(sizeof(Edge_set_f)));

        edges->pts = malloc(edges_in->num_edges * sizeof(Edge_point_f*));
        assert(edges->pts);

        *(edges->pts) = malloc(edges_in->total_num_pts * sizeof(Edge_point_f));
        assert(*(edges->pts));

        edges->num_pts = malloc(edges_in->num_edges * sizeof(uint32_t));
        assert(edges->num_pts);
    }
    else
    {
        edges = *edges_out;

        edges->pts = realloc(edges->pts, 
                edges_in->num_edges * sizeof(Edge_point_f*));
        if (edges_in->num_edges > 0) assert(edges->pts);

        *(edges->pts) = realloc(*(edges->pts), 
                edges_in->total_num_pts * sizeof(Edge_point_f));
        if (edges_in->total_num_pts > 0) assert(*(edges->pts));

        edges->num_pts = realloc(edges->num_pts, 
                edges_in->num_edges * sizeof(uint32_t));
        if (edges_in->num_edges > 0) assert(edges->num_pts);
    }

    pts = *(edges->pts);
    num_pts = edges->num_pts;

    edges->total_num_pts = 0;
    edges->num_edges = 0;
    edges->num_rows = edges_in->num_rows;
    edges->num_cols = edges_in->num_cols;

    for (e = 0; e < edges_in->num_edges; e++)
    {
        *num_pts = 0;

        for (pt = 0; pt < edges_in->num_pts[ e ]; pt++)
        {
            if (sample_bernoulli_pmf_d(p))
            {
                pts->row  = edges_in->pts[ e ][ pt ].row;
                pts->col  = edges_in->pts[ e ][ pt ].col;
                pts->dcol = edges_in->pts[ e ][ pt ].dcol;
                pts->drow = edges_in->pts[ e ][ pt ].drow;
                pts->mag  = edges_in->pts[ e ][ pt ].mag;

                (*num_pts)++;
                edges->total_num_pts++;

                pts++;
            }
        }

        if (*num_pts != 0)
        {
            num_pts++;
            edges->num_edges++;
        }
    }

    edges->pts = realloc(edges->pts, edges->num_edges * sizeof(Edge_point_f*));
    if (edges->num_edges > 0) assert(edges->pts);

    *(edges->pts) = realloc(*(edges->pts), 
            edges->total_num_pts * sizeof(Edge_point_f));
    if (edges->total_num_pts > 0) assert(*(edges->pts));

    edges->num_pts = realloc(edges->num_pts, 
            edges->num_edges * sizeof(uint32_t));
    if (edges->num_edges > 0) assert(edges->num_pts);


    pts = *(edges->pts);
    for (e = 1; e < edges->num_edges; e++)
    {
        pts += edges->num_pts[ e-1 ];
        edges->pts[ e ] = pts;
    }

    if (*edges_out == edges_in)
    {
        free_edge_set_f(*edges_out);
    }

    *edges_out = edges;
}

Error* color_edge_set_f
(
    Image_f**         img_out,
    const Image_f*    img_in,
    const Edge_set_f* edges, 
    const Pixel_f*    pxl
)
{
    uint32_t e;
    Error* err;

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

    for (e = 0; e < edges->num_edges; e++)
    {
        if ((err = color_edge_points_f(img_out, *img_out, edges->pts[ e ],
                edges->num_pts[ e ], pxl)))
        {
            return err;
        }
    }

    return NULL;
}

Error* randomly_color_edge_set_f
(
    Image_f**         img_out,
    const Image_f*    img_in,
    const Edge_set_f* edges
)
{
    uint32_t e;
    Error* err;

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

    for (e = 0; e < edges->num_edges; e++)
    {
        if ((err = randomly_color_edge_points_f(img_out, *img_out, 
                        edges->pts[ e ], edges->num_pts[ e ])))
        {
            return err;
        }
    }

    return NULL;
}

Error* color_edge_points_f
(
    Image_f**           img_out,
    const Image_f*      img_in,
    const Edge_point_f* pts, 
    uint32_t            num_pts,
    const Pixel_f*      pxl
)
{
    uint32_t pt;
    uint32_t row, col;
    Image_f* img;

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

    for (pt = 0; pt < num_pts; pt++)
    {
        row = pts[ pt ].row;
        col = pts[ pt ].col;

        if (row >= img_in->num_rows || col >= img_in->num_cols)
        {
            if (*img_out != img_in)
            {
                free_image_f(*img_out);
                *img_out = 0;
            }
            return JWSC_EARG("Edge points outside image to color in");
        }

        img->pxls[ row ][ col ] = *pxl;
    }

    return NULL;
}

Error* randomly_color_edge_points_f
(
    Image_f**           img_out,
    const Image_f*      img_in,
    const Edge_point_f* pts, 
    uint32_t            num_pts
)
{
    Pixel_f pxl;

    pxl.r = (float)rand() / (float)RAND_MAX;
    pxl.g = (float)rand() / (float)RAND_MAX;
    pxl.b = (float)rand() / (float)RAND_MAX;

    return color_edge_points_f(img_out, img_in, pts, num_pts, &pxl);
}

Error* read_edge_set_f
(
    Edge_set_f** edges_out,
    const char*  fname
)
{
    uint32_t i, j;
    FILE*    fp = NULL;
    Error*   e = NULL;

    Edge_set_f*  edges;
    Edge_point_f* pt;

    if (*edges_out != NULL)
    {
        free_edge_set_f(*edges_out);
        *edges_out = NULL;
    }

    if ((fp = fopen(fname, "r")) == NULL)
    {
        e = JWSC_EIO(strerror(errno));
        goto error;
    }

    assert(*edges_out = malloc(sizeof(Edge_set_f)));
    edges = *edges_out;

    if (fscanf(fp, "num_rows=%u\n", &(edges->num_rows)) != 1 ||
        fscanf(fp, "num_cols=%u\n", &(edges->num_cols)) != 1 ||
        fscanf(fp, "num_edges=%u\n", &(edges->num_edges)) != 1 ||
        fscanf(fp, "total_num_pts=%u\n", &(edges->total_num_pts)) != 1)
    {
        e = JWSC_EARG("File not properly formatted to read edge points");
        goto error;
    }

    assert(edges->pts = malloc(edges->num_edges * sizeof(Edge_point_f*)));
    assert(*(edges->pts) = malloc(edges->total_num_pts * sizeof(Edge_point_f)));
    assert(edges->num_pts = malloc(edges->num_edges * sizeof(uint32_t)));

    pt = *(edges->pts);

    for (i = 0; i < edges->num_edges; i++)
    {
        if (fscanf(fp, "num_pts=%u\n", &(edges->num_pts[ i ])) != 1)
        {
            e = JWSC_EARG("File not properly formatted to read edge points");
            goto error;
        }

        edges->pts[ i ] = pt;

        for (j = 0; j < edges->num_pts[ i ]; j++)
        {
            if (fscanf(fp, "col=%u row=%u dcol=%f drow=%f mag=%f\n", 
                        &(pt->col), &(pt->row), &(pt->dcol), &(pt->drow),
                        &(pt->mag)) != 5)
            {
                e = JWSC_EARG("Edge point file improperly formatted");
                goto error;
            }

            pt++;
        }
    }

    if (fclose(fp) != 0)
    {
        e = JWSC_EIO(strerror(errno));
        goto error;
    }

    return NULL;

error:
    if (fp) fclose(fp);
    free_edge_set_f(*edges_out); 
    *edges_out = NULL;
    return e;
}

Error* write_edge_set_f
(
    const Edge_set_f* edges,
    const char*       fname
)
{
    uint32_t i, j;
    FILE*    fp;

    const Edge_point_f* pt;

    if ((fp = fopen(fname, "w")) == NULL)
    {
        return JWSC_EIO(strerror(errno));
    }

    fprintf(fp, "num_rows=%u\n", edges->num_rows);
    fprintf(fp, "num_cols=%u\n", edges->num_cols);
    fprintf(fp, "num_edges=%u\n", edges->num_edges);
    fprintf(fp, "total_num_pts=%u\n", edges->total_num_pts);

    for (i = 0; i < edges->num_edges; i++)
    {
        fprintf(fp, "num_pts=%u\n", edges->num_pts[ i ]);

        for (j = 0; j < edges->num_pts[ i ]; j++)
        {
            pt = &(edges->pts[ i ][ j ]);

            fprintf(fp, "col=%u row=%u", pt->col, pt->row);
            fprintf(fp, " dcol=%f drow=%f", pt->dcol, pt->drow);
            fprintf(fp, " mag=%f\n", pt->mag);
        }
    }

    if (fclose(fp) != 0)
    {
        return JWSC_EIO(strerror(errno));
    }

    return NULL;
}