spacetelescope / stsci.imagemanip

STScI image manipulation.
Other
0 stars 1 forks source link

The C code in the package is likely incorrect on many levels #12

Closed mcara closed 5 years ago

mcara commented 5 years ago

I started working on finding an appropriate replacement function/algorithm for this package as discussed in https://github.com/spacetelescope/stsci.imagemanip/issues/9 IMO this package is so incorrect, that drizzlepac would have been better off setting LFLTFILE to all 1 instead of using stsci.imagemanip. However, I am not an expert in STIS flat field files and so I am invoking higher powers here: @philhodge @ibusko - the experts in STIS flat fields, based on references found online.

From the STIS book:

To create the single combined flat field file, calstis first expands the large scale sensitivity flat (LFLTFILE) to full format, using bilinear interpolation. The pixel-to-pixel flat, delta flat, and expanded low order flat are then multiplied together.

stsci.imagemanip is used by the drizzlepac to perform this "expansion" of LFLTFILE.

I believe the major issue with stsci.imagemanip code is that it seem that the code uses indices incorrectly for numpy arrays in the sense that it incorrectly "flattens" 2D array. The code may work correctly if arrays were using Fortran order (maybe, but my experiments failed to fix it by playing with array ordering). I am also not too sure the interpolation is truly "bi"-linear in the sense that it may not work well in the presence of cross-terms.

The test module for this package is useless and it does not help at all in understanding how STIS LFLTFILE files should be expanded (i.e., what is the intention) nor does it actually test that package works as intended.


Small synthetic test

from stsci.imagemanip import interp2d
import numpy as np

# generate a "binned" flat field:
y, x = np.indices((10, 12), dtype=np.float32)
f = 0.5 * x + 0.75 * y

# compute oversampled flat field based on the theoretical model:
y, x = np.indices((20, 24), dtype=np.float32)
x /= 2.0
y /= 2.0 
x -= 0.25
y -= 0.25 
g = 0.5 * x + 0.75 * y  # "truth"

# oversample using stsci.imagemanip:
h = interp2d.expand2d(f, (20, 24))

# compare the results:
print(g)
print(h)

Results are:

>>> print(g)  # reference
[[-0.312 -0.062  0.188  0.438  0.688  0.938  1.188  1.438  1.688  1.938  2.188  2.438  2.688  2.938  3.188  3.438  3.688  3.938  4.188  4.438  4.688  4.938  5.188  5.438]
 [ 0.062  0.312  0.562  0.812  1.062  1.312  1.562  1.812  2.062  2.312  2.562  2.812  3.062  3.312  3.562  3.812  4.062  4.312  4.562  4.812  5.062  5.312  5.562  5.812]
 [ 0.438  0.688  0.938  1.188  1.438  1.688  1.938  2.188  2.438  2.688  2.938  3.188  3.438  3.688  3.938  4.188  4.438  4.688  4.938  5.188  5.438  5.688  5.938  6.188]
 [ 0.812  1.062  1.312  1.562  1.812  2.062  2.312  2.562  2.812  3.062  3.312  3.562  3.812  4.062  4.312  4.562  4.812  5.062  5.312  5.562  5.812  6.062  6.312  6.562]
 [ 1.188  1.438  1.688  1.938  2.188  2.438  2.688  2.938  3.188  3.438  3.688  3.938  4.188  4.438  4.688  4.938  5.188  5.438  5.688  5.938  6.188  6.438  6.688  6.938]
 [ 1.562  1.812  2.062  2.312  2.562  2.812  3.062  3.312  3.562  3.812  4.062  4.312  4.562  4.812  5.062  5.312  5.562  5.812  6.062  6.312  6.562  6.812  7.062  7.312]
 [ 1.938  2.188  2.438  2.688  2.938  3.188  3.438  3.688  3.938  4.188  4.438  4.688  4.938  5.188  5.438  5.688  5.938  6.188  6.438  6.688  6.938  7.188  7.438  7.688]
 [ 2.312  2.562  2.812  3.062  3.312  3.562  3.812  4.062  4.312  4.562  4.812  5.062  5.312  5.562  5.812  6.062  6.312  6.562  6.812  7.062  7.312  7.562  7.812  8.062]
 [ 2.688  2.938  3.188  3.438  3.688  3.938  4.188  4.438  4.688  4.938  5.188  5.438  5.688  5.938  6.188  6.438  6.688  6.938  7.188  7.438  7.688  7.938  8.188  8.438]
 [ 3.062  3.312  3.562  3.812  4.062  4.312  4.562  4.812  5.062  5.312  5.562  5.812  6.062  6.312  6.562  6.812  7.062  7.312  7.562  7.812  8.062  8.312  8.562  8.812]
 [ 3.438  3.688  3.938  4.188  4.438  4.688  4.938  5.188  5.438  5.688  5.938  6.188  6.438  6.688  6.938  7.188  7.438  7.688  7.938  8.188  8.438  8.688  8.938  9.188]
 [ 3.812  4.062  4.312  4.562  4.812  5.062  5.312  5.562  5.812  6.062  6.312  6.562  6.812  7.062  7.312  7.562  7.812  8.062  8.312  8.562  8.812  9.062  9.312  9.562]
 [ 4.188  4.438  4.688  4.938  5.188  5.438  5.688  5.938  6.188  6.438  6.688  6.938  7.188  7.438  7.688  7.938  8.188  8.438  8.688  8.938  9.188  9.438  9.688  9.938]
 [ 4.562  4.812  5.062  5.312  5.562  5.812  6.062  6.312  6.562  6.812  7.062  7.312  7.562  7.812  8.062  8.312  8.562  8.812  9.062  9.312  9.562  9.812 10.062 10.312]
 [ 4.938  5.188  5.438  5.688  5.938  6.188  6.438  6.688  6.938  7.188  7.438  7.688  7.938  8.188  8.438  8.688  8.938  9.188  9.438  9.688  9.938 10.188 10.438 10.688]
 [ 5.312  5.562  5.812  6.062  6.312  6.562  6.812  7.062  7.312  7.562  7.812  8.062  8.312  8.562  8.812  9.062  9.312  9.562  9.812 10.062 10.312 10.562 10.812 11.062]
 [ 5.688  5.938  6.188  6.438  6.688  6.938  7.188  7.438  7.688  7.938  8.188  8.438  8.688  8.938  9.188  9.438  9.688  9.938 10.188 10.438 10.688 10.938 11.188 11.438]
 [ 6.062  6.312  6.562  6.812  7.062  7.312  7.562  7.812  8.062  8.312  8.562  8.812  9.062  9.312  9.562  9.812 10.062 10.312 10.562 10.812 11.062 11.312 11.562 11.812]
 [ 6.438  6.688  6.938  7.188  7.438  7.688  7.938  8.188  8.438  8.688  8.938  9.188  9.438  9.688  9.938 10.188 10.438 10.688 10.938 11.188 11.438 11.688 11.938 12.188]
 [ 6.812  7.062  7.312  7.562  7.812  8.062  8.312  8.562  8.812  9.062  9.312  9.562  9.812 10.062 10.312 10.562 10.812 11.062 11.312 11.562 11.812 12.062 12.312 12.562]]

>>> print(h)  # imagemanip
[[-1.375 -1.125 -0.875 -0.297  0.609  1.188  1.438  1.688  1.938  2.188  2.438  2.688  2.938  3.188  3.438  3.688  3.938  4.188  4.438  4.688  1.125  1.375  1.625  1.547]
 [ 1.141  1.062  1.312  1.562  1.812  2.062  2.312  2.562  2.812  3.062  3.312  3.562  3.812  4.062  4.312  4.562  3.625  3.875  4.125  3.391  1.672  0.938  1.188  1.438]
 [ 1.688  1.938  2.188  2.438  2.688  2.938  3.188  3.438  3.688  3.938  4.188  4.438  4.812  5.062  5.312  4.578  2.859  2.125  2.375  2.297  1.891  1.812  2.062  2.312]
 [ 2.562  2.812  3.062  3.312  3.562  3.812  4.062  4.312  4.688  4.938  5.188  5.109  4.703  4.625  4.875  4.141  2.422  1.688  1.938  2.188  2.438  2.688  2.938  3.188]
 [ 3.438  3.688  3.938  4.188  4.562  4.812  5.062  5.312  5.562  5.812  6.062  5.328  3.609  2.875  3.125  3.047  2.641  2.562  2.812  3.062  3.312  3.562  3.812  4.062]
 [ 4.438  4.688  4.938  5.188  5.438  5.688  5.938  5.859  5.453  5.375  5.625  4.891  3.172  2.438  2.688  2.938  3.188  3.438  3.688  3.938  4.312  4.562  4.812  5.062]
 [ 5.312  5.562  5.812  6.062  6.312  6.562  6.812  6.078  4.359  3.625  3.875  3.797  3.391  3.312  3.562  3.812  4.188  4.438  4.688  4.938  5.188  5.438  5.688  5.938]
 [ 6.188  6.438  6.688  6.609  6.203  6.125  6.375  5.641  3.922  3.188  3.438  3.688  4.062  4.312  4.562  4.812  5.062  5.312  5.562  5.812  6.062  6.312  6.562  6.812]
 [ 7.062  7.312  7.562  6.828  5.109  4.375  4.625  4.875  3.938  4.188  4.438  4.688  4.938  5.188  5.438  5.688  5.938  6.188  6.438  6.688  6.938  7.188  7.438  7.359]
 [ 6.953  6.875  7.125  7.375  3.812  4.062  4.312  4.562  4.812  5.062  5.312  5.562  5.812  6.062  6.312  6.562  6.812  7.062  7.312  7.562  7.812  8.062  8.312  8.562]
 [ 3.688  3.938  4.188  4.438  4.688  4.938  5.188  5.438  5.688  5.938  6.188  6.438  6.688  6.938  7.188  7.438  7.688  7.938  8.188  8.438  4.875  5.125  5.375  5.297]
 [ 4.891  4.812  5.062  5.312  5.562  5.812  6.062  6.312  6.562  6.812  7.062  7.312  7.562  7.812  8.062  8.312  7.375  7.625  7.875  7.141  5.422  4.688  4.938  5.188]
 [ 5.438  5.688  5.938  6.188  6.438  6.688  6.938  7.188  7.438  7.688  7.938  8.188  8.562  8.812  9.062  8.328  6.609  5.875  6.125  6.047  5.641  5.562  5.812  6.062]
 [ 6.312  6.562  6.812  7.062  7.312  7.562  7.812  8.062  8.438  8.688  8.938  8.859  8.453  8.375  8.625  7.891  6.172  5.438  5.688  5.938  6.188  6.438  6.688  6.938]
 [ 7.188  7.438  7.688  7.938  8.312  8.562  8.812  9.062  9.312  9.562  9.812  9.078  7.359  6.625  6.875  6.797  6.391  6.312  6.562  6.812  7.062  7.312  7.562  7.812]
 [ 8.188  8.438  8.688  8.938  9.188  9.438  9.688  9.609  9.203  9.125  9.375  8.641  6.922  6.188  6.438  6.688  6.938  7.188  7.438  7.688  8.062  8.312  8.562  8.812]
 [ 9.062  9.312  9.562  9.812 10.062 10.312 10.562  9.828  8.109  7.375  7.625  7.547  7.141  7.062  7.312  7.562  7.938  8.188  8.438  8.688  8.938  9.188  9.438  9.688]
 [ 9.938 10.188 10.438 10.359  9.953  9.875 10.125  9.391  7.672  6.938  7.188  7.438  7.812  8.062  8.312  8.562  8.812  9.062  9.312  9.562  9.812 10.062 10.312 10.562]
 [10.812 11.062 11.312 10.578  8.859  8.125  8.375  8.625  7.688  7.938  8.188  8.438  8.688  8.938  9.188  9.438  9.688  9.938 10.188 10.438 10.688 10.938 11.188 11.109]
 [10.703 10.625 10.875 11.125  7.562  7.812  8.062  8.312  8.562  8.812  9.062  9.312  9.562  9.812 10.062 10.312 10.562 10.812 11.062 11.641 12.547 13.125 13.375 13.625]]

Possibly the only instance when this code works is when input arrays are square and maybe this is good enough for STIS flats, ... don't know.

A Quick Fix:

The issue of incorrect indices could be fixed in the caller like this:

interp2d.expand2d(f.ravel().reshape((12, 10)), (24, 20)).ravel().reshape((20, 24))

I can't believe this was the intended way of calling this function.

mcara commented 5 years ago

Now I believe I understand the intention of the code better and most likely the single issue with the code is that it uses wrong array order indexing as confirmed by the "quick fix" presented above. An alternative algorithm entirely in Python/numpy was implemented in https://github.com/spacetelescope/drizzlepac/pull/227 thus making this package unnecessary: drizzlepac was the only package that was using stsci.imagemanip.