aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Makefile.am3
-rw-r--r--src/Makefile.in8
-rw-r--r--src/cubeit.c304
3 files changed, 296 insertions, 19 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index 007aff1c..58f2ff95 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -60,7 +60,8 @@ facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
image.c diffraction.c sfac.c geometry.c
facetron_LDADD = @LIBS@
-cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c
+cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c diffraction.c \
+ sfac.c render.c filters.c image.c
cubeit_LDADD = @LIBS@
INCLUDES = "-I$(top_srcdir)/data"
diff --git a/src/Makefile.in b/src/Makefile.in
index 4fe63294..0a1db948 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -65,7 +65,9 @@ am_compare_hkl_OBJECTS = compare_hkl.$(OBJEXT) sfac.$(OBJEXT) \
compare_hkl_OBJECTS = $(am_compare_hkl_OBJECTS)
compare_hkl_DEPENDENCIES =
am_cubeit_OBJECTS = cubeit.$(OBJEXT) cell.$(OBJEXT) \
- hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT)
+ hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT) \
+ diffraction.$(OBJEXT) sfac.$(OBJEXT) render.$(OBJEXT) \
+ filters.$(OBJEXT) image.$(OBJEXT)
cubeit_OBJECTS = $(am_cubeit_OBJECTS)
cubeit_DEPENDENCIES =
am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) \
@@ -291,7 +293,9 @@ facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
image.c diffraction.c sfac.c geometry.c
facetron_LDADD = @LIBS@
-cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c
+cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c diffraction.c \
+ sfac.c render.c filters.c image.c
+
cubeit_LDADD = @LIBS@
INCLUDES = "-I$(top_srcdir)/data"
all: all-am
diff --git a/src/cubeit.c b/src/cubeit.c
index ba7472aa..d08d30aa 100644
--- a/src/cubeit.c
+++ b/src/cubeit.c
@@ -21,10 +21,14 @@
#include <unistd.h>
#include <getopt.h>
#include <errno.h>
+#include <assert.h>
+#include <png.h>
#include "image.h"
#include "cell.h"
#include "hdf5-file.h"
+#include "diffraction.h"
+#include "render.h"
#define MAX_HITS (1024)
@@ -40,6 +44,10 @@ static void show_help(const char *s)
"\n"
" -i, --input=<filename> Specify the name of the input stream.\n"
" Can be '-' for stdin.\n"
+" -g. --geometry=<file> Get detector geometry from file.\n"
+" -x, --prefix=<p> Prefix filenames from input file with <p>.\n"
+" --basename Remove the directory parts of the filenames.\n"
+" --no-check-prefix Don't attempt to correct the --prefix.\n"
);
}
@@ -113,6 +121,180 @@ static int find_chunk(FILE *fh, UnitCell **cell, char **filename)
}
+static void interpolate_linear(double *vals, double v,
+ int xs, int ys, int zs,
+ int xv, int yv, double zv)
+{
+ int k;
+ double val1, val2;
+ float f, c;
+
+ c = (zv+0.5)*(float)zs;
+ k = (int)c;
+ f = c - (float)k;
+ assert(f >= 0.0);
+ assert(k < zs);
+
+ val1 = v*(1.0-f);
+ val2 = v*f;
+
+ vals[xs*ys*k + xs*yv + xv] = val1;
+ vals[xs*ys*(k+1) + xs*yv + xv] = val2;
+}
+
+
+static void interpolate_bilinear(double *vals, double v,
+ int xs, int ys, int zs,
+ int xv, double yv, double zv)
+{
+ int k;
+ double val1, val2;
+ float f, c;
+
+ c = (yv+0.5)*(float)ys;
+ k = (int)c;
+ f = c - (float)k;
+ assert(f >= 0.0);
+ assert(k < ys);
+
+ val1 = v*(1.0-f);
+ val2 = v*f;
+
+ interpolate_linear(vals, val1, xs, ys, zs, xv, k, zv);
+ interpolate_linear(vals, val2, xs, ys, zs, xv, k+1, zv);
+}
+
+
+static void interpolate_onto_grid(double *vals, double v,
+ int xs, int ys, int zs,
+ double xv, double yv, double zv)
+{
+ int k;
+ double val1, val2;
+ float f, c;
+
+ c = (xv+0.5)*(float)xs;
+ k = (int)c;
+ f = c - (float)k;
+ assert(f >= 0.0);
+ assert(k < xs);
+
+ val1 = v*(1.0-f);
+ val2 = v*f;
+
+ interpolate_bilinear(vals, val1, xs, ys, zs, k, yv, zv);
+ interpolate_bilinear(vals, val2, xs, ys, zs, k+1, yv, zv);
+}
+
+
+static void write_slice(const char *filename, double *vals, int z,
+ int xs, int ys, int zs)
+{
+#ifdef HAVE_LIBPNG
+ FILE *fh;
+ png_structp png_ptr;
+ png_infop info_ptr;
+ png_bytep *row_pointers;
+ int x, y;
+ float max = 0.0;
+ int w, h;
+
+ w = xs;
+ h = ys;
+
+ if ( max <= 6 ) { max = 10; }
+
+ for ( y=0; y<h; y++ ) {
+ for ( x=0; x<w; x++ ) {
+
+ float val;
+
+ val = vals[xs*ys*z + xs*y + x];
+
+ if ( val > max ) max = val;
+
+ }
+ }
+
+ fh = fopen(filename, "wb");
+ if ( !fh ) {
+ ERROR("Couldn't open output file.\n");
+ return;
+ }
+ png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING,
+ NULL, NULL, NULL);
+ if ( !png_ptr ) {
+ ERROR("Couldn't create PNG write structure.\n");
+ fclose(fh);
+ return;
+ }
+ info_ptr = png_create_info_struct(png_ptr);
+ if ( !info_ptr ) {
+ png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
+ ERROR("Couldn't create PNG info structure.\n");
+ fclose(fh);
+ return;
+ }
+ if ( setjmp(png_jmpbuf(png_ptr)) ) {
+ png_destroy_write_struct(&png_ptr, &info_ptr);
+ fclose(fh);
+ ERROR("PNG write failed.\n");
+ return;
+ }
+ png_init_io(png_ptr, fh);
+
+ png_set_IHDR(png_ptr, info_ptr, w, h, 8,
+ PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
+ PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
+
+ row_pointers = malloc(h*sizeof(png_bytep *));
+
+ /* Write the image data */
+ max /= 1.0;
+ if ( max <= 6 ) { max = 10; }
+
+ for ( y=0; y<h; y++ ) {
+
+ row_pointers[y] = malloc(w*3);
+
+ for ( x=0; x<w; x++ ) {
+
+ float r, g, b;
+ float val;
+
+ val = vals[xs*ys*z + xs*y + x];
+
+ render_scale(val, max, SCALE_COLOUR, &r, &g, &b);
+ row_pointers[y][3*x] = (png_byte)255*r;
+ row_pointers[y][3*x+1] = (png_byte)255*g;
+ row_pointers[y][3*x+2] = (png_byte)255*b;
+
+ }
+ }
+
+ for ( y=0; y<h/2+1; y++ ) {
+ png_bytep scratch;
+ scratch = row_pointers[y];
+ row_pointers[y] = row_pointers[h-y-1];
+ row_pointers[h-y-1] = scratch;
+ }
+
+ png_set_rows(png_ptr, info_ptr, row_pointers);
+ png_write_png(png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL);
+
+ png_destroy_write_struct(&png_ptr, &info_ptr);
+ for ( y=0; y<h; y++ ) {
+ free(row_pointers[y]);
+ }
+ free(row_pointers);
+ fclose(fh);
+
+#else
+ STATUS("No PNG support.\n");
+#endif
+}
+
+
int main(int argc, char *argv[])
{
int c;
@@ -122,16 +304,31 @@ int main(int argc, char *argv[])
char *filename;
unsigned int angles[180];
int i;
+ char *prefix = NULL;
+ char *geomfile = NULL;
+ int config_basename = 0;
+ int config_checkprefix = 1;
+ struct detector *det;
+ int config_angles = 0;
+ const signed int ht = 3;
+ const signed int kt = 4;
+ const signed int lt = 5;
+ const int gs = 16;
+ double *vals;
/* Long options */
const struct option longopts[] = {
{"help", 0, NULL, 'h'},
{"input", 1, NULL, 'i'},
+ {"geometry", 1, NULL, 'g'},
+ {"prefix", 1, NULL, 'x'},
+ {"basename", 0, &config_basename, 1},
+ {"no-check-prefix", 0, &config_checkprefix, 0},
{0, 0, NULL, 0}
};
/* Short options */
- while ((c = getopt_long(argc, argv, "hi:", longopts, NULL)) != -1) {
+ while ((c = getopt_long(argc, argv, "hi:g:x:", longopts, NULL)) != -1) {
switch (c) {
case 'h' :
@@ -142,6 +339,14 @@ int main(int argc, char *argv[])
infile = strdup(optarg);
break;
+ case 'g' :
+ geomfile = strdup(optarg);
+ break;
+
+ case 'x' :
+ prefix = strdup(optarg);
+ break;
+
case 0 :
break;
@@ -162,23 +367,42 @@ int main(int argc, char *argv[])
return ENOENT;
}
+ if ( prefix == NULL ) {
+ prefix = strdup("");
+ } else {
+ if ( config_checkprefix ) {
+ prefix = check_prefix(prefix);
+ }
+ }
+
+ det = get_detector_geometry(geomfile);
+ if ( det == NULL ) {
+ ERROR("Failed to read detector geometry from '%s'\n", geomfile);
+ return 1;
+ }
+ free(geomfile);
+
/* Initialise histogram */
for ( i=0; i<180; i++ ) angles[i] = 0;
+ /* Initialise shape transform array */
+ vals = calloc(gs*gs*gs, sizeof(double));
+
/* Loop over all successfully indexed patterns */
while ( find_chunk(fh, &cell, &filename) == 0 ) {
struct image image;
struct hdfile *hdfile;
double ang;
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
unsigned int bin;
+ struct rvec q;
+ int x, y;
- //STATUS("Processing '%s'\n", filename);
+ STATUS("Processing '%s'\n", filename);
- #if 0
hdfile = hdfile_open(filename);
if ( hdfile == NULL ) {
return ENOENT;
@@ -188,20 +412,68 @@ int main(int argc, char *argv[])
}
hdf5_read(hdfile, &image, 0);
- #endif
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
- ang = angle_between(csx, csy, csz, 0.0, 0.0, 1.0); /* 0->pi */
- ang = rad2deg(ang); /* 0->180 deg */
- bin = rint(ang);
- angles[bin]++;
+ cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by,
+ &bz, &cx, &cy, &cz);
+
+ for ( x=0; x<image.width; x++ ) {
+ for ( y=0; y<image.height; y++ ) {
+
+ double hd, kd, ld;
+ signed int h, k, l;
+ double dh, dk, dl;
+
+ q = get_q(&image, x, y, 1, NULL, 1.0/image.lambda);
+
+ hd = q.u * ax + q.v * ay + q.w * az;
+ kd = q.u * bx + q.v * by + q.w * bz;
+ ld = q.u * cx + q.v * cy + q.w * cz;
+
+ h = (signed int)rint(hd);
+ k = (signed int)rint(kd);
+ l = (signed int)rint(ld);
+
+ if ( !((h==ht) && (k==kt) && (l==lt)) ) continue;
+
+ dh = hd - h;
+ dk = kd - k;
+ dl = ld - l;
+
+ double v = image.data[x+image.width*y];
+
+ interpolate_onto_grid(vals, v, gs, gs, gs, dh, dk, dl);
+
+ }
+ }
+
+ if ( config_angles ) {
+
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ ang = angle_between(csx, csy, csz, 0.0, 0.0, 1.0);
+ ang = rad2deg(ang); /* 0->180 deg */
+ bin = rint(ang);
+ angles[bin]++;
+
+ }
}
- for ( i=0; i<180; i++ ) {
- STATUS("%i %i\n", i, angles[i]);
+ for ( i=0; i<gs; i++ ) {
+ char line[64];
+ snprintf(line, 63, "slice-%i.png", i);
+ write_slice(line, vals, i, gs, gs, gs);
+ }
+
+ if ( config_angles ) {
+ for ( i=0; i<180; i++ ) {
+ STATUS("%i %i\n", i, angles[i]);
+ }
}
return 0;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 17:52:54 +0000