Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Oleksii Turkot
CrystFEL
Commits
8a89c359
Commit
8a89c359
authored
Jan 05, 2010
by
Thomas White
Browse files
Introduce process_images
parent
9a0a1c22
Changes
11
Hide whitespace changes
Inline
Side-by-side
.gitignore
View file @
8a89c359
...
...
@@ -13,4 +13,5 @@ src/pattern_sim
src/process_hkl
src/get_hkl
src/hdfsee
src/process_images
*~
Makefile.am
View file @
8a89c359
...
...
@@ -2,5 +2,5 @@ EXTRA_DIST = configure src/cell.h src/hdf5-file.h src/image.h src/relrod.h \
src/utils.h src/diffraction.h src/detector.h src/ewald.h
\
src/sfac.h src/intensities.h src/reflections.h src/list_tmp.h
\
src/statistics.h src/displaywindow.h src/render.h src/hdfsee.h
\
data/displaywindow.ui
data/displaywindow.ui
src/dirax.h
SUBDIRS
=
src data
configure.ac
View file @
8a89c359
...
...
@@ -31,6 +31,6 @@ AC_MSG_WARN([GTK not found. hdfsee will not be built.]))
AM_CONDITIONAL([HAVE_GTK], test x$havegtk = xtrue)
CFLAGS="$CFLAGS $HDF5_CFLAGS $GTK_CFLAGS"
LIBS="$LIBS $HDF5_LIBS -lm -lz -lgsl -lgslcblas $GTK_LIBS -lgthread-2.0"
LIBS="$LIBS $HDF5_LIBS -lm -lz -lgsl -lgslcblas $GTK_LIBS -lgthread-2.0
-lutil
"
AC_OUTPUT(Makefile src/Makefile data/Makefile)
src/Makefile.am
View file @
8a89c359
bin_PROGRAMS
=
pattern_sim process_hkl get_hkl
bin_PROGRAMS
=
pattern_sim process_hkl get_hkl
process_images
if
HAVE_GTK
bin_PROGRAMS
+=
hdfsee
...
...
@@ -16,6 +16,9 @@ process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \
reflections.c
process_hkl_LDADD
=
@LIBS@
process_images_SOURCES
=
process_images.c hdf5-file.c utils.c dirax.c cell.c
process_images_LDADD
=
@LIBS@
if
HAVE_GTK
hdfsee_SOURCES
=
hdfsee.c displaywindow.c render.c hdf5-file.c utils.c
hdfsee_LDADD
=
@LIBS@
...
...
src/cell.c
View file @
8a89c359
...
...
@@ -146,6 +146,36 @@ void cell_set_cartesian(UnitCell *cell,
}
void
cell_set_cartesian_x
(
UnitCell
*
cell
,
double
ax
,
double
bx
,
double
cx
)
{
if
(
!
cell
)
return
;
cell
->
ax
=
ax
;
cell
->
bx
=
bx
;
cell
->
cx
=
cx
;
cell_update_crystallographic
(
cell
);
}
void
cell_set_cartesian_y
(
UnitCell
*
cell
,
double
ay
,
double
by
,
double
cy
)
{
if
(
!
cell
)
return
;
cell
->
ay
=
ay
;
cell
->
by
=
by
;
cell
->
cy
=
cy
;
cell_update_crystallographic
(
cell
);
}
void
cell_set_cartesian_z
(
UnitCell
*
cell
,
double
az
,
double
bz
,
double
cz
)
{
if
(
!
cell
)
return
;
cell
->
az
=
az
;
cell
->
bz
=
bz
;
cell
->
cz
=
cz
;
cell_update_crystallographic
(
cell
);
}
UnitCell
*
cell_new_from_parameters
(
double
a
,
double
b
,
double
c
,
double
alpha
,
double
beta
,
double
gamma
)
{
...
...
src/cell.h
View file @
8a89c359
...
...
@@ -60,6 +60,13 @@ extern void cell_get_cartesian(UnitCell *cell,
double
*
bx
,
double
*
by
,
double
*
bz
,
double
*
cx
,
double
*
cy
,
double
*
cz
);
extern
void
cell_set_cartesian_x
(
UnitCell
*
cell
,
double
ax
,
double
bx
,
double
cx
);
extern
void
cell_set_cartesian_y
(
UnitCell
*
cell
,
double
ay
,
double
by
,
double
cy
);
extern
void
cell_set_cartesian_z
(
UnitCell
*
cell
,
double
az
,
double
bz
,
double
cz
);
extern
void
cell_get_reciprocal
(
UnitCell
*
cell
,
double
*
asx
,
double
*
asy
,
double
*
asz
,
double
*
bsx
,
double
*
bsy
,
double
*
bsz
,
...
...
src/dirax.c
0 → 100644
View file @
8a89c359
/*
* dirax.c
*
* Invoke the DirAx auto-indexing program
*
* (c) 2006-2009 Thomas White <taw@physics.org>
*
* Part of CrystFEL - crystallography with a FEL
*
*/
#ifdef HAVE_CONFIG_H
#include
<config.h>
#endif
#include
<glib.h>
#include
<stdlib.h>
#include
<stdio.h>
#include
<math.h>
#include
<string.h>
#include
<pty.h>
#include
<unistd.h>
#include
<sys/wait.h>
#include
<fcntl.h>
#include
<assert.h>
#include
<sys/ioctl.h>
#include
<termio.h>
#include
<sgtty.h>
#include
"image.h"
#include
"dirax.h"
typedef
enum
{
DIRAX_INPUT_NONE
,
DIRAX_INPUT_LINE
,
DIRAX_INPUT_PROMPT
}
DirAxInputType
;
static
void
dirax_parseline
(
const
char
*
line
,
struct
image
*
image
)
{
int
i
,
rf
;
char
*
copy
;
copy
=
strdup
(
line
);
for
(
i
=
0
;
i
<
strlen
(
copy
);
i
++
)
{
if
(
copy
[
i
]
==
'\r'
)
copy
[
i
]
=
'r'
;
if
(
copy
[
i
]
==
'\n'
)
copy
[
i
]
=
'\0'
;
}
STATUS
(
"DirAx: %s
\n
"
,
copy
);
free
(
copy
);
if
(
strstr
(
line
,
"reflections from file"
)
)
{
ERROR
(
"DirAx can't understand this data."
);
return
;
}
/* Is this the first line of a unit cell specification? */
rf
=
0
;
i
=
0
;
while
(
(
i
<
strlen
(
line
))
&&
((
line
[
i
]
==
'R'
)
||
(
line
[
i
]
==
'D'
)
||
(
line
[
i
]
==
' '
))
)
{
if
(
line
[
i
]
==
'R'
)
rf
=
1
;
if
(
(
line
[
i
]
==
'D'
)
&&
rf
)
{
image
->
dirax_read_cell
=
1
;
if
(
image
->
cell
)
{
free
(
image
->
cell
);
}
image
->
cell
=
cell_new
();
return
;
}
i
++
;
}
/* Parse unit cell vectors as appropriate */
if
(
image
->
dirax_read_cell
==
1
)
{
/* First row of unit cell values */
float
x1
,
x2
,
x3
;
sscanf
(
line
,
"%f %f %f"
,
&
x1
,
&
x2
,
&
x3
);
cell_set_cartesian_x
(
image
->
cell
,
x1
*
1e10
,
x2
*
1e10
,
x3
*
1e10
);
image
->
dirax_read_cell
++
;
return
;
}
else
if
(
image
->
dirax_read_cell
==
2
)
{
/* First row of unit cell values */
float
y1
,
y2
,
y3
;
sscanf
(
line
,
"%f %f %f"
,
&
y1
,
&
y2
,
&
y3
);
cell_set_cartesian_y
(
image
->
cell
,
y1
*
1e10
,
y2
*
1e10
,
y3
*
1e10
);
image
->
dirax_read_cell
++
;
return
;
}
else
if
(
image
->
dirax_read_cell
==
3
)
{
/* First row of unit cell values */
float
z1
,
z2
,
z3
;
sscanf
(
line
,
"%f %f %f"
,
&
z1
,
&
z2
,
&
z3
);
cell_set_cartesian_z
(
image
->
cell
,
z1
*
1e10
,
z2
*
1e10
,
z3
*
1e10
);
STATUS
(
"Read a reciprocal unit cell from DirAx
\n
"
);
/* FIXME: Do something */
image
->
dirax_read_cell
=
0
;
return
;
}
image
->
dirax_read_cell
=
0
;
}
static
void
dirax_sendline
(
const
char
*
line
,
struct
image
*
image
)
{
char
*
copy
;
int
i
;
write
(
image
->
dirax_pty
,
line
,
strlen
(
line
));
copy
=
strdup
(
line
);
for
(
i
=
0
;
i
<
strlen
(
copy
);
i
++
)
{
if
(
copy
[
i
]
==
'\r'
)
copy
[
i
]
=
'\0'
;
if
(
copy
[
i
]
==
'\n'
)
copy
[
i
]
=
'\0'
;
}
STATUS
(
"To DirAx: '%s'
\n
"
,
copy
);
free
(
copy
);
}
static
void
dirax_send_next
(
struct
image
*
image
)
{
switch
(
image
->
dirax_step
)
{
case
1
:
{
dirax_sendline
(
"
\\
echo off
\n
"
,
image
);
image
->
dirax_step
++
;
break
;
}
case
2
:
{
dirax_sendline
(
"read xfel.drx
\n
"
,
image
);
image
->
dirax_step
++
;
break
;
}
case
3
:
{
dirax_sendline
(
"dmax 10
\n
"
,
image
);
image
->
dirax_step
++
;
break
;
}
case
4
:
{
dirax_sendline
(
"indexfit 2
\n
"
,
image
);
image
->
dirax_step
++
;
break
;
}
case
5
:
{
dirax_sendline
(
"levelfit 200
\n
"
,
image
);
image
->
dirax_step
++
;
break
;
}
case
6
:
{
dirax_sendline
(
"go
\n
"
,
image
);
image
->
dirax_step
++
;
break
;
}
case
7
:
{
dirax_sendline
(
"cell
\n
"
,
image
);
image
->
dirax_step
++
;
break
;
}
default:
{
image
->
dirax_step
=
0
;
STATUS
(
"DirAx is idle
\n
"
);
}
}
}
static
gboolean
dirax_readable
(
GIOChannel
*
dirax
,
GIOCondition
condition
,
struct
image
*
image
)
{
int
rval
;
rval
=
read
(
image
->
dirax_pty
,
image
->
dirax_rbuffer
+
image
->
dirax_rbufpos
,
image
->
dirax_rbuflen
-
image
->
dirax_rbufpos
);
if
(
(
rval
==
-
1
)
||
(
rval
==
0
)
)
{
ERROR
(
"Lost connection to DirAx
\n
"
);
waitpid
(
image
->
dirax_pid
,
NULL
,
0
);
g_io_channel_shutdown
(
image
->
dirax
,
FALSE
,
NULL
);
image
->
dirax
=
NULL
;
return
FALSE
;
}
else
{
int
no_string
=
0
;
image
->
dirax_rbufpos
+=
rval
;
assert
(
image
->
dirax_rbufpos
<=
image
->
dirax_rbuflen
);
while
(
(
!
no_string
)
&&
(
image
->
dirax_rbufpos
>
0
)
)
{
int
i
;
int
block_ready
=
0
;
DirAxInputType
type
=
DIRAX_INPUT_NONE
;
/* See if there's a full line in the buffer yet */
for
(
i
=
0
;
i
<
image
->
dirax_rbufpos
-
1
;
i
++
)
{
/* Means the last value looked at is rbufpos-2 */
/* Is there a prompt in the buffer? */
if
(
i
+
7
<=
image
->
dirax_rbufpos
)
{
if
(
(
strncmp
(
image
->
dirax_rbuffer
+
i
,
"Dirax> "
,
7
)
==
0
)
||
(
strncmp
(
image
->
dirax_rbuffer
+
i
,
"PROMPT:"
,
7
)
==
0
)
)
{
block_ready
=
1
;
type
=
DIRAX_INPUT_PROMPT
;
break
;
}
}
if
(
(
image
->
dirax_rbuffer
[
i
]
==
'\r'
)
&&
(
image
->
dirax_rbuffer
[
i
+
1
]
==
'\n'
)
)
{
block_ready
=
1
;
type
=
DIRAX_INPUT_LINE
;
break
;
}
}
if
(
block_ready
)
{
unsigned
int
new_rbuflen
;
unsigned
int
endbit_length
;
switch
(
type
)
{
case
DIRAX_INPUT_LINE
:
{
char
*
block_buffer
=
NULL
;
block_buffer
=
malloc
(
i
+
1
);
memcpy
(
block_buffer
,
image
->
dirax_rbuffer
,
i
);
block_buffer
[
i
]
=
'\0'
;
if
(
block_buffer
[
0
]
==
'\r'
)
{
memmove
(
block_buffer
,
block_buffer
+
1
,
i
);
}
dirax_parseline
(
block_buffer
,
image
);
free
(
block_buffer
);
endbit_length
=
i
+
2
;
break
;
}
case
DIRAX_INPUT_PROMPT
:
{
dirax_send_next
(
image
);
endbit_length
=
i
+
7
;
break
;
}
default
:
{
ERROR
(
" Unrecognised input mode (this never happens!)
\n
"
);
abort
();
}
}
/* Now the block's been parsed, it should be
* forgotten about */
memmove
(
image
->
dirax_rbuffer
,
image
->
dirax_rbuffer
+
endbit_length
,
image
->
dirax_rbuflen
-
endbit_length
);
/* Subtract the number of bytes removed */
image
->
dirax_rbufpos
=
image
->
dirax_rbufpos
-
endbit_length
;
new_rbuflen
=
image
->
dirax_rbuflen
-
endbit_length
;
if
(
new_rbuflen
==
0
)
{
new_rbuflen
=
256
;
}
image
->
dirax_rbuffer
=
realloc
(
image
->
dirax_rbuffer
,
new_rbuflen
);
image
->
dirax_rbuflen
=
new_rbuflen
;
}
else
{
if
(
image
->
dirax_rbufpos
==
image
->
dirax_rbuflen
)
{
/* More buffer space is needed */
image
->
dirax_rbuffer
=
realloc
(
image
->
dirax_rbuffer
,
image
->
dirax_rbuflen
+
256
);
image
->
dirax_rbuflen
=
image
->
dirax_rbuflen
+
256
;
/* The new space gets used at the next
* read, shortly... */
}
no_string
=
1
;
}
}
}
return
TRUE
;
}
static
int
map_position
(
struct
image
*
image
,
double
x
,
double
y
,
double
*
rx
,
double
*
ry
,
double
*
rz
)
{
/* "Input" space */
double
d
;
/* Angular description of reflection */
double
theta
,
psi
,
k
;
x
-=
image
->
x_centre
;
y
-=
image
->
y_centre
;
k
=
1
.
0
/
image
->
lambda
;
if
(
image
->
fmode
==
FORMULATION_CLEN
)
{
/* Convert pixels to metres */
x
/=
image
->
resolution
;
y
/=
image
->
resolution
;
x
=
x
*
k
/
image
->
camera_len
;
y
=
y
*
k
/
image
->
camera_len
;
x
/=
image
->
resolution
;
y
/=
image
->
resolution
;
/* Convert pixels to metres */
d
=
sqrt
((
x
*
x
)
+
(
y
*
y
));
theta
=
atan2
(
d
,
image
->
camera_len
);
}
else
if
(
image
->
fmode
==
FORMULATION_PIXELSIZE
)
{
/* Convert pixels to metres^-1 */
x
=
x
*
image
->
pixel_size
;
y
=
y
*
image
->
pixel_size
;
x
*=
image
->
pixel_size
;
y
*=
image
->
pixel_size
;
/* Convert pixels to metres^-1 */
d
=
sqrt
((
x
*
x
)
+
(
y
*
y
));
theta
=
atan2
(
d
,
k
);
}
else
{
ERROR
(
"Unrecognised formulation mode in mapping_scale.
\n
"
);
return
-
1
;
}
psi
=
atan2
(
y
,
x
);
*
rx
=
k
*
sin
(
theta
)
*
cos
(
psi
);
*
ry
=
k
*
sin
(
theta
)
*
sin
(
psi
);
*
rz
=
k
-
k
*
cos
(
theta
);
return
0
;
}
void
index_pattern
(
struct
image
*
image
)
{
FILE
*
fh
;
unsigned
int
opts
;
int
saved_stderr
;
int
x
,
y
;
GMainLoop
*
ml
;
fh
=
fopen
(
"xfel.drx"
,
"w"
);
if
(
!
fh
)
{
ERROR
(
"Couldn't open temporary file xfel.drx
\n
"
);
return
;
}
fprintf
(
fh
,
"%f
\n
"
,
0
.
5
);
/* Lie about the wavelength. */
for
(
x
=
0
;
x
<
image
->
width
;
x
++
)
{
for
(
y
=
0
;
y
<
image
->
height
;
y
++
)
{
int
val
;
val
=
image
->
data
[
y
+
image
->
height
*
(
image
->
width
-
1
-
x
)];
if
(
val
>
1000
)
{
double
rx
,
ry
,
rz
;
/* Map and record reflection */
map_position
(
image
,
x
,
y
,
&
rx
,
&
ry
,
&
rz
);
fprintf
(
fh
,
"%10f %10f %10f %8f
\n
"
,
rx
/
1e10
,
ry
/
1e10
,
rz
/
1e10
,
1
.
0
);
}
}
}
fclose
(
fh
);
saved_stderr
=
dup
(
STDERR_FILENO
);
image
->
dirax_pid
=
forkpty
(
&
image
->
dirax_pty
,
NULL
,
NULL
,
NULL
);
if
(
image
->
dirax_pid
==
-
1
)
{
ERROR
(
"Failed to fork for DirAx
\n
"
);
return
;
}
if
(
image
->
dirax_pid
==
0
)
{
/* Child process: invoke DirAx */
struct
termios
t
;
/* Turn echo off */
tcgetattr
(
STDIN_FILENO
,
&
t
);
t
.
c_lflag
&=
~
(
ECHO
|
ECHOE
|
ECHOK
|
ECHONL
);
tcsetattr
(
STDIN_FILENO
,
TCSANOW
,
&
t
);
/* Reconnect stderr */
dup2
(
saved_stderr
,
STDERR_FILENO
);
execlp
(
"dirax"
,
""
,
(
char
*
)
NULL
);
ERROR
(
"Failed to invoke DirAx.
\n
"
);
_exit
(
0
);
}
image
->
dirax_rbuffer
=
malloc
(
256
);
image
->
dirax_rbuflen
=
256
;
image
->
dirax_rbufpos
=
0
;
/* Set non-blocking */
opts
=
fcntl
(
image
->
dirax_pty
,
F_GETFL
);
fcntl
(
image
->
dirax_pty
,
F_SETFL
,
opts
|
O_NONBLOCK
);
image
->
dirax_step
=
1
;
/* This starts the "initialisation" procedure */
image
->
dirax_read_cell
=
0
;
image
->
dirax
=
g_io_channel_unix_new
(
image
->
dirax_pty
);
g_io_add_watch
(
image
->
dirax
,
G_IO_IN
|
G_IO_HUP
,
(
GIOFunc
)
dirax_readable
,
image
);
ml
=
g_main_loop_new
(
NULL
,
FALSE
);
g_main_loop_run
(
ml
);
return
;
}
src/dirax.h
0 → 100644
View file @
8a89c359
/*
* dirax.h
*
* Invoke the DirAx auto-indexing program
*
* (c) 2006-2009 Thomas White <taw@physics.org>
*
* Part of CrystFEL - crystallography with a FEL
*
*/
#ifndef DIRAX_H
#define DIRAX_H
#ifdef HAVE_CONFIG_H
#include
<config.h>
#endif
extern
void
index_pattern
(
struct
image
*
image
);
#endif
/* DIRAX_H */
src/hdf5-file.c
View file @
8a89c359
...
...
@@ -246,8 +246,14 @@ int hdf5_read(struct hdfile *f, struct image *image)
image
->
data
=
buf
;
image
->
height
=
f
->
nx
;
image
->
width
=
f
->
ny
;
/* FIXME: The following are basically made up... */
image
->
x_centre
=
image
->
width
/
2
;
image
->
y_centre
=
image
->
height
/
2
;
image
->
lambda
=
ph_en_to_lambda
(
J_to_eV
(
1793
));
image
->
fmode
=
FORMULATION_CLEN
;
image
->
camera_len
=
75.0e-3
;
/* 75 mm camera length */
image
->
resolution
=
13333
.
3
;
/* 75 micron pixel size */
return
0
;
}
...
...
src/image.h
View file @
8a89c359
...
...
@@ -19,8 +19,10 @@
#include
<stdint.h>
#include
<complex.h>
#include
<glib.h>
#include
"utils.h"
#include
"cell.h"
/* How is the scaling of the image described? */
...
...
@@ -92,6 +94,19 @@ struct image {
ImageFeatureList
*
features
;
/* "Experimental" features */
ImageFeatureList
*
rflist
;
/* "Predicted" features */
UnitCell
*
cell
;
/* DirAx auto-indexing low-level stuff */
GIOChannel
*
dirax
;
int
dirax_pty
;
pid_t
dirax_pid
;
char
*
dirax_rbuffer
;
int
dirax_rbufpos
;
int
dirax_rbuflen
;
/* DirAx auto-indexing high-level stuff */
int
dirax_step
;
int
dirax_read_cell
;
};
/* An opaque type representing a list of images */
...
...
src/process_images.c
0 → 100644
View file @
8a89c359
/*
* process_images.c
*
* Find hits, index patterns, output hkl+intensity etc.
*
* (c) 2006-2009 Thomas White <taw@physics.org>
*
* Part of CrystFEL - crystallography with a FEL
*
*/