Geometry fixes

[jadball]

Addresses #366

Highlights:

dty units of microns now enforced in dataset.py

Microns are overwhelmingly our choice for .par files and I think other solutions introduce significant potentials for confusion.
I solve this by checking the BLISS axis units for dtymotor and converting from mm to microns if needed:

def import_motors_from_master(self):
    dty_is_mm = False
            (...)
            # read from h5:
            dty = hin[scan]["instrument/positioners"][self.dtymotor]
            # get the unit for dty:
            dty_unit = hin[scan]["instrument/positioners"][self.dtymotor].attrs['units']
            if dty_unit == 'mm':
                dty_is_mm = True
            (...)
    # ensure dty is in microns - this is for compatability with geometry
    if dty_is_mm:
        self.dty = self.dty * 1000
        logging.info("Converted dty from mm to um")
    logging.info("imported omega/dty")

Introduced a new geometric parameter ybeam

This is the nominal y-value of the beam in the lab frame. Usually 0 for NSCOPE, 14 mm for TDXRD.
Old approach:

def sample_to_lab_sincos(sx, sy, y0, dty, sinomega, cosomega):
    # first, rotate sx and sy by omega about its origin (rotation axis)
    sxr = sx * cosomega - sy * sinomega
    syr = sx * sinomega + sy * cosomega

    # then translate about the rotated frame
    lx = sxr
    ly = syr + dty - y0

    return lx, ly

New approach:

def sample_to_lab_sincos(sx, sy, y0, ybeam, dty, sinomega, cosomega):
    # first, rotate sx and sy by omega about its origin (rotation axis)
    sxr = sx * cosomega - sy * sinomega
    syr = sx * sinomega + sy * cosomega

    # then translate about the rotated frame
    lx = sxr
    ly = syr + dty + (ybeam - y0)

    return lx, ly

I have updated the rest of geometry.py to match, including tests and Jupyter notebooks.
This affects less than you might think, because ybeam often cancels out and we rarely use the lab frame in computation.

Updated geometry.py function names and docstrings for clarity

Decoupled dtyi (discretisation of dty) from "step space" (si, sj)

New dty to dtyi conversion function now requires ymin but is conceptually simpler and does not induce a sign flip any more:

def dty_to_dtyi(dty, ystep, ymin):
    """
    We take continuous dty space and discretise it
    We do this by counting the number of ysteps away from ymin
    ymin should probably be ds.ymin
    This should be equivalent to determining the closest bin in ds.ybincens
    But this is neater as it accounts for values outside the ds.ybincens range

    dtyi = (dty - ymin) / ystep
    """
    dtyi = np.round((dty - ymin) / ystep).astype(int)
    return dtyi

I've updated the point_by_point.py code accordingly to include ymin and ybeam where needed.
ybeam is basically only needed when considering shifts in tomographic reconstructions.

Introduced new standard way to generate grids in (si, sj) to reconstruct on

This is geometry.step_grid_from_ybincens which I discuss here: #366 (comment)

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[jadball]

Caveat here: we don't store the units for dty in sparse files...

[jonwright]

A couple of quick comments:

• I don't understand the difference between ybeam and y0. Why not just say y0=14 for TDXRD? It seems to add a redundant parameter.

• the units thing should be for a new version and not in this release. Old datasets have no dty units and I used mm for distance and pixel size a few times in the past. Forcing microns now seems to break the fitting I did with Ryan.

[jadball]

• I don't understand the difference between ybeam and y0. Why not just say y0=14 for TDXRD? It seems to add a redundant parameter.

You're right - I now understand this better. If we define the lab frame via the beam (e.g ly = 0 when a grain is in the beam rather than ly = ybeam) then it all resolves without issue.

The final sticking point is the shift calculation for the tomographic reconstruction, which I think I've nailed now. New commits incoming!

[jonwright]

To modify the dty numbers is a breaking change. It should not go into the next release.

[jadball]

New highlights:

Updated geometry.py function names and docstrings for clarity

Decoupled dtyi (discretisation of dty) from "step space" (si, sj)

New dty to dtyi conversion function now requires ymin but is conceptually simpler and does not induce a sign flip any more:

def dty_to_dtyi(dty, ystep, ymin):
    """
    We take continuous dty space and discretise it
    We do this by counting the number of ysteps away from ymin
    ymin should probably be ds.ymin
    This should be equivalent to determining the closest bin in ds.ybincens
    But this is neater as it accounts for values outside the ds.ybincens range

    dtyi = (dty - ymin) / ystep
    """
    dtyi = np.round((dty - ymin) / ystep).astype(int)
    return dtyi

I've updated the point_by_point.py code accordingly

Introduced new standard way to generate grids in (si, sj) to reconstruct on

This is geometry.step_grid_from_ybincens which I discuss here: #366 (comment)

Accounting for shift and pad properly in tomographic reconstructions:

def sino_shift_and_pad(y0, ny, ymin, ystep):
    """
    Determine the difference in sinogram pixels between
    the central dtyi value of the sinogram and the dtyi value of the rotation axis (from y0)
    To do this, we convert y0 into dtyi without rounding
    Also determine the minimum pad to get the whole sample in the frame
    Which should be double the shift
    """
    # get the middle row of the sinogram in pixel space
    ymid_px = ny/2
    # get the value of y0 in pixel space
    # this is the same as dty_to_dtyi without rounding
    y0_px = (y0 - ymin)/ystep
    shift = ymid_px - y0_px
    pad = np.ceil(np.abs(shift) * 2).astype(int) + 1
    return shift, pad

[jadball]

@jonwright ready when you are if you're happy to merge :)