lamino-reco/lamino.c

#include <stdlib.h>
#include <ufo-0/ufo/ufo.h>

typedef struct {
    gchar *radios;
    gchar *darks;
    gchar *flats;
    gchar *output;
    gint width;
    gint height;
    guint num_radios;
    guint num_darks;
    gdouble theta;
    gdouble tau;
    gdouble psi;
    gdouble px;
    gdouble py;
    gdouble v_origin[3];
    guint v_size[3];
} Params;

static void
check (GError *error)
{
    if (error != NULL) {
        if (error->message != NULL)
            g_printerr ("%s", error->message);

        exit (1);
    }
}

static UfoTaskNode *
make_task (UfoPluginManager *pm, const gchar *name)
{
    GError *error = NULL;
    UfoTaskNode *task;

    task = ufo_plugin_manager_get_task (pm, name, &error);
    check (error);

    return task;
}

static guint
next_power_of_two (guint32 x)
{
    --x;
    x |= x >> 1;
    x |= x >> 2;
    x |= x >> 4;
    x |= x >> 8;
    x |= x >> 16;
    return x + 1;
}

static void
run_reconstruction (Params *params)
{
    UfoPluginManager *pm;
    UfoTaskGraph *graph;
    UfoBaseScheduler *sched;
    UfoTaskNode *radio_reader;
    UfoTaskNode *dark_reader;
    UfoTaskNode *flat_reader;
    UfoTaskNode *ffc;
    UfoTaskNode *pad;
    UfoTaskNode *ramp;
    UfoTaskNode *fft1;
    UfoTaskNode *fft2;
    UfoTaskNode *ifft;
    UfoTaskNode *conv;
    UfoTaskNode *reco;
    UfoTaskNode *writer;
    guint padded_width;
    guint padded_height;
    guint xl, xr, yt, yb;
    GError *error = NULL;

    pm = ufo_plugin_manager_new (NULL);
    graph = UFO_TASK_GRAPH (ufo_task_graph_new ());

    radio_reader = make_task (pm, "reader");
    flat_reader = make_task (pm, "reader");
    dark_reader = make_task (pm, "reader");
    ffc = make_task (pm, "flat-field-correction");
    pad = make_task (pm, "padding-2d");
    ramp = make_task (pm, "lamino-ramp");
    fft1 = make_task (pm, "fft");
    fft2 = make_task (pm, "fft");
    ifft = make_task (pm, "ifft");
    conv = make_task (pm, "lamino-conv");
    reco = make_task (pm, "lamino-bp");
    writer = make_task (pm, "writer");

    g_object_set (radio_reader,
                  "path", params->radios,
                  "end", params->num_radios - 1,
                  NULL);

    g_object_set (flat_reader, "path", params->flats, NULL);
    g_object_set (dark_reader, "path", params->darks, NULL);

    g_object_set (ffc, "num-darks", params->num_darks, NULL);
    g_object_set (fft1, "dimensions", 2, NULL);
    g_object_set (fft2, "dimensions", 2, NULL);
    g_object_set (ifft, "dimensions", 2, NULL);

    padded_width = next_power_of_two ((guint32) params->width);
    padded_height = next_power_of_two ((guint32) params->height);

    xl = (padded_width - params->width) / 2;
    xr = xl;
    yt = (padded_height - params->height) / 2;
    yb = yt;

    g_object_set (pad,
                  "xl", xl,
                  "xr", xr,
                  "yt", yt,
                  "yb", yb,
                  "mode", "brep",
                  NULL);

    g_object_set (ramp,
                  "width", padded_width,
                  "height", padded_height,
                  "theta", G_PI * 2 / params->theta,
                  "tau", params->tau,
                  NULL);

    g_object_set (reco,
                  "angle-step", (G_PI * 2.0) / params->num_radios,
                  "theta", G_PI * 2 / params->theta,
                  "psi", params->psi,
                  "proj-ox", params->px + xl,
                  "proj-oy", params->py,
                  "vol-ox", params->v_origin[0],
                  "vol-oy", params->v_origin[1],
                  "vol-oz", params->v_origin[2],
                  "vol-sx", params->v_size[0],
                  "vol-sy", params->v_size[1],
                  "vol-sz", params->v_size[2],
                  NULL);

    g_object_set (writer,
                  "filename", "/data/visitor/ma2183/id19/volume-%i.tif",
                  NULL);

    ufo_task_graph_connect_nodes_full (graph, radio_reader, ffc, 0);
    ufo_task_graph_connect_nodes_full (graph, dark_reader, ffc, 1);
    ufo_task_graph_connect_nodes_full (graph, flat_reader, ffc, 2);
    ufo_task_graph_connect_nodes (graph, ffc, pad);
    ufo_task_graph_connect_nodes (graph, pad, fft1);
    ufo_task_graph_connect_nodes (graph, ramp, fft2);
    ufo_task_graph_connect_nodes_full (graph, fft1, conv, 0);
    ufo_task_graph_connect_nodes_full (graph, fft2, conv, 1);
    ufo_task_graph_connect_nodes (graph, conv, ifft);
    ufo_task_graph_connect_nodes (graph, ifft, reco);
    ufo_task_graph_connect_nodes (graph, reco, writer);

    sched = UFO_BASE_SCHEDULER (ufo_scheduler_new (NULL, NULL));
    ufo_base_scheduler_run (sched, graph, &error);
    check (error);

    g_object_unref (pm);
    g_object_unref (graph);
    g_object_unref (sched);

    g_object_unref (radio_reader);
    g_object_unref (flat_reader);
    g_object_unref (dark_reader);
    g_object_unref (ffc);
    g_object_unref (pad);
    g_object_unref (fft1);
    g_object_unref (fft2);
    g_object_unref (conv);
    g_object_unref (reco);
    g_object_unref (writer);
}

static void
parse_params (Params *params, int argc, char **argv)
{
    GOptionContext *context;
    GError *error = NULL;

    GOptionEntry entries[] = {
        { "width", 0, 0, G_OPTION_ARG_INT, &params->width, "Width", "" },
        { "height", 0, 0, G_OPTION_ARG_INT, &params->height, "Height", "" },
        { "num-radios", 0, 0, G_OPTION_ARG_INT, &params->num_radios, "Number of radios", "" },
        { "num-darks", 0, 0, G_OPTION_ARG_INT, &params->num_darks, "Number of darks", "" },
        { "radios", 0, 0, G_OPTION_ARG_STRING, &params->radios, "Radios", "" },
        { "darks", 0, 0, G_OPTION_ARG_STRING, &params->darks, "Darks", "" },
        { "flats", 0, 0, G_OPTION_ARG_STRING, &params->flats, "Flats", "" },
        { "output", 0, 0, G_OPTION_ARG_STRING, &params->output, "Output", "" },
        { "theta", 0, 0, G_OPTION_ARG_DOUBLE, &params->theta, "Tilt (theta)", "" },
        { "tau", 0, 0, G_OPTION_ARG_DOUBLE, &params->theta, "Pixel size (theta)", "" },
        { "psi", 0, 0, G_OPTION_ARG_DOUBLE, &params->theta, "Misalignment (psi)", "" },
        { "px", 0, 0, G_OPTION_ARG_DOUBLE, &params->px, "X coordinate of axis", "" },
        { "py", 0, 0, G_OPTION_ARG_DOUBLE, &params->py, "Y coordinate of axis", "" },
        { "vx", 0, 0, G_OPTION_ARG_DOUBLE, &params->v_origin[0], "X coordinate of box origin", "" },
        { "vy", 0, 0, G_OPTION_ARG_DOUBLE, &params->v_origin[1], "Y coordinate of box origin", "" },
        { "vz", 0, 0, G_OPTION_ARG_DOUBLE, &params->v_origin[2], "Z coordinate of box origin", "" },
        { "vw", 0, 0, G_OPTION_ARG_INT, &params->v_size[0], "Width of box", "" },
        { "vh", 0, 0, G_OPTION_ARG_INT, &params->v_size[1], "Height of box", "" },
        { "vd", 0, 0, G_OPTION_ARG_INT, &params->v_size[2], "Depth of box", "" },
        { NULL }
    };

    context = g_option_context_new ("- test tree model performance");
    g_option_context_add_main_entries (context, entries, NULL);

    if (!g_option_context_parse (context, &argc, &argv, &error)) {
        g_print ("option parsing failed: %s\n", error->message);
        exit (1);
    }

    if (params->width == 0 || params->height == 0 || params->num_radios == 0 ||
        !params->radios || !params->darks || !params->flats ||
        params->theta == G_MAXDOUBLE ||
        params->px == G_MAXDOUBLE || params->py == G_MAXDOUBLE) {
            g_printerr ("Error: Parameters missing.\n\n%s",
                        g_option_context_get_help (context, TRUE, NULL));

            exit (1);
    }

    g_option_context_free (context);
}

static void
check_params (Params *params)
{
}

int
main (int argc, char **argv)
{
    Params params = {
        .radios = NULL,
        .darks = NULL,
        .flats = NULL,
        .output = "volume-%i.tif",
        .width = 0,
        .height = 0,
        .num_radios = 0,
        .num_darks = 1,
        .theta = G_MAXDOUBLE,
        .tau = 1.0,
        .psi = 0.0,
        .px = G_MAXDOUBLE,
        .py = G_MAXDOUBLE,
        .v_size = {256, 256, 256},
        .v_origin = {0.0, 0.0, 0.0},
    };

#if !(GLIB_CHECK_VERSION (2, 36, 0))
    g_type_init ();
#endif

    parse_params (&params, argc, argv);
    check_params (&params);
    run_reconstruction (&params);

    return 0;
}