spacetelescope / drizzle

A package for combining dithered images into a single image
https://spacetelescope-drizzle.readthedocs.io/en/latest/
Other
51 stars 26 forks source link

Basic examples no longer behaving as expected #132

Closed dholzber closed 9 months ago

dholzber commented 10 months ago

We have been utilizing the lower level functionality of drizzle. Several of the basic examples that are used for unit testing started failing once we updated to 1.14.4. My guess is that this is due to the changes that were made in 1.14.0, as that's when the unexpected behavior started occurring.

Here's a script with 2 of the examples that now produce unexpected results, along with the output and what we expect the results to be (which I had previously computed according to the equations in the DrizzlePac handbook):

import drizzle.drizzle 
import numpy as np 

print( "VERSION:", drizzle.__version__ )

print( "Test 1" )
inPic = np.ones( ( 2, 2 ), dtype=np.float32 )
pixFrac = 1.0
scale = 1
inWeight = np.ones( inPic.shape, dtype=np.float32 )
pixMap = np.array( [ # identity
   [ [ 0., 0. ],
     [ 1., 0. ] ],
   [ [ 0., 1. ],
     [ 1., 1. ] ],
], dtype=np.float32 )

outShape = ( inPic.shape[0]*scale, inPic.shape[1]*scale )
outPic = np.zeros( outShape, dtype=np.float32 )
outCounts = np.zeros( outShape, dtype=np.int32 )
outWeight = np.zeros( outShape, dtype=np.float32 )

drizzle.cdrizzle.tdriz( inPic, inWeight, pixMap, outPic, outWeight, 
                        outCounts, pixfrac=pixFrac, scale=scale,
                        kernel="square", in_units="cps", 
                        expscale=1.0, wtscale=1.0 )

expected = np.ones( ( outShape ), dtype=np.float32 )

print( "Expected: " )
print( expected )

print( "Actual: " )
print( outPic )

print()
print( "Test 2" )
inPic = np.zeros( ( 2, 2 ), dtype=np.float32 )
pixFrac = 0.5
scale = 4
inWeight = np.ones( inPic.shape, dtype=np.float32 )
pixMap = np.array( [ # perfectly aligned with the new image
   [ [ 1.5, 1.5 ],
     [ 5.5, 1.5 ] ],
   [ [ 1.5, 5.5 ],
     [ 5.5, 5.5 ] ],
], dtype=np.float32 )

outShape = ( inPic.shape[0]*scale, inPic.shape[1]*scale )
outPic = np.ones( outShape, dtype=np.float32 )
outCounts = np.zeros( outShape, dtype=np.int32 )
outWeight = np.ones( outShape, dtype=np.float32 )

drizzle.cdrizzle.tdriz( inPic, inWeight, pixMap, outPic, outWeight, 
                        outCounts, pixfrac=pixFrac, scale=scale,
                        kernel="square", in_units="cps", 
                        expscale=1.0, wtscale=1.0 )

# per DrizzlePac equations on p. 13, consider xy = 0,0:
#
# a_xy = .25 since i_xy will overlap partially w/ 4 pixels in I_xy.
# w_xy = 1
# W_xy = 1
# i_xy = 0
# I_xy = 1
#
# W'_xy = a_xy * w_xy + W_xy
#       = .25 * 1 + 1
#       = 1.25
# I'_xy = ( a_xy * i_xy * w_xy + I_xy * W_xy ) / W'_xy
#       = ( 0 + 1 ) / 1.25
#       = 0.8
#
# so all pixels overlapping w/ the original shrunken pixels should be 0.8
expected = [[1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,],
            [1,  0.8, 0.8, 1.,  1.,  0.8, 0.8, 1., ],
            [1,  0.8, 0.8, 1.,  1.,  0.8, 0.8, 1., ],
            [1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,],
            [1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,],
            [1,  0.8, 0.8, 1.,  1.,  0.8, 0.8, 1., ],
            [1,  0.8, 0.8, 1.,  1.,  0.8, 0.8, 1., ],
            [1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,]]
expected = np.array( expected, dtype=np.float32 )

print( "Expected: " )
print( expected )

print( "Actual: " )
print( outPic )

Previous output:

VERSION: 1.13.6 # Version 1.13.7 is the last version that produced this output
Test 1
Expected: 
[[1. 1.]
 [1. 1.]]
Actual: 
[[1. 1.]
 [1. 1.]]

Test 2
Expected: 
[[1.  1.  1.  1.  1.  1.  1.  1. ]
 [1.  0.8 0.8 1.  1.  0.8 0.8 1. ]
 [1.  0.8 0.8 1.  1.  0.8 0.8 1. ]
 [1.  1.  1.  1.  1.  1.  1.  1. ]
 [1.  1.  1.  1.  1.  1.  1.  1. ]
 [1.  0.8 0.8 1.  1.  0.8 0.8 1. ]
 [1.  0.8 0.8 1.  1.  0.8 0.8 1. ]
 [1.  1.  1.  1.  1.  1.  1.  1. ]]
Actual: 
[[1.  1.  1.  1.  1.  1.  1.  1. ]
 [1.  0.8 0.8 1.  1.  0.8 0.8 1. ]
 [1.  0.8 0.8 1.  1.  0.8 0.8 1. ]
 [1.  1.  1.  1.  1.  1.  1.  1. ]
 [1.  1.  1.  1.  1.  1.  1.  1. ]
 [1.  0.8 0.8 1.  1.  0.8 0.8 1. ]
 [1.  0.8 0.8 1.  1.  0.8 0.8 1. ]
 [1.  1.  1.  1.  1.  1.  1.  1. ]]

Latest output:

VERSION: 1.14.4 # Same output since 1.14.0 
Test 1
Expected: 
[[1. 1.]
 [1. 1.]]
Actual: 
[[0. 0.]
 [0. 0.]]

Test 2
Expected: 
[[1.  1.  1.  1.  1.  1.  1.  1. ]
 [1.  0.8 0.8 1.  1.  0.8 0.8 1. ]
 [1.  0.8 0.8 1.  1.  0.8 0.8 1. ]
 [1.  1.  1.  1.  1.  1.  1.  1. ]
 [1.  1.  1.  1.  1.  1.  1.  1. ]
 [1.  0.8 0.8 1.  1.  0.8 0.8 1. ]
 [1.  0.8 0.8 1.  1.  0.8 0.8 1. ]
 [1.  1.  1.  1.  1.  1.  1.  1. ]]
Actual: 
[[1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1.]]

Please let me know if you require any additional information/clarification. Thank you! :)

mcara commented 9 months ago

Thank you for reporting this. Fixed via #137