aboutsummaryrefslogtreecommitdiff
path: root/src/cflip.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/cflip.c')
-rw-r--r--src/cflip.c910
1 files changed, 910 insertions, 0 deletions
diff --git a/src/cflip.c b/src/cflip.c
new file mode 100644
index 0000000..2b3dfd6
--- /dev/null
+++ b/src/cflip.c
@@ -0,0 +1,910 @@
+/*
+ * cflip.c
+ *
+ * Charge flipping iteration algorithms
+ *
+ * (c) 2006-2007 Thomas White <taw27@cam.ac.uk>
+ * (c) 2008 Alex Eggeman <ase25@cam.ac.uk>
+ *
+ * synth2d - two-dimensional Fourier synthesis
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <gtk/gtk.h>
+#include <fftw3.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+
+#include "main.h"
+#include "data.h"
+#include "displaywindow.h"
+#include "symmetry.h"
+#include "gtk-symmetry.h"
+#include "reflist.h"
+#include "dpsynth.h"
+
+static int cf_window_open = 0;
+
+typedef struct struct_cflipdialog {
+
+ /* Dialog box bits */
+ GtkWidget *cf_window;
+ GtkWidget *cf_go;
+ GtkWidget *cf_stop;
+ GtkWidget *cf_next;
+ GtkWidget *cf_reset;
+ GtkWidget *cf_iterations;
+
+ /* Algorithm guts */
+ ReflectionList *cf_reflections; /* List of 'modulus constraints' */
+ ReflectionList *cf_listout; /* List of structure factors returned at the end of an iteration */
+ Symmetry cf_symmetry; /* Symmetry of the algorithm */
+ fftw_complex *cf_drive_old;
+ fftw_complex *cf_drive;
+ fftw_complex *cf_in; /* Input realspace array */
+ fftw_complex *cf_out; /* Array on the other end of the FFT of 'in' */
+ unsigned int *cf_support;
+ unsigned int cf_tangent; /*Coehlo Constraints*/
+ double *intensities; /*Array of pixel intensities for auto-thresholding*/
+ double *amplitudes; /*Array of reflection intiensities for auto-thresholding*/
+ fftw_plan cf_plan_i2d;
+ fftw_plan cf_plan_d2i;
+ unsigned int width; /* Width of the array in pixels */
+ unsigned int height; /* Height of the array in pixels */
+ double delta;
+ //DPSynthWindow *dpsynth;
+
+ /* Tracking */
+
+ unsigned int running;
+ unsigned int n_iterations; /* Number of iterations performed so far */
+ unsigned int period; /* Number of iterations per symmetry application */
+
+} CFlipDialog;
+
+int compare_doubles (const void *xp,const void *yp) {
+
+ const double *x, *y;
+ int z;
+
+ x=xp;
+ y=yp;
+
+
+ if (*x > *y) {
+ z = 1;
+ } else if (*x < *y) {
+ z = -1;
+ } else {
+ z = 0;
+ }
+ return z;
+}
+
+static int cflip_tangent (ReflectionList *listout, double emin, double gmin) {
+
+ ReflectionList *strongest_reflections;
+ TripletList *triplet_list;
+ unsigned int i;
+ unsigned int p;
+ unsigned int q;
+ unsigned int n_total_triplets;
+ double E1;
+ double E2;
+ double E3;
+ double G;
+ double Th;
+ double Bh;
+ double *Mh;
+ double maxMh;
+ double *phase_tangent;
+ double alpha_tan, phi_a, phi_tan;
+ unsigned int r, t;
+ signed int ph, pk, pl, pph, qh, qk, ql, qph, rh, rk, rl, rph;
+
+ strongest_reflections = reflist_new_parent(listout);
+ /* Initialise triplet list */
+ triplet_list = malloc(sizeof(TripletList));
+ triplet_list->n_triplets = 0;
+
+ for ( i=0; i<listout->n_reflections; i++ ) {
+ listout->refs[i].parent_index = i;
+ }
+
+ for ( i=1; i<listout->n_reflections; i++ ) {
+ if ( (listout->refs[i].amplitude >= emin) && (listout->refs[i].h != 0) && (listout->refs[i].k != 0) ) {
+ reflist_addref_parent(strongest_reflections, listout->refs[i].h, listout->refs[i].k, listout->refs[i].l,
+ listout->refs[i].parent_index);
+ }
+ }
+
+ Mh = malloc(strongest_reflections->n_reflections*sizeof(double));
+
+ phase_tangent = malloc(strongest_reflections->n_reflections*sizeof(double));
+
+ /* Iterate doubly over the reflections */
+ n_total_triplets = 0;
+ for ( p=1; p<strongest_reflections->n_reflections; p++ ) {
+ for ( q=1; q<strongest_reflections->n_reflections; q++ ) {
+
+ ph = strongest_reflections->refs[p].h;
+ pk = strongest_reflections->refs[p].k;
+ pl = strongest_reflections->refs[p].l;
+ pph = strongest_reflections->refs[p].phase_known;
+ qh = strongest_reflections->refs[q].h;
+ qk = strongest_reflections->refs[q].k;
+ ql = strongest_reflections->refs[q].l;
+ qph = strongest_reflections->refs[q].phase_known;
+
+ assert(ph || pk || pl);
+ assert(qh || qk || ql);
+
+ /* Determine the third reflection in the triplet */
+ rh = -ph-qh;
+ rk = -pk-qk;
+ if ( pl != ql ) {
+ printf("Coelho: pl != ql: %3i %3i\n", pl, ql);
+ continue;
+ }
+ //assert(pl == ql); /* This failure mode should have been weeded out ages ago */
+ rl = pl;
+
+ /* See if the third reflection is available */
+ if ( (r = reflist_inlist(strongest_reflections, rh, rk, rl)) ) {
+
+ /* Trust me on this ... */
+ E1 = strongest_reflections->parent->refs[strongest_reflections->refs[p].parent_index].amplitude;
+ E2 = strongest_reflections->parent->refs[strongest_reflections->refs[q].parent_index].amplitude;
+ E3 = strongest_reflections->parent->refs[strongest_reflections->refs[r].parent_index].amplitude;
+
+ G = E1 * E2 * E3;
+
+ rph = strongest_reflections->parent->refs[strongest_reflections->refs[r].parent_index].phase_known;
+
+ /* Convincing enough? */
+ if ( G > gmin ) {
+ triplet_list->triplets[triplet_list->n_triplets].p.h = ph;
+ triplet_list->triplets[triplet_list->n_triplets].p.k = pk;
+ triplet_list->triplets[triplet_list->n_triplets].p.l = pl;
+ triplet_list->triplets[triplet_list->n_triplets].q.h = qh;
+ triplet_list->triplets[triplet_list->n_triplets].q.k = qk;
+ triplet_list->triplets[triplet_list->n_triplets].q.l = ql;
+ triplet_list->triplets[triplet_list->n_triplets].G = G;
+ triplet_list->n_triplets++;
+ if ( triplet_list->n_triplets == MAX_TRIPLETS ) {
+ printf("DM: Too many triplets!\n");
+ return NULL;
+ }
+
+ }
+ n_total_triplets++;
+
+ }
+
+ }
+ }
+ maxMh = 0;
+
+ for ( i=0; i<strongest_reflections->n_reflections; i++ ) {
+
+ signed int h = strongest_reflections->refs[i].h;
+ signed int k = strongest_reflections->refs[i].k;
+
+ for ( t=0; t<triplet_list->n_triplets; t++ ) {
+
+ signed int ph = triplet_list->triplets[t].p.h;
+ signed int pk = triplet_list->triplets[t].p.k;
+
+
+ /* Does this triplet concern this reflection? */
+ if ( (ph == h) && (pk == k) ) {
+
+ signed int qh = triplet_list->triplets[t].q.h;
+ signed int qk = triplet_list->triplets[t].q.k;
+ signed int rh = -ph-qh;
+ signed int rk = -pk-qk;
+ double G = triplet_list->triplets[t].G;
+
+ /* Are the other reflections in the triplet still in the list? */
+ if ( reflist_inlist_2d(listout, qh, qk) && reflist_inlist_2d(listout, rh, rk) ) {
+ Th += G * sin(( qph - rph));
+ Bh += G * cos(( qph - rph));
+ }
+
+
+ }
+
+ }
+ phase_tangent[i] = atan( Th / Bh );
+ Mh[i] = sqrt( (Th*Th) / (Bh*Bh) );
+ if ( Mh[i] > maxMh ) {
+ maxMh = Mh[i];
+ }
+
+ }
+
+ for ( i=0; i<strongest_reflections->n_reflections; i++ ) {
+
+ alpha_tan = Mh[i] / maxMh;
+
+ phi_a = strongest_reflections->parent->refs[strongest_reflections->refs[i].parent_index].phase_known;
+
+ phi_tan = phi_a + alpha_tan*(phase_tangent[i]-phi_a);
+
+ listout->refs[strongest_reflections->refs[i].parent_index].phase_known = phi_tan;
+
+ }
+
+ free(strongest_reflections);
+ free(triplet_list);
+
+ return 0;
+
+}
+
+void cflip_mask(CFlipDialog *cflip_dialog, fftw_complex *in) {
+
+
+ unsigned int width = cflip_dialog->width;
+ unsigned int height = cflip_dialog->height;
+ unsigned int i;
+ fftw_complex *in_new = malloc(width*height*sizeof(fftw_complex));
+ bzero(in_new, width*height*sizeof(fftw_complex));
+ ReflectionList *reflections;
+
+ reflections = cflip_dialog->cf_reflections;
+
+ for ( i=1; i<reflections->n_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 CF iteration */
+void cflip_iteration(CFlipDialog *cflip_dialog) {
+
+ double re, im;
+ double am, ph, max, i_max;
+ signed int h, k;
+ unsigned int j, indexj;
+ unsigned int x, y;
+ double fzero;
+ double emin = 0;
+ double gmin = 0;
+ double mean, neg;
+ double coehlo_thresh;
+ unsigned int mean_n;
+ double meansq;
+ double sig;
+ double delta;
+ unsigned int i;
+ unsigned int width = cflip_dialog->width;
+ unsigned int height = cflip_dialog->height;
+ ReflectionList *reflections = cflip_dialog->cf_reflections;
+ ReflectionList *listout = cflip_dialog->cf_listout;
+
+ cflip_dialog->intensities = malloc(width*height*sizeof(double));
+ cflip_dialog->amplitudes = malloc(reflections->n_reflections*sizeof(double));
+
+// listout = malloc(sizeof(ReflectionList));
+ listout = NULL;
+ assert(cflip_dialog->cf_drive_old != NULL);
+ assert(cflip_dialog->cf_drive != NULL);
+ assert(cflip_dialog->cf_in != NULL);
+ assert(cflip_dialog->cf_out != NULL);
+
+ if ( !cflip_dialog->cf_tangent ) {
+
+ /* Transform back to reciprocal space and scale. This time, the driving function is
+ being transformed. */
+ fftw_execute(cflip_dialog->cf_plan_d2i);
+ for ( x=0; x<width; x++ ) {
+ for ( y=0; y<((height/2)+1); y++ ) {
+ cflip_dialog->cf_in[y+((height/2)+1)*x][0] = cflip_dialog->cf_in[y+((height/2)+1)*x][0] / (width*height); /* Re */
+ cflip_dialog->cf_in[y+((height/2)+1)*x][1] = cflip_dialog->cf_in[y+((height/2)+1)*x][1] / (width*height); /* Im */
+ }
+ }
+
+ /* Impose observed intensities */
+ for ( i=0; i<reflections->n_reflections; i++ ) {
+
+ h = reflections->refs[i].h;
+ k = reflections->refs[i].k;
+
+ if ( h < 0 ) { h = width+h; }
+ if ( k < 0 ) { k = height+k; }
+
+ re = cflip_dialog->cf_in[k + height*h][0];
+ im = cflip_dialog->cf_in[k + height*h][1];
+ ph = atan2(im, re);
+ am = reflections->refs[i].amplitude;
+
+ cflip_dialog->cf_in[k + height*h][0] = am*cos(ph);
+ cflip_dialog->cf_in[k + height*h][1] = am*sin(ph);
+
+ }
+
+ fzero = cflip_dialog->cf_in[0][0];
+ cflip_mask(cflip_dialog, cflip_dialog->cf_in);
+
+ /* Impose symmetry */
+ if ( cflip_dialog->n_iterations % cflip_dialog->period == 0 ) {
+
+ listout = reflist_new_from_array(cflip_dialog->cf_in, width, height);
+ symmetry_symmetrise(listout, cflip_dialog->cf_symmetry, SYMFLAG_PHASES_KNOWN);
+ reflist_fill_array(cflip_dialog->cf_in, listout, width, height);
+ //dpsynth_update(cflip_dialog->dpsynth, listout);
+ reflist_free(listout);
+
+
+ } else {
+
+ listout = reflist_new_from_array(cflip_dialog->cf_in, width, height);
+ symmetry_symmetrise(listout, cflip_dialog->cf_symmetry & SYMMETRY_FRIEDEL, SYMFLAG_PHASES_KNOWN);
+ reflist_fill_array(cflip_dialog->cf_in, listout, width, height);
+ //dpsynth_update(cflip_dialog->dpsynth, listout);
+ reflist_free(listout);
+
+ }
+
+
+
+ /* Transform to real space */
+ fftw_execute(cflip_dialog->cf_plan_i2d);
+
+ memcpy(cflip_dialog->cf_out, cflip_dialog->cf_drive, width*height*sizeof(fftw_complex));
+ displaywindow_switchview();
+ mean = 0; mean_n = 0; meansq = 0; sig = 0; delta = 0;
+ for ( y=0; y<height; y++ ) {
+ for ( x=0; x<width; x++ ) {
+ re = cflip_dialog->cf_drive[y+height*x][0];
+ im = cflip_dialog->cf_drive[y+height*x][1];
+ am = sqrt(re*re + im*im);
+ mean += am;
+ mean_n += 1;
+ meansq += am * am;
+ }
+ }
+
+ sig = sqrt((meansq/mean_n)-((mean/mean_n) * (mean/mean_n)));
+ delta = sig * cflip_dialog->delta;
+
+ /* Determine the support */
+ for ( y=0; y<height; y++ ) {
+ for ( x=0; x<width; x++ ) {
+ if ( cflip_dialog->cf_drive[y+height*x][0] > delta ) {
+ cflip_dialog->cf_support[y+height*x] = 1;
+ } else {
+ cflip_dialog->cf_support[y+height*x] = 0;
+ }
+ }
+ }
+
+
+ /* Impose non-negativity, sharpen peaks and apply feedback */
+
+ /*mean = 0; neg = 0;
+
+ for ( y=0; y<height; y++ ) {
+ for ( x=0; x<width; x++ ) {
+ re = cflip_dialog->cf_drive[y+height*x][0];
+ im = cflip_dialog->cf_drive[y+height*x][1];
+ am = sqrt(re*re + im*im);
+
+ mean += am;
+ if ( re<0 ) {
+ neg += am;
+ }
+
+ }
+ }
+
+ mean = mean /(width * height);
+
+ double entropy;
+ entropy = 0;
+ for ( y=0; y<height; y++ ) {
+ for ( x=0; x<width; x++ ) {
+ entropy += (am/mean) * log (am/mean);
+ }
+ }
+ entropy = -entropy;
+ printf("step %i; Entropy : %f, Negativity : %f \n", cflip_dialog->n_iterations, entropy, neg);*/
+
+
+ for ( y=0; y<height; y++ ) {
+ for ( x=0; x<width; x++ ) {
+
+ if ( !cflip_dialog->cf_support[y+height*x] ) {
+
+ cflip_dialog->cf_drive[y+height*x][0] = -cflip_dialog->cf_drive[y+height*x][0];
+ cflip_dialog->cf_drive[y+height*x][1] = -cflip_dialog->cf_drive[y+height*x][1];
+
+ }
+
+ }
+ }
+
+ } else { /*Rework the Algorithm from here to incorporate the Coehlo constraints*/
+
+ /* Transform back to reciprocal space and scale. This time, the driving function is being transformed. */
+ fftw_execute(cflip_dialog->cf_plan_d2i);
+ for ( x=0; x<width; x++ ) {
+ for ( y=0; y<((height/2)+1); y++ ) {
+ cflip_dialog->cf_in[y+((height/2)+1)*x][0] = cflip_dialog->cf_in[y+((height/2)+1)*x][0] / (width*height); /* Re */
+ cflip_dialog->cf_in[y+((height/2)+1)*x][1] = cflip_dialog->cf_in[y+((height/2)+1)*x][1] / (width*height); /* Im */
+ }
+ }
+
+ /* Do Tangent Formula Wizardry */
+ cflip_mask(cflip_dialog, cflip_dialog->cf_in);
+ listout = reflist_new_from_array(cflip_dialog->cf_in, width, height);
+
+ i_max = 0;
+
+ for ( i=0; i<listout->n_reflections; i++ ) {
+
+ cflip_dialog->amplitudes[i]=listout->refs[i].amplitude;
+
+ if ( listout->refs[i].amplitude > i_max ) {
+ i_max = listout->refs[i].amplitude;
+ }
+ }
+
+ emin = 0.5 * i_max;
+ gmin = 0.2 * i_max;
+
+ cflip_tangent (listout, emin, gmin);
+
+ i_max = 0; j=0; max = 0; indexj=0;
+
+ for ( i=0; i<reflections->n_reflections; i++ ) {
+
+ cflip_dialog->amplitudes[i]=reflections->refs[i].amplitude;
+ if ( reflections->refs[i].amplitude > i_max ) {
+ i_max = reflections->refs[i].amplitude;
+ }
+ }
+
+ qsort(cflip_dialog->amplitudes, listout->n_reflections, sizeof(double), compare_doubles);
+
+ j = (int)(listout->n_reflections / 2);
+ j = (int)(listout->n_reflections / 4);
+ max = cflip_dialog->amplitudes[j];
+
+ /* Impose observed intensities */
+ //unsigned int n_sub = 0;
+ for ( i=0; i<listout->n_reflections; i++ ) {
+
+ h = listout->refs[i].h;
+ k = listout->refs[i].k;
+
+ indexj = reflist_inlist(reflections, h, k, 0);
+
+ if ( h < 0 ) { h = width+h; }
+ if ( k < 0 ) { k = height+k; }
+
+ ph = listout->refs[i].phase_known;
+ am = reflections->refs[indexj].amplitude;
+ if (am < max ) {
+ am = 0;
+ //n_sub += 1;
+ }
+
+ //printf("number of Zero reflections is %i.\n", n_sub);
+ cflip_dialog->cf_in[k + height*h][0] = am*cos(ph);
+ cflip_dialog->cf_in[k + height*h][1] = am*sin(ph);
+
+ }
+
+ fzero = cflip_dialog->cf_in[0][0];
+ //free(listout);
+
+ /* Impose symmetry */
+ if ( cflip_dialog->n_iterations % cflip_dialog->period == 0 ) {
+
+ listout = reflist_new_from_array(cflip_dialog->cf_in, width, height);
+ symmetry_symmetrise(listout, cflip_dialog->cf_symmetry, SYMFLAG_PHASES_KNOWN);
+ reflist_fill_array(cflip_dialog->cf_in, listout, width, height);
+ //dpsynth_update(cflip_dialog->dpsynth, listout);
+ reflist_free(listout);
+
+
+ } else {
+
+ listout = reflist_new_from_array(cflip_dialog->cf_in, width, height);
+ symmetry_symmetrise(listout, cflip_dialog->cf_symmetry & SYMMETRY_FRIEDEL, SYMFLAG_PHASES_KNOWN);
+ reflist_fill_array(cflip_dialog->cf_in, listout, width, height);
+ //dpsynth_update(cflip_dialog->dpsynth, listout);
+ reflist_free(listout);
+
+ }
+
+ /* Transform to real space */
+ fftw_execute(cflip_dialog->cf_plan_i2d);
+
+ max = 0;
+
+ memcpy(cflip_dialog->cf_out, cflip_dialog->cf_drive, width*height*sizeof(fftw_complex));
+ displaywindow_switchview();
+
+ /* Find the Threshold Real-Space Intensity */
+ mean = 0; mean_n = 0; sig = 0; delta = 0;
+ for ( y=0; y<height; y++ ) {
+ for ( x=0; x<width; x++ ) {
+ re = cflip_dialog->cf_drive[y+height*x][0];
+ im = cflip_dialog->cf_drive[y+height*x][1];
+ am = sqrt((re*re) + (im*im));
+ if ( am > max ) {
+ max = am;
+ }
+
+ }
+ }
+
+ for ( y=0; y<height; y++ ) {
+ for ( x=0; x<width; x++ ) {
+ re = cflip_dialog->cf_drive[y+height*x][0];
+ im = cflip_dialog->cf_drive[y+height*x][1];
+ am = sqrt(re*re + im*im) / max;
+ mean += am;
+ cflip_dialog->intensities[y+height*x]=am;
+ }
+ }
+
+ qsort(cflip_dialog->intensities, mean_n, sizeof(double), compare_doubles);
+
+ coehlo_thresh = 0;
+ i=0;
+ while ( coehlo_thresh < (0.6 * mean) ) {
+
+ coehlo_thresh += cflip_dialog->intensities[i];
+ delta = cflip_dialog->intensities[i];
+ i += 1;
+ }
+
+ /* Determine the support */
+ for ( y=0; y<height; y++ ) {
+ for ( x=0; x<width; x++ ) {
+ if ( cflip_dialog->cf_drive[y+height*x][0] > delta ) {
+ cflip_dialog->cf_support[y+height*x] = 1;
+ } else {
+ cflip_dialog->cf_support[y+height*x] = 0;
+ }
+ }
+ }
+
+ /* Perform the Charge Flip */
+ for ( y=0; y<height; y++ ) {
+ for ( x=0; x<width; x++ ) {
+
+ if ( !cflip_dialog->cf_support[y+height*x] ) {
+
+ cflip_dialog->cf_drive[y+height*x][0] = -cflip_dialog->cf_drive[y+height*x][0];
+ cflip_dialog->cf_drive[y+height*x][1] = -cflip_dialog->cf_drive[y+height*x][1];
+
+ } else {
+ //cflip_dialog->cf_drive[y+height*x][0] = delta + sqrt((cflip_dialog->cf_drive[y+height*x][0])-delta);
+ cflip_dialog->cf_drive[y+height*x][0] = cflip_dialog->cf_drive[y+height*x][0];
+ cflip_dialog->cf_drive[y+height*x][1] = cflip_dialog->cf_drive[y+height*x][1];
+
+ }
+ }
+ }
+ }
+free(cflip_dialog->intensities);
+free(cflip_dialog->amplitudes);
+}
+
+static gint cflip_threshsup_deltachanged(GtkWidget *widget, CFlipDialog *cflip_dialog) {
+
+ cflip_dialog->delta = gtk_range_get_value(GTK_RANGE(widget));
+ return 0;
+}
+
+void cflip_finalise(CFlipDialog *cflip_dialog) {
+
+ free(cflip_dialog->cf_drive); cflip_dialog->cf_drive = NULL;
+ free(cflip_dialog->cf_drive_old); cflip_dialog->cf_drive_old = NULL;
+ free(cflip_dialog->cf_support); cflip_dialog->cf_support = NULL;
+ free(cflip_dialog->cf_in); cflip_dialog->cf_in = NULL;
+ reflist_free(cflip_dialog->cf_reflections); cflip_dialog->cf_reflections = NULL;
+ /* Don't free cf_out, since it is now backing up the display */
+
+ fftw_destroy_plan(cflip_dialog->cf_plan_i2d);
+ fftw_destroy_plan(cflip_dialog->cf_plan_d2i);
+
+}
+
+
+static gint cflip_window_periodchanged(GtkWidget *cfw_period, CFlipDialog *cflip_dialog) {
+
+ const char *str;
+
+ str = gtk_entry_get_text(GTK_ENTRY(cfw_period));
+ cflip_dialog->period = atoi(str);
+ return 0;
+}
+
+static void cflip_window_close(GtkWidget *widget, gint arg1, CFlipDialog *cflip_dialog) {
+
+ if ( cflip_dialog->running ) {
+ cflip_dialog->running = 0;
+ }
+ cf_window_open = 0;
+
+ cflip_finalise(cflip_dialog);
+
+ gtk_widget_destroy(widget);
+
+}
+
+static void cflip_update_iterations(CFlipDialog *cflip_dialog) {
+
+ char *text = malloc(40);
+ snprintf(text, 39, "Iterations performed so far: %i", cflip_dialog->n_iterations);
+ gtk_label_set_text(GTK_LABEL(cflip_dialog->cf_iterations), text);
+ free(text);
+
+}
+
+void cflip_reset(CFlipDialog *cflip_dialog) {
+
+ reflist_free(cflip_dialog->cf_reflections);
+ cflip_dialog->cf_reflections = reflist_copy(main_reflist());
+
+ unsigned int width = cflip_dialog->width;
+ unsigned int height = cflip_dialog->height;
+ cflip_dialog->n_iterations = 0;
+ cflip_update_iterations(cflip_dialog);
+ ReflectionList *reflections = cflip_dialog->cf_reflections;
+
+ unsigned int i;
+
+ for ( i=0; i<reflections->n_reflections; i++ ) {
+
+ double am, ph;
+ signed int h, k;
+
+ h = reflections->refs[i].h;
+ k = reflections->refs[i].k;
+
+ if ( h < 0 ) { h = width+h; }
+ if ( k < 0 ) { k = height+k; }
+
+ am = reflections->refs[i].amplitude;
+ ph = (((double)random())/RAND_MAX) * (2*M_PI);
+
+ cflip_dialog->cf_in[k + height*h][0] = am*cos(ph);
+ cflip_dialog->cf_in[k + height*h][1] = am*sin(ph);
+
+ }
+ symmetry_symmetrise_array(cflip_dialog->cf_in, width, height, cflip_dialog->cf_symmetry & SYMMETRY_FRIEDEL);
+ fftw_execute(cflip_dialog->cf_plan_i2d);
+
+ memcpy(cflip_dialog->cf_out, cflip_dialog->cf_drive, width*height*sizeof(fftw_complex));
+ displaywindow_switchview();
+
+ //dpsynth_update(cflip_dialog->dpsynth, reflections);
+
+}
+
+void cflip_initialise(CFlipDialog *cflip_dialog) {
+
+ ReflectionList *reflections = reflist_copy(main_reflist());
+ Symmetry symmetry;
+ cflip_dialog->width = data_width();
+ cflip_dialog->height = data_height();
+ cflip_dialog->n_iterations = 0;
+ cflip_dialog->period = 1;
+ unsigned int width = cflip_dialog->width;
+ unsigned int height = cflip_dialog->height;
+ symmetry = cflip_dialog->cf_symmetry;
+
+ cflip_dialog->cf_drive_old = fftw_malloc(width * height * sizeof(fftw_complex));
+ cflip_dialog->cf_drive = fftw_malloc(width * height * sizeof(fftw_complex));
+ cflip_dialog->cf_support = fftw_malloc(width * height * sizeof(unsigned int));
+ cflip_dialog->cf_in = fftw_malloc(width * height * sizeof(fftw_complex));
+ memset(cflip_dialog->cf_in, 0, width * height * sizeof(fftw_complex));
+ cflip_dialog->cf_out = fftw_malloc(width * height * sizeof(fftw_complex));
+
+ cflip_dialog->cf_plan_i2d = fftw_plan_dft_2d(width, height, cflip_dialog->cf_in, cflip_dialog->cf_drive, FFTW_BACKWARD, FFTW_MEASURE | FFTW_PRESERVE_INPUT);
+ cflip_dialog->cf_plan_d2i = fftw_plan_dft_2d(width, height, cflip_dialog->cf_drive, cflip_dialog->cf_in, FFTW_FORWARD, FFTW_MEASURE | FFTW_PRESERVE_INPUT);
+
+ cflip_dialog->cf_reflections = reflections;
+
+ displaywindow_set_realspace(cflip_dialog->cf_out, DWR_GSF);
+ //cflip_dialog->dpsynth = dpsynth_open(reflections, "Current Diffraction Pattern", 1);
+
+ cflip_reset(cflip_dialog);
+
+}
+
+static gint cflip_window_go(GtkWidget *widget, CFlipDialog *cflip_dialog) {
+
+
+ if ( cflip_dialog->running ) {
+ return 0;
+ }
+ cflip_dialog->running = 1;
+ while ( cflip_dialog->running ) {
+ cflip_dialog->n_iterations++;
+ cflip_update_iterations(cflip_dialog);
+ cflip_iteration(cflip_dialog);
+ while ( gtk_events_pending() ) gtk_main_iteration();
+ }
+ return 0;
+
+}
+
+static gint cflip_window_next(GtkWidget *widget, CFlipDialog *cflip_dialog) {
+
+ cflip_dialog->n_iterations++;
+ cflip_update_iterations(cflip_dialog);
+ cflip_iteration(cflip_dialog);
+
+ return 0;
+}
+
+static gint cflip_window_stop(GtkWidget *widget, CFlipDialog *cflip_dialog) {
+
+ cflip_dialog->running = 0;
+ return 0;
+
+}
+
+static gint cflip_window_reset(GtkWidget *widget, CFlipDialog *cflip_dialog) {
+
+ cflip_dialog->running = 0;
+ cflip_reset(cflip_dialog);
+ return 0;
+}
+
+static gint cflip_symmetry_changed(GtkWidget *widget, CFlipDialog *cflip_dialog) {
+
+ cflip_dialog->cf_symmetry = gtk_symmetry_get_symmetry(GTK_SYMMETRY(widget));
+ return 0;
+}
+
+GtkWidget *cfw_tangent;
+void cflip_tangent_toggled(GtkWidget *widget, CFlipDialog *cflip_dialog) {
+
+ cflip_dialog->cf_tangent = 1-cflip_dialog->cf_tangent;
+
+}
+
+void cflip_dialog_open() {
+
+ GtkWidget *cfw_window;
+ GtkWidget *cfw_period;
+ GtkWidget *cfw_period_hbox;
+ GtkWidget *cfw_threshsup_delta;
+ GtkWidget *cfw_threshsup_delta_label;
+ GtkWidget *cfw_threshsup_hbox;
+ GtkWidget *vbox;
+ GtkWidget *hbox;
+ GtkWidget *sym_define;
+ GtkWidget *cfw_support_label;
+ GtkWidget *cfw_process_label;
+ GtkWidget *cfw_period_label;
+ GtkWidget *cfw_gostop_hbox;
+ GtkWidget *cfw_tangent;
+ char *text;
+ CFlipDialog *cflip_dialog;
+
+ cflip_dialog = malloc(sizeof(CFlipDialog));
+
+ if ( cf_window_open ) {
+ return;
+ }
+ cf_window_open = 1;
+
+ cflip_dialog->running = 0;
+ cflip_dialog->cf_tangent = 0;
+ cflip_dialog->n_iterations = 0;
+ cflip_dialog->delta = 0;
+
+ cfw_window = gtk_dialog_new_with_buttons("Charge Flipping", 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(cfw_window)->vbox), GTK_WIDGET(hbox), FALSE, FALSE, 7);
+ gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5);
+
+ cflip_dialog->cf_iterations = gtk_label_new("");
+ text = malloc(40);
+ snprintf(text, 39, "Iterations performed so far: %i", cflip_dialog->n_iterations);
+ gtk_label_set_text(GTK_LABEL(cflip_dialog->cf_iterations), text);
+ free(text);
+ gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cflip_dialog->cf_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(cflip_symmetry_changed), cflip_dialog);
+ cflip_dialog->cf_symmetry = PLANEGROUP_P1;
+
+ cfw_period_label = gtk_label_new("Iterations to symmetrize:");
+ gtk_misc_set_alignment(GTK_MISC(cfw_period_label), 1, 0.5);
+ cfw_period = gtk_entry_new();
+ gtk_entry_set_width_chars(GTK_ENTRY(cfw_period), 5);
+ gtk_entry_set_text(GTK_ENTRY(cfw_period), "1");
+ cfw_period_hbox = gtk_hbox_new(FALSE, 0);
+ gtk_box_pack_start(GTK_BOX(cfw_period_hbox), GTK_WIDGET(cfw_period_label), FALSE, FALSE, 5);
+ gtk_box_pack_start(GTK_BOX(cfw_period_hbox), GTK_WIDGET(cfw_period), FALSE, FALSE, 5);
+ gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cfw_period_hbox), FALSE, FALSE, 5);
+ g_signal_connect(G_OBJECT(cfw_period), "activate", G_CALLBACK(cflip_window_periodchanged), cflip_dialog);
+
+ cfw_support_label = gtk_label_new("");
+ gtk_label_set_markup(GTK_LABEL(cfw_support_label), "<span weight=\"bold\">Threshold</span>");
+ gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cfw_support_label), FALSE, FALSE, 10);
+
+ cfw_tangent = gtk_check_button_new_with_label("Coehlo Constraints");
+ gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cfw_tangent), FALSE, FALSE, 5);
+ if ( cflip_dialog->cf_tangent ) {
+ gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(cfw_tangent), TRUE);
+ } else {
+ gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(cfw_tangent), FALSE);
+ }
+
+ cfw_process_label = gtk_label_new("Maximum unflipped intensity");
+ gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cfw_process_label), FALSE, FALSE, 5);
+
+ cfw_threshsup_hbox = gtk_hbox_new(FALSE, 0);
+ gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cfw_threshsup_hbox), FALSE, FALSE, 5);
+ cfw_threshsup_delta_label = gtk_label_new("δ:");
+ gtk_box_pack_start(GTK_BOX(cfw_threshsup_hbox), GTK_WIDGET(cfw_threshsup_delta_label), FALSE, FALSE, 5);
+ cfw_threshsup_delta = gtk_hscale_new_with_range(-5, 5, 0.1);
+ gtk_scale_set_value_pos(GTK_SCALE(cfw_threshsup_delta), GTK_POS_RIGHT);
+ gtk_box_pack_start(GTK_BOX(cfw_threshsup_hbox), GTK_WIDGET(cfw_threshsup_delta), TRUE, TRUE, 5);
+ gtk_range_set_value(GTK_RANGE(cfw_threshsup_delta), cflip_dialog->delta);
+
+ gtk_widget_set_sensitive(cfw_threshsup_delta, TRUE);
+ gtk_widget_set_sensitive(cfw_threshsup_delta_label, TRUE);
+
+ cfw_gostop_hbox = gtk_hbox_new(FALSE, 0);
+ cflip_dialog->cf_go = gtk_button_new_with_label("Run");
+ gtk_button_set_image(GTK_BUTTON(cflip_dialog->cf_go), gtk_image_new_from_stock(GTK_STOCK_MEDIA_PLAY, GTK_ICON_SIZE_BUTTON));
+ cflip_dialog->cf_next = gtk_button_new_with_label("Next");
+ gtk_button_set_image(GTK_BUTTON(cflip_dialog->cf_next), gtk_image_new_from_stock(GTK_STOCK_MEDIA_NEXT, GTK_ICON_SIZE_BUTTON));
+ cflip_dialog->cf_stop = gtk_button_new_from_stock(GTK_STOCK_MEDIA_STOP);
+ cflip_dialog->cf_reset = gtk_button_new_with_label("Reset");
+ gtk_button_set_image(GTK_BUTTON(cflip_dialog->cf_reset), gtk_image_new_from_stock(GTK_STOCK_MEDIA_PREVIOUS, GTK_ICON_SIZE_BUTTON));
+ gtk_box_pack_start(GTK_BOX(cfw_gostop_hbox), GTK_WIDGET(cflip_dialog->cf_reset), FALSE, FALSE, 5);
+ gtk_box_pack_start(GTK_BOX(cfw_gostop_hbox), GTK_WIDGET(cflip_dialog->cf_stop), FALSE, FALSE, 5);
+ gtk_box_pack_start(GTK_BOX(cfw_gostop_hbox), GTK_WIDGET(cflip_dialog->cf_next), FALSE, FALSE, 5);
+ gtk_box_pack_start(GTK_BOX(cfw_gostop_hbox), GTK_WIDGET(cflip_dialog->cf_go), FALSE, FALSE, 5);
+ gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cfw_gostop_hbox), FALSE, FALSE, 5);
+ g_signal_connect(G_OBJECT(cflip_dialog->cf_go), "clicked", G_CALLBACK(cflip_window_go), cflip_dialog);
+ g_signal_connect(G_OBJECT(cflip_dialog->cf_stop), "clicked", G_CALLBACK(cflip_window_stop), cflip_dialog);
+ g_signal_connect(G_OBJECT(cflip_dialog->cf_next), "clicked", G_CALLBACK(cflip_window_next), cflip_dialog);
+ g_signal_connect(G_OBJECT(cflip_dialog->cf_reset), "clicked", G_CALLBACK(cflip_window_reset), cflip_dialog);
+ g_signal_connect(G_OBJECT(cfw_tangent), "toggled", G_CALLBACK(cflip_tangent_toggled), cflip_dialog);
+ g_signal_connect(G_OBJECT(cfw_window), "response", G_CALLBACK(cflip_window_close), cflip_dialog);
+ g_signal_connect(G_OBJECT(cfw_threshsup_delta), "value-changed", G_CALLBACK(cflip_threshsup_deltachanged), cflip_dialog);
+ g_signal_connect(G_OBJECT(displaywindow_gtkwindow()), "delete-event", G_CALLBACK(cflip_window_stop), cflip_dialog);
+
+ cflip_initialise(cflip_dialog);
+ displaywindow_forceview(DWV_REALSPACE);
+
+ gtk_widget_show_all(cfw_window);
+
+}