NSLS-II-CSX / csxtools

Useful python tools for CSX (23-ID)
http://nsls-ii-csx.github.io/csxtools
Other
4 stars 13 forks source link

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

Open leofang opened 4 years ago

leofang commented 4 years ago

@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 commented 4 years ago

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 commented 4 years ago

@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 commented 4 years ago

@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 commented 4 years ago

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

ambarb commented 4 years ago

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 commented 4 years ago

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 https://github.com/NSLS-II-CSX/csxtools/issues/68#issue-558382278.

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 commented 4 years ago

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 commented 4 years ago

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

mrakitin commented 4 years ago

Looks like the current is older: image

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 commented 4 years ago

oops. i added the line you need

ambarb commented 4 years ago

!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 commented 4 years ago

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 commented 4 years ago

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

danielballan commented 4 years ago

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 commented 4 years ago

@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 https://github.com/NSLS-II-CSX/csxtools/blob/81a64d76074bafff012afcb23119b50ccd3aaa50/csxtools/utils.py#L125-L127 the _crop() function was invoked and hence led to the reported error, which can be fixed by the patch (https://github.com/NSLS-II-CSX/csxtools/issues/68#issuecomment-580991283). 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 commented 4 years ago

@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 commented 4 years ago

@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 commented 4 years ago

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: https://github.com/NSLS-II-CSX/csxtools/blob/81a64d76074bafff012afcb23119b50ccd3aaa50/csxtools/utils.py#L212-L216 where the input argument roi is already of the form [x, y, x+w, y+h] based on this block: https://github.com/NSLS-II-CSX/csxtools/blob/81a64d76074bafff012afcb23119b50ccd3aaa50/csxtools/utils.py#L63-L69 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": https://github.com/NSLS-II-CSX/csxtools/blob/81a64d76074bafff012afcb23119b50ccd3aaa50/csxtools/utils.py#L24 https://github.com/NSLS-II-CSX/csxtools/blob/81a64d76074bafff012afcb23119b50ccd3aaa50/csxtools/utils.py#L203 but in get_fastccd_flatfield the rotation (due to np.rot90) is "ccw": https://github.com/NSLS-II-CSX/csxtools/blob/81a64d76074bafff012afcb23119b50ccd3aaa50/csxtools/utils.py#L273 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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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.

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 .

Manual_Flatfiled_Correction_VS_csxtools_function_p13 Manual_Flatfiled_Correction_VS_csxtools_function_p15

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 commented 4 years ago

rel: #73