Position to data coords

[tmichela]

Find pixels coordinates in data array from the positions on an assembled image.

e.g. visualization of pixel coordinates from ring positions given by pyfai:

regular grid around the centre:

[takluyver]

If I'm following correctly, these 'x' and 'y' things should really be 'fast scan' and 'slow scan'. E.g. for AGIPD, the fast scan direction, which you've called 'x' here, is typically along the y axis.

(I know I sound like a broken record about this, but geometry code is hard to follow at the best of times, so I really want to avoid making it harder with half-right names)

[tmichela]

These are the x and y coordinates in a (2D) data array. I can call them fast and slow scan if you prefer but I found it more confusing as it is not the same as what we call ss/fs usually in this code base. And in this case x will always be fast scan and y slow scan.

[takluyver]

For new methods like this, where we can be consistent across all detectors, how would you feel about saying the input is always in metres, and if you've got e.g. data in pixel units, you apply the unit like data * geom.pixel_size?

[takluyver]

The sub-array excluding invalid coordinates feels a little bit janky. What about a boolean array which is True for valid coordinates and False for invalid ones? I think that makes it easy to match up the input with the output if needed.

[takluyver]

I think it should be possible - and not much more complex - to do this 'properly', preserving subpixel offsets and rotations and so on. Something like (pseudocode):

for mod in modules:
    for t, tile in enumerate(mod):
        ss_coord, fs_coord = find_data_coords_for_tile(data, tile)
        on_tile = (0 < ss_coord) & (ss_coord < tile.ss_pixels) & (0 < fs_coord) & (fs_coord < tile.fs_pixels)
        ss_coord, fs_coord = tile_coords_to_module(ss_coord[on_tile], fs_coord[on_tile], t)
        # Record these results to result arrays

I'm assuming the numbers of tiles per detector are small enough to get away with just checking each tile one by one. So far it's a max of 128 tiles (LPD-1M), IIRC. If that's too slow, we'd need a bit more complexity to find the relevant tiles efficiently, but hopefully not much more.

[tmichela]

Hum, the issue I can see is if one relies on our (snapped) assembled data to find coordinates, using the subpixels/rotated positions here might lead to wrong coordinate mapping :/

[JamesWrigley]

FYI, I did something vaguely similar for MID recently. They were drawing small ROIs on an assembled image, so to avoid loading/assembling all the modules I wrote a function that took in the ROI bounds and created a mask in array coordinates:

# Helper function to create a mask in an unassembled (module, y, x)
# shape from the ROI bounds in an assembled image.
def generate_mask(geom, xmin, xmax, ymin, ymax):
    # Create a LUT, a flat array the same size of a data array where the elements are the assembled-array indices of the same pixel.
    lut_data = np.empty((16, 512, 128))
    lut_data.flat[:] = np.arange(lut_data.size)    
    assembled, _ = geom.position_modules(lut_data)
    lut = np.argsort(assembled.flat)[:lut_data.size].astype(np.uint64)

    # Create an assembled mask image from ROI coordinates
    assembled_mask = np.full_like(assembled, False, dtype=bool)
    assembled_mask[ymin:ymax, xmin:xmax] = True
    assembled_mask[np.isnan(assembled)] = False

    # Create a mask in array coordinates and use the LUT to assign the masked elements
    mask = np.full((16, 512, 128), False, dtype=bool)
    mask.flat[assembled_mask.flat[lut]] = True

    # Find the modules used in the mask
    module_numbers = np.any(mask, axis=(1, 2)).nonzero()[0]
    
    return mask[module_numbers], module_numbers

This was specific to the AGIPD and a rectangular ROI, but the same principle should work for any geometry and mask. I wonder if it would be simpler to add this thing as a generic disassemble() function, so all positions_to_data_coords() would have to do is convert the assembled positions into an assembled mask and pass it to disassemble()?