HOWTO: Support for New Detectors to dxtbx
-----------------------------------------

A new framework has been implemented within dxtbx to make it more straight-
forward to add support for new detector types and beamlines. In essence all
that is needed is to implement a Python class which extends the Format class
to add some specific details about this detector and the associated beamline /
experimental environment.

In particular there are two groups of things which need to be
implemented - a static method named "understand" which will take a
look at the image and return True if it understands it, and a bunch of
class methods which need to override the construction of the
description.

"understand" Static Method
--------------------------

This method is the key to how the whole framework operates - you write code
which looks at the image to decide whether it is right for this class. If it is
not you must return False - i.e. if you are making a custom class for a given
detector serial number and it is given an image from a different detector.

Ideally your implementation will inherit from a similar Format class
and just apply further customizations.  Your implementation will be
chosen to read the image if it is the most customized, i.e. it derives
from the longest chain of ancestors, all of which claim to understand
the image.

Class Methods
-------------

The class methods need to use the built in factories to construct descriptions
of the experimental apparatus from the image, namely the goniometer, detector,
beam and scan. In many cases the "simple" model will be the best which is
often trivial. In other cases it may be more complex but will hopefully
correspond to an already existing factory method.

Adding Support for RIGAKU SATURN
--------------------------------

The simplest way to demonstrate how to do something is to go ahead and do it -
which is what we will do here. The Rigaku Saturn detectors record the images
as a comprehensive text header in SMV format followed by a straightforward
unsigned short matrix containing the pixel intensities. At this stage we will
focus on making use of the header information. However before we do anything,
obtain an example of the image you are trying to support - you will need this.

Fortunately there is already a parser available for SMV format images -
FormatSMV - so we will make life a little easier by starting with that. Also
let's call the class something clear like FormatSMVRigakuSaturn.py:

    #!/usr/bin/env python
    # FormatSMVRigakuSaturn.py
    #   Copyright (C) 2011 Diamond Light Source, Graeme Winter
    #
    #   This code is distributed under the BSD license, a copy of which is
    #   included in the root directory of this package.
    #
    # An implementation of the SMV image reader for Rigaku Saturn images.
    # Inherits from FormatSMV.

    from FormatSMV import FormatSMV

    class FormatSMVRigakuSaturn(FormatSMV):
        '''A class for reading SMV format Rigaku Saturn images, and correctly
        constructing a model for the experiment from this.'''

(if you are coding this in your $HOME/.dxtbx directory, you may need to add

    import os
    import sys
    assert (os.environ['XIA2_ROOT'] in sys.path)

    FormatRoot = os.path.append(os.environ['XIA2_ROOT'], 'dxtbx', 'format')

    if not FormatRoot in sys.path:
        sys.path.append(FormatRoot)

to the start of the class file. However I would personally just do all of the
coding in the dxtbx directory as that will keep things simple. Anyway, we
digress. The first thing which needs to be added is a static method named
"understand" which will say whether we can make sense of this image and in
essence whether this is the tool for the job. This cannot make use of any
of the class methods as it is static - however if you were to implement an
image header parser as a static method, that can be used here. Here we're in
luck - the FormatSMV already has a static method which will return the contents
of the image header as a dictionary.

As this is a demonstration case, the check here is rather extensive looking
for every header keyword we will want. In general you would probably not
be quite this thorough.

        @staticmethod
        def understand(image_file):
            '''Check to see if this looks like a Rigaku Saturn SMV format
            image, i.e. we can make sense of it. Essentially that will be if
            it contains all of the keys we are looking for.'''

            size, header = FormatSMV.get_smv_header(image_file)

            wanted_header_items = [
                'DETECTOR_NUMBER', 'DETECTOR_NAMES',
                'CRYSTAL_GONIO_NUM_VALUES', 'CRYSTAL_GONIO_NAMES',
                'CRYSTAL_GONIO_UNITS', 'CRYSTAL_GONIO_VALUES',
                'DTREK_DATE_TIME',
                'ROTATION', 'ROTATION_AXIS_NAME', 'ROTATION_VECTOR',
                'SOURCE_VECTORS', 'SOURCE_WAVELENGTH',
                'SOURCE_POLARZ', 'DIM', 'SIZE1', 'SIZE2',
                ]

            for header_item in wanted_header_items:
                if not header_item in header:
                    return False

            detector_prefix = header['DETECTOR_NAMES'].split()[0].strip()

            more_wanted_header_items = [
                'DETECTOR_DIMENSIONS', 'DETECTOR_SIZE', 'DETECTOR_VECTORS',
            'GONIO_NAMES', 'GONIO_UNITS', 'GONIO_VALUES', 'GONIO_VECTORS',
                'SPATIAL_BEAM_POSITION'
                ]

            for header_item in more_wanted_header_items:
                if not '%s%s' % (detector_prefix, header_item) in header:
                    return False

            return True

And - just to make sure our coding is going along the right lines, let's
add a __main__ method to give it a tiny workout:

    if __name__ == '__main__':

        import sys

        for arg in sys.argv[1:]:
            print FormatSMVRigakuSaturn.understand(arg)

Clearly when you run this and give it an image from the right kind of detector
you would expect this to return True. If it does not, then you have some
debugging to do. Let's assume it returns True - it does for me. It is then time
to implement the methods to do the rest of the work.

Constructor
-----------

The motif for this is that the constructor is passed the image file, and that
that is the only opportunity to pass information in to the system.

        def __init__(self, image_file):
            '''Initialise the image structure from the given file, including a
            proper model of the experiment. Easy from Rigaku Saturn images as
            they contain everything pretty much we need...'''

            assert(self.understand(image_file))

            FormatSMV.__init__(self, image_file)

            return

This is very straightforward - just call the superclass constructor which
will do most of the work we need. Keep in mind though that it will also
call the extra methods included below... which do turn out to be rather more
complex than anticipated for implementation but are therefore rather
comprehensive.

Goniometer
----------

This one is made slightly simpler by the expedient that the answer is given
in the file - however if it were not it should be straightforward to compute
the true rotation axis and also the fixed rotation, from the information
given. This would probably closely model the code for the detector (which
did need all of the sums doing properly)

        def _xgoniometer(self):
            '''Initialize the structure for the goniometer - this will need to
            correctly compose the axes given in the image header. In this case
            this is made rather straightforward as the image header has the
            calculated rotation axis stored in it. We could work from the
            rest of the header and construct a goniometer model.'''

            axis = tuple(map(float, self._header_dictionary[
                'ROTATION_VECTOR'].split()))

            return self._xgoniometer_factory.KnownAxis(axis)

Detector
--------

This worked out much more complex to do properly in a general manner, as
much as anything as it is in principle attached to six axes - three rotations
and three translations. In addition to this the detector description needs a
whole load of stuff from the header which needs to be fetched out and unpacked.
At this moment it's good at least that the reference frame is consistent.

        def _xdetector(self):
            '''Return a model for the detector, allowing for two-theta offsets
            and the detector position. This will be rather more complex...'''

            detector_name = self._header_dictionary[
                'DETECTOR_NAMES'].split()[0].strip()

            detector_axes = map(float, self._header_dictionary[
                '%sDETECTOR_VECTORS' % detector_name].split())

            detector_fast = matrix.col(tuple(detector_axes[:3]))
            detector_slow = matrix.col(tuple(detector_axes[3:]))

            beam_pixels = map(float, self._header_dictionary[
                '%sSPATIAL_DISTORTION_INFO' % detector_name].split()[:2])
            pixel_size = map(float, self._header_dictionary[
                '%sSPATIAL_DISTORTION_INFO' % detector_name].split()[2:])
            image_size = map(int, self._header_dictionary[
                '%sDETECTOR_DIMENSIONS' % detector_name].split())

            detector_origin = - (
                beam_pixels[0] * pixel_size[0] * detector_fast + \
                beam_pixels[1] * pixel_size[1] * detector_slow)

Given the position of the detector and it's orientation, now get the list of
rotations and translations which are applied to this to determine where
it will actually be in the laboratory frame. N.B. we compute the list of R's
and T's then apply them in reverse.

            gonio_axes = map(float, self._header_dictionary[
                '%sGONIO_VECTORS' % detector_name].split())
            gonio_values = map(float, self._header_dictionary[
                '%sGONIO_VALUES' % detector_name].split())
            gonio_units = self._header_dictionary[
                '%sGONIO_UNITS' % detector_name].split()
            gonio_num_axes = int(self._header_dictionary[
                '%sGONIO_NUM_VALUES' % detector_name])

            rotations = []
            translations = []

            for j, unit in enumerate(gonio_units):
                axis = matrix.col(gonio_axes[3 * j:3 * (j + 1)])
                if unit == 'deg':
                    rotations.append(axis.axis_and_angle_as_r3_rotation_matrix(
                        gonio_values[j], deg = True))
                    translations.append(matrix.col((0.0, 0.0, 0.0)))
                elif unit == 'mm':
                    rotations.append(matrix.sqr((1.0, 0.0, 0.0,
                                                 0.0, 1.0, 0.0,
                                                 0.0, 0.0, 1.0)))
                    translations.append(gonio_values[j] * axis)
                else:
                    raise RuntimeError, 'unknown axis unit %s' % unit

            rotations.reverse()
            translations.reverse()

            for j in range(gonio_num_axes):
                detector_fast = rotations[j] * detector_fast
                detector_slow = rotations[j] * detector_slow
                detector_origin = rotations[j] * detector_origin
                detector_origin = translations[j] + detector_origin

            overload = int(self._header_dictionary['SATURATED_VALUE'])

            return self._xdetector_factory.Complex(
                detector_origin.elems, detector_fast.elems,
                detector_slow.elems, pixel_size, image_size, overload)

Beam
----

This is relatively simple as everything we need is already in place...

        def _xbeam(self):
            '''Return a simple model for the beam.'''

            beam_direction = map(float, self._header_dictionary[
                'SOURCE_VECTORS'].split()[:3])

            polarization = map(float, self._header_dictionary[
                'SOURCE_POLARZ'].split())

            p_fraction = polarization[0]
            p_plane = polarization[1:]

            wavelength = float(self._header_dictionary['SCAN_WAVELENGTH'])

            return self._xbeam_factory.Complex(
                beam_direction, p_fraction, p_plane, wavelength)

Scan
----

The scan is also relatively simple - however there is no epoch recorded for
each image individually, which is untidy. A better way to do this may
be to find the first image in the scan and add the cumulative exposure time
to it.

        def _xscan(self):
            '''Return the scan information for this image.'''

            rotation = map(float, self._header_dictionary['ROTATION'].split())

            format = self._xscan_factory.Format('SMV')
            epoch = time.mktime(time.strptime(self._header_dictionary[
                'DTREK_DATE_TIME'], '%d-%b-%Y %H:%M:%S'))

            exposure_time = rotation[3]
            osc_start = rotation[0]
            osc_range = rotation[1]

            return self._xscan_factory.Single(
                self._image_file, format, exposure_time,
                osc_start, osc_range, epoch)
