roi argument for get_fastccd_images() is broken ... plus some flatfield dependency discussion

[leofang]

@wen-hu needs the flatfield support to process some ptychography datasets, but I found an error with the latest csxtools (0.1.15):

>>> from csxtools.utils import get_fastccd_flatfield, get_fastccd_images, get_images_to_4D
>>> import numpy as np
>>> from databroker import Broker
>>>
>>> db = Broker.named('csx')
>>> dark0_f = 127832; dark1_f = dark0_f + 1; dark2_f = dark0_f +2
>>> flat_im = get_fastccd_flatfield(db[127831], dark=(db[dark0_f], db[dark1_f], db[dark2_f]))
/home/leofang/conda_envs/ptycho_exp_Jan2020/lib/python3.6/site-packages/csxtools/utils.py:190: UserWarning: Images and get_images are deprecated. Use Header.data('fccd_image') instead.
  images = header.db.get_images(header, tag)

/home/leofang/conda_envs/ptycho_exp_Jan2020/lib/python3.6/site-packages/csxtools/utils.py:271: RuntimeWarning: invalid value encountered in greater
  flat[flat > limits[1]] = np.nan
Flatfield correction removed 79554 pixels (8.29 %)
>>> 
>>> flat_im.shape
(1000, 960)
>>> silcerator = get_fastccd_images(db[127849], dark_headers=(db[127846], db[127847], db[127848]), flat=flat_im, roi=[0,0,1000,960])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/leofang/conda_envs/ptycho_exp_Jan2020/lib/python3.6/site-packages/csxtools/utils.py", line 94, in get_fastccd_images
    b = get_images_to_3D(bgnd_events, dtype=np.uint16)
  File "/home/leofang/conda_envs/ptycho_exp_Jan2020/lib/python3.6/site-packages/csxtools/utils.py", line 175, in get_images_to_3D
    im = np.vstack([np.asarray(im, dtype=dtype) for im in images])
  File "/home/leofang/conda_envs/ptycho_exp_Jan2020/lib/python3.6/site-packages/csxtools/utils.py", line 175, in <listcomp>
    im = np.vstack([np.asarray(im, dtype=dtype) for im in images])
  File "/home/leofang/conda_envs/ptycho_exp_Jan2020/lib/python3.6/site-packages/slicerator/__init__.py", line 474, in <genexpr>
    return (self._get(i) for i in range(len(self)))
  File "/home/leofang/conda_envs/ptycho_exp_Jan2020/lib/python3.6/site-packages/slicerator/__init__.py", line 461, in _get
    return self._proc_func(*(copy(a[key]) for a in self._ancestors))
  File "/home/leofang/conda_envs/ptycho_exp_Jan2020/lib/python3.6/site-packages/slicerator/__init__.py", line 685, in proc_func
    return func(*(x + args), **kwargs)
  File "/home/leofang/conda_envs/ptycho_exp_Jan2020/lib/python3.6/site-packages/csxtools/utils.py", line 209, in _crop_images
    return _crop(image, roi)
  File "/home/leofang/conda_envs/ptycho_exp_Jan2020/lib/python3.6/site-packages/csxtools/utils.py", line 216, in _crop
    return image.T[roi[1]:roi[3], roi[0]:roi[2]].T
AttributeError: 'ImageStack' object has no attribute 'T'

Seems like image in _crop is assumed to be a numpy array but it actually isn't (or is not handled correctly by @pipeline). Perhaps we can hot-fix it by not doing a rotate90 or np.rot90 to avoid transpose (T) operation to the flatfield? But not sure if this would break existing code, or if existing code just never worked 🤣

Any thoughts? @danielballan @mrakitin

[leofang]

This will make the obtained silcerator work with get_images_to_3D/get_images_to_4D (note the casting). Should I create a PR?

diff --git a/csxtools/utils.py b/csxtools/utils.py
index 62d3b98..1fcf686 100644
--- a/csxtools/utils.py
+++ b/csxtools/utils.py
@@ -213,6 +213,8 @@ def _crop(image, roi):
     image_shape = image.shape
     # Assuming ROI is specified in the "rotated" (correct) orientation
     roi = [image_shape[-2]-roi[3], roi[0], image_shape[-1]-roi[1], roi[2]]
+    # need to perform an "unsafe" cast here, or fastccd.correct_images would be cranky!
+    image = np.asarray(image, dtype=np.uint16)
     return image.T[roi[1]:roi[3], roi[0]:roi[2]].T

[ambarb]

@leofang I know @stuwilkins wanted to rotate the images at particular times. Maybe you can speak with him.

I haven't had issues with the flatfield myself. What version of csxtool are you and @wen-hu using?

[ambarb]

@leofang. Sorry, I did not see the version that you provided.

I am using 0.1.13 csxtools and my code implementation of the flat field works correctly.
Can you point me to the version that you are speaking about? I would like to test my notebooks agains this version.

Also, there is discussion with @danielballan @mrakitin and @jklynch to put data streams that contact a fully correct flatfield.

[ambarb]

note this is this discussion of the flatfield stream is to make things play nicely with xi-cam

[ambarb]

Disregard location of version. I see that master is at 0.1.15 and just has not been pushed to the facilities current standard conda version for analysis at CSX. I will clone and update so I can test my notebooks.

[leofang]

Oh, Andi, sorry for my long silence! I just found the notification from my GitHub email sea 😅

Yes, I used 0.1.15 to test it. Please try it with the reproducer from #68 (comment).

Also, please let me know how you made it work with 0.1.13, so that I can investigate further. It could be that I was misled by the documentation and used the API wrong.

Thanks!

[ambarb]

I as on vacation part of last week. For 0.1.13, you can get the flatfield to generally work with

from csxtools.utils import get_fastccd_images,get_images_to_3D, get_images_to_4D, get_fastccd_flatfield,fccd_mask
from csxtools.ipynb import image_stack_to_movie, show_image_stack
from csxtools.image import stackmean, images_mean, images_sum

bgnd1 = db[117332+3]
bgnd2 = db[117332+2]
bgnd8 = db[117332+1]
flat = db[117332]
ff = get_fastccd_flatfield(flat, (bgnd8, bgnd2, bgnd1))

n=0
bgnd2 = None
bgnd1 = None


scans_no     = [117493 for i in range(117493, 117499+1,2)]  #Escans for 300K for 4 Ls
scans_bgnd8 = [x-1 for x in scans_no]
                

bgnd8 = db[scans_bgnd8[n],]#db[111531,]
scan_no = scans_no[n]  # random sat thing
h = db[scan_no]

#gains, most to least sensiteve
images = get_fastccd_images(h, (bgnd8[0], bgnd2, bgnd1), flat=ff)
#### FOR A SCAN  ###
stack = get_images_to_4D(images)

YOu should be able to just copy and paste to test in the conda env for CSX analysis from 2019-1.2 ("current" on nsls-ii notebook jupyternub)

There are options to the function likeflat=fccd_mask(), limits=(0.8, 1.2), but they are not required, and the fccd_mask() is no longer valid unless someone has defined a new one.

[ambarb]

There is also a pre-dated issue #69 for the whole story for anyone who is interested.

[mrakitin]

Looks like the current is older:

Also, I get an error for undefined h. What should it be?

/opt/conda_envs/analysis-2018-2.1/lib/python3.6/site-packages/csxtools/utils.py:181: UserWarning: Images and get_images are deprecated. Use Header.data('fccd_image') instead.
  images = header.db.get_images(header, tag)
/opt/conda_envs/analysis-2018-2.1/lib/python3.6/site-packages/csxtools/utils.py:262: RuntimeWarning: invalid value encountered in greater
  flat[flat > limits[1]] = np.nan
Flatfield correction removed 74169 pixels (7.73 %)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-c8e2e6e3fb0a> in <module>()
     24 
     25 #gains, most to least sensiteve
---> 26 images = get_fastccd_images(h, (bgnd8[0], bgnd2, bgnd1), flat=ff)
     27 #### FOR A SCAN  ###
     28 stack = get_images_to_4D(images)

NameError: name 'h' is not defined

[ambarb]

oops. i added the line you need

[ambarb]

!which python is cute. Thanks. I was just importing the packages and looking at the version numbers in the bash session. I forgot to confirm with the notebook itself.

[mrakitin]

Yes, the fastest way to get your path. I sometimes use:

import os
print(os.__file__)

but it's longer to type.

As for the example, it works for me now, but returns this (I assume it's OK):

/opt/conda_envs/analysis-2018-2.1/lib/python3.6/site-packages/csxtools/utils.py:181: UserWarning: Images and get_images are deprecated. Use Header.data('fccd_image') instead.
  images = header.db.get_images(header, tag)
/opt/conda_envs/analysis-2018-2.1/lib/python3.6/site-packages/csxtools/utils.py:262: RuntimeWarning: invalid value encountered in greater
  flat[flat > limits[1]] = np.nan
Flatfield correction removed 74169 pixels (7.73 %)
Missing dark image for gain setting 2
Missing dark image for gain setting 1

[ambarb]

YEs, that is correct. There was a issue for that, but it was closed after some changes before christmas.

[danielballan]

If there is confusion about which Python you are currently in, use import sys; sys.executable which is always correct. Running !which python will often work but it is easy to come up with scenarios where the first python on the PATH is not the currently-running python.

[leofang]

@ambarb @mrakitin Thanks for the discussion. I think the problem in my reproducer was that I also set roi in get_fastccd_images(). Due to this piece of code

csxtools/csxtools/utils.py

Lines 125 to 127 in 81a64d7

	# Crop Flatfield image
	if flat is not None and roi is not None:
	flat = _crop(flat, roi)

the _crop() function was invoked and hence led to the reported error, which can be fixed by the patch (#68 (comment)). This bug has been there since at least 0.1.13, so you can try it out yourselves.

If roi were not set, we would not see the error.

I am still a bit concerned with the rotation of the flat field in get_fastccd_flatfield(), though. Will need to think about it and determine if there's any issue.

[ambarb]

@leofang thanks for explaining your use case with the roi. _crop was added as a way to deal with memory issues of large data sets. After lazy loading was added, this was no longer a problem. Does the cropping save that much time?

As for the rotation, I will have a look, but the output from get_fastccd_flatfield() is the raw results as seem from the camera's perspective. We have the camera laying on it's side in the chamber. The idea was to use the flat filed output as is, apply it in get_fastccd_images() and then once those images were corrected, to rotate the images so that it matches what your brain expect for the final step before outputting normalized images.

The reason this work flow makes sense to me is that however we configure the camera to run, all the image normalization is done on raw images so you don't have to worry about the cropping, extra added columns, addition of overscan, frame store versus non-frame store, binning of columns to decrease readout time.

[ambarb]

@leofang To keep our existing workflow described, do you think we need to think about refactoring the roi argument usage so that we are more future proof to the camera configurations?

[leofang]

Thanks for the detailed explanation, Andi! I think as long as the roi argument is clearly documented we should be fine.

I know almost nothing about the experimental setup, and by inspecting the code currently I see two problems:

1. The _crop function makes no sense to me:
   

   csxtools/csxtools/utils.py

   Lines 212 to 216 in 81a64d7

   	def _crop(image, roi):
   	image_shape = image.shape
   	# Assuming ROI is specified in the "rotated" (correct) orientation
   	roi = [image_shape[-2]-roi[3], roi[0], image_shape[-1]-roi[1], roi[2]]
   	return image.T[roi[1]:roi[3], roi[0]:roi[2]].T

   
   where the input argument roi is already of the form [x, y, x+w, y+h] based on this block:
   

   csxtools/csxtools/utils.py

   Lines 63 to 69 in 81a64d7

   	# Now lets sort out the ROI
   	if roi is not None:
   	roi = list(roi)
   	# Convert ROI to start:stop from start:size
   	roi[2] = roi[0] + roi[2]
   	roi[3] = roi[1] + roi[3]
   	logger.info("Computing with ROI of %s", str(roi))

   
   So I have no idea what line 215 is about. The entire _crop() function should be rewritten and thoroughly tested to handle the case when both flat and roi are present so as to close this issue.
2. Regarding the rotation, I feel the rotations are inconsistent. For example, in get_fastccd_images the rotation is "cw":
   

   csxtools/csxtools/utils.py

   Line 24 in 81a64d7

   - Rotation (returned images are rotated 90 deg cw)

   
   

   csxtools/csxtools/utils.py

   Line 203 in 81a64d7

   image = rotate90(image, 'cw')

   
   but in get_fastccd_flatfield the rotation (due to np.rot90) is "ccw":
   

   csxtools/csxtools/utils.py

   Line 273 in 81a64d7

   flat = np.rot90(flat)

   
   So, is it expected that we pass in a counterclockwise flatfield and get the ccd images rotated clockwise? Could it be a bug due to the difference in facing upstream or downstream?

[ambarb]

I see your point for the rotation. I am doing some image math to make sure that this is brought out by some analysis. Almost finished with it and will update here.

[ambarb]

I suggest we make a new issue report for the roi and leave it as an enhancement. I think whomever wants to tackle this will have to also keep in mind the "future proof-ness" of it and the plans to add metadata for concatenating the normalized images to be done by wrapper function or optional argument for get_fastccd_images()

Any objections to this approach?

[leofang]

I see your point for the rotation. I am doing some image math to make sure that this is brought out by some analysis. Almost finished with it and will update here.

Sounds good. Let me know if it is indeed a bug.

I suggest we make a new issue report for the roi and leave it as an enhancement.

I suppose we could just rename the title of this issue to indicate this is an roi-related problem? Without specifying roi it seems to be fine after all (apart from the possible rotation issue).

[ambarb]

I don't think there is a bug for the flatfield usage at large. Images E and F are 100% manual flatflield correction to the data used to generate the flatfield using a flatfield that was either rotated cw or ccw.

• Image A: cw rotation is the improper rotation direction if I were to do this by hand. it does not match the lab frame

• Image B: ccw rotation is the correct rotation direction if I were to do this by hand. it does match the lab frame - what we defined as the camera's perspective.

• Image C is that processed data from get_fastccd_images(flat=None)

• image D is C but with flat=get_fastccd_flatfield().

If the code in the repo is correct, D should match F. H shows the difference between D and F which is very small. G and F have vmin=-1, vmax=1 . Tested both version 0.1.13 and 0.1.15 .

Now, I am not sure if we want to chase down where the correction for the inconsistency that you point out is. But we should make sure that the unit test for the flat field is in place and that it is always used whenever someone touches get_fastccd_images() and any other functions that the get_fastccd_flatfield() is dependant on?

[leofang]

rel: #73