Implement cached reconstruction

lamino.c

@@ -14,8 +14,6 @@ typedef struct {
 } CommandMap;
 
 
-
-
 static void
 print (Params *params, gchar **argv)
 {
@@ -97,7 +95,7 @@ static void
 start_shell (Params *params)
 {
     static CommandMap map[] = {
-        { "run", run_simple_reconstruction },
+        { "run", run_cached_reconstruction },
         { "print", print },
         { "set", set },
         { "quit", quit },

reco.c

@@ -1,8 +1,28 @@
 #include <stdlib.h>
+#include <string.h>
 #include <ufo/ufo.h>
 #include "reco.h"
 #include "io.h"
 
+typedef struct {
+    UfoPluginManager *pm;
+    UfoTaskGraph *graph;
+    UfoTaskNode *pad;
+    UfoTaskNode *ramp;
+    UfoTaskNode *fft1;
+    UfoTaskNode *fft2;
+    UfoTaskNode *ifft;
+    UfoTaskNode *conv;
+    UfoTaskNode *reco;
+    UfoTaskNode *writer;
+
+    /* interactive */
+    guint current;
+    UfoTaskNode *output;
+    UfoTaskNode *input;
+    Params *params;
+} LaminoData;
+
 
 static void
 print_dots (gpointer user)
@@ -60,59 +80,37 @@ params_okay (Params *params)
            params->py < G_MAXDOUBLE;
 }
 
-void
-run_simple_reconstruction (Params *params, gchar **argv)
+gboolean
+with_flat_field_correction (Params *params)
+{
+    return params->darks != NULL && params->flats != NULL;
+}
+
+static LaminoData *
+reco_graph_new (Params *params)
 {
+    LaminoData *data;
     guint xl, xr, yt, yb;
     guint padded_width;
     guint padded_height;
     gdouble theta_rad;
     gdouble angle_step;
     guint fwidth;
-    gdouble time = 0.0;
-    UfoPluginManager *pm = NULL;
-    UfoTaskGraph *graph = NULL;
-    UfoBaseScheduler *sched = NULL;
-    UfoTaskNode *radio_reader = NULL;
-    UfoTaskNode *dark_reader = NULL;
-    UfoTaskNode *flat_reader = NULL;
-    UfoTaskNode *ffc = NULL;
-    UfoTaskNode *pad = NULL;
-    UfoTaskNode *ramp = NULL;
-    UfoTaskNode *fft1 = NULL;
-    UfoTaskNode *fft2 = NULL;
-    UfoTaskNode *ifft = NULL;
-    UfoTaskNode *conv = NULL;
-    UfoTaskNode *reco = NULL;
-    UfoTaskNode *writer = NULL;
-    GError *error = NULL;
 
-    if (!params_okay (params)) {
-        err ("Parameters missing.\n");
-        return;
-    }
+    data = g_malloc0 (sizeof (LaminoData));
 
-    pm = ufo_plugin_manager_new (NULL);
-    graph = UFO_TASK_GRAPH (ufo_task_graph_new ());
+    data->params = params;
+    data->pm = ufo_plugin_manager_new (NULL);
+    data->graph = UFO_TASK_GRAPH (ufo_task_graph_new ());
 
-    radio_reader = make_task (pm, "reader");
-    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 (fft1, "dimensions", 2, NULL);
-    g_object_set (fft2, "dimensions", 2, NULL);
-    g_object_set (ifft, "dimensions", 2, NULL);
+    data->pad = make_task (data->pm, "padding-2d");
+    data->ramp = make_task (data->pm, "lamino-ramp");
+    data->fft1 = make_task (data->pm, "fft");
+    data->fft2 = make_task (data->pm, "fft");
+    data->ifft = make_task (data->pm, "ifft");
+    data->conv = make_task (data->pm, "lamino-conv");
+    data->reco = make_task (data->pm, "lamino-bp");
+    data->writer = make_task (data->pm, "writer");
 
     fwidth = ((guint) params->px) * 2;
     padded_width = next_power_of_two (fwidth) * 2;
@@ -126,26 +124,30 @@ run_simple_reconstruction (Params *params, gchar **argv)
     angle_step = (G_PI * 2.0) / params->num_radios;
     theta_rad = params->theta / 360. * G_PI * 2;
 
-    info ("Axis: x=%.1f  y=%.1f  variation=%.1f\n",
+    info ("Axis     x=%.1f  y=%.1f  variation=%.1f\n",
           params->px, params->py, params->px_variation);
 
-    info ("Lamino: theta=%.3f  tau=%.3f  step=%.5f  fwidth=%d\n",
+    info ("Lamino   theta=%.3f  tau=%.3f  step=%.5f  fwidth=%d\n",
           theta_rad, params->tau, angle_step, fwidth);
 
-    info ("Padding: size=[%d %d]  (xl=%d xr=%d yt=%d yb=%d)\n",
+    info ("Padding  size=[%d %d]  (xl=%d xr=%d yt=%d yb=%d)\n",
           padded_width, padded_height, xl, xr, yt, yb);
 
-    info ("Volume: origin=[%.1f %.1f %.1f]  size=[%d %d %d]\n",
+    info ("Volume   origin=[%.1f %.1f %.1f]  size=[%d %d %d]\n",
           params->v_origin[0], params->v_origin[1], params->v_origin[2],
           params->v_size[0], params->v_size[1], params->v_size[2]);
 
-    g_object_set (pad,
+    g_object_set (data->fft1, "dimensions", 2, NULL);
+    g_object_set (data->fft2, "dimensions", 2, NULL);
+    g_object_set (data->ifft, "dimensions", 2, NULL);
+
+    g_object_set (data->pad,
                   "xl", xl, "xr", xr,
                   "yt", yt, "yb", yb,
                   "mode", "brep",
                   NULL);
 
-    g_object_set (ramp,
+    g_object_set (data->ramp,
                   "width", padded_width,
                   "height", padded_height,
                   "theta", theta_rad,
@@ -153,7 +155,7 @@ run_simple_reconstruction (Params *params, gchar **argv)
                   "fwidth", fwidth,
                   NULL);
 
-    g_object_set (reco,
+    g_object_set (data->reco,
                   "angle-step", angle_step,
                   "theta", theta_rad,
                   "psi", params->psi,
@@ -168,62 +170,228 @@ run_simple_reconstruction (Params *params, gchar **argv)
                   "vol-sz", params->v_size[2],
                   NULL);
 
-    g_object_set (writer,
+    g_object_set (data->writer,
                   "filename", params->output,
                   NULL);
 
-    if (params->darks != NULL && params->flats != NULL) {
-        flat_reader = make_task (pm, "reader");
-        dark_reader = make_task (pm, "reader");
-        ffc = make_task (pm, "flat-field-correction");
+    ufo_task_graph_connect_nodes (data->graph, data->pad, data->fft1);
+    ufo_task_graph_connect_nodes (data->graph, data->ramp, data->fft2);
+    ufo_task_graph_connect_nodes_full (data->graph, data->fft1, data->conv, 0);
+    ufo_task_graph_connect_nodes_full (data->graph, data->fft2, data->conv, 1);
+    ufo_task_graph_connect_nodes (data->graph, data->conv, data->ifft);
+    ufo_task_graph_connect_nodes (data->graph, data->ifft, data->reco);
+    ufo_task_graph_connect_nodes (data->graph, data->reco, data->writer);
+
+    return data;
+}
+
+static void
+reco_graph_free (LaminoData *data)
+{
+    g_object_unref (data->pm);
+    g_object_unref (data->graph);
+
+    g_object_unref (data->pad);
+    g_object_unref (data->fft1);
+    g_object_unref (data->fft2);
+    /* g_object_unref (conv); */
+    g_object_unref (data->reco);
+    g_object_unref (data->writer);
+}
+
+static void
+update_reader (UfoTaskNode *reader, Params *params)
+{
+    g_object_set (reader,
+                  "path", params->radios,
+                  "end", params->num_radios - 1,
+                  NULL);
+}
+
+static void
+show_time (UfoBaseScheduler *sched)
+{
+    gdouble time = 0.0;
+
+    g_object_get (sched, "time", &time, NULL);
+    g_print ("\n");
+    info ("Finished in %3.3fs\n", time);
+}
+
+void
+run_simple_reconstruction (Params *params, gchar **argv)
+{
+    LaminoData *data;
+    UfoBaseScheduler *sched = NULL;
+    UfoTaskNode *radio_reader = NULL;
+    UfoTaskNode *dark_reader = NULL;
+    UfoTaskNode *flat_reader = NULL;
+    UfoTaskNode *ffc = NULL;
+    GError *error = NULL;
+
+    if (!params_okay (params)) {
+        err ("Parameters missing.\n");
+        return;
+    }
+
+    data = reco_graph_new (params);
+    g_signal_connect (data->reco, "processed", G_CALLBACK (print_dots), NULL);
+
+    /* Set up input */
+    radio_reader = make_task (data->pm, "reader");
+    update_reader (radio_reader, params);
+
+    if (with_flat_field_correction (params)) {
+        flat_reader = make_task (data->pm, "reader");
+        dark_reader = make_task (data->pm, "reader");
+        ffc = make_task (data->pm, "flat-field-correction");
 
         g_object_set (ffc, "dark-scale", params->dark_scale, NULL);
         g_object_set (flat_reader, "path", params->flats, NULL);
         g_object_set (dark_reader, "path", params->darks, 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_full (data->graph, radio_reader, ffc, 0);
+        ufo_task_graph_connect_nodes_full (data->graph, dark_reader, ffc, 1);
+        ufo_task_graph_connect_nodes_full (data->graph, flat_reader, ffc, 2);
+        ufo_task_graph_connect_nodes (data->graph, ffc, data->pad);
     }
     else {
         warn ("> No flat/dark field correction\n");
-        ufo_task_graph_connect_nodes (graph, radio_reader, pad);
+        ufo_task_graph_connect_nodes (data->graph, radio_reader, data->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);
-
-    g_signal_connect (reco, "processed", G_CALLBACK (print_dots), NULL);
-
+    /* Run reconstruction */
+    info ("Processing ");
     sched = UFO_BASE_SCHEDULER (ufo_scheduler_new (NULL, NULL));
-    ufo_base_scheduler_run (sched, graph, &error);
+    ufo_base_scheduler_run (sched, data->graph, &error);
     check (error);
+    show_time (sched);
 
-    g_object_get (sched, "time", &time, NULL);
-    g_print ("\n");
-    info("Finished in %3.3fs\n", time);
+    /* Clean up */
+    reco_graph_free (data);
 
-    g_object_unref (pm);
-    g_object_unref (graph);
     g_object_unref (sched);
-
     g_object_unref (radio_reader);
-    g_object_unref (pad);
-    g_object_unref (fft1);
-    g_object_unref (fft2);
-    /* g_object_unref (conv); */
-    g_object_unref (reco);
-    g_object_unref (writer);
 
-    if (params->darks != NULL && params->flats != NULL) {
+    if (with_flat_field_correction (params)) {
         g_object_unref (ffc);
         g_object_unref (flat_reader);
         g_object_unref (dark_reader);
     }
 }
+
+static void
+move_data (UfoTaskNode *reader, LaminoData *data)
+{
+    UfoBuffer *buffer;
+    gfloat *mem;
+    gsize size;
+    static int n = 0;
+
+    /* FIXME: we miss one projection */
+    /* We will block here, even though reader sent the "generated" signal */
+    if (data->current == 0) {
+        data->current++;
+        return;
+    }
+
+    buffer = ufo_output_task_get_output_buffer (UFO_OUTPUT_TASK (data->output));
+    mem = ufo_buffer_get_host_array (buffer, NULL);
+    size = ufo_buffer_get_size (buffer);
+    memcpy (data->params->cache + (data->current - 1) * size, mem, size);
+    ufo_output_task_release_output_buffer (UFO_OUTPUT_TASK (data->output), buffer);
+    data->current++;
+
+    if ((n++ % 10) == 0)
+        g_print (".");
+}
+
+static gpointer
+input_thread (LaminoData *data)
+{
+    UfoBuffer *buffer;
+    UfoRequisition req;
+    guint current;
+
+    current = 0;
+    req.n_dims = 2;
+    req.dims[0] = data->params->width;
+    req.dims[1] = data->params->height;
+    buffer = ufo_buffer_new (&req, NULL);
+
+    while (current < data->params->num_radios) {
+        gfloat *mem;
+        gsize size;
+
+        mem = ufo_buffer_get_host_array (buffer, NULL);
+        size = ufo_buffer_get_size (buffer);
+        memcpy (mem, data->params->cache + current * size, size);
+        ufo_input_task_release_input_buffer (UFO_INPUT_TASK (data->input), buffer);
+        current++;
+
+        buffer = ufo_input_task_get_input_buffer (UFO_INPUT_TASK (data->input));
+
+        if ((current % 10) == 0)
+            g_print (".");
+    }
+
+    ufo_input_task_stop (UFO_INPUT_TASK (data->input));
+    return NULL;
+}
+
+void
+run_cached_reconstruction (Params *params, gchar **argv)
+{
+    LaminoData *data;
+    GThread *thread;
+    UfoBaseScheduler *sched = NULL;
+    GError *error = NULL;
+
+    data = reco_graph_new (params);
+
+    if (params->cache == NULL) {
+        UfoTaskGraph *read_graph;
+        UfoTaskNode *radio_reader;
+        UfoBaseScheduler *sched;
+        GError *error = NULL;
+
+        params->cache = g_malloc (params->num_radios * params->width * params->height * sizeof (gfloat));
+
+        /* fill the cache and reconstruct simultaneously */
+        radio_reader = make_task (data->pm, "reader");
+        update_reader (radio_reader, params);
+        data->output = UFO_TASK_NODE (ufo_output_task_new (2));
+
+        read_graph = UFO_TASK_GRAPH (ufo_task_graph_new ());
+        ufo_task_graph_connect_nodes (read_graph, radio_reader, data->output);
+        g_signal_connect (radio_reader, "generated", G_CALLBACK (move_data), data);
+        data->current = 0;
+
+        sched = UFO_BASE_SCHEDULER (ufo_scheduler_new (NULL, NULL));
+        info ("Reading ");
+        ufo_base_scheduler_run (sched, read_graph, &error);
+        check (error);
+        g_print ("\n");
+
+        g_object_unref (radio_reader);
+        g_object_unref (read_graph);
+        g_object_unref (data->output);
+    }
+
+    data->current = 0;
+    data->input = UFO_TASK_NODE (ufo_input_task_new ());
+    ufo_task_graph_connect_nodes (data->graph, data->input, data->pad);
+
+    /* Run reconstruction */
+    info ("Processing ");
+    thread = g_thread_new (NULL, (GThreadFunc) input_thread, data);
+
+    sched = UFO_BASE_SCHEDULER (ufo_scheduler_new (NULL, NULL));
+    ufo_base_scheduler_run (sched, data->graph, &error);
+    check (error);
+    show_time (sched);
+
+    g_thread_unref (thread);
+    g_object_unref (data->input);
+    reco_graph_free (data);
+}

reco.h

@@ -24,11 +24,13 @@ typedef struct {
     guint v_size[3];
     GOptionEntry *entries;
 
-    gfloat *cache;
+    gchar *cache;
 } Params;
 
 gboolean    params_okay                 (Params *params);
 void        run_simple_reconstruction   (Params *params,
                                          gchar **argv);
+void        run_cached_reconstruction   (Params *params,
+                                         gchar **argv);
 
 #endif