yayahjb / cbflib

CBFlib repository cloned from SF CBFlib repository as of 1 Dec 15
8 stars 20 forks source link

Full CBF with still images causes some surprising problems #6

Closed graeme-winter closed 7 years ago

graeme-winter commented 7 years ago

Not necessarily an issue report, more a case where we should have a think about the correct behaviour.

imgCIF header:

loop_
_diffrn_scan_axis.scan_id
_diffrn_scan_axis.axis_id
_diffrn_scan_axis.angle_start
_diffrn_scan_axis.angle_range
_diffrn_scan_axis.angle_increment
_diffrn_scan_axis.displacement_start
_diffrn_scan_axis.displacement_range
_diffrn_scan_axis.displacement_increment
SCAN1 GON_OMEGA  30.0000 0.0000 0.0000 0.0 0.0 0.0
SCAN1 GON_CHI    0.0000 0.0000 0.0000 0.0 0.0 0.0
SCAN1 GON_PHI    0.0000 0.0000 0.0000 0.0 0.0 0.0
SCAN1 DET_Z      0.0 0.0 0.0 209.35 0.0 0.0
SCAN1 DET_Y      0.0 0.0 0.0 0.0 0.0 0.0
SCAN1 DET_X      0.0 0.0 0.0 0.0 0.0 0.0

Code:

import pycbf, sys
h = pycbf.cbf_handle_struct()
h.read_widefile(sys.argv[1], pycbf.MSG_DIGEST)
g = h.construct_goniometer()
angles = tuple(g.get_rotation_range())

I would have expected angles to be 30, 30 but I get

dials.python bugreport.py grid_full_cbf_0005.cbf
Traceback (most recent call last):
  File "bugreport.py", line 5, in <module>
    angles = tuple(g.get_rotation_range())
  File "/Users/graeme/svn/cctbx/build/lib/pycbf.py", line 725, in get_rotation_range
    return _pycbf.cbf_positioner_struct_get_rotation_range(self)
Exception: CBFlib Error(s): CBF_NOTFOUND 

presumably because cbf_simple looks for axes which are changing - perhaps we should have Python bindings to a lower level, or handle this differently in cbflib? Example image attached. grid_full_cbf_0005.zip

graeme-winter commented 7 years ago

OK, inspecting the code is helpful - this is it:

        for (axis = 0; axis < goniometer->axes; axis++)
            if (goniometer->axis [axis].type == CBF_ROTATION_AXIS)
                if (goniometer->axis [axis].increment)
                {
                    if (start)
                        *start = goniometer->axis [axis].start;
                    if (increment)
                        *increment = goniometer->axis [axis].increment;
                    return 0;
                }
        return CBF_NOTFOUND;

I am going to guess that the fact that it checks increment here is the issue. I see the logic - this will keep going until it finds one which has increment != 0 then return else complain.

I also note that this gives the right behaviour:

import pycbf, sys
h = pycbf.cbf_handle_struct()
h.read_widefile(sys.argv[1], pycbf.MSG_DIGEST)
g = h.construct_goniometer()
x = g.rotate_vector(0.0, 1, 0, 0)
y = g.rotate_vector(0.0, 0, 1, 0)
z = g.rotate_vector(0.0, 0, 0, 1)

print x
print y
print z

based on the values above, giving

[1.0, 0.0, 0.0]
[0.0, 0.8660254037844387, 0.49999999999999994]
[0.0, -0.49999999999999994, 0.8660254037844387]

i.e. a 30 degree rotation about x...

graeme-winter commented 7 years ago

So, adding the following loop to cbf_simple.c

    /* Second pass just look for first rotation axis */ 

        for (axis = 0; axis < goniometer->axes; axis++) {
      if (goniometer->axis [axis].type == CBF_ROTATION_AXIS) {
        if (start) *start = goniometer->axis [axis].start;
        if (increment) *increment = goniometer->axis [axis].increment;
        return 0;
      }
    }

to cbf_get_rotation_range just before returning CBF_NOTFOUND gives me exactly the behaviour I would expect and is well defined...

graeme-winter commented 7 years ago

Proposed patch

diff --git a/src/cbf_simple.c b/src/cbf_simple.c
index fe3c712..f00148c 100644
--- a/src/cbf_simple.c
+++ b/src/cbf_simple.c
@@ -6830,7 +6830,8 @@ extern "C" {
             return CBF_ARGUMENT;

-        /* Currently just return the range of the first rotation axis */
+        /* Currently just return the range of the first rotation axis with 
+      non-zero increment */

         for (axis = 0; axis < goniometer->axes; axis++)

@@ -6848,7 +6849,17 @@ extern "C" {

                     return 0;
                 }
-
+   
+   /* Second pass just look for first rotation axis */ 
+   
+        for (axis = 0; axis < goniometer->axes; axis++) {
+     if (goniometer->axis [axis].type == CBF_ROTATION_AXIS) {
+       if (start) *start = goniometer->axis [axis].start;
+       if (increment) *increment = goniometer->axis [axis].increment;
+       return 0;
+     }
+   }
+   
         return CBF_NOTFOUND;
     }
yayahjb commented 7 years ago

Having some sort of a fix seems reasonable, but I am puzzled by the data saying that none of the axes are moving, but wanting the code to say that one of them is moving by 30 degrees. Shouldn't we say what the data says -- that the range is 0, i.e., that this is a still shot taken at the static position, omega=30, chi=phi=0?

graeme-winter commented 7 years ago

None of the axes are moving, but the first rotation axis has a static setting of 30 degrees, so get_rotation_range (with the patch above) gives 30 for start, 0 for range.

graeme-winter commented 7 years ago

Getting all of the angles and reading this "properly" would give you phi=chi=0 omega=30 which when combined with the axis definitions of phi, chi and omega would give you the rotation matrix which is produced when you ask the gonio to rotate. The confusion here is when using cbf_simple via pycbf, which only allows for a singe start and range

yayahjb commented 7 years ago

Dear Graeme,

On reflection, I agree with your interpretation and your patch, since for a still, specifying any axis as the rotation axis is as good as specifying any other axis as the rotation axis. I've merged your pull request into master.

Regards, Herbert

On Thu, Apr 20, 2017 at 6:47 AM, Graeme Winter notifications@github.com wrote:

Getting all of the angles and reading this "properly" would give you phi=chi=0 omega=30 which when combined with the axis definitions of phi, chi and omega would give you the rotation matrix which is produced when you ask the gonio to rotate. The confusion here is when using cbf_simple via pycbf, which only allows for a singe start and range

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/yayahjb/cbflib/issues/6#issuecomment-295679023, or mute the thread https://github.com/notifications/unsubscribe-auth/AEPiAYlCgj1tA8_wivvlM3uoT2LOEUAmks5rxzfSgaJpZM4NCkJN .