/* * gsf.c * * Gerchberg-Saxton-Fienup iteration algorithms * * (c) 2006-2007 Thomas White * (c) 2008 Alex Eggeman * * synth2d - two-dimensional Fourier synthesis * */ #ifdef HAVE_CONFIG_H #include #endif #include #include #include #include #include #include #include "main.h" #include "data.h" #include "displaywindow.h" #include "symmetry.h" #include "gtk-symmetry.h" #include "reflist.h" GtkWidget *gsfw_iterations = NULL; static int gsf_window_open = 0; static int gsf_running = 0; static int gsf_n_iterations = 0; static int gsf_centroid = 0; static int gsf_atomicity = 0; static int gsf_zerofzero = 1; static int gsf_fracsup = 0; static int gsf_threshsup = 1; static double gsf_fracsup_phi = 0.15; static double gsf_hio_beta = 0.15; static double gsf_threshsup_delta = 0; /* i.e. all positive values are in the support */ Symmetry gsf_symmetry = PLANEGROUP_P1; fftw_complex *gsf_drive = NULL; fftw_complex *gsf_drive_old = NULL; unsigned int *gsf_support = NULL; fftw_complex *gsf_in = NULL; fftw_complex *gsf_out = NULL; fftw_plan gsf_plan_i2d; fftw_plan gsf_plan_d2i; ReflectionList *gsf_reflections = NULL; static int gsf_antistag = 0; static void gsf_pin_centroid(fftw_complex *out, fftw_complex *drive, unsigned int width, unsigned int height) { unsigned int x, y; double xc = 0; double yc = 0; double total_int = 0; signed int xs; signed int ys; fftw_complex *out_copy; for ( x=0; xn_reflections; i++ ) { signed int h = reflections->refs[i].h; signed int k = reflections->refs[i].k; if ( h < 0 ) { h = width+h; } if ( k < 0 ) { k = height+k; } in_new[k + height*h][0] = in[k + height*h][0]; in_new[k + height*h][1] = in[k + height*h][1]; } memcpy(in, in_new, width*height*sizeof(fftw_complex)); free(in_new); } /* Perform a single GSF iteration */ void gsf_iteration(int gsf_n_iterations) { unsigned int x, y; double mean; double top; unsigned int mean_n; double fzero; unsigned int i; unsigned int width = data_width(); unsigned int height = data_height(); assert(gsf_drive_old != NULL); assert(gsf_drive != NULL); assert(gsf_in != NULL); assert(gsf_out != NULL); /* Copy this structure if required */ memcpy(gsf_drive_old, gsf_drive, width*height*sizeof(fftw_complex)); /* Transform back to reciprocal space and scale. This time, the driving function is being transformed. */ fftw_execute(gsf_plan_d2i); for ( x=0; xn_reflections; i++ ) { double re, im; double am, ph; signed int h, k; h = gsf_reflections->refs[i].h; k = gsf_reflections->refs[i].k; if ( h < 0 ) { h = width+h; } if ( k < 0 ) { k = height+k; } re = gsf_in[k + height*h][0]; im = gsf_in[k + height*h][1]; ph = atan2(im, re); am = gsf_reflections->refs[i].amplitude; gsf_in[k + height*h][0] = am*cos(ph); gsf_in[k + height*h][1] = am*sin(ph); } fzero = gsf_in[0][0]; gsf_mask(gsf_reflections, gsf_in, width, height); if ( !gsf_zerofzero ) { gsf_in[0][0] = fzero; } /* Impose Symmetry */ symmetry_symmetrise_array(gsf_in, width, height, gsf_symmetry); /* Transform to real space */ fftw_execute(gsf_plan_i2d); if ( gsf_centroid ) { gsf_pin_centroid(gsf_drive, gsf_drive_old, width, height); } memcpy(gsf_out, gsf_drive, width*height*sizeof(fftw_complex)); displaywindow_switchview(); /* Determine the support */ if ( gsf_threshsup && !gsf_fracsup ) { /* Threshold support */ for ( y=0; y gsf_threshsup_delta ) { gsf_support[y+height*x] = 1; } else { gsf_support[y+height*x] = 0; } if ( gsf_antistag && ((x+y)>width) ) { gsf_support[y+height*x] = 0; } } } } else if ( gsf_fracsup ) { /* 'Fractional' support */ double max; unsigned int max_x, max_y; int n, target, i; n = height*width; target = n * gsf_fracsup_phi; i = 0; while ( i < target ) { max_x = 0; max_y = 0; max = -1e10; for ( y=0; y 0) ) { if ( ( am > max) && (!gsf_support[y+height*x]) ) { max_x = x; max_y = y; max = am; } } } } gsf_support[max_y+height*max_x] = 1; i++; } } else { /* No options selected for support - consider everything valid */ for ( y=0; y 0 ) { mean += am; mean_n++; if ( am > top ) top = am; } } } mean = mean / mean_n; for ( y=0; yn_reflections; i++ ) { double am, ph; signed int h, k; h = gsf_reflections->refs[i].h; k = gsf_reflections->refs[i].k; if ( h < 0 ) { h = width+h; } if ( k < 0 ) { k = height+k; } am = gsf_reflections->refs[i].amplitude; ph = (((double)random())/RAND_MAX) * (2*M_PI); gsf_in[k + height*h][0] = am*cos(ph); gsf_in[k + height*h][1] = am*sin(ph); } fftw_execute(gsf_plan_i2d); memcpy(gsf_out, gsf_drive, width*height*sizeof(fftw_complex)); displaywindow_switchview(); } void gsf_initialise(ReflectionList *reflections) { gsf_drive_old = fftw_malloc(data_width()*data_height()*sizeof(fftw_complex)); gsf_drive = fftw_malloc(data_width()*data_height()*sizeof(fftw_complex)); gsf_support = fftw_malloc(data_width()*data_height()*sizeof(unsigned int)); gsf_in = fftw_malloc(data_width()*data_height()*sizeof(fftw_complex)); memset(gsf_in, 0, data_width()*data_height()*sizeof(fftw_complex)); if ( !gsf_out ) gsf_out = fftw_malloc(data_width()*data_height()*sizeof(fftw_complex)); gsf_plan_i2d = fftw_plan_dft_2d(data_width(), data_height(), gsf_in, gsf_drive, FFTW_BACKWARD, FFTW_MEASURE | FFTW_PRESERVE_INPUT); gsf_plan_d2i = fftw_plan_dft_2d(data_width(), data_height(), gsf_drive, gsf_in, FFTW_FORWARD, FFTW_MEASURE | FFTW_PRESERVE_INPUT); gsf_reflections = malloc(sizeof(ReflectionList)); memcpy(gsf_reflections, reflections, sizeof(ReflectionList)); displaywindow_set_realspace(gsf_out, DWR_GSF); gsf_reset(reflections); } static gint gsf_window_go(GtkWidget *gsf_go) { if ( gsf_running ) { return 0; } gsf_running = 1; while ( gsf_running ) { gsf_n_iterations++; gsf_update_iterations(); gsf_iteration(gsf_n_iterations); while ( gtk_events_pending() ) gtk_main_iteration(); } return 0; } static gint gsf_window_stop(GtkWidget *widget, gpointer data) { gsf_running = 0; return 0; } static gint gsf_window_reset(GtkWidget *widget) { gsf_running = 0; main_gsf_reset(); return 0; } static gint gsf_symmetry_changed(GtkWidget *widget, gpointer data) { gsf_symmetry = gtk_symmetry_get_symmetry(GTK_SYMMETRY(widget)); return 0; } void gsf_dialog_open() { GtkWidget *gsfw_window; GtkWidget *gsfw_go; GtkWidget *gsfw_stop; GtkWidget *gsfw_reset; GtkWidget *gsfw_gostop_hbox; GtkWidget *gsfw_centroid; GtkWidget *gsfw_atomicity; GtkWidget *gsfw_zerofzero; GtkWidget *gsfw_antistag; GtkWidget *gsfw_fracsup; GtkWidget *gsfw_fracsup_phi; GtkWidget *gsfw_fracsup_hbox; GtkWidget *gsfw_hio_beta; GtkWidget *gsfw_hio_hbox; GtkWidget *gsfw_threshsup; GtkWidget *gsfw_threshsup_delta; GtkWidget *gsfw_threshsup_hbox; GtkWidget *vbox; GtkWidget *hbox; GtkWidget *sym_define; GtkWidget *gsfw_support_label; GtkWidget *gsfw_process_label; char *text; if ( gsf_window_open ) { return; } gsf_window_open = 1; gsfw_window = gtk_dialog_new_with_buttons("Projection Iteration", GTK_WINDOW(displaywindow_gtkwindow()), GTK_DIALOG_DESTROY_WITH_PARENT, GTK_STOCK_CLOSE, GTK_RESPONSE_CLOSE, NULL); vbox = gtk_vbox_new(FALSE, 0); hbox = gtk_hbox_new(TRUE, 0); gtk_box_pack_start(GTK_BOX(GTK_DIALOG(gsfw_window)->vbox), GTK_WIDGET(hbox), FALSE, FALSE, 7); gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5); gsfw_iterations = gtk_label_new(""); text = malloc(40); snprintf(text, 39, "Iterations performed so far: %i", gsf_n_iterations); gtk_label_set_text(GTK_LABEL(gsfw_iterations), text); free(text); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_iterations), FALSE, FALSE, 5); sym_define = gtk_symmetry_new(2, 2, TRUE); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(sym_define), FALSE, FALSE, 5); g_signal_connect(G_OBJECT(sym_define), "changed", G_CALLBACK(gsf_symmetry_changed), NULL); gsf_symmetry = PLANEGROUP_P1; gsfw_centroid = gtk_check_button_new_with_label("Pin centroid"); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_centroid), FALSE, FALSE, 5); if ( gsf_centroid ) { gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(gsfw_centroid), TRUE); } else { gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(gsfw_centroid), FALSE); } gsfw_atomicity = gtk_check_button_new_with_label("Atomicity"); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_atomicity), FALSE, FALSE, 5); if ( gsf_atomicity ) { gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(gsfw_atomicity), TRUE); } else { gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(gsfw_atomicity), FALSE); } gsfw_zerofzero = gtk_check_button_new_with_label("Zero F(00)"); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_zerofzero), FALSE, FALSE, 5); if ( gsf_zerofzero ) { gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(gsfw_zerofzero), TRUE); } else { gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(gsfw_zerofzero), FALSE); } gsfw_antistag = gtk_check_button_new_with_label("Anti-stagnation"); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_antistag), FALSE, FALSE, 5); if ( gsf_antistag ) { gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(gsfw_antistag), TRUE); } else { gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(gsfw_antistag), FALSE); } gsfw_support_label = gtk_label_new(""); gtk_label_set_markup(GTK_LABEL(gsfw_support_label), "Dynamic Support"); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_support_label), FALSE, FALSE, 10); gsfw_fracsup = gtk_radio_button_new_with_label(NULL, "Fraction of most positive points"); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_fracsup), FALSE, FALSE, 5); if ( gsf_fracsup ) { gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(gsfw_fracsup), TRUE); } else { gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(gsfw_fracsup), FALSE); } gsfw_fracsup_hbox = gtk_hbox_new(FALSE, 0); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_fracsup_hbox), FALSE, FALSE, 5); gsfw_fracsup_phi_label = gtk_label_new("Φ:"); gtk_box_pack_start(GTK_BOX(gsfw_fracsup_hbox), GTK_WIDGET(gsfw_fracsup_phi_label), FALSE, FALSE, 5); gsfw_fracsup_phi = gtk_hscale_new_with_range(0, 1, 0.05); gtk_scale_set_value_pos(GTK_SCALE(gsfw_fracsup_phi), GTK_POS_RIGHT); gtk_box_pack_start(GTK_BOX(gsfw_fracsup_hbox), GTK_WIDGET(gsfw_fracsup_phi), TRUE, TRUE, 5); gtk_range_set_value(GTK_RANGE(gsfw_fracsup_phi), gsf_hio_beta); if ( gsf_fracsup) { gtk_widget_set_sensitive(gsfw_fracsup_phi, TRUE); gtk_widget_set_sensitive(gsfw_fracsup_phi_label, TRUE); } else { gtk_widget_set_sensitive(gsfw_fracsup_phi, FALSE); gtk_widget_set_sensitive(gsfw_fracsup_phi_label, FALSE); } gsfw_threshsup = gtk_radio_button_new_with_label_from_widget(GTK_RADIO_BUTTON(gsfw_fracsup), "Points above threshold"); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_threshsup), FALSE, FALSE, 5); if ( gsf_threshsup ) { gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(gsfw_threshsup), TRUE); } else { gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(gsfw_threshsup), FALSE); } gsfw_threshsup_hbox = gtk_hbox_new(FALSE, 0); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_threshsup_hbox), FALSE, FALSE, 5); gsfw_threshsup_delta_label = gtk_label_new("δ:"); gtk_box_pack_start(GTK_BOX(gsfw_threshsup_hbox), GTK_WIDGET(gsfw_threshsup_delta_label), FALSE, FALSE, 5); gsfw_threshsup_delta = gtk_hscale_new_with_range(-100, 100, 1); gtk_scale_set_value_pos(GTK_SCALE(gsfw_threshsup_delta), GTK_POS_RIGHT); gtk_box_pack_start(GTK_BOX(gsfw_threshsup_hbox), GTK_WIDGET(gsfw_threshsup_delta), TRUE, TRUE, 5); gtk_range_set_value(GTK_RANGE(gsfw_threshsup_delta), gsf_threshsup_delta); if ( gsf_threshsup ) { gtk_widget_set_sensitive(gsfw_threshsup_delta, TRUE); gtk_widget_set_sensitive(gsfw_threshsup_delta_label, TRUE); } else { gtk_widget_set_sensitive(gsfw_threshsup_delta, FALSE); gtk_widget_set_sensitive(gsfw_threshsup_delta_label, FALSE); } gsfw_process_label = gtk_label_new(""); gtk_label_set_markup(GTK_LABEL(gsfw_process_label), "Points Outside Support"); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_process_label), FALSE, FALSE, 10); gsfw_hio_hbox = gtk_hbox_new(FALSE, 0); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_hio_hbox), FALSE, FALSE, 5); gsfw_hio_beta_label = gtk_label_new("β:"); gtk_box_pack_start(GTK_BOX(gsfw_hio_hbox), GTK_WIDGET(gsfw_hio_beta_label), FALSE, FALSE, 5); gsfw_hio_beta = gtk_hscale_new_with_range(0, 1, 0.05); gtk_scale_set_value_pos(GTK_SCALE(gsfw_hio_beta), GTK_POS_RIGHT); gtk_box_pack_start(GTK_BOX(gsfw_hio_hbox), GTK_WIDGET(gsfw_hio_beta), TRUE, TRUE, 5); gtk_range_set_value(GTK_RANGE(gsfw_hio_beta), gsf_hio_beta); gsfw_gostop_hbox = gtk_hbox_new(FALSE, 0); gsfw_go = gtk_button_new_with_label("Run"); gtk_button_set_image(GTK_BUTTON(gsfw_go), gtk_image_new_from_stock(GTK_STOCK_MEDIA_PLAY, GTK_ICON_SIZE_BUTTON)); gsfw_stop = gtk_button_new_from_stock(GTK_STOCK_MEDIA_STOP); gsfw_reset = gtk_button_new_with_label("Reset"); gtk_button_set_image(GTK_BUTTON(gsfw_reset), gtk_image_new_from_stock(GTK_STOCK_MEDIA_PREVIOUS, GTK_ICON_SIZE_BUTTON)); gtk_box_pack_start(GTK_BOX(gsfw_gostop_hbox), GTK_WIDGET(gsfw_reset), FALSE, FALSE, 5); gtk_box_pack_start(GTK_BOX(gsfw_gostop_hbox), GTK_WIDGET(gsfw_stop), FALSE, FALSE, 5); gtk_box_pack_start(GTK_BOX(gsfw_gostop_hbox), GTK_WIDGET(gsfw_go), FALSE, FALSE, 5); gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(gsfw_gostop_hbox), FALSE, FALSE, 5); g_signal_connect(G_OBJECT(gsfw_go), "clicked", G_CALLBACK(gsf_window_go), NULL); g_signal_connect(G_OBJECT(gsfw_stop), "clicked", G_CALLBACK(gsf_window_stop), NULL); g_signal_connect(G_OBJECT(gsfw_reset), "clicked", G_CALLBACK(gsf_window_reset), NULL); g_signal_connect(G_OBJECT(gsfw_window), "response", G_CALLBACK(gsf_window_close), NULL); g_signal_connect(G_OBJECT(gsfw_window), "delete-event", G_CALLBACK(gsf_window_stop), NULL); g_signal_connect(G_OBJECT(gsfw_centroid), "toggled", G_CALLBACK(gsf_centroid_toggled), NULL); g_signal_connect(G_OBJECT(gsfw_atomicity), "toggled", G_CALLBACK(gsf_atomicity_toggled), NULL); g_signal_connect(G_OBJECT(gsfw_zerofzero), "toggled", G_CALLBACK(gsf_zerofzero_toggled), NULL); g_signal_connect(G_OBJECT(gsfw_fracsup), "toggled", G_CALLBACK(gsf_fracsup_toggled), gsfw_fracsup_phi); g_signal_connect(G_OBJECT(gsfw_fracsup_phi), "value-changed", G_CALLBACK(gsf_fracsup_phichanged), NULL); g_signal_connect(G_OBJECT(gsfw_threshsup), "toggled", G_CALLBACK(gsf_threshsup_toggled), gsfw_threshsup_delta); g_signal_connect(G_OBJECT(gsfw_threshsup_delta), "value-changed", G_CALLBACK(gsf_threshsup_deltachanged), NULL); g_signal_connect(G_OBJECT(gsfw_hio_beta), "value-changed", G_CALLBACK(gsf_hio_betachanged), NULL); g_signal_connect(G_OBJECT(gsfw_antistag), "toggled", G_CALLBACK(gsf_antistag_toggled), NULL); g_signal_connect(G_OBJECT(displaywindow_gtkwindow()), "delete-event", G_CALLBACK(gsf_window_stop), NULL); main_gsf_initialise(); displaywindow_forceview(DWV_REALSPACE); gtk_widget_show_all(gsfw_window); }