surface.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 <config.h>

#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <inttypes.h>
#include <string.h>

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

#include <jwsc/base/error.h>
#include <jwsc/base/option.h>
#include <jwsc/matblock/matblock.h>
#include <jwsc/imgblock/imgblock.h>
#include <jwsc/imgblock/imgblock_io.h>
#include <jwsc/imgblock/imgblock_util.h>
#include <jwsc/imgblock/surface.h>


#define DEFAULT_BEGIN_THRESH    0.06f
#define DEFAULT_END_THRESH      0.025f
#define DEFAULT_SIGMA           1.0f
#define DEFAULT_SAMPLE_RATE     1.0f
#define DEFAULT_GRAD_IMG_WEIGHT 1.0f
#define DEFAULT_GRAD_ROW_WEIGHT 1.0f
#define DEFAULT_GRAD_COL_WEIGHT 1.0f
#define DEFAULT_RES_IMG         1.0f
#define DEFAULT_RES_ROW         1.0f
#define DEFAULT_RES_COL         1.0f

#define NUM_OPTS_NO_ARG   3
#define NUM_OPTS_WITH_ARG 9


static Option_no_arg opts_no_arg[NUM_OPTS_NO_ARG];

static Option_with_arg opts_with_arg[NUM_OPTS_WITH_ARG];


static const char* pts_fname;

static const char* map_fname;

static const char* colored_fname;

static float begin_thresh;

static float end_thresh;

static float sigma;

static float sample_rate;

static float grad_img_weight;

static float grad_row_weight;

static float grad_col_weight;

static uint8_t images_named;

static float res_img;

static float res_row;

static float res_col;


static void print_usage(void)
{
    fprintf(stderr, "usage: surface OPTIONS [img-%%d num_imgs | img-1 ... img-n]\n");
    fprintf(stderr, "where OPTIONS is one or more of the following:\n");
    print_options(stderr, 25, NUM_OPTS_NO_ARG, opts_no_arg, NUM_OPTS_WITH_ARG,
            opts_with_arg);
}


static Error* process_help_opt(void)
{
    print_usage();
    exit(EXIT_SUCCESS);
    return NULL;
}


static Error* process_version_opt(void)
{
    fprintf(stderr, "%s\n", SURFACE_PACKAGE_STRING);
    exit(EXIT_SUCCESS);
    return NULL;
}


static Error* process_named_opt(void)
{
    images_named = 1;
    return NULL;
}


static Error* process_colored_out_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'colored-out' requires an argument");
    }
    colored_fname = arg;
    return NULL;
}


static Error* process_map_out_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'map-out' requires an argument");
    }
    map_fname = arg;
    return NULL;
}


static Error* process_pts_out_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'pts-out' requires an argument");
    }
    pts_fname = arg;
    return NULL;
}


static Error* process_begin_thresh_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'begin-thresh' requires an argument");
    }
    if (sscanf(arg, "%f", &begin_thresh) != 1 || begin_thresh < 0)
    {
        return JWSC_EARG("Option 'begin-thresh' must be >= 0");
    }
    return NULL;
}


static Error* process_end_thresh_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'end-thresh' requires an argument");
    }
    if (sscanf(arg, "%f", &end_thresh) != 1 || end_thresh < 0)
    {
        return JWSC_EARG("Option 'end-thresh' must be >= 0");
    }
    return NULL;
}


static Error* process_sigma_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'sigma' requires an argument");
    }
    if (sscanf(arg, "%f", &sigma) != 1 || sigma <= 0)
    {
        return JWSC_EARG("Option 'sigma' must be > 0");
    }
    return NULL;
}


static Error* process_sample_rate_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'sample-rate' requires an argument");
    }
    if (sscanf(arg, "%f", &sample_rate) != 1 || sample_rate <= 0 ||
            sample_rate > 1.0)
    {
        return JWSC_EARG("Option 'sample-rate' must be in (0, 1]");
    }
    return NULL;
}


static Error* process_grad_weight_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'grad-weight' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f", &grad_img_weight, &grad_row_weight, 
                &grad_col_weight) != 3)
    {
        return JWSC_EARG("Option 'grad-weight' has format img,row,col");
    }
    if (grad_img_weight < 0 || grad_row_weight < 0 || grad_col_weight < 0)
    {
        return JWSC_EARG("Option 'grad-weight' must have img,row,col >= 0");
    }
    if ((grad_img_weight + grad_row_weight + grad_col_weight) <= 1.0e-8)
    {
        return JWSC_EARG("Option 'grad-weight' must have img+row+col > 0");
    }
    return NULL;
}


static Error* process_resolution_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'resolution' requires an argument");
    }
    if (sscanf(arg, "%f,%f,%f", &res_img, &res_row, &res_col) != 3)
    {
        return JWSC_EARG("Option 'resolution' has format img,row,col");
    }
    if (res_img <= 0 || res_img > 1 || res_row <= 0 || res_row > 1 ||
            res_col <= 0 || res_col > 1)
    {
        return JWSC_EARG("Option 'resolution' must be in (0,1]");
    }
    return NULL;
}





int main(int argc, const char** argv)
{
    uint32_t num_rows, num_cols;
    uint32_t num_pts;
    uint32_t i, p;
    int32_t  num_imgs;
    int      argi;
    int      arg;

    char*       chars;
    char**      names;
    const char* input_fname;

    Pixel_f pxl;

    Imgblock_f*      img_in  = NULL;
    Imgblock_f*      img_out = NULL;
    Matblock_f*      m       = NULL;
    Error*           e       = NULL;
    Surface_point_f* pts     = NULL;
    Surface_point_f* pt      = NULL;

    /* Init the options. */
    init_option_no_arg(&(opts_no_arg[0]), "help", 'h',
            "Prints program usage.", process_help_opt);
    init_option_no_arg(&(opts_no_arg[1]), "version", 'v',
            "Prints program version.", process_version_opt);
    init_option_no_arg(&(opts_no_arg[2]), "named", 'n',
            "Indicates that the images to read are fully named and listed in sequence on the command line. The default is to read the images as numbered by a printf-formatted name.", process_named_opt);
    init_option_with_arg(&(opts_with_arg[0]), "colored-out", 'c',
            "Original images with surface colored. Use a printf formatted name to output the images in a numbered sequence.", process_colored_out_opt);
    init_option_with_arg(&(opts_with_arg[1]), "map-out", 'm',
            "Bi-level image surface map. Use a printf formatted name to output the images in a numbered sequence.", process_map_out_opt);
    init_option_with_arg(&(opts_with_arg[2]), "pts-out", 'p',
            "ASCII formatted surface points.", process_pts_out_opt);
    init_option_with_arg(&(opts_with_arg[3]), "begin-thresh", 'b',
            "Hysterisis starting threshold.", process_begin_thresh_opt);
    init_option_with_arg(&(opts_with_arg[4]), "end-thresh", 'e',
            "Hysterisis stopping threshold.", process_end_thresh_opt);
    init_option_with_arg(&(opts_with_arg[5]), "sigma", 's',
            "Gaussian smoothing sigma.", process_sigma_opt);
    init_option_with_arg(&(opts_with_arg[6]), "sample-rate", 'r',
            "Surface point sampling rate.", process_sample_rate_opt);
    init_option_with_arg(&(opts_with_arg[7]), "grad-weight", 'w',
            "Weights for each direction of the gradient. Format for the argument is 'img,row,col'. They all must be non-negative and sum to > 0.", 
            process_grad_weight_opt);
    init_option_with_arg(&(opts_with_arg[8]), "resolution", 'w',
            "Down-sampling resolution for the input image stack. Format for the argument is 'img,row,col'. They all must be in (0,1].",
            process_resolution_opt);

    /* Init default option values. */
    map_fname       = NULL;
    colored_fname   = NULL;
    pts_fname       = NULL;
    begin_thresh    = DEFAULT_BEGIN_THRESH;
    end_thresh      = DEFAULT_END_THRESH;
    sigma           = DEFAULT_SIGMA;
    sample_rate     = DEFAULT_SAMPLE_RATE;
    grad_img_weight = DEFAULT_GRAD_IMG_WEIGHT;
    grad_row_weight = DEFAULT_GRAD_ROW_WEIGHT;
    grad_col_weight = DEFAULT_GRAD_COL_WEIGHT;
    res_img         = DEFAULT_RES_IMG;
    res_row         = DEFAULT_RES_ROW;
    res_col         = DEFAULT_RES_COL;
    images_named    = 0;

    /* Check for the help option first. */
    process_option_no_arg(argc, argv, &argi, &(opts_no_arg[0]));

    if ((e = process_options(argc, argv, &argi, NUM_OPTS_NO_ARG, opts_no_arg,
                    NUM_OPTS_WITH_ARG, opts_with_arg)) != NULL)
    {
        print_error_msg_exit("surface", e->msg);
    }

    /* Either read the image block as a printf numbered sequence, or a list of
     * names in argv. */
    if (images_named)
    {
        if ((argc - argi) < 1)
        {
            print_error_msg_exit("surface", "No input images");
        }

        num_imgs = argc - argi;
        chars = calloc(num_imgs*256, sizeof(char));
        names = malloc(num_imgs*sizeof(char*));
        arg = argi;
        for (i = 0; i < num_imgs; i++)
        {
            names[ i ] = &(chars[ i*256 ]);
            strncpy(names[ i ], argv[ arg++ ], 255);
        }
        if ((e = read_imgblock_named_f(&img_in, (const char**)names, num_imgs)) 
                != NULL)
        {
            print_error_msg_exit("surface", e->msg);
        }
        free(chars);
        free(names);
    }
    else
    {
        if ((argc - argi) < 2)
        {
            print_error_msg_exit("surface", "No input image format and count");
        }

        input_fname = argv[ argi ];
        if (sscanf(argv[ argi+1 ], "%d", &num_imgs) != 1 || num_imgs <= 0)
        {
            print_error_msg_exit("surface", "No input image count");
        }

        if ((e = read_imgblock_numbered_f(&img_in, input_fname, num_imgs)) 
                != NULL)
        {
            print_error_msg_exit("surface", e->msg);
        }
    }


    /* Down sample the image data if requested. */
    assert(down_sample_imgblock_f(&img_in, img_in, res_img, res_row, res_col) 
            == NULL);
    num_imgs = img_in->num_imgs;
    num_rows = img_in->num_rows;
    num_cols = img_in->num_cols;


    /* Detect surface and sample them if requested. */
    create_matblock_from_imgblock_f(&m, img_in, 1.0f, 1.0f, 1.0f);
    if ((e = detect_matblock_surface_points_f(&pts, &num_pts, m, sigma, 
                    grad_img_weight, grad_row_weight, grad_col_weight, 
                    begin_thresh, end_thresh)) != NULL)
    {
        print_error_msg_exit("surface", e->msg);
        return EXIT_FAILURE;
    }

    fprintf(stdout, "Number of detected surface points: %d\n", num_pts);
    sample_surface_points_f(&pts, &num_pts, pts, num_pts, sample_rate);
    fprintf(stdout, "Number of sampled surface points: %d\n", num_pts);

    /* Write a surface points if requested. */
    if (pts_fname)
    {
        if ((e = write_surface_points_f(pts, num_pts, pts_fname)) != NULL)
        {
            print_error_msg("surface", NULL);
            print_error_msg(pts_fname, e->msg);
            free_error(e);
        }
    }

    /* Write a surface map image block if requested. */
    if (map_fname)
    {
        create_zero_matblock_f(&m, num_imgs, num_rows, num_cols);

        for (p = 0; p < num_pts; p++)
        {
            pt = &(pts[ p ]);
            m->elts[ pt->mat ][ pt->row ][ pt->col ] = 1.0f;
        }

        create_imgblock_from_matblock_f(&img_out, m);

        if ((e = write_imgblock_numbered_f(img_out, map_fname)) != NULL)
        {
            print_error_msg("surface", e->msg);
            free_error(e);
        }
    }

    /* Write the image block with the surface colored on it if requested. */
    if (colored_fname)
    {
        pxl.r = 1.0f;
        pxl.g = 0;
        pxl.b = 0;

        color_surface_points_f(&img_out, img_in, pts, num_pts, &pxl);

        if ((e = write_imgblock_numbered_f(img_out, colored_fname)) != NULL)
        {
            print_error_msg("surface", e->msg);
            free_error(e);
        }
    }

    free(pts);
    free_imgblock_f(img_in);
    free_imgblock_f(img_out);
    free_matblock_f(m);

    return EXIT_SUCCESS;
}