aboutsummaryrefslogtreecommitdiff
path: root/src/itrans-stat.c
blob: 7099c8d4fe01adb8d5f844410564cab372ac4005 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
/*
 * itrans-stat.c
 *
 * Peak detection by iterative statistical analysis and processing
 *
 * (c) 2007 Gordon Ball <gfb21@cam.ac.uk>
 *	    Thomas White <taw27@cam.ac.uk>
 *
 *  dtr - Diffraction Tomography Reconstruction
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <math.h>
#include <gsl/gsl_matrix.h>

#include "reflections.h"
#include "control.h"
#include "imagedisplay.h"
#include "utils.h"

/* 
 * Create a gsl_matrix for performing image operations on
 * from a raw image and control context
 * Converting to gsl_matrix because many of the required operations
 * can be done as matrices to save effort
 * Renormalises matrix to 0->1
 */
static gsl_matrix *itrans_peaksearch_stat_createImageMatrix(uint16_t *image, int width, int height) {

	gsl_matrix *raw;
	int i, j;

	raw = gsl_matrix_alloc(width, height);
	for ( i=0; i<raw->size1; i++ ) {
		for ( j=0; j<raw->size2; j++ ) {
			gsl_matrix_set(raw, i, j, (double)image[i+j*width]);
		}
	}
	matrix_renormalise(raw);
	
	return raw;

}

/* Find and return the mean value of the matrix */
static double itrans_peaksearch_stat_image_mean(gsl_matrix *m) {

	int i, j;
	double mean = 0.;

	for ( i=0; i<m->size1; i++ ) {
		for ( j=0; j<m->size2; j++ ) {
			mean += gsl_matrix_get(m,i,j);
		}
	}

	return mean / (m->size1 * m->size2);

}

/*
 * Return the standard deviation, sigma, of the matrix values
 * \sqrt(\sum((x-mean)^2)/n)
 */
static double itrans_peaksearch_stat_image_sigma(gsl_matrix *m, double mean) {

	int i,j;
	double diff2 = 0;

	for ( i=0; i<m->size1; i++ ) {
		for ( j=0; j<m->size2; j++ ) {
			diff2 += (gsl_matrix_get(m,i,j)-mean)*(gsl_matrix_get(m,i,j)-mean);
		}
	}
	
	return sqrt(diff2/(m->size1 * m->size2));

}

/*
 * Filter all pixels with value < mean + k*sigma to 0
 * Set all matching pixels to 1 d
 */
static void itrans_peaksearch_stat_sigma_filter(gsl_matrix *m, double k) {

	double mean,sigma;
	int i,j;

	mean = itrans_peaksearch_stat_image_mean(m);
	sigma = itrans_peaksearch_stat_image_sigma(m,mean);

	for ( i=0; i<m->size1; i++ ) {
		for ( j=0; j<m->size2; j++ ) {
			if (gsl_matrix_get(m,i,j) >= mean+k*sigma) {
				gsl_matrix_set(m, i, j, 1.);
			} else {
				gsl_matrix_set(m, i, j, 0.);
			}
		}
	}
	
}

/*
 * Calculate the mean within a circle centred (x,y) of radius r
 *
 * TODO: Use a mask instead of calculating valid points
 */
static double itrans_peaksearch_stat_circle_mean(gsl_matrix *m, int x, int y, int r, gsl_matrix *mask) {

	double mean = 0.;
	int i, j;
	int n = 0;

	for ( i=x-r; i<=x+r; i++ ) {
		for ( j=y-r; j<=y+r; j++ ) {
			//printf("cm: ij=(%d,%d) mask=(%d,%d)\n",i,j,i-x+r,j-y+r);
			if ( gsl_matrix_get(mask,i-x+r,j-y+r)>0. ) {
				mean += gsl_matrix_get(m, i, j);
				//printf("cm: (%d,%d) mean=%lf val=%lf\n",i,j,mean,gsl_matrix_get(m,i,j));
				n++;
			}
		}
	}
	
	//printf("cm: (%d,%d) summean=%lf n=%d\n",x,y,mean,n);
	return mean/n;
	
}

/*
 * Calculate sigma within a circle centred (x,y) of radius r
 *
 * TODO: Use a mask instead of calculating valid points
 */
static double itrans_peaksearch_stat_circle_sigma(gsl_matrix *m, int x, int y, int r, gsl_matrix *mask, double mean) {

	double diff2 = 0.;
	int i, j;
	int n = 0;

	for ( i=x-r; i<=x+r; i++ ) {
		for ( j=y-r; j<=y+r; j++ ) {
			if ( gsl_matrix_get(mask, i-x+r, j-y+r) > 0 ) {
				diff2 += (gsl_matrix_get(m,i,j)-mean)*(gsl_matrix_get(m,i,j)-mean);
				n++;
			}
		}
	}
	
	return sqrt(diff2/n);

}

/*
 * Calculate a circular mask to save recalculating it every time
 */
static gsl_matrix *itrans_peaksearch_stat_circle_mask(int r) {

	gsl_matrix *m;
	int i,j;

	m = gsl_matrix_calloc(2*r+1, 2*r+1);
	for ( i=0; i<2*r+1; i++ ) {
		for ( j=0; j<2*r+1; j++ ) {
			if ( sqrt((r-i)*(r-i)+(r-j)*(r-j)) <= r ) {
				 gsl_matrix_set(m, i, j, 1.);
			}
			
		}
	}
	
	return m;
	
}

/*
 * Variation on the above filter where instead of using
 * sigma and mean for the whole image, it is found for a local
 * circle of radius r pixels.
 * The central point is also calculated as an approximately gaussian
 * average over local pixels rather than just a single pixel.
 * This takes a long time - 20-30 seconds for a 1024^2 image and r=10,
 * obviously large values of r will make it even slower.
 * The appropriate r value needs to be considered carefully - it needs to
 * be somewhat smaller than the average inter-reflection seperation.
 * 10 pixels worked well for the sapphire.mrc images, but this might
 * not be the case for other images. Images for which this would be very
 * large probably need to be resampled to a lower resolution.
 *
 * TODO: Pass a mask to the ancilliary functions instead of having
 * them calculate several hundred million sqrts
 * 
 * self-referencing problem being dealt with - output being written onto the array before the next point it computed
 * problem carried over from the OO version where a new object was created by each operation
 */
static void itrans_peaksearch_stat_local_sigma_filter(gsl_matrix *m, int r, double k) {

	double mean,sigma;
	double local;
	gsl_matrix *mask;
	gsl_matrix *new;
	int i,j;
	//int interval;

	mask = itrans_peaksearch_stat_circle_mask(r);
	new = gsl_matrix_alloc(m->size1, m->size2);

	//interval = (m->size1-r)/20;
	//printf("lsf: starting loop\n");
	//printf("lsf: ");
	//for (i=r;i<m->size1-r;i++) {
	//	for (j=r;j<m->size2-r;j++) {

	for ( i=0; i<m->size1; i++ ) {
		for ( j=0; j<m->size2; j++ ) {
			if ( ((i >= r) && (i < m->size1-r)) && ((j >= r) && (j < m->size2-r)) ) {
			
				//printf("lsf: evaluating (%d,%d)\n",i,j);
				mean = itrans_peaksearch_stat_circle_mean(m, i, j, r, mask);
				//printf("lsf: mean=%lf",mean);
				sigma = itrans_peaksearch_stat_circle_sigma(m, i, j, r, mask, mean);
				//printf(" sigma=%lf",sigma);
				local = gsl_matrix_get(m, i, j);
				local += gsl_matrix_get(m, i+1, j) + gsl_matrix_get(m, i-1, j) + gsl_matrix_get(m, i, j+1) + gsl_matrix_get(m, i, j-1);
				local += .5*(gsl_matrix_get(m, i+1, j+1) + gsl_matrix_get(m, i-1, j+1) + gsl_matrix_get(m, i+1, j-1) + gsl_matrix_get(m, i-1, j-1));
				local /= 7.;
				//printf(" local=%lf\n",local);
				if ( local > mean+k*sigma ) {
					gsl_matrix_set(new, i, j, 1.);
				} else {
					gsl_matrix_set(new, i, j, 0.);
				}
				
			} else {
				gsl_matrix_set(new, i, j, 0.);
			}
		}
		//if (i % interval == 0) printf(".");
	}
	
	//printf("done\n");
	gsl_matrix_memcpy(m, new);
	gsl_matrix_free(new);

}

/*
 * Apply an arbitary kernel to the image - each point takes the value
 * of the sum of the kernel convolved with the surrounding pixels.
 * The kernel needs to be a (2n+1)^2 array (centred on (n,n)). It needs
 * to be square and should probably be normalised.
 * Don't do daft things like rectangular kernels or kernels larger
 * than the image - this doesn't check such things and will cry.
 * 
 * Also suffers from self-reference problem
 */
static void itrans_peaksearch_stat_apply_kernel(gsl_matrix *m, gsl_matrix *kernel) {

	int size;
	int half;
	gsl_matrix *l;
	gsl_matrix_view lv;
	gsl_matrix *new;
	double val;
	int i, j, x, y;

	size = kernel->size1;
	half = (size-1)/2;
	new = gsl_matrix_calloc(m->size1, m->size2);

	for ( i=0; i<m->size1; i++ ) {
		for ( j=0; j<m->size2; j++ ) {
			if ( ((i >= half) && (i < m->size1-half)) && ((j >= half) && (j < m->size2-half)) ) {
			
				lv = gsl_matrix_submatrix(m, i-half, j-half, size, size);
				l = &lv.matrix;
				val = 0.;
				for ( x=0; x<size; x++ ) {
					for ( y=0; y<size; y++ ) {
						val += gsl_matrix_get(l, x, y)*gsl_matrix_get(kernel, x, y);
					}
				}
				//gsl_matrix_free(l);
				gsl_matrix_set(new,i,j,val);
				
			}
		}
	}
	
	gsl_matrix_memcpy(m,new);
	gsl_matrix_free(new);
	
}

/* Generate the simplist possible kernel - a flat one */
static gsl_matrix *itrans_peaksearch_stat_generate_flat_kernel(int half) {

	gsl_matrix *k;

	k = gsl_matrix_alloc(2*half+1,2*half+1);
	gsl_matrix_set_all(k,1./((2*half+1)*(2*half+1)));

	return k;

}

/* Expands or contracts a gsl_matrix by copying the columns to a new one */
static gsl_matrix *itrans_peaksearch_stat_matrix_expand(gsl_matrix *m, int oldsize, int newsize) {

	gsl_matrix *new;
	int j;
	
	//printf("me: %d->%d\n",oldsize,newsize);
	
	new = gsl_matrix_calloc(2, newsize);
	//printf("me: calloc(2,%d)\n",newsize);
	for ( j=0; j<oldsize; j++) {
		if ( j < newsize ) {
			//printf("me: copying col %d\n",j);
			gsl_matrix_set(new, 0, j, gsl_matrix_get(m, 0, j));
			gsl_matrix_set(new, 1, j, gsl_matrix_get(m, 1, j));
		}
	}

	//printf("me: freeing old matrix\n");
	gsl_matrix_free(m);

	//printf("me: new s1=%d s2=%d\n",new->size1,new->size2);
	//printf("me: m s1=%d s2=%d\n",m->size1,m->size2);
	return new;
	
}

/*
 * Stack-based flood-fill iteration routine 
 * have to return a pointer to com each time because if the matrix size has to be changed then we need to keep
 * track of the location of the resized instance
 */
static gsl_matrix *itrans_peaksearch_stat_do_ff(int i, int j, int* mask, double threshold, gsl_matrix *m, gsl_matrix *com, int *com_n, int *com_size) {

	if ( (i >= 0) && (i < m->size1) ) {
		if ( (j >= 0) && (j < m->size2) ) {
			if ( mask[i+j*m->size1] == 0 ) {
				if ( gsl_matrix_get(m, i, j) > threshold ) {
				
					//printf("dff: found valid point (%d,%d)\n",i,j);
					gsl_matrix_set(com, 0, *com_n, i);
					gsl_matrix_set(com, 1, *com_n, j);
					*com_n = *com_n + 1;
					if ( *com_n == *com_size ) {
						com = itrans_peaksearch_stat_matrix_expand(com, *com_size, *com_size*2);
						*com_size *= 2;
					}
					mask[i+j*m->size1] = 1;
					com = itrans_peaksearch_stat_do_ff(i+1, j, mask,threshold, m, com, com_n, com_size);
					com = itrans_peaksearch_stat_do_ff(i-1, j, mask,threshold, m, com, com_n, com_size);
					com = itrans_peaksearch_stat_do_ff(i, j+1, mask,threshold, m, com, com_n, com_size);
					com = itrans_peaksearch_stat_do_ff(i, j-1, mask,threshold, m, com, com_n, com_size);
					
				}
			}
		}
	}
	
	return com;

}

/*
 * Find points by floodfilling all pixels above threshold
 * Implements a stack-based flood-filling method which may
 * cause stack overflows if the blocks are extremely large -
 * dependent on implementation (default 4MiB stack for unix?)
 * Implements a crude variable-sized array method which hopefully works
 * Returns a gsl_matrix with x coordinates in row 0 and y
 * coordinates in row 1, which should be of the right length
 * Variable count is set to the number of points found
 */
static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, double threshold, int *count) {

	int *mask;
	int size, com_size, i, j, k, n;
	int com_n;
	gsl_matrix *p;
	gsl_matrix *com;
	double com_x, com_y;

	mask = calloc(m->size1*m->size2, sizeof(int));
	size = 32;
	n = 0;
	p = gsl_matrix_calloc(2, size);

	//printf("ff: starting loop\n");
	for ( i=0; i<m->size1; i++ ) {
		for ( j=0; j<m->size2; j++ ) {
			if ( gsl_matrix_get(m, i, j) > threshold ) {
				if ( mask[i+j*m->size1] == 0 ) {
				
					//printf("ff: found starting point (%d,%d)\n",i,j);
					com_size = 32;
					com_n = 0;
					com_x = com_y = 0.;					
					com = gsl_matrix_calloc(2, com_size); //this is going to hold the points found for this location
					//printf("ff: starting floodfill stack\n");
					com = itrans_peaksearch_stat_do_ff(i, j, mask, threshold, m, com, &com_n, &com_size);
					//printf("ff: ended floodfill stack\n");
					for ( k=0; k<com_n; k++ ) {
						com_x += gsl_matrix_get(com, 0, k);
						com_y += gsl_matrix_get(com, 1, k);
					}
					com_x /= com_n;
					com_y /= com_n;
					//printf("ff: point CoM (%lf,%lf)\n",com_x,com_y);
					gsl_matrix_set(p, 0, n, com_x);
					gsl_matrix_set(p, 1, n, com_y);
					n++;
					if ( n == size ) {
						p = itrans_peaksearch_stat_matrix_expand(p, size, size*2);
						size *= 2;
					}
					
				}
			}
		}
	}
	//printf("ff: ending loop, found %d\n",n);

	*count = n;
	//printf("pcheck s1=%d s2=%d\n",p->size1,p->size2);
	p = itrans_peaksearch_stat_matrix_expand(p, size, n);
	//printf("pcheck s1=%d s2=%d\n",p->size1,p->size2);
	
	return p;
	
}

/* Implements the iteration based automatic method
 * returns a gsl_matrix formatted as described in flood-fill */
static gsl_matrix *itrans_peaksearch_stat_iterate(gsl_matrix *m, unsigned int *count) {

	int old;
	int cur;
	double k;
	double mean,sigma;
	gsl_matrix *p;
	gsl_matrix *kernel;
	
	old = m->size1*m->size2;
	
	//printf("Iterate: starting\n");
	//printf("Iterate: generating kernel\n");
	kernel = itrans_peaksearch_stat_generate_flat_kernel(1);
	//printf("Iterate: performing local_sigma_filter\n");
	itrans_peaksearch_stat_local_sigma_filter(m, 10, 1.);
	
	//printf("Iterate: starting loop\n");
	while ( 1 ) {
	
		//printf("Iterate: smoothing");
		itrans_peaksearch_stat_apply_kernel(m, kernel);
		//printf(" (1)");
		itrans_peaksearch_stat_apply_kernel(m, kernel);
		//printf(" (2)\n");
		mean = itrans_peaksearch_stat_image_mean(m);
		sigma = itrans_peaksearch_stat_image_sigma(m,mean);
		//printf("Iterate: mean=%lf sigma=%lf\n",mean,sigma);
		k = (0.5-mean)/sigma;
		//printf("Iterate: applying sigma_filter(%lf)\n",k);
		itrans_peaksearch_stat_sigma_filter(m,k);
		//printf("Iterate: floodfilling\n");
		p = itrans_peaksearch_stat_floodfill(m, 0, &cur);
		//printf("Iterate: %d points found\n", cur);
		if ( old < 1.05*cur ) break;
		old = cur;
		
	}
	
	gsl_matrix_free(kernel); 
	//printf("Iterate: finished\n");
	*count = cur;
	
	return p;

}

unsigned int itrans_peaksearch_stat(ImageRecord imagerecord, ControlContext *ctx) {

	unsigned int n_reflections;	
	gsl_matrix *m;
	gsl_matrix *p;
	int i;
	double px,py;
	uint16_t *image = imagerecord.image;

	m = itrans_peaksearch_stat_createImageMatrix(image, imagerecord.width, imagerecord.height);
	p = itrans_peaksearch_stat_iterate(m, &n_reflections);
	
	for ( i=0; i<n_reflections; i++ ) {
	
		px = gsl_matrix_get(p,0,i);
		py = gsl_matrix_get(p,1,i);
		reflection_add_from_dp(ctx, (px-imagerecord.x_centre), (py-imagerecord.y_centre), imagerecord, 1.0);
		
	}
	gsl_matrix_free(m);
	gsl_matrix_free(p);
	
	return n_reflections;
	
}