123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397 |
- #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)
- {
- static int n = 0;
- if ((n++ % 10) == 0)
- g_print (".");
- }
- 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;
- }
- gboolean
- params_okay (Params *params)
- {
- return params->width != 0 &&
- params->height != 0 &&
- params->num_radios != 0 &&
- params->radios != NULL &&
- params->theta < G_MAXDOUBLE &&
- params->px < G_MAXDOUBLE &&
- params->py < G_MAXDOUBLE;
- }
- 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;
- data = g_malloc0 (sizeof (LaminoData));
- data->params = params;
- data->pm = ufo_plugin_manager_new (NULL);
- data->graph = UFO_TASK_GRAPH (ufo_task_graph_new ());
- 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;
- padded_height = next_power_of_two ((guint32) params->height + 1);
- xl = params->px - params->width / 2;
- xr = padded_width - params->width - xl;
- yt = params->py - params->height / 2;
- yb = padded_height - params->height - yt;
- 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",
- params->px, params->py, params->px_variation);
- 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",
- padded_width, padded_height, xl, xr, yt, yb);
- 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 (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 (data->ramp,
- "width", padded_width,
- "height", padded_height,
- "theta", theta_rad,
- "tau", params->tau,
- "fwidth", fwidth,
- NULL);
- g_object_set (data->reco,
- "angle-step", angle_step,
- "theta", theta_rad,
- "psi", params->psi,
- "proj-ox", params->px,
- "proj-oy", params->py,
- "proj-ox-variation", params->px_variation,
- "vol-ox", params->v_size[0] / 2 + params->v_origin[0],
- "vol-oy", params->v_size[1] / 2 + params->v_origin[1],
- "vol-oz", params->v_size[2] / 2 + 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 (data->writer,
- "filename", params->output,
- NULL);
- 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 (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 (data->graph, radio_reader, data->pad);
- }
- /* Run reconstruction */
- info ("Processing ");
- sched = UFO_BASE_SCHEDULER (ufo_scheduler_new (NULL, NULL));
- ufo_base_scheduler_run (sched, data->graph, &error);
- check (error);
- show_time (sched);
- /* Clean up */
- reco_graph_free (data);
- g_object_unref (sched);
- g_object_unref (radio_reader);
- 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);
- }
|