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
|
/*
* mapping.c
*
* 3D Mapping
*
* (c) 2007 Thomas White <taw27@cam.ac.uk>
*
* dtr - Diffraction Tomography Reconstruction
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "reflections.h"
#include "control.h"
#include "imagedisplay.h"
#include "itrans.h"
ReflectionContext *mapping_create(ControlContext *ctx) {
ReflectionContext *rctx;
ImageDisplay *sum_id;
int16_t max_val;
int max_x, max_y;
int twidth, theight;
int x, y, i;
int16_t *image_total;
/* Create reflection context */
rctx = reflection_init();
/* Find centre of image stack
* Determine maximum size of image to accommodate, and allocate memory */
twidth = 0; theight = 0;
for ( i=0; i<ctx->n_images; i++ ) {
if ( ctx->images[i].width > twidth ) twidth = ctx->images[i].width;
if ( ctx->images[i].height > theight ) theight = ctx->images[i].height;
}
image_total = malloc(twidth * theight * sizeof(int16_t));
memset(image_total, 0, twidth * theight * sizeof(int16_t));
/* Add the image stack together */
for ( i=0; i<ctx->n_images; i++ ) {
int xoffs, yoffs;
xoffs = (twidth - ctx->images[i].width)/2;
yoffs = (theight - ctx->images[i].height)/2;
for ( y=0; y<theight; y++ ) {
for ( x=0; x<twidth; x++ ) {
image_total[x + twidth*y] += ctx->images[i].image[x + twidth*y]/ctx->n_images;
}
}
}
/* Locate the highest point */
max_val = 0; max_x = 0; max_y = 0;
for ( y=0; y<theight; y++ ) {
for ( x=0; x<twidth; x++ ) {
if ( image_total[x + twidth*y] > max_val ) {
max_val = image_total[x + twidth*y];
max_x = x; max_y = y;
}
}
}
ctx->x_centre = max_x;
ctx->y_centre = max_y;
/* Display */
sum_id = imagedisplay_open(image_total, twidth, theight, "Sum of All Images");
imagedisplay_mark_point(sum_id, ctx->x_centre, ctx->y_centre);
imagedisplay_add_tilt_axis(sum_id, ctx, ctx->omega);
/* Pass all images through itrans
* (let itrans add the reflections to rctx for now) */
ctx->reflectionctx = rctx;
for ( i=0; i<ctx->n_images; i++ ) {
ctx->images[i].x_centre = ctx->x_centre;
ctx->images[i].y_centre = ctx->y_centre;
itrans_process_image(ctx->images[i].image, ctx, ctx->images[i].tilt);
}
return rctx;
}
|
