pattern_sim: Remove image->data (CPU only)

2 files changed, 89 insertions, 57 deletions

diff --git a/src/diffraction.c b/src/diffraction.c
index 5a936809..93674f52 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -361,31 +361,30 @@ static double molecule_factor(const double *intensities, const double *phases,
 }
 
 
-
-static void diffraction_at_k(struct image *image, const double *intensities,
-                             const double *phases, const unsigned char *flags,
-                             UnitCell *cell, GradientMethod m,
-                             const SymOpList *sym, double k,
-                             double ax, double ay, double az,
-                             double bx, double by, double bz,
-                             double cx, double cy, double cz,
-                             double *lut_a, double *lut_b, double *lut_c,
-                             double weight)
+static void diffraction_panel(struct image *image, const double *intensities,
+                              const double *phases, const unsigned char *flags,
+                              UnitCell *cell, GradientMethod m,
+                              const SymOpList *sym, double k,
+                              double ax, double ay, double az,
+                              double bx, double by, double bz,
+                              double cx, double cy, double cz,
+                              double *lut_a, double *lut_b, double *lut_c,
+                              int pn, double weight)
 {
-	unsigned int fs, ss;
+	int fs, ss;
 	const int nxs = 4;
 	const int nys = 4;
+	struct panel *p = &image->det->panels[pn];
 
 	weight /= nxs*nys;
 
-	for ( fs=0; fs<image->width; fs++ ) {
-	for ( ss=0; ss<image->height; ss++ ) {
+	for ( ss=0; ss<p->h; ss++ ) {
+	for ( fs=0; fs<p->w; fs++ ) {
 
 		int idx;
 		double f_lattice, I_lattice;
 		double I_molecule;
 		struct rvec q;
-		double twotheta;
 		int xs, ys;
 		float xo, yo;
 
@@ -395,7 +394,7 @@ static void diffraction_at_k(struct image *image, const double *intensities,
 			xo = (1.0/nxs) * xs;
 			yo = (1.0/nys) * ys;
 
-			q = get_q(image, fs+xo, ss+yo, &twotheta, k);
+			q = get_q_for_panel(p, fs+xo, ss+yo, NULL, k);
 
 			f_lattice = lattice_factor(q, ax, ay, az,
 					           bx, by, bz,
@@ -411,14 +410,33 @@ static void diffraction_at_k(struct image *image, const double *intensities,
 
 			I_lattice = pow(f_lattice, 2.0);
 
-			idx = fs + image->width*ss;
-			image->data[idx] += I_lattice * I_molecule * weight;
-			image->twotheta[idx] = twotheta;
+			idx = fs + p->w*ss;
+			image->dp[pn][idx] += I_lattice * I_molecule * weight;
 
 		}
 		}
 	}
-	progress_bar(fs, image->width-1, "Calculating diffraction");
+	progress_bar(ss, p->h-1, "Calculating diffraction");
+	}
+}
+
+
+static void diffraction_at_k(struct image *image, const double *intensities,
+                             const double *phases, const unsigned char *flags,
+                             UnitCell *cell, GradientMethod m,
+                             const SymOpList *sym, double k,
+                             double ax, double ay, double az,
+                             double bx, double by, double bz,
+                             double cx, double cy, double cz,
+                             double *lut_a, double *lut_b, double *lut_c,
+                             double weight)
+{
+	int i;
+
+	for ( i=0; i<image->det->n_panels; i++ ) {
+		diffraction_panel(image, intensities, phases, flags, cell, m,
+		                  sym, k, ax, ay, az, bx, by, bz, cx, cy, cz,
+		                  lut_a, lut_b, lut_c, i, weight);
 	}
 }
 
@@ -710,12 +728,6 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 
 	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 
-	/* Allocate (and zero) the "diffraction array" */
-	image->data = calloc(image->width * image->height, sizeof(float));
-
-	/* Needed later for Lorentz calculation */
-	image->twotheta = malloc(image->width * image->height * sizeof(double));
-
 	lut_a = get_sinc_lut(na, no_fringes);
 	lut_b = get_sinc_lut(nb, no_fringes);
 	lut_c = get_sinc_lut(nc, no_fringes);
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index bd97ead7..f724057e 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -249,8 +249,7 @@ int main(int argc, char *argv[])
 	int c;
 	struct image image;
 	struct gpu_context *gctx = NULL;
-	struct image *powder;
-	float *powder_data;
+	struct image powder;
 	char *intfile = NULL;
 	double *intensities;
 	char *rval;
@@ -293,6 +292,7 @@ int main(int argc, char *argv[])
 	double bandwidth = 0.01;
 	double photon_energy = 9000.0;
 	struct beam_params beam;
+	int i;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -667,7 +667,21 @@ int main(int argc, char *argv[])
 	/* Initialise stuff */
 	image.filename = NULL;
 	image.features = NULL;
-	image.flags = NULL;
+	image.bad = NULL;
+	image.dp = malloc(image.det->n_panels*sizeof(float *));
+	if ( image.dp == NULL ) {
+		ERROR("Failed to allocate data\n");
+		return 1;
+	}
+	for ( i=0; i<image.det->n_panels; i++ ) {
+		image.dp[i] = calloc(image.det->panels[i].w *
+		                     image.det->panels[i].h,
+		                     sizeof(float));
+		if ( image.dp[i] == NULL ) {
+			ERROR("Failed to allocate data\n");
+			return 1;
+		}
+	}
 
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 	if ( config_random ) {
@@ -679,12 +693,23 @@ int main(int argc, char *argv[])
 		gsl_rng_set(rng, seed);
 	}
 
-	powder = calloc(1, sizeof(struct image));
-	powder->width = image.width;
-	powder->height = image.height;
-	powder->det = image.det;
-	powder_data = calloc(image.width*image.height, sizeof(float));
-	powder->data = powder_data;
+	powder.width = image.width;
+	powder.height = image.height;
+	powder.det = image.det;
+	powder.dp = malloc(image.det->n_panels*sizeof(float *));
+	if ( powder.dp == NULL ) {
+		ERROR("Failed to allocate powder data\n");
+		return 1;
+	}
+	for ( i=0; i<image.det->n_panels; i++ ) {
+		powder.dp[i] = calloc(image.det->panels[i].w *
+		                      image.det->panels[i].h,
+		                      sizeof(float));
+		if ( powder.dp[i] == NULL ) {
+			ERROR("Failed to allocate powder data\n");
+			return 1;
+		}
+	}
 
 	/* Splurge a few useful numbers */
 	STATUS("Simulation parameters:\n");
@@ -732,6 +757,7 @@ int main(int argc, char *argv[])
 		int na, nb, nc;
 		double a, b, c, d;
 		UnitCell *cell;
+		int err = 0;
 
 		if ( random_size ) {
 
@@ -825,10 +851,6 @@ int main(int argc, char *argv[])
 
 		}
 
-		/* Ensure no residual information */
-		image.data = NULL;
-		image.twotheta = NULL;
-
 		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
 		STATUS("Particle size = %i x %i x %i"
 		       " ( = %5.2f x %5.2f x %5.2f nm)\n",
@@ -836,8 +858,6 @@ int main(int argc, char *argv[])
 
 		if ( config_gpu ) {
 
-			int err;
-
 			if ( gctx == NULL ) {
 				gctx = setup_gpu(config_nosfac,
 				                 intensities, flags, sym_str,
@@ -845,13 +865,12 @@ int main(int argc, char *argv[])
 			}
 			err = get_diffraction_gpu(gctx, &image, na, nb, nc,
 			                          cell, no_fringes);
-			if ( err ) image.data = NULL;
 
 		} else {
 			get_diffraction(&image, na, nb, nc, intensities, phases,
 			                flags, cell, grad, sym, no_fringes);
 		}
-		if ( image.data == NULL ) {
+		if ( err ) {
 			ERROR("Diffraction calculation failed.\n");
 			done = 1;
 			goto skip;
@@ -862,18 +881,23 @@ int main(int argc, char *argv[])
 
 		if ( powder_fn != NULL ) {
 
-			int x, y, w;
+			int pn;
 
-			w = image.width;
+			for ( pn=0; pn<image.det->n_panels; pn++ ) {
+
+				size_t w, i;
+
+				w = image.det->panels[pn].w
+				  * image.det->panels[pn].h;
+
+				for ( i=0; i<w; i++ ) {
+					powder.dp[pn][i] += (double)image.dp[pn][i];
+				}
 
-			for ( x=0; x<image.width; x++ ) {
-			for ( y=0; y<image.height; y++ ) {
-				powder->data[x+w*y] += (double)image.data[x+w*y];
-			}
 			}
 
 			if ( !(ndone % 10) ) {
-				hdf5_write_image(powder_fn, powder, NULL);
+				hdf5_write_image(powder_fn, &powder, NULL);
 			}
 		}
 
@@ -895,10 +919,6 @@ int main(int argc, char *argv[])
 
 		}
 
-		/* Clean up */
-		free(image.data);
-		free(image.twotheta);
-
 		cell_free(cell);
 
 skip:
@@ -909,18 +929,18 @@ skip:
 	} while ( !done );
 
 	if ( powder_fn != NULL ) {
-		hdf5_write_image(powder_fn, powder, NULL);
+		hdf5_write_image(powder_fn, &powder, NULL);
 	}
 
 	if ( gctx != NULL ) {
 		cleanup_gpu(gctx);
 	}
 
-	free(image.det->panels);
+	for ( i=0; i<image.det->n_panels; i++ ) {
+		free(powder.dp[i]);
+	}
 	free(intfile);
 	free(image.det);
-	free(powder->data);
-	free(powder);
 	cell_free(input_cell);
 	free(intensities);
 	free(outfile);