Initial commit

Makefile

@@ -0,0 +1,16 @@
+PKGS=ufo
+CFLAGS+=$(shell pkg-config --cflags $(PKGS))
+LDFLAGS+=$(shell pkg-config --libs $(PKGS))
+
+SRC=lamino.c
+BIN=lamino
+
+all: $(BIN)
+
+.PHONY: clean
+
+$(BIN): $(SRC)
+	@$(CC) $(CFLAGS) -o $@ $< $(LDFLAGS)
+
+clean:
+	rm -f $(BIN)

README.md

@@ -0,0 +1,27 @@
+# lamino-c
+
+Laminographic reconstruction with the UFO framework and Anton Myagotin's
+lamino-filters.
+
+## Parameters
+
+Parameter     | Meaning
+------------- | ---------------------------
+`width`       | Projecion width in pixels
+`height`      | Projecion height in pixels
+`num-radios`  | Number of projections
+`radios`      | Location of radios
+`darks`       | Location of darks
+`flats`       | Location of flats
+`output`      | Output filename
+`theta`       | Laminographic tilt angle
+`px`          | X coordinate of projection axis
+`py`          | Y coordinate of projection axis
+`vx`          | X coordinate of box origin
+`vy`          | Y coordinate of box origin
+`vz`          | Z coordinate of box origin
+`vw`          | Width of box
+`vh`          | Height of box
+`vd`          | Depth of box
+`tau`         | Pixel size
+`psi`         | Misalignment

lamino.c

@@ -0,0 +1,267 @@
+#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;
+}