From b850d47cd31139f211eb088fe325cca4f89798e2 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Wed, 23 Oct 2024 16:59:53 -0400 Subject: [PATCH] API re-design: simplify interface, remove I/O, make instrument agnostic (#134) API re-design: simplify interface, remove I/O, make instrument agnostic --- .coveragerc | 37 +- .flake8 | 11 +- .ruff.toml | 211 ++++ CHANGES.rst | 12 + README.rst | 175 +--- docs/conf.py | 25 +- docs/drizzle/api.rst | 6 +- docs/drizzle/user.rst | 2 + docs/exts/numfig.py | 110 +++ docs/rtd_environment.yaml | 1 + drizzle/__init__.py | 3 +- drizzle/calc_pixmap.py | 52 - drizzle/doblot.py | 80 -- drizzle/dodrizzle.py | 189 ---- drizzle/drizzle.py | 569 ----------- drizzle/resample.py | 702 ++++++++++++++ drizzle/tests/test_cdrizzle.py | 23 +- drizzle/tests/test_drizzle.py | 834 ---------------- drizzle/tests/test_file_io.py | 173 ---- drizzle/tests/test_overlap_calc.py | 4 +- drizzle/tests/test_pixmap.py | 76 -- drizzle/tests/test_resample.py | 1437 ++++++++++++++++++++++++++++ drizzle/tests/test_utils.py | 193 ++++ drizzle/util.py | 256 +---- drizzle/utils.py | 239 +++++ pyproject.toml | 131 ++- setup.py | 2 +- src/cdrizzleapi.c | 64 +- src/cdrizzlebox.c | 127 +-- src/cdrizzlemap.c | 11 +- src/cdrizzleutil.c | 1 - src/cdrizzleutil.h | 1 - src/tests/utest_cdrizzle.c | 2 - 33 files changed, 3162 insertions(+), 2597 deletions(-) create mode 100644 .ruff.toml create mode 100644 docs/exts/numfig.py delete mode 100644 drizzle/calc_pixmap.py delete mode 100644 drizzle/doblot.py delete mode 100644 drizzle/dodrizzle.py delete mode 100644 drizzle/drizzle.py create mode 100644 drizzle/resample.py delete mode 100644 drizzle/tests/test_drizzle.py delete mode 100644 drizzle/tests/test_file_io.py delete mode 100644 drizzle/tests/test_pixmap.py create mode 100644 drizzle/tests/test_resample.py create mode 100644 drizzle/tests/test_utils.py create mode 100644 drizzle/utils.py diff --git a/.coveragerc b/.coveragerc index 330eedf4..644d5542 100644 --- a/.coveragerc +++ b/.coveragerc @@ -1,21 +1,32 @@ [run] -source = {packagename} -omit = */tests/* +source = drizzle + +omit = + drizzle/__init__* + drizzle/**/conftest.py + drizzle/**/setup* + drizzle/**/tests/* + drizzle/version* + drizzle/_version.py + */drizzle/__init__* + */drizzle/**/conftest.py + */drizzle/**/setup* + */drizzle/**/tests/* [report] exclude_lines = - # Have to re-enable the standard pragma - pragma: no cover + # Have to re-enable the standard pragma + pragma: no cover - # Don't complain about packages we have installed - except ImportError + # Don't complain about packages we have installed + except ImportError - # Don't complain if tests don't hit assertions - raise AssertionError - raise NotImplementedError + # Don't complain if tests don't hit defensive assertion code: + raise AssertionError + raise NotImplementedError - # Don't complain about script hooks - def main\(.*\): + # Don't complain about script hooks + 'def main(.*):' - # Ignore branches that don't pertain to this version of Python - pragma: py{ignore_python_version} + # Ignore branches that don't pertain to this version of Python + pragma: py{ignore_python_version} diff --git a/.flake8 b/.flake8 index 2afc4819..2e302c1e 100644 --- a/.flake8 +++ b/.flake8 @@ -2,9 +2,12 @@ max-line-length = 100 select = F, W, E, C ignore = - E231, # Missing whitespace after ',', ';', or ':' - W503, # Line break occurred before a binary operator - W504, # Line break occurred after a binary operator + # Line break occurred before a binary operator + W503, + # Line break occurred after a binary operator + W504, exclude = docs, - .eggs + build + .eggs, + .tox, diff --git a/.ruff.toml b/.ruff.toml new file mode 100644 index 00000000..8f35bfaa --- /dev/null +++ b/.ruff.toml @@ -0,0 +1,211 @@ +extend = "pyproject.toml" +lint.ignore = [ + # flake8-annotations (ANN) : static typing + "ANN001", # Function argument without type annotation + "ANN003", # `**kwargs` without type annotation + "ANN202", # Private function without return type annotation + "ANN401", # Use of `Any` type + + # flake8-unused-arguments (ARG) + "ARG001", # unused-function-argument + "ARG002", # unused-method-argument + "ARG003", # unused-class-method-argument + "ARG004", # unused-static-method-argument + "ARG005", # unused-lambda-argument + + # flake8-bugbear (B) + "B006", # MutableArgumentDefault + "B007", # UnusedLoopControlVariable + "B023", # FunctionUsesLoopVariable + "B028", # No-explicit-stacklevel + "B904", # RaiseWithoutFromInsideExcept + "B905", # ZipWithoutExplicitStrict # requires 3.10+ + + # flake8-blind-except (BLE) + "BLE001", # blind-except + + # mccabe (C90) : code complexity + # TODO: configure maximum allowed complexity. + "C901", # McCabeComplexity + + # pydocstyle (D) + # Missing Docstrings + "D100", # undocumented-public-module + "D101", # undocumented-public-class + "D103", # undocumented-public-function + "D104", # undocumented-public-package + "D205", # blank-line-after-summary + # Quotes Issues + "D300", # triple-single-quotes + "D301", # escape-sequence-in-docstring + # Docstring Content Issues + "D403", # first-line-capitalized + "D404", # docstring-starts-with-this + "D401", # non-imperative-mood. + "D414", # empty-docstring-section + "D419", # docstring is empty + + # flake8-datetimez (DTZ) + "DTZ001", # call-datetime-without-tzinfo + "DTZ005", # call-datetime-now-without-tzinfo + + # pycodestyle (E, W) + "E501", # line-too-long + "E721", # type-comparison + "E731", # lambda-assignment + + # flake8-errmsg (EM) : nicer error tracebacks + "EM101", # raw-string-in-exception + "EM102", # f-string-in-exception + "EM103", # dot-format-in-exception + + # eradicate (ERA) + # NOTE: be careful that developer notes are kept. + "ERA001", # commented-out-code + + # flake8-executable (EXE) + "EXE002", # shebang-missing-executable-file + + # Pyflakes (F) + "F841", # unused-variable + + # flake8-boolean-trap (FBT) : boolean flags should be kwargs, not args + # NOTE: a good thing to fix, but changes API. + "FBT001", # boolean-positional-arg-in-function-definition + "FBT002", # boolean-default-value-in-function-definition + "FBT003", # boolean-positional-value-in-function-call + + # flake8-fixme (FIX) + "FIX001", # Line contains FIXME. this should be fixed or at least FIXME replaced with TODO + "FIX004", # Line contains HACK. replace HACK with NOTE. + + # pep8-naming (N) + # NOTE: some of these can/should be fixed, but this changes the API. + "N801", # invalid-class-name + "N802", # invalid-function-name + "N803", # invalid-argument-name + "N804", # invalid-first-argument-name-for-class-method + "N805", # invalid-first-argument-name-for-method + "N807", # dunder-function-name + "N813", # camelcase-imported-as-lowercase + "N815", # mixed-case-variable-in-class-scope + "N816", # mixed-case-variable-in-global-scope + "N818", # error-suffix-on-exception-name + + # NumPy-specific rules (NPY) + "NPY002", # Replace legacy `np.random.rand` call with `np.random.Generator` (2023-05-03) + + # Perflint (PERF) + "PERF203", # `try`-`except` within a loop incurs performance overhead + "PERF401", # Use a list comprehension to create a transformed list + + # Pylint (PLC, PLE, PLR, PLW) + "PLE0101", # return-in-init + "PLR0124", # Name compared with itself + "PLR0402", # ConsiderUsingFromImport + "PLR0912", # too-many-branches + "PLR0913", # too-many-args + "PLR0915", # too-many-statements + "PLR1714", # Consider merging multiple comparisons + "PLR2004", # MagicValueComparison + "PLR5501", # collapsible-else-if + "PLW0603", # global-statement + "PLW2901", # redefined-loop-name + + # flake8-pytest-style (PT) + "PT001", # pytest-fixture-incorrect-parentheses-style + "PT003", # pytest-extraneous-scope-function + "PT004", # pytest-missing-fixture-name-underscore + "PT006", # pytest-parametrize-names-wrong-type + "PT007", # pytest-parametrize-values-wrong-type + "PT011", # pytest-raises-too-broad + "PT012", # pytest-raises-with-multiple-statements + "PT017", # pytest-assert-in-exceptinstead + "PT018", # pytest-composite-assertion + "PT022", # pytest-useless-yield-fixture + "PT023", # pytest-incorrect-mark-parentheses-style + + # flake8-use-pathlib (PTH) + "PTH100", # os-path-abspath + "PTH101", # os-chmod + "PTH102", # os-mkdir + "PTH103", # os-makedirs + "PTH104", # os-rename + "PTH107", # os-remove + "PTH108", # os-unlink + "PTH109", # os-getcwd + "PTH110", # os-path-exists + "PTH111", # os-path-expanduser + "PTH112", # os-path-isdir + "PTH113", # os-path-isfile + "PTH118", # os-path-join + "PTH119", # os-path-basename + "PTH120", # os-path-dirname + "PTH122", # os-path-splitext + "PTH123", # builtin-open + "PTH202", # `os.path.getsize` should be replaced by `Path.stat().st_size` + "PTH207", # Replace `glob` with `Path.glob` or `Path.rglob` + + # flake8-return (RET) + "RET501", # unnecessary-return-none + "RET502", # implicit-return-value + "RET503", # implicit-return + "RET504", # unnecessary-assign + "RET507", # superfluous-else-continue + + # flake8-raise (RSE) + "RSE102", # unnecessary-paren-on-raise-exception + + # Ruff-specific rules (RUF) + "RUF001", # ambiguous-unicode-character-string + "RUF002", # ambiguous-unicode-character-docstring + "RUF010", # use conversion in f-string + "RUF012", # Mutable class attributes should be annotated with `typing.ClassVar` + + # flake8-bandit (S) + "S101", # Use of `assert` detected + "S105", # hardcoded-password-string + "S110", # try-except-pass + "S112", # try-except-continue + "S301", # suspicious-pickle-usage + "S307", # Use of possibly insecure function; consider using `ast.literal_eval` + "S311", # suspicious-non-cryptographic-randomness + "S324", # hashlib-insecure-hash-function + "S506", # UnsafeYAMLLoad + "S310", # Suspicious-url-open-usage + "S603", # `subprocess` call: check for execution of untrusted input + "S607", # Starting a process with a partial executable path + + # flake8-simplify (SIM) + "SIM102", # NestedIfStatements + "SIM105", # UseContextlibSuppress + "SIM108", # UseTernaryOperator + "SIM114", # if-with-same-arms + "SIM115", # OpenFileWithContextHandler + "SIM117", # MultipleWithStatements + "SIM118", # KeyInDict + "SIM201", # NegateEqualOp + "SIM300", # yoda condition + + # flake8-print (T20) + "T201", # PrintUsed + + # flake8-todos (TD) + "TD001", # Invalid TODO tag + "TD003", # Missing issue link on the line following this TODO + "TD004", # Missing colon in TODO + "TD007", # Missing space after colon in TODO + + # tryceratops (TRY) + "TRY002", # raise-vanilla-class + "TRY003", # raise-vanilla-args + "TRY004", # prefer-type-error + "TRY201", # verbose-raise + "TRY301", # raise-within-try + + # flake8-quotes (Q) + "Q000", # use double quotes +] +lint.unfixable = [ + "E711" # NoneComparison. Hard to fix b/c numpy has it's own None. +] diff --git a/CHANGES.rst b/CHANGES.rst index dded7ae5..046e02ec 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -18,6 +18,18 @@ Release Notes - build wheels with Numpy 2.0 release candidate [#149] +2.0.0 (Unreleased) +================== + +- Backward incompatible major re-design of API. [#134] + +- Bug fix: exposure time was undefined when in_units were not cps. [#134] + +- BUG FIX: ``cdrizzle.tdriz`` signature was out of date. [#134] + +- Removed support for the ``'tophat'`` kernel. [#134] + + 1.15.1 (2024-03-05) =================== diff --git a/README.rst b/README.rst index cb80daa2..d4e821f9 100644 --- a/README.rst +++ b/README.rst @@ -137,177 +137,8 @@ missing information. The blot methods perform the inverse operation of drizzle. That is, blotting performs the inverse mapping to transform the dithered median image back into the coordinate system of the original input image. Blotting is primarily used -for identifying cosmic rays in the original image. Like the original drizzle -task, blot requires the user to provide the world coordinate system (WCS) -transformations as inputs. +for identifying cosmic rays in the original image. Blot requires the user +to provide the world coordinate system (WCS)-based transformations in the +form of a pixel map array as input. .. [Driz2012] Gonzaga, S., Hack, W., Fruchter, A., Mack, J., eds. 2012, The DrizzlePac Handbook. (Baltimore, STScI) - - -The Drizzle Library -------------------- - -The Drizzle library is object-oriented and you use it by first creating an object of -the ``Drizzle`` class. To create a new Drizzle output image, supply an Astropy -WCS object representing the coordinate system of the output image. -The other parameters are the linear pixel dimension described in the previous -section, the drizzle kernel used, how each input image is scaled (by exposure -time or time squared), and the pixel value set in the output image where the -input images do not overlap. - -After creating a Drizzle object, you add one or more images by calling the -``add_fits_file`` method. The arguments are the name of the FITS file containing -the input image and optionally the name of a FITS file containing the pixel -weighting. Both file names can be followed by an extension name or number in -square brackets. Optionally you can pass the name of the header keywords -containing the exposure time and units. Two units are understood: counts and -cps (counts per second). - -The following function is a demonstration of how you can create a new output -image:: - - def drizzle_demo_one(reference, outfile, infiles): - """ - First demonstration of drizzle - - Parameters - ========== - reference - A file containing the wcs of the output image - - outfile - The name of the output image - - infiles - The names of the input images to be combined - """ - # Get the WCS for the output image - hdulist = fits.open(reference) - reference_wcs = wcs.WCS(hdulist[1].header) - - # Initialize the output with the WCS - driz = drizzle.drizzle.Drizzle(outwcs=reference_wcs) - - # Combine the input images into on drizzle image - for infile in infiles: - driz.add_fits_file(infile) - - # Write the drizzled image out - driz.write(outfile) - -Optionally you can supply the input and weight images as Numpy arrays by using -the ``add_image`` method. If you use this method, you must supply the extra -information that would otherwise be read from the FITS image: The WCS -of the input image, the exposure time, and image units. - -Here is an example of how you would call ``add_image``:: - - def drizzle_demo_two(reference, outfile, infiles): - """ - Demonstration of drizzle with add image. - - Parameters - ========== - reference - A file containing the wcs of the output image. - - outfile - The name of the output image. - - infiles - The names of the input images to be combined. - """ - # Get the WCS for the output image - reflist = fits.open(reference) - reference_wcs = wcs.WCS(reflist[1].header) - - # Initialize the output with the WCS - driz = drizzle.drizzle.Drizzle(outwcs=reference_wcs) - - # Combine the input images into on drizzle image - for infile in infiles: - # Open the file and read the image and wcs - # This is a contrived example, we would not do this - # unless the data came from another source - # than a FITS file - imlist = fits.open(reference) - image = imlist[1].data - image_wcs = wcs.WCS(imlist[1].header) - driz.add_image(image, image_wcs) - - # Write the drizzled image out - driz.write(outfile) - -After combining all the input images, you write the output image into a FITS -file with the ``write`` method. You must pass the name of the output image and -optionally the units. You can also supply a set of header cards to be added -to the primary header of the output FITS file. - -You can also add more images to an existing Drizzle output file by creating -a new Drizzle object and passing the existing output file name as the new -object is created. In that case the output WCS and all -other parameters are read from the file. - -Here is a demonstration of adding additional input images to a drizzled image:: - - def drizzle_demo_three(outfile, infiles): - """ - Demonstration of drizzle and adding to an existing output. - - Parameters - ========== - outfile - Name of output image that new files will be appended to. - - infiles - The names of the input images to be added. - """ - # Re-open the output file - driz = drizzle.drizzle.Drizzle(infile=outfile) - - # Add the input images to the existing output image - for infile in infiles: - driz.add_fits_file(infile) - - # Write the modified drizzled image out - driz.write(outfile) - -You can use the methods ``blot_fits_file`` and ``blot_image`` to transform the drizzled -output image into another WCS. Most usually this is the -coordinates of one of the input images and is used to identify cosmic rays or -other defects. The two methods ``blot_fits_file`` and ``blot_image`` allow you to -retrieve the WCS from the FITS file header or input it directly. -The optional parameter ``interp`` allows you to selct the method used to resample -the pixels on the new grid, and ``sincscl`` is used to scale the sinc function if one -of the sinc interpolation methods is used. This function demonstrates how both -methods are called:: - - def drizzle_demo_four(outfile, blotfile): - """ - Demonstration of blot methods. - - Parameters - ========== - outfile - Name of output image that will be converted. - - blotfile - Name of image containing wcs to be transformed to. - """ - # Open drizzle using the output file - # Transform it to another coordinate system - driz = drizzle.drizzle.Drizzle(infile=outfile) - driz.blot_fits_file(blotfile) - driz.write(outfile) - - # Read the WCS and transform using it instead - # This is a contrived example - blotlist = fits.open(blotfile) - blot_wcs = wcs.WCS(blotlist[1].header) - driz = drizzle.drizzle.Drizzle(infile=outfile) - driz.blot_image(blot_wcs) - driz.write(outfile) - -The lower level function ``dodrizzle`` is present for backwards compatibility with -the existing STScI DrizzlePac code and should not be used unless you are also -concerned with this compatibility. diff --git a/docs/conf.py b/docs/conf.py index 5c49cbbc..6382db5d 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -29,25 +29,12 @@ import sys from pathlib import Path - -# Get configuration information from pyproject.toml -# try: -# import tomllib -# except ImportError: -# import toml as tomllib -# -# sys.path.insert(1, '..') -# -# setup_cfg = tomllib.load( -# path.join(path.dirname(__file__), '..', 'pyproject.toml') -# )['project'] - if sys.version_info < (3, 11): import tomli as tomllib else: import tomllib -#sys.path.insert(0, os.path.abspath('../')) +sys.path.insert(0, str(Path("exts/").resolve())) def find_mod_objs_patched(*args, **kwargs): @@ -73,12 +60,22 @@ def setup(app): extensions = [ + 'numfig', 'sphinx_automodapi.automodapi', 'sphinx.ext.autodoc', 'sphinx.ext.autosummary', 'sphinx.ext.doctest', 'sphinx.ext.intersphinx', + 'pytest_doctestplus.sphinx.doctestplus', + 'sphinx.ext.todo', + 'sphinx.ext.inheritance_diagram', + 'sphinx.ext.viewcode', + 'sphinx.ext.napoleon', + 'sphinx_automodapi.automodsumm', + 'sphinx_automodapi.autodoc_enhancements', + 'sphinx_automodapi.smart_resolver', ] + # -- General configuration ---------------------------------------------------- master_doc = 'index' diff --git a/docs/drizzle/api.rst b/docs/drizzle/api.rst index 8eeb16a8..9a60ef92 100644 --- a/docs/drizzle/api.rst +++ b/docs/drizzle/api.rst @@ -1,7 +1,5 @@ API Documentation ================= -.. automodapi:: drizzle.drizzle -.. automodapi:: drizzle.dodrizzle -.. automodapi:: drizzle.doblot -.. automodapi:: drizzle.calc_pixmap +.. automodapi:: drizzle.resample +.. automodapi:: drizzle.utils diff --git a/docs/drizzle/user.rst b/docs/drizzle/user.rst index a6210d3d..90291d6e 100644 --- a/docs/drizzle/user.rst +++ b/docs/drizzle/user.rst @@ -1 +1,3 @@ +.. _main-user-doc: + .. include:: ../../README.rst diff --git a/docs/exts/numfig.py b/docs/exts/numfig.py new file mode 100644 index 00000000..dfd871e8 --- /dev/null +++ b/docs/exts/numfig.py @@ -0,0 +1,110 @@ +from docutils.nodes import figure, caption, Text, reference, raw, SkipNode +from sphinx.roles import XRefRole + + +# Element classes + +class page_ref(reference): + pass + + +class num_ref(reference): + pass + + +# Visit/depart functions + +def skip_page_ref(self, node): + raise SkipNode + + +def latex_visit_page_ref(self, node): + self.body.append("\\pageref{%s:%s}" % (node['refdoc'], node['reftarget'])) + raise SkipNode + + +def latex_visit_num_ref(self, node): + fields = node['reftarget'].split('#') + if len(fields) > 1: + label, target = fields + ref_link = '%s:%s' % (node['refdoc'], target) + latex = "\\hyperref[%s]{%s \\ref*{%s}}" % (ref_link, label, ref_link) + self.body.append(latex) + else: + self.body.append('\\ref{%s:%s}' % (node['refdoc'], fields[0])) + + raise SkipNode + + +def doctree_read(app, doctree): + # first generate figure numbers for each figure + env = app.builder.env + figid_docname_map = getattr(env, 'figid_docname_map', {}) + + for figure_info in doctree.traverse(figure): + for id in figure_info['ids']: + figid_docname_map[id] = env.docname + + env.figid_docname_map = figid_docname_map + + +def doctree_resolved(app, doctree, docname): + i = 1 + figids = {} + for figure_info in doctree.traverse(figure): + if app.builder.name != 'latex' and app.config.number_figures: + for cap in figure_info.traverse(caption): + cap[0] = Text("%s %d: %s" % (app.config.figure_caption_prefix, i, cap[0])) + + for id in figure_info['ids']: + figids[id] = i + + i += 1 + + # replace numfig nodes with links + if app.builder.name != 'latex': + for ref_info in doctree.traverse(num_ref): + if '#' in ref_info['reftarget']: + label, target = ref_info['reftarget'].split('#') + labelfmt = label + " %d" + else: + labelfmt = '%d' + target = ref_info['reftarget'] + + if target not in figids: + continue + + if app.builder.name == 'html': + target_doc = app.builder.env.figid_docname_map[target] + link = "%s#%s" % (app.builder.get_relative_uri(docname, target_doc), + target) + html = '%s' % (link, labelfmt % (figids[target])) + ref_info.replace_self(raw(html, html, format='html')) + else: + ref_info.replace_self(Text(labelfmt % (figids[target]))) + + +def clean_env(app): + app.builder.env.i = 1 + app.builder.env.figid_docname_map = {} + + +def setup(app): + app.add_config_value('number_figures', True, True) + app.add_config_value('figure_caption_prefix', "Figure", True) + + app.add_node(page_ref, + text=(skip_page_ref, None), + html=(skip_page_ref, None), + latex=(latex_visit_page_ref, None)) + + app.add_role('page', XRefRole(nodeclass=page_ref)) + + app.add_node(num_ref, + latex=(latex_visit_num_ref, None)) + + app.add_role('num', XRefRole(nodeclass=num_ref)) + + app.connect("builder-inited", clean_env) + app.connect('doctree-read', doctree_read) + app.connect('doctree-resolved', doctree_resolved) diff --git a/docs/rtd_environment.yaml b/docs/rtd_environment.yaml index 161fef38..d1129f22 100644 --- a/docs/rtd_environment.yaml +++ b/docs/rtd_environment.yaml @@ -6,4 +6,5 @@ dependencies: - python=3.11 - pip - graphviz + - pytest-doctestplus - sphinx_rtd_theme>1.2.0 diff --git a/drizzle/__init__.py b/drizzle/__init__.py index 9418e47b..c4e910fe 100644 --- a/drizzle/__init__.py +++ b/drizzle/__init__.py @@ -3,9 +3,10 @@ """ from importlib.metadata import PackageNotFoundError, version +from drizzle import cdrizzle, resample, utils # noqa: F401 try: __version__ = version(__name__) except PackageNotFoundError: # package is not installed - __version__ = 'unknown' + __version__ = '' diff --git a/drizzle/calc_pixmap.py b/drizzle/calc_pixmap.py deleted file mode 100644 index 911c5464..00000000 --- a/drizzle/calc_pixmap.py +++ /dev/null @@ -1,52 +0,0 @@ -import numpy as np - - -def calc_pixmap(first_wcs, second_wcs): - """ - Calculate a mapping between the pixels of two images. - - Parameters - ---------- - - first_wcs : wcs - A WCS object representing the coordinate system you are - converting from - - seond_wcs : wcs - A WCS object representing the coordinate system you are - converting to - - Returns - ------- - - A three dimensional array representing the transformation between - the two. The last dimension is of length two and contains the x and - y coordinates of a pixel center, repectively. The other two coordinates - correspond to the two coordinates of the image the first WCS is from. - """ - - first_naxis1, first_naxis2 = first_wcs.pixel_shape - - # We add one to the pixel co-ordinates before the transformation and subtract - # it afterwards because wcs co-ordinates are one based, while pixel co-ordinates - # are zero based, The result is the final values in pixmap give a tranformation - # between the pixel co-ordinates in the first image to pixel co-ordinates in the - # co-ordinate system of the second. - - one = np.ones(2, dtype='float64') - idxmap = np.indices((first_naxis1, first_naxis2), dtype='float64') - idxmap = idxmap.transpose() + one - - idxmap = idxmap.reshape(first_naxis2 * first_naxis1, 2) - - worldmap = first_wcs.all_pix2world(idxmap, 1) - - if second_wcs.sip is None: - pixmap = second_wcs.wcs_world2pix(worldmap, 1) - else: - pixmap = second_wcs.all_world2pix(worldmap, 1) - - pixmap = pixmap.reshape(first_naxis2, first_naxis1, 2) - pixmap = pixmap - one - - return pixmap diff --git a/drizzle/doblot.py b/drizzle/doblot.py deleted file mode 100644 index 818ebedd..00000000 --- a/drizzle/doblot.py +++ /dev/null @@ -1,80 +0,0 @@ -""" -STScI Python compatable blot module -""" -import numpy as np - -from drizzle import calc_pixmap -from drizzle import cdrizzle - - -def doblot(source, source_wcs, blot_wcs, exptime, coeffs=True, - interp='poly5', sinscl=1.0, stepsize=10, wcsmap=None): - """ - Low level routine for performing the 'blot' operation. - - Create a single blotted image from a single source image. The - interface is compatible with STScI code. All distortion information - is assumed to be included in the WCS specification of the 'output' - blotted image given in 'blot_wcs'. - - Parameters - ---------- - - source : 2d array - Input numpy array of the source image in units of 'cps'. - - source_wcs : wcs - The source image WCS. - - blot_wcs : wcs - The blotted image WCS. The WCS that the source image will be - resampled to. - - exptime : float - The exposure time of the input image. - - interp : str, optional - The type of interpolation used in the blotting. The - possible values are "nearest" (nearest neighbor interpolation), - "linear" (bilinear interpolation), "poly3" (cubic polynomial - interpolation), "poly5" (quintic polynomial interpolation), - "sinc" (sinc interpolation), "lan3" (3rd order Lanczos - interpolation), and "lan5" (5th order Lanczos interpolation). - - sincscl : float, optional - The scaling factor for sinc interpolation. - - Returns - ------- - - A 2d numpy array with the blotted image - - Other Parameters - ---------------- - - coeffs : bool, optional - Not used. Only kept for backwards compatibility. - - stepsize : float, optional - Was used when input to output mapping was computed - internally. Is no longer used and only here for backwards compatibility. - - wcsmap : function, optional - Was used when input to output mapping was computed - internally. Is no longer used and only here for backwards compatibility. - """ - _outsci = np.zeros(blot_wcs.pixel_shape[::-1], dtype=np.float32) - - # compute the undistorted 'natural' plate scale - blot_wcs.sip = None - blot_wcs.cpdis1 = None - blot_wcs.cpdis2 = None - blot_wcs.det2im = None - - pixmap = calc_pixmap.calc_pixmap(blot_wcs, source_wcs) - pix_ratio = source_wcs.pscale / blot_wcs.pscale - - cdrizzle.tblot(source, pixmap, _outsci, scale=pix_ratio, kscale=1.0, - interp=interp, exptime=exptime, misval=0.0, sinscl=sinscl) - - return _outsci diff --git a/drizzle/dodrizzle.py b/drizzle/dodrizzle.py deleted file mode 100644 index 7e19c19a..00000000 --- a/drizzle/dodrizzle.py +++ /dev/null @@ -1,189 +0,0 @@ -""" -STScI Python compatable drizzle module -""" -import numpy as np - -from . import util -from . import calc_pixmap -from . import cdrizzle - - -def dodrizzle(insci, input_wcs, inwht, - output_wcs, outsci, outwht, outcon, - expin, in_units, wt_scl, - wcslin_pscale=1.0, uniqid=1, - xmin=0, xmax=0, ymin=0, ymax=0, - pixfrac=1.0, kernel='square', fillval="INDEF"): - """ - Low level routine for performing 'drizzle' operation.on one image. - - The interface is compatible with STScI code. All images are Python - ndarrays, instead of filenames. File handling (input and output) is - performed by the calling routine. - - Parameters - ---------- - - insci : 2d array - A 2d numpy array containing the input image to be drizzled. - it is an error to not supply an image. - - input_wcs : 2d array - The world coordinate system of the input image. - - inwht : 2d array - A 2d numpy array containing the pixel by pixel weighting. - Must have the same dimensions as insci. If none is supplied, - the weghting is set to one. - - output_wcs : wcs - The world coordinate system of the output image. - - outsci : 2d array - A 2d numpy array containing the output image produced by - drizzling. On the first call it should be set to zero. - Subsequent calls it will hold the intermediate results - - outwht : 2d array - A 2d numpy array containing the output counts. On the first - call it should be set to zero. On subsequent calls it will - hold the intermediate results. - - outcon : 2d or 3d array, optional - A 2d or 3d numpy array holding a bitmap of which image was an input - for each output pixel. Should be integer zero on first call. - Subsequent calls hold intermediate results. - - expin : float - The exposure time of the input image, a positive number. The - exposure time is used to scale the image if the units are counts. - - in_units : str - The units of the input image. The units can either be "counts" - or "cps" (counts per second.) - - wt_scl : float - A scaling factor applied to the pixel by pixel weighting. - - wcslin_pscale : float, optional - The pixel scale of the input image. Conceptually, this is the - linear dimension of a side of a pixel in the input image, but it - is not limited to this and can be set to change how the drizzling - algorithm operates. - - uniqid : int, optional - The id number of the input image. Should be one the first time - this function is called and incremented by one on each subsequent - call. - - xmin : float, optional - This and the following three parameters set a bounding rectangle - on the input image. Only pixels on the input image inside this - rectangle will have their flux added to the output image. Xmin - sets the minimum value of the x dimension. The x dimension is the - dimension that varies quickest on the image. If the value is zero, - no minimum will be set in the x dimension. All four parameters are - zero based, counting starts at zero. - - xmax : float, optional - Sets the maximum value of the x dimension on the bounding box - of the input image. If the value is zero, no maximum will - be set in the x dimension, the full x dimension of the output - image is the bounding box. - - ymin : float, optional - Sets the minimum value in the y dimension on the bounding box. The - y dimension varies less rapidly than the x and represents the line - index on the input image. If the value is zero, no minimum will be - set in the y dimension. - - ymax : float, optional - Sets the maximum value in the y dimension. If the value is zero, no - maximum will be set in the y dimension, the full x dimension - of the output image is the bounding box. - - pixfrac : float, optional - The fraction of a pixel that the pixel flux is confined to. The - default value of 1 has the pixel flux evenly spread across the image. - A value of 0.5 confines it to half a pixel in the linear dimension, - so the flux is confined to a quarter of the pixel area when the square - kernel is used. - - kernel: str, optional - The name of the kernel used to combine the input. The choice of - kernel controls the distribution of flux over the kernel. The kernel - names are: "square", "turbo", "point", "gaussian", "lanczos2", - and "lanczos3". - - .. warning:: - The "gaussian" and "lanczos2/3" kernels **DO NOT** conserve flux. - - fillval: str, optional - The value a pixel is set to in the output if the input image does - not overlap it. The default value of INDEF does not set a value. - - Returns - ------- - A tuple with three values: a version string, the number of pixels - on the input image that do not overlap the output image, and the - number of complete lines on the input image that do not overlap the - output input image. - - """ - - # Ensure that the fillval parameter gets properly interpreted - # for use with tdriz - if util.is_blank(fillval): - fillval = 'INDEF' - else: - fillval = str(fillval) - - if in_units == 'cps': - expscale = 1.0 - else: - expscale = expin - - # Add input weight image if it was not passed in - - if (insci.dtype > np.float32): - insci = insci.astype(np.float32) - - if inwht is None: - inwht = np.ones_like(insci) - - # Compute what plane of the context image this input would - # correspond to: - planeid = int((uniqid - 1) / 32) - - # Check if the context image has this many planes - if outcon.ndim == 3: - nplanes = outcon.shape[0] - elif outcon.ndim == 2: - nplanes = 1 - else: - nplanes = 0 - - if nplanes <= planeid: - raise IndexError("Not enough planes in drizzle context image") - - # Alias context image to the requested plane if 3d - if outcon.ndim == 3: - outcon = outcon[planeid] - - pix_ratio = output_wcs.pscale / wcslin_pscale - - # Compute the mapping between the input and output pixel coordinates - pixmap = calc_pixmap.calc_pixmap(input_wcs, output_wcs) - - # - # Call 'drizzle' to perform image combination - # This call to 'cdriz.tdriz' uses the new C syntax - # - _vers, nmiss, nskip = cdrizzle.tdriz( - insci, inwht, pixmap, outsci, outwht, outcon, - uniqid=uniqid, xmin=xmin, xmax=xmax, - ymin=ymin, ymax=ymax, scale=pix_ratio, pixfrac=pixfrac, - kernel=kernel, in_units=in_units, expscale=expscale, - wtscale=wt_scl, fillstr=fillval) - - return _vers, nmiss, nskip diff --git a/drizzle/drizzle.py b/drizzle/drizzle.py deleted file mode 100644 index 37c8c3fd..00000000 --- a/drizzle/drizzle.py +++ /dev/null @@ -1,569 +0,0 @@ -""" -The `drizzle` module defines the `Drizzle` class, for combining input -images into a single output image using the drizzle algorithm. -""" -import os -import os.path - -import numpy as np -from astropy import wcs -from astropy.io import fits - -from . import util -from . import doblot -from . import dodrizzle - - -class Drizzle(object): - """ - Combine images using the drizzle algorithm - """ - def __init__(self, infile="", outwcs=None, - wt_scl="exptime", pixfrac=1.0, kernel="square", - fillval="INDEF"): - """ - Create a new Drizzle output object and set the drizzle parameters. - - All parameters are optional, but either infile or outwcs must be supplied. - If infile initializes the object from a file written after a - previous run of drizzle. Results from the previous run will be combined - with new results. The value passed in outwcs will be ignored. If infile is - not set, outwcs will be used to initilize a new run of drizzle. - - Parameters - ---------- - - infile : str, optional - A fits file containing results from a previous run. The three - extensions SCI, WHT, and CTX contain the combined image, total counts - and image id bitmap, repectively. The WCS of the combined image is - also read from the SCI extension. - - outwcs : wcs, optional - The world coordinate system (WCS) of the combined image. This - parameter must be present if no input file is given and is ignored if - one is. - - wt_scl : str, optional - How each input image should be scaled. The choices are `exptime` - which scales each image by its exposure time, `expsq` which scales - each image by the exposure time squared, or an empty string, which - allows each input image to be scaled individually. - - pixfrac : float, optional - The fraction of a pixel that the pixel flux is confined to. The - default value of 1 has the pixel flux evenly spread across the image. - A value of 0.5 confines it to half a pixel in the linear dimension, - so the flux is confined to a quarter of the pixel area when the square - kernel is used. - - kernel : str, optional - The name of the kernel used to combine the inputs. The choice of - kernel controls the distribution of flux over the kernel. The kernel - names are: "square", "turbo", "point", "gaussian", "lanczos2", - and "lanczos3". - - .. warning:: - The "gaussian" and "lanczos2/3" kernels **DO NOT** conserve flux. - - fillval : str, otional - The value a pixel is set to in the output if the input image does - not overlap it. The default value of INDEF does not set a value. - """ - - # Initialize the object fields - - self.outsci = None - self.outwht = None - self.outcon = None - - self.outexptime = 0.0 - self.uniqid = 0 - - self.outwcs = outwcs - self.wt_scl = wt_scl - self.kernel = kernel - self.fillval = fillval - self.pixfrac = float(pixfrac) - - self.sciext = "SCI" - self.whtext = "WHT" - self.ctxext = "CTX" - - out_units = "cps" - - if not util.is_blank(infile): - if os.path.exists(infile): - handle = fits.open(infile) - - # Read parameters from image header - self.outexptime = util.get_keyword(handle, "DRIZEXPT", default=0.0) - self.uniqid = util.get_keyword(handle, "NDRIZIM", default=0) - - self.sciext = util.get_keyword(handle, "DRIZOUDA", default="SCI") - self.whtext = util.get_keyword(handle, "DRIZOUWE", default="WHT") - self.ctxext = util.get_keyword(handle, "DRIZOUCO", default="CTX") - - self.wt_scl = util.get_keyword(handle, "DRIZWTSC", default=wt_scl) - self.kernel = util.get_keyword(handle, "DRIZKERN", default=kernel) - self.fillval = util.get_keyword(handle, "DRIZFVAL", default=fillval) - self.pixfrac = float(util.get_keyword(handle, - "DRIZPIXF", default=pixfrac)) - - out_units = util.get_keyword(handle, "DRIZOUUN", default="cps") - - try: - hdu = handle[self.sciext] - self.outsci = hdu.data.copy().astype(np.float32) - self.outwcs = wcs.WCS(hdu.header, fobj=handle) - except KeyError: - pass - - try: - hdu = handle[self.whtext] - self.outwht = hdu.data.copy().astype(np.float32) - except KeyError: - pass - - try: - hdu = handle[self.ctxext] - self.outcon = hdu.data.copy().astype(np.int32) - if self.outcon.ndim == 2: - self.outcon = np.reshape(self.outcon, (1, - self.outcon.shape[0], - self.outcon.shape[1])) - - elif self.outcon.ndim == 3: - pass - - else: - msg = ("Drizzle context image has wrong dimensions: " + - infile) - raise ValueError(msg) - - except KeyError: - pass - - handle.close() - - # Check field values - - if self.outwcs: - util.set_pscale(self.outwcs) - else: - raise ValueError("Either an existing file or wcs must be supplied to Drizzle") - - if util.is_blank(self.wt_scl): - self.wt_scl = '' - elif self.wt_scl != "exptime" and self.wt_scl != "expsq": - raise ValueError("Illegal value for wt_scl: %s" % out_units) - - if out_units == "counts": - np.divide(self.outsci, self.outexptime, self.outsci) - elif out_units != "cps": - raise ValueError("Illegal value for wt_scl: %s" % out_units) - - # Initialize images if not read from a file - outwcs_naxis1, outwcs_naxis2 = self.outwcs.pixel_shape - if self.outsci is None: - self.outsci = np.zeros(self.outwcs.pixel_shape[::-1], - dtype=np.float32) - - if self.outwht is None: - self.outwht = np.zeros(self.outwcs.pixel_shape[::-1], - dtype=np.float32) - if self.outcon is None: - self.outcon = np.zeros((1, outwcs_naxis2, outwcs_naxis1), - dtype=np.int32) - - def add_fits_file(self, infile, inweight="", - xmin=0, xmax=0, ymin=0, ymax=0, - unitkey="", expkey="", wt_scl=1.0): - """ - Combine a fits file with the output drizzled image. - - Parameters - ---------- - - infile : str - The name of the fits file, possibly including an extension. - - inweight : str, otional - The name of a file containing a pixel by pixel weighting - of the input data. If it is not set, an array will be generated - where all values are set to one. - - xmin : float, otional - This and the following three parameters set a bounding rectangle - on the output image. Only pixels on the output image inside this - rectangle will have their flux updated. Xmin sets the minimum value - of the x dimension. The x dimension is the dimension that varies - quickest on the image. If the value is zero or less, no minimum will - be set in the x dimension. All four parameters are zero based, - counting starts at zero. - - xmax : float, otional - Sets the maximum value of the x dimension on the bounding box - of the ouput image. If the value is zero or less, no maximum will - be set in the x dimension. - - ymin : float, optional - Sets the minimum value in the y dimension on the bounding box. The - y dimension varies less rapidly than the x and represents the line - index on the output image. If the value is zero or less, no minimum - will be set in the y dimension. - - ymax : float, optional - Sets the maximum value in the y dimension. If the value is zero or - less, no maximum will be set in the y dimension. - - unitkey : string, optional - The name of the header keyword containing the image units. The - units can either be "counts" or "cps" (counts per second.) If it is - left blank, the value is assumed to be "cps." If the value is counts, - before using the input image it is scaled by dividing it by the - exposure time. - - expkey : string, optional - The name of the header keyword containing the exposure time. The - exposure time is used to scale the image if the units are counts and - to scale the image weighting if the drizzle was initialized with - wt_scl equal to "exptime" or "expsq." If the value of this parameter - is blank, the exposure time is set to one, implying no scaling. - - wt_scl : float, optional - If drizzle was initialized with wt_scl left blank, this value will - set a scaling factor for the pixel weighting. If drizzle was - initialized with wt_scl set to "exptime" or "expsq", the exposure time - will be used to set the weight scaling and the value of this parameter - will be ignored. - """ - - insci = None - inwht = None - - if not util.is_blank(infile): - fileroot, extn = util.parse_filename(infile) - - if os.path.exists(fileroot): - handle = fits.open(fileroot) - hdu = util.get_extn(handle, extn=extn) - - if hdu is not None: - insci = hdu.data - inwcs = wcs.WCS(header=hdu.header) - insci = hdu.data.copy() - handle.close() - - if insci is None: - raise ValueError("Drizzle cannot find input file: %s" % infile) - - if not util.is_blank(inweight): - fileroot, extn = util.parse_filename(inweight) - - if os.path.exists(fileroot): - handle = fits.open(fileroot) - hdu = util.get_extn(handle, extn=extn) - - if hdu is not None: - inwht = hdu.data.copy() - handle.close() - - in_units = util.get_keyword(fileroot, unitkey, "cps") - expin = util.get_keyword(fileroot, expkey, 1.0) - - self.add_image(insci, inwcs, inwht=inwht, - xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, - expin=expin, in_units=in_units, wt_scl=wt_scl) - - def add_image(self, insci, inwcs, inwht=None, - xmin=0, xmax=0, ymin=0, ymax=0, - expin=1.0, in_units="cps", wt_scl=1.0): - """ - Combine an input image with the output drizzled image. - - Instead of reading the parameters from a fits file, you can set - them by calling this lower level method. `Add_fits_file` calls - this method after doing its setup. - - Parameters - ---------- - - insci : array - A 2d numpy array containing the input image to be drizzled. - it is an error to not supply an image. - - inwcs : wcs - The world coordinate system of the input image. This is - used to convert the pixels to the output coordinate system. - - inwht : array, optional - A 2d numpy array containing the pixel by pixel weighting. - Must have the same dimenstions as insci. If none is supplied, - the weghting is set to one. - - xmin : float, optional - This and the following three parameters set a bounding rectangle - on the output image. Only pixels on the output image inside this - rectangle will have their flux updated. Xmin sets the minimum value - of the x dimension. The x dimension is the dimension that varies - quickest on the image. If the value is zero or less, no minimum will - be set in the x dimension. All four parameters are zero based, - counting starts at zero. - - xmax : float, optional - Sets the maximum value of the x dimension on the bounding box - of the ouput image. If the value is zero or less, no maximum will - be set in the x dimension. - - ymin : float, optional - Sets the minimum value in the y dimension on the bounding box. The - y dimension varies less rapidly than the x and represents the line - index on the output image. If the value is zero or less, no minimum - will be set in the y dimension. - - ymax : float, optional - Sets the maximum value in the y dimension. If the value is zero or - less, no maximum will be set in the y dimension. - - expin : float, optional - The exposure time of the input image, a positive number. The - exposure time is used to scale the image if the units are counts and - to scale the image weighting if the drizzle was initialized with - wt_scl equal to "exptime" or "expsq." - - in_units : str, optional - The units of the input image. The units can either be "counts" - or "cps" (counts per second.) If the value is counts, before using - the input image it is scaled by dividing it by the exposure time. - - wt_scl : float, optional - If drizzle was initialized with wt_scl left blank, this value will - set a scaling factor for the pixel weighting. If drizzle was - initialized with wt_scl set to "exptime" or "expsq", the exposure time - will be used to set the weight scaling and the value of this parameter - will be ignored. - """ - - insci = insci.astype(np.float32) - util.set_pscale(inwcs) - - if inwht is None: - inwht = np.ones(insci.shape, dtype=insci.dtype) - else: - inwht = inwht.astype(np.float32) - - if self.wt_scl == "exptime": - wt_scl = expin - elif self.wt_scl == "expsq": - wt_scl = expin * expin - - self.increment_id() - self.outexptime += expin - - dodrizzle.dodrizzle(insci, inwcs, inwht, self.outwcs, - self.outsci, self.outwht, self.outcon, - expin, in_units, wt_scl, - wcslin_pscale=inwcs.pscale, uniqid=self.uniqid, - xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, - pixfrac=self.pixfrac, kernel=self.kernel, - fillval=self.fillval) - - def blot_fits_file(self, infile, interp='poly5', sinscl=1.0): - """ - Resample the output using another image's world coordinate system. - - Parameters - ---------- - - infile : str - The name of the fits file containing the world coordinate - system that the output file will be resampled to. The name may - possibly include an extension. - - interp : str, optional - The type of interpolation used in the resampling. The - possible values are "nearest" (nearest neighbor interpolation), - "linear" (bilinear interpolation), "poly3" (cubic polynomial - interpolation), "poly5" (quintic polynomial interpolation), - "sinc" (sinc interpolation), "lan3" (3rd order Lanczos - interpolation), and "lan5" (5th order Lanczos interpolation). - - sincscl : float, optional - The scaling factor for sinc interpolation. - """ - blotwcs = None - - fileroot, extn = util.parse_filename(infile) - - if os.path.exists(fileroot): - handle = fits.open(fileroot) - hdu = util.get_extn(handle, extn=extn) - - if hdu is not None: - blotwcs = wcs.WCS(header=hdu.header) - handle.close() - - if not blotwcs: - raise ValueError("Drizzle did not get a blot reference image") - - self.blot_image(blotwcs, interp=interp, sinscl=sinscl) - - def blot_image(self, blotwcs, interp='poly5', sinscl=1.0): - """ - Resample the output image using an input world coordinate system. - - Parameters - ---------- - - blotwcs : wcs - The world coordinate system to resample on. - - interp : str, optional - The type of interpolation used in the resampling. The - possible values are "nearest" (nearest neighbor interpolation), - "linear" (bilinear interpolation), "poly3" (cubic polynomial - interpolation), "poly5" (quintic polynomial interpolation), - "sinc" (sinc interpolation), "lan3" (3rd order Lanczos - interpolation), and "lan5" (5th order Lanczos interpolation). - - sincscl : float, optional - The scaling factor for sinc interpolation. - """ - - util.set_pscale(blotwcs) - self.outsci = doblot.doblot(self.outsci, self.outwcs, blotwcs, - 1.0, interp=interp, sinscl=sinscl) - - self.outwcs = blotwcs - - def increment_id(self): - """ - Increment the id count and add a plane to the context image if needed - - Drizzle tracks which input images contribute to the output image - by setting a bit in the corresponding pixel in the context image. - The uniqid indicates which bit. So it must be incremented each time - a new image is added. Each plane in the context image can hold 32 bits, - so after each 32 images, a new plane is added to the context. - """ - - # Compute what plane of the context image this input would - # correspond to: - planeid = int(self.uniqid / 32) - - # Add a new plane to the context image if planeid overflows - - if self.outcon.shape[0] == planeid: - plane = np.zeros_like(self.outcon[0]) - self.outcon = np.append(self.outcon, [plane], axis=0) - - # Increment the id - self.uniqid += 1 - - def write(self, outfile, out_units="cps", outheader=None): - """ - Write the output from a set of drizzled images to a file. - - The output file will contain three extensions. The "SCI" extension - contains the resulting image. The "WHT" extension contains the - combined weights. The "CTX" extension is a bit map. The nth bit - is set to one if the nth input image contributed non-zero flux - to the output image. The "CTX" image is three dimensionsional - to account for the possibility that there are more than 32 input - images. - - Parameters - ---------- - - outfile : str - The name of the output file. If the file already exists, - the old file is deleted after writing the new file. - - out_units : str, optional - The units of the output image, either `counts` or `cps` - (counts per second.) If the units are counts, the resulting - image will be multiplied by the computed exposure time. - - outheader : header, optional - A fits header containing cards to be added to the primary - header of the output image. - """ - - if out_units != "counts" and out_units != "cps": - raise ValueError("Illegal value for out_units: %s" % str(out_units)) - - # Write the WCS to the output image - - handle = self.outwcs.to_fits() - phdu = handle[0] - - # Write the class fields to the primary header - phdu.header['DRIZOUDA'] = \ - (self.sciext, 'Drizzle, output data image') - phdu.header['DRIZOUWE'] = \ - (self.whtext, 'Drizzle, output weighting image') - phdu.header['DRIZOUCO'] = \ - (self.ctxext, 'Drizzle, output context image') - phdu.header['DRIZWTSC'] = \ - (self.wt_scl, 'Drizzle, weighting factor for input image') - phdu.header['DRIZKERN'] = \ - (self.kernel, 'Drizzle, form of weight distribution kernel') - phdu.header['DRIZPIXF'] = \ - (self.pixfrac, 'Drizzle, linear size of drop') - phdu.header['DRIZFVAL'] = \ - (self.fillval, 'Drizzle, fill value for zero weight output pix') - phdu.header['DRIZOUUN'] = \ - (out_units, 'Drizzle, units of output image - counts or cps') - - # Update header keyword NDRIZIM to keep track of how many images have - # been combined in this product so far - phdu.header['NDRIZIM'] = (self.uniqid, 'Drizzle, number of images') - - # Update header of output image with exptime used to scale the output data - # if out_units is not counts, this will simply be a value of 1.0 - # the keyword 'exptime' will always contain the total exposure time - # of all input image regardless of the output units - - phdu.header['EXPTIME'] = \ - (self.outexptime, 'Drizzle, total exposure time') - - outexptime = 1.0 - if out_units == 'counts': - np.multiply(self.outsci, self.outexptime, self.outsci) - outexptime = self.outexptime - phdu.header['DRIZEXPT'] = (outexptime, 'Drizzle, exposure time scaling factor') - - # Copy the optional header to the primary header - - if outheader: - phdu.header.extend(outheader, unique=True) - - # Add three extensions containing, the drizzled output image, - # the total counts, and the context bitmap, in that order - - extheader = self.outwcs.to_header() - - ehdu = fits.ImageHDU() - ehdu.data = self.outsci - ehdu.header['EXTNAME'] = (self.sciext, 'Extension name') - ehdu.header['EXTVER'] = (1, 'Extension version') - ehdu.header.extend(extheader, unique=True) - handle.append(ehdu) - - whdu = fits.ImageHDU() - whdu.data = self.outwht - whdu.header['EXTNAME'] = (self.whtext, 'Extension name') - whdu.header['EXTVER'] = (1, 'Extension version') - whdu.header.extend(extheader, unique=True) - handle.append(whdu) - - xhdu = fits.ImageHDU() - xhdu.data = self.outcon - xhdu.header['EXTNAME'] = (self.ctxext, 'Extension name') - xhdu.header['EXTVER'] = (1, 'Extension version') - xhdu.header.extend(extheader, unique=True) - handle.append(xhdu) - - handle.writeto(outfile, overwrite=True) - handle.close() diff --git a/drizzle/resample.py b/drizzle/resample.py new file mode 100644 index 00000000..1fcd5af4 --- /dev/null +++ b/drizzle/resample.py @@ -0,0 +1,702 @@ +""" +The `drizzle` module defines the `Drizzle` class, for combining input +images into a single output image using the drizzle algorithm. +""" +import numpy as np + +from drizzle import cdrizzle + +__all__ = ["blot_image", "Drizzle"] + +SUPPORTED_DRIZZLE_KERNELS = [ + "square", + "gaussian", + "point", + "turbo", + "lanczos2", + "lanczos3", +] + +CTX_PLANE_BITS = 32 + + +class Drizzle: + """ + A class for managing resampling and co-adding of multiple images onto a + common output grid. The main method of this class is :py:meth:`add_image`. + The main functionality of this class is to resample and co-add multiple + images onto one output image using the "drizzle" algorithm described in + `Fruchter and Hook, PASP 2002 `_. + In the simplest terms, it redistributes flux + from input pixels to one or more output pixels based on the chosen kernel, + supplied weights, and input-to-output coordinate transformations as defined + by the ``pixmap`` argument. For more details, see :ref:`main-user-doc`. + + This class keeps track of the total exposure time of all co-added images + and also of which input images have contributed to an output (resampled) + pixel. This is accomplished via *context image*. + + Main outputs of :py:meth:`add_image` can be accessed as class properties + ``out_img``, ``out_wht``, ``out_ctx``, and ``exptime``. + + .. warning:: + Output arrays (``out_img``, ``out_wht``, and ``out_ctx``) can be + pre-allocated by the caller and be passed to the initializer or the + class initializer can allocate these arrays based on other input + parameters such as ``output_shape``. If caller-supplied output arrays + have the correct type (`numpy.float32` for ``out_img`` and ``out_wht`` + and `numpy.int32` for the ``out_ctx`` array) and if ``out_ctx`` is + large enough not to need to be resized, these arrays will be used as is + and may be modified by the :py:meth:`add_image` method. If not, + a copy of these arrays will be made when converting to the expected + type (or expanding the context array). + + Output Science Image + -------------------- + + Output science image is obtained by adding input pixel fluxes according to + equations (4) and (5) in + `Fruchter and Hook, PASP 2002 `_. + The weights and coefficients in those equations will depend on the chosen + kernel, input image weights, and pixel overlaps computed from ``pixmap``. + + Output Weight Image + ------------------- + + Output weight image stores the total weight of output science pixels + according to equation (4) in + `Fruchter and Hook, PASP 2002 `_. + It depends on the chosen kernel, input image weights, and pixel overlaps + computed from ``pixmap``. + + Output Context Image + -------------------- + + Each pixel in the context image is a bit field that encodes + information about which input image has contributed to the corresponding + pixel in the resampled data array. Context image uses 32 bit integers to + encode this information and hence it can keep track of only 32 input images. + First bit corresponds to the first input image, second bit corrsponds to the + second input image, and so on. We call this (0-indexed) order "context ID" + which is represented by the ``ctx_id`` parameter/property. If the number of + input images exceeds 32, then it is necessary to have multiple context + images ("planes") to hold information about all input images with the first + plane encoding which of the first 32 images contributed to the output data + pixel, second plane representing next 32 input images (number 33-64), etc. + For this reason, context array is either a 2D array (if the total number + of resampled images is less than 33) of the type `numpy.int32` and shape + ``(ny, nx)`` or a a 3D array of shape ``(np, ny, nx)`` where ``nx`` and + ``ny`` are dimensions of image's data. ``np`` is the number of "planes" + equal to ``(number of input images - 1) // 32 + 1``. If a bit at position + ``k`` in a pixel with coordinates ``(p, y, x)`` is 0 then input image number + ``32 * p + k`` (0-indexed) did not contribute to the output data pixel + with array coordinates ``(y, x)`` and if that bit is 1 then input image + number ``32 * p + k`` did contribute to the pixel ``(y, x)`` in the + resampled image. + + As an example, let's assume we have 8 input images. Then, when ``out_ctx`` + pixel values are displayed using binary representation (and decimal in + parenthesis), one could see values like this:: + + 00000001 (1) - only first input image contributed to this output pixel; + 00000010 (2) - 2nd input image contributed; + 00000100 (4) - 3rd input image contributed; + 10000000 (128) - 8th input image contributed; + 10000100 (132=128+4) - 3rd and 8th input images contributed; + 11001101 (205=1+4+8+64+128) - input images 1, 3, 4, 7, 8 have contributed + to this output pixel. + + In order to test if a specific input image contributed to an output pixel, + one needs to use bitwise operations. Using the example above, to test + whether input images number 4 and 5 have contributed to the output pixel + whose corresponding ``out_ctx`` value is 205 (11001101 in binary form) we + can do the following: + + >>> bool(205 & (1 << (5 - 1))) # (205 & 16) = 0 (== 0 => False): did NOT contribute + False + >>> bool(205 & (1 << (4 - 1))) # (205 & 8) = 8 (!= 0 => True): did contribute + True + + In general, to get a list of all input images that have contributed to an + output resampled pixel with image coordinates ``(x, y)``, and given a + context array ``ctx``, one can do something like this: + + .. doctest-skip:: + + >>> import numpy as np + >>> np.flatnonzero([v & (1 << k) for v in ctx[:, y, x] for k in range(32)]) + + For convenience, this functionality was implemented in the + :py:func:`~drizzle.utils.decode_context` function. + + References + ---------- + A full description of the drizzling algorithm can be found in + `Fruchter and Hook, PASP 2002 `_. + + Examples + -------- + .. highlight:: python + .. code-block:: python + + # wcs1 - WCS of the input image usually with distortions (to be resampled) + # wcs2 - WCS of the output image without distortions + + import numpy as np + from drizzle.resample import Drizzle + from drizzle.utils import calc_pixmap + + # simulate some data and a pixel map: + data = np.ones((240, 570)) + pixmap = calc_pixmap(wcs1, wcs2) + # or simulate a mapping from input image to output image frame: + # y, x = np.indices((240, 570), dtype=np.float64) + # pixmap = np.dstack([x, y]) + + # initialize Drizzle object + d = Drizzle(out_shape=(240, 570)) + d.add_image(data, exptime=15, pixmap=pixmap) + + # access outputs: + d.out_img + d.out_ctx + d.out_wht + + """ + + def __init__(self, kernel="square", fillval=None, out_shape=None, + out_img=None, out_wht=None, out_ctx=None, exptime=0.0, + begin_ctx_id=0, max_ctx_id=None, disable_ctx=False): + """ + kernel: str, optional + The name of the kernel used to combine the input. The choice of + kernel controls the distribution of flux over the kernel. The kernel + names are: "square", "gaussian", "point", "turbo", + "lanczos2", and "lanczos3". The square kernel is the default. + + .. warning:: + The "gaussian" and "lanczos2/3" kernels **DO NOT** + conserve flux. + + out_shape : tuple, None, optional + Shape (`numpy` order ``(Ny, Nx)``) of the output images (context + image will have a third dimension of size proportional to the number + of input images). This parameter is helpful when neither + ``out_img``, ``out_wht``, nor ``out_ctx`` images are provided. + + fillval: float, None, str, optional + The value of output pixels that did not have contributions from + input images' pixels. When ``fillval`` is either `None` or + ``"INDEF"`` and ``out_img`` is provided, the values of ``out_img`` + will not be modified. When ``fillval`` is either `None` or + ``"INDEF"`` and ``out_img`` is **not provided**, the values of + ``out_img`` will be initialized to `numpy.nan`. If ``fillval`` + is a string that can be converted to a number, then the output + pixels with no contributions from input images will be set to this + ``fillval`` value. + + out_img : 2D array of float32, None, optional + A 2D numpy array containing the output image produced by + drizzling. On the first call the array values should be set to zero. + Subsequent calls it will hold the intermediate results. + + out_wht : 2D array of float32, None, optional + A 2D numpy array containing the output counts. On the first + call it should be set to zero. On subsequent calls it will + hold the intermediate results. + + out_ctx : 2D or 3D array of int32, None, optional + A 2D or 3D numpy array holding a bitmap of which image was an input + for each output pixel. Should be integer zero on first call. + Subsequent calls hold intermediate results. This parameter is + ignored when ``disable_ctx`` is `True`. + + exptime : float, optional + Exposure time of previously resampled images when provided via + parameters ``out_img``, ``out_wht``, ``out_ctx``. + + begin_ctx_id : int, optional + The context ID number (0-based) of the first image that will be + resampled (using `add_image`). Subsequent images will be asigned + consecutively increasing ID numbers. This parameter is ignored + when ``disable_ctx`` is `True`. + + max_ctx_id : int, None, optional + The largest integer context ID that is *expected* to be used for + an input image. When it is a non-negative number and ``out_ctx`` is + `None`, it allows to pre-allocate the necessary array for the output + context image. If the actual number of input images that will be + resampled will exceed initial allocation for the context image, + additional context planes will be added as needed (context array + will "grow" in the third dimention as new input images are added.) + The default value of `None` is equivalent to setting ``max_ctx_id`` + equal to ``begin_ctx_id``. This parameter is ignored either when + ``out_ctx`` is provided or when ``disable_ctx`` is `True`. + + disable_ctx : bool, optional + Indicates to not create a context image. If ``disable_ctx`` is set + to `True`, parameters ``out_ctx``, ``begin_ctx_id``, and + ``max_ctx_id`` will be ignored. + + """ + self._disable_ctx = disable_ctx + + if disable_ctx: + self._ctx_id = None + self._max_ctx_id = None + else: + if begin_ctx_id < 0: + raise ValueError("Invalid context image ID") + self._ctx_id = begin_ctx_id # the ID of the *last* image to be resampled + if max_ctx_id is None: + max_ctx_id = begin_ctx_id + elif max_ctx_id < begin_ctx_id: + raise ValueError("'max_ctx_id' cannot be smaller than 'begin_ctx_id'.") + self._max_ctx_id = max_ctx_id + + if exptime < 0.0: + raise ValueError("Exposure time must be non-negative.") + + if (exptime > 0.0 and out_img is None and out_ctx is None and out_wht is None): + raise ValueError( + "Exposure time must be 0.0 for the first resampling " + "(when no ouput resampled images have been provided)." + ) + + if ( + exptime == 0.0 and + ( + (out_ctx is not None and np.sum(out_ctx) > 0) or + (out_wht is not None and np.sum(out_wht) > 0) + ) + ): + raise ValueError( + "Inconsistent exposure time and context and/or weight images: " + "Exposure time cannot be 0 when context and/or weight arrays " + "are non-zero." + ) + + self._texptime = exptime + + if kernel.lower() not in SUPPORTED_DRIZZLE_KERNELS: + raise ValueError(f"Kernel '{kernel}' is not supported.") + self._kernel = kernel + + if fillval is None: + fillval = "INDEF" + + elif isinstance(fillval, str): + fillval = fillval.strip() + if fillval.upper() in ["", "INDEF"]: + fillval = "INDEF" + else: + float(fillval) + fillval = str(fillval) + + else: + fillval = str(fillval) + + if out_img is None and fillval == "INDEF": + fillval = "NaN" + + self._fillval = fillval + + # shapes will collect user specified 'out_shape' and shapes of + # out_* arrays (if provided) in order to check all shapes are the same. + shapes = set() + + if out_img is not None: + out_img = np.asarray(out_img, dtype=np.float32) + shapes.add(out_img.shape) + + if out_wht is not None: + out_wht = np.asarray(out_wht, dtype=np.float32) + shapes.add(out_wht.shape) + + if out_ctx is not None: + out_ctx = np.asarray(out_ctx, dtype=np.int32) + if out_ctx.ndim == 2: + out_ctx = out_ctx[None, :, :] + elif out_ctx.ndim != 3: + raise ValueError("'out_ctx' must be either a 2D or 3D array.") + shapes.add(out_ctx.shape[1:]) + + if out_shape is not None: + shapes.add(tuple(out_shape)) + + if len(shapes) == 1: + self._out_shape = shapes.pop() + self._alloc_output_arrays( + out_shape=self._out_shape, + max_ctx_id=max_ctx_id, + out_img=out_img, + out_wht=out_wht, + out_ctx=out_ctx, + ) + elif len(shapes) > 1: + raise ValueError( + "Inconsistent data shapes specified: 'out_shape' and/or " + "out_img, out_wht, out_ctx have different shapes." + ) + else: + self._out_shape = None + self._out_img = None + self._out_wht = None + self._out_ctx = None + + @property + def fillval(self): + """Fill value for output pixels without contributions from input images.""" + return self._fillval + + @property + def kernel(self): + """Resampling kernel.""" + return self._kernel + + @property + def ctx_id(self): + """Context image "ID" (0-based ) of the next image to be resampled.""" + return self._ctx_id + + @property + def out_img(self): + """Output resampled image.""" + return self._out_img + + @property + def out_wht(self): + """Output weight image.""" + return self._out_wht + + @property + def out_ctx(self): + """Output "context" image.""" + return self._out_ctx + + @property + def total_exptime(self): + """Total exposure time of all resampled images.""" + return self._texptime + + def _alloc_output_arrays(self, out_shape, max_ctx_id, out_img, out_wht, + out_ctx): + # allocate arrays as needed: + if out_wht is None: + self._out_wht = np.zeros(out_shape, dtype=np.float32) + else: + self._out_wht = out_wht + + if self._disable_ctx: + self._out_ctx = None + else: + if out_ctx is None: + n_ctx_planes = max_ctx_id // CTX_PLANE_BITS + 1 + ctx_shape = (n_ctx_planes, ) + out_shape + self._out_ctx = np.zeros(ctx_shape, dtype=np.int32) + else: + self._out_ctx = out_ctx + + if not (out_wht is None and out_ctx is None): + # check that input data make sense: weight of pixels with + # non-zero context values must be different from zero: + if np.any( + np.bitwise_xor( + self._out_wht > 0.0, + np.sum(self._out_ctx, axis=0) > 0 + ) + ): + raise ValueError( + "Inconsistent values of supplied 'out_wht' and " + "'out_ctx' arrays. Pixels with non-zero context " + "values must have positive weights and vice-versa." + ) + + if out_img is None: + self._out_img = np.zeros(out_shape, dtype=np.float32) + else: + self._out_img = out_img + + def _increment_ctx_id(self): + """ + Returns a pair of the *current* plane number and bit number in that + plane and increments context image ID + (after computing the return value). + """ + if self._disable_ctx: + return None, 0 + + self._plane_no = self._ctx_id // CTX_PLANE_BITS + depth = self._out_ctx.shape[0] + + if self._plane_no >= depth: + # Add a new plane to the context image if planeid overflows + plane = np.zeros((1, ) + self._out_shape, np.int32) + self._out_ctx = np.append(self._out_ctx, plane, axis=0) + + plane_info = (self._plane_no, self._ctx_id % CTX_PLANE_BITS) + # increment ID for the *next* image to be added: + self._ctx_id += 1 + + return plane_info + + def add_image(self, data, exptime, pixmap, scale=1.0, + weight_map=None, wht_scale=1.0, pixfrac=1.0, in_units='cps', + xmin=None, xmax=None, ymin=None, ymax=None): + """ + Resample and add an image to the cumulative output image. Also, update + output total weight image and context images. + + Parameters + ---------- + data : 2D numpy.ndarray + A 2D numpy array containing the input image to be drizzled. + + exptime : float + The exposure time of the input image, a positive number. The + exposure time is used to scale the image if the units are counts. + + pixmap : 3D array + A mapping from input image (``data``) coordinates to resampled + (``out_img``) coordinates. ``pixmap`` must be an array of shape + ``(Ny, Nx, 2)`` where ``(Ny, Nx)`` is the shape of the input image. + ``pixmap[..., 0]`` forms a 2D array of X-coordinates of input + pixels in the ouput frame and ``pixmap[..., 1]`` forms a 2D array of + Y-coordinates of input pixels in the ouput coordinate frame. + + scale : float, optional + The pixel scale of the input image. Conceptually, this is the + linear dimension of a side of a pixel in the input image, but it + is not limited to this and can be set to change how the drizzling + algorithm operates. + + weight_map : 2D array, None, optional + A 2D numpy array containing the pixel by pixel weighting. + Must have the same dimensions as ``data``. + + When ``weight_map`` is `None`, the weight of input data pixels will + be assumed to be 1. + + wht_scale : float + A scaling factor applied to the pixel by pixel weighting. + + pixfrac : float, optional + The fraction of a pixel that the pixel flux is confined to. The + default value of 1 has the pixel flux evenly spread across the image. + A value of 0.5 confines it to half a pixel in the linear dimension, + so the flux is confined to a quarter of the pixel area when the square + kernel is used. + + in_units : str + The units of the input image. The units can either be "counts" + or "cps" (counts per second.) + + xmin : float, optional + This and the following three parameters set a bounding rectangle + on the input image. Only pixels on the input image inside this + rectangle will have their flux added to the output image. Xmin + sets the minimum value of the x dimension. The x dimension is the + dimension that varies quickest on the image. If the value is zero, + no minimum will be set in the x dimension. All four parameters are + zero based, counting starts at zero. + + xmax : float, optional + Sets the maximum value of the x dimension on the bounding box + of the input image. If the value is zero, no maximum will + be set in the x dimension, the full x dimension of the output + image is the bounding box. + + ymin : float, optional + Sets the minimum value in the y dimension on the bounding box. The + y dimension varies less rapidly than the x and represents the line + index on the input image. If the value is zero, no minimum will be + set in the y dimension. + + ymax : float, optional + Sets the maximum value in the y dimension. If the value is zero, no + maximum will be set in the y dimension, the full x dimension + of the output image is the bounding box. + + Returns + ------- + nskip : float + The number of lines from the box defined by + ``((xmin, xmax), (ymin, ymax))`` in the input image that were + ignored and did not contribute to the output image. + + nmiss : float + The number of pixels from the box defined by + ``((xmin, xmax), (ymin, ymax))`` in the input image that were + ignored and did not contribute to the output image. + + """ + # this enables initializer to not need output image shape at all and + # set output image shape based on output coordinates from the pixmap. + # + if self._out_shape is None: + pmap_xmin = int(np.floor(np.nanmin(pixmap[:, :, 0]))) + pmap_xmax = int(np.ceil(np.nanmax(pixmap[:, :, 0]))) + pmap_ymin = int(np.floor(np.nanmin(pixmap[:, :, 1]))) + pmap_ymax = int(np.ceil(np.nanmax(pixmap[:, :, 1]))) + pixmap = pixmap.copy() + pixmap[:, :, 0] -= pmap_xmin + pixmap[:, :, 1] -= pmap_ymin + self._out_shape = ( + pmap_xmax - pmap_xmin + 1, + pmap_ymax - pmap_ymin + 1 + ) + + self._alloc_output_arrays( + out_shape=self._out_shape, + max_ctx_id=max(self._max_ctx_id, self._ctx_id), + out_img=None, + out_wht=None, + out_ctx=None, + ) + + plane_no, id_in_plane = self._increment_ctx_id() + + if exptime <= 0.0: + raise ValueError("'exptime' *must* be a strictly positive number.") + + # Ensure that the fillval parameter gets properly interpreted + # for use with tdriz + if in_units == 'cps': + expscale = 1.0 + else: + expscale = exptime + + self._texptime += exptime + + data = np.asarray(data, dtype=np.float32) + pixmap = np.asarray(pixmap, dtype=np.float64) + in_ymax, in_xmax = data.shape + + if pixmap.shape[:2] != data.shape: + raise ValueError( + "'pixmap' shape is not consistent with 'data' shape." + ) + + if xmin is None or xmin < 0: + xmin = 0 + + if ymin is None or ymin < 0: + ymin = 0 + + if xmax is None or xmax > in_xmax - 1: + xmax = in_xmax - 1 + + if ymax is None or ymax > in_ymax - 1: + ymax = in_ymax - 1 + + if weight_map is not None: + weight_map = np.asarray(weight_map, dtype=np.float32) + else: # TODO: this should not be needed after C code modifications + weight_map = np.ones_like(data) + + pixmap = np.asarray(pixmap, dtype=np.float64) + + if self._disable_ctx: + ctx_plane = None + else: + if self._out_ctx.ndim == 2: + raise AssertionError("Context image is expected to be 3D") + ctx_plane = self._out_ctx[plane_no] + + # TODO: probably tdriz should be modified to not return version. + # we should not have git, Python, C, ... versions + + # TODO: While drizzle code in cdrizzlebox.c supports weight_map=None, + # cdrizzleapi.c does not. It should be modified to support this + # for performance reasons. + + _vers, nmiss, nskip = cdrizzle.tdriz( + input=data, + weights=weight_map, + pixmap=pixmap, + output=self._out_img, + counts=self._out_wht, + context=ctx_plane, + uniqid=id_in_plane + 1, + xmin=xmin, + xmax=xmax, + ymin=ymin, + ymax=ymax, + scale=scale, # scales image intensity. usually equal to pixel scale + pixfrac=pixfrac, + kernel=self._kernel, + in_units=in_units, + expscale=expscale, + wtscale=wht_scale, + fillstr=self._fillval, + ) + self._cversion = _vers # TODO: probably not needed + + return nmiss, nskip + + +def blot_image(data, pixmap, pix_ratio, exptime, output_pixel_shape, + interp='poly5', sinscl=1.0): + """ + Resample the ``data`` input image onto an output grid defined by + the ``pixmap`` array. ``blot_image`` performs resampling using one of + the several interpolation algorithms and, unlike the "drizzle" algorithm + with 'square', 'turbo', and 'point' kernels, this resampling is not + flux-conserving. + + This method works best for with well sampled images and thus it is + typically used to resample the output of :py:class:`Drizzle` back to the + coordinate grids of input images of :py:meth:`Drizzle.add_image`. + The output of :py:class:`Drizzle` are usually well sampled images especially + if it was created from a set of dithered images. + + Parameters + ---------- + data : 2D array + Input numpy array of the source image in units of 'cps'. + + pixmap : 3D array + A mapping from input image (``data``) coordinates to resampled + (``out_img``) coordinates. ``pixmap`` must be an array of shape + ``(Ny, Nx, 2)`` where ``(Ny, Nx)`` is the shape of the input image. + ``pixmap[..., 0]`` forms a 2D array of X-coordinates of input + pixels in the ouput frame and ``pixmap[..., 1]`` forms a 2D array of + Y-coordinates of input pixels in the ouput coordinate frame. + + output_pixel_shape : tuple of int + A tuple of two integer numbers indicating the dimensions of the output + image ``(Nx, Ny)``. + + pix_ratio : float + Ratio of the input image pixel scale to the ouput image pixel scale. + + exptime : float + The exposure time of the input image. + + interp : str, optional + The type of interpolation used in the resampling. The + possible values are: + + - "nearest" (nearest neighbor interpolation); + - "linear" (bilinear interpolation); + - "poly3" (cubic polynomial interpolation); + - "poly5" (quintic polynomial interpolation); + - "sinc" (sinc interpolation); + - "lan3" (3rd order Lanczos interpolation); and + - "lan5" (5th order Lanczos interpolation). + + sincscl : float, optional + The scaling factor for "sinc" interpolation. + + Returns + ------- + out_img : 2D numpy.ndarray + A 2D numpy array containing the resampled image data. + + """ + out_img = np.zeros(output_pixel_shape[::-1], dtype=np.float32) + + cdrizzle.tblot(data, pixmap, out_img, scale=pix_ratio, kscale=1.0, + interp=interp, exptime=exptime, misval=0.0, sinscl=sinscl) + + return out_img diff --git a/drizzle/tests/test_cdrizzle.py b/drizzle/tests/test_cdrizzle.py index c24f30d9..a9c7e1df 100644 --- a/drizzle/tests/test_cdrizzle.py +++ b/drizzle/tests/test_cdrizzle.py @@ -9,16 +9,21 @@ def test_cdrizzle(): """ size = 100 - data = np.zeros((size,size), dtype='float32') - weights = np.ones((size,size), dtype='float32') + data = np.zeros((size, size), dtype='float32') + weights = np.ones((size, size), dtype='float32') - pixmap = np.indices((size,size), dtype='float64') + pixmap = np.indices((size, size), dtype='float64') pixmap = pixmap.transpose() - output_data = np.zeros((size,size), dtype='float32') - output_counts = np.zeros((size,size), dtype='float32') - output_context = np.zeros((size,size), dtype='int32') + output_data = np.zeros((size, size), dtype='float32') + output_counts = np.zeros((size, size), dtype='float32') + output_context = np.zeros((size, size), dtype='int32') - cdrizzle.test_cdrizzle(data, weights, pixmap, - output_data, output_counts, - output_context) + cdrizzle.test_cdrizzle( + data, + weights, + pixmap, + output_data, + output_counts, + output_context, + ) diff --git a/drizzle/tests/test_drizzle.py b/drizzle/tests/test_drizzle.py deleted file mode 100644 index 54c07977..00000000 --- a/drizzle/tests/test_drizzle.py +++ /dev/null @@ -1,834 +0,0 @@ -import math -import os -import pytest - -import numpy as np -from astropy import wcs -from astropy.io import fits - -from drizzle import drizzle, cdrizzle - -TEST_DIR = os.path.abspath(os.path.dirname(__file__)) -DATA_DIR = os.path.join(TEST_DIR, 'data') -ok = False - - -def bound_image(image): - """ - Compute region where image is non-zero - """ - coords = np.nonzero(image) - ymin = coords[0].min() - ymax = coords[0].max() - xmin = coords[1].min() - xmax = coords[1].max() - return (ymin, ymax, xmin, xmax) - - -def centroid(image, size, center): - """ - Compute the centroid of a rectangular area - """ - ylo = int(center[0] - size / 2) - yhi = min(ylo + size, image.shape[0]) - xlo = int(center[1] - size / 2) - xhi = min(xlo + size, image.shape[1]) - - yx1 = np.mgrid[ylo:yhi, xlo:xhi, 1:2] - center = (yx1[..., 0] * image[ylo:yhi, xlo:xhi]).sum( - axis=(1, 2), - dtype=np.float64, - ) - - if center[2] == 0.0: - return None - - center[0] /= center[2] - center[1] /= center[2] - return center - - -def centroid_close(list_of_centroids, size, point): - """ - Find if any centroid is close to a point - """ - for i in range(len(list_of_centroids) - 1, -1, -1): - if (abs(list_of_centroids[i][0] - point[0]) < int(size / 2) and - abs(list_of_centroids[i][1] - point[1]) < int(size / 2)): - return 1 - - return 0 - - -def centroid_compare(centroid): - return centroid[1] - - -def centroid_distances(image1, image2, amp, size): - """ - Compute a list of centroids and the distances between them in two images - """ - distances = [] - list_of_centroids = centroid_list(image2, amp, size) - for center2 in list_of_centroids: - center1 = centroid(image1, size, center2) - if center1 is None: - continue - - disty = center2[0] - center1[0] - distx = center2[1] - center1[1] - dist = math.sqrt(disty * disty + distx * distx) - dflux = abs(center2[2] - center1[2]) - distances.append([dist, dflux, center1, center2]) - - distances.sort(key=centroid_compare) - return distances - - -def centroid_list(image, amp, size): - """ - Find the next centroid - """ - list_of_centroids = [] - points = np.transpose(np.nonzero(image > amp)) - for point in points: - if not centroid_close(list_of_centroids, size, point): - center = centroid(image, size, point) - list_of_centroids.append(center) - - return list_of_centroids - - -def centroid_statistics(title, fname, image1, image2, amp, size): - """ - write centroid statistics to compare differences btw two images - """ - stats = ("minimum", "median", "maximum") - images = (None, None, image1, image2) - im_type = ("", "", "test", "reference") - - diff = [] - distances = centroid_distances(image1, image2, amp, size) - indexes = (0, int(len(distances) / 2), len(distances) - 1) - fd = open(fname, 'w') - fd.write("*** %s ***\n" % title) - - if len(distances) == 0: - diff = [0.0, 0.0, 0.0] - fd.write("No matches!!\n") - - elif len(distances) == 1: - diff = [distances[0][0], distances[0][0], distances[0][0]] - - fd.write("1 match\n") - fd.write("distance = %f flux difference = %f\n" % - (distances[0][0], distances[0][1])) - - for j in range(2, 4): - ylo = int(distances[0][j][0]) - 1 - yhi = int(distances[0][j][0]) + 2 - xlo = int(distances[0][j][1]) - 1 - xhi = int(distances[0][j][1]) + 2 - subimage = images[j][ylo:yhi,xlo:xhi] - fd.write("\n%s image centroid = (%f,%f) image flux = %f\n" % - (im_type[j], distances[0][j][0], distances[0][j][1], - distances[0][j][2])) - fd.write(str(subimage) + "\n") - - else: - fd.write("%d matches\n" % len(distances)) - - for k in range(0,3): - i = indexes[k] - diff.append(distances[i][0]) - fd.write("\n%s distance = %f flux difference = %f\n" % - (stats[k],distances[i][0], distances[i][1])) - - for j in range(2, 4): - ylo = int(distances[i][j][0]) - 1 - yhi = int(distances[i][j][0]) + 2 - xlo = int(distances[i][j][1]) - 1 - xhi = int(distances[i][j][1]) + 2 - subimage = images[j][ylo:yhi,xlo:xhi] - fd.write("\n%s %s image centroid = (%f,%f) image flux = %f\n" % - (stats[k], im_type[j], distances[i][j][0], - distances[i][j][1], distances[i][j][2])) - fd.write(str(subimage) + "\n") - - fd.close() - return tuple(diff) - - -def make_point_image(input_image, point, value): - """ - Create an image with a single point set - """ - output_image = np.zeros(input_image.shape, dtype=input_image.dtype) - output_image[point] = value - return output_image - - -def make_grid_image(input_image, spacing, value): - """ - Create an image with points on a grid set - """ - output_image = np.zeros(input_image.shape, dtype=input_image.dtype) - - shape = output_image.shape - half_space = int(spacing / 2) - for y in range(half_space, shape[0], spacing): - for x in range(half_space, shape[1], spacing): - output_image[y,x] = value - - return output_image - - -def print_wcs(title, wcs): - """ - Print the wcs header cards - """ - print("=== %s ===" % title) - print(wcs.to_header_string()) - - -def read_image(filename): - """ - Read the image from a fits file - """ - path = os.path.join(DATA_DIR, filename) - hdu = fits.open(path) - - image = hdu[1].data - hdu.close() - return image - - -def read_wcs(filename): - """ - Read the wcs of a fits file - """ - path = os.path.join(DATA_DIR, filename) - hdu = fits.open(path) - the_wcs = wcs.WCS(hdu[1].header) - hdu.close() - return the_wcs - - -def test_square_with_point(tmpdir): - """ - Test do_driz square kernel with point - """ - output = str(tmpdir.join('output_square_point.fits')) - output_difference = str(tmpdir.join('difference_square_point.txt')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') - output_template = os.path.join(DATA_DIR, 'reference_square_point.fits') - - insci = read_image(input_file) - inwcs = read_wcs(input_file) - insci = make_point_image(insci, (500, 200), 100.0) - inwht = np.ones(insci.shape,dtype=insci.dtype) - output_wcs = read_wcs(output_template) - - driz = drizzle.Drizzle(outwcs=output_wcs, wt_scl="") - driz.add_image(insci, inwcs, inwht=inwht) - - if ok: - driz.write(output_template) - else: - driz.write(output) - template_data = read_image(output_template) - - min_diff, med_diff, max_diff = centroid_statistics("square with point", - output_difference, - driz.outsci, - template_data, 20.0, 8) - - assert med_diff < 1.0e-6 - assert max_diff < 1.0e-5 - - -@pytest.mark.parametrize( - 'kernel', ['square', 'point', 'turbo', 'gaussian', 'lanczos3'], -) -def test_zero_input_weight(kernel): - """ - Test do_driz square kernel with grid - """ - # initialize input: - insci = np.ones((200, 400), dtype=np.float32) - inwht = np.ones((200, 400), dtype=np.float32) - inwht[:, 150:155] = 0 - - # initialize output: - outsci = np.zeros((210, 410), dtype=np.float32) - outwht = np.zeros((210, 410), dtype=np.float32) - outctx = np.zeros((210, 410), dtype=np.int32) - - # define coordinate mapping: - pixmap = np.moveaxis(np.mgrid[1:201, 1:401][::-1], 0, -1) - - # resample: - cdrizzle.tdriz( - insci, inwht, pixmap, - outsci, outwht, outctx, - uniqid=1, - xmin=0, xmax=400, - ymin=0, ymax=200, - pixfrac=1, - kernel=kernel, - in_units='cps', - expscale=1, - wtscale=1, - fillstr='INDEF', - ) - - # check that no pixel with 0 weight has any counts: - assert np.sum(np.abs(outsci[(outwht == 0)])) == 0.0 - - -def test_square_with_grid(tmpdir): - """ - Test do_driz square kernel with grid - """ - output = str(tmpdir.join('output_square_grid.fits')) - output_difference = str(tmpdir.join('difference_square_grid.txt')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') - output_template = os.path.join(DATA_DIR, 'reference_square_grid.fits') - - insci = read_image(input_file) - inwcs = read_wcs(input_file) - insci = make_grid_image(insci, 64, 100.0) - - inwht = np.ones(insci.shape,dtype=insci.dtype) - output_wcs = read_wcs(output_template) - - driz = drizzle.Drizzle(outwcs=output_wcs, wt_scl="") - driz.add_image(insci, inwcs, inwht=inwht) - - if ok: - driz.write(output_template) - else: - driz.write(output) - template_data = read_image(output_template) - - min_diff, med_diff, max_diff = centroid_statistics("square with grid", - output_difference, - driz.outsci, - template_data, 20.0, 8) - assert med_diff < 1.0e-6 - assert max_diff < 1.0e-5 - - -def test_turbo_with_grid(tmpdir): - """ - Test do_driz turbo kernel with grid - """ - output = str(tmpdir.join('output_turbo_grid.fits')) - output_difference = str(tmpdir.join('difference_turbo_grid.txt')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') - output_template = os.path.join(DATA_DIR, 'reference_turbo_grid.fits') - - insci = read_image(input_file) - inwcs = read_wcs(input_file) - insci = make_grid_image(insci, 64, 100.0) - inwht = np.ones(insci.shape,dtype=insci.dtype) - output_wcs = read_wcs(output_template) - - driz = drizzle.Drizzle(outwcs=output_wcs, wt_scl="", kernel='turbo') - driz.add_image(insci, inwcs, inwht=inwht) - - if ok: - driz.write(output_template) - else: - driz.write(output) - template_data = read_image(output_template) - - min_diff, med_diff, max_diff = centroid_statistics("turbo with grid", - output_difference, - driz.outsci, - template_data, 20.0, 8) - - assert med_diff < 1.0e-6 - assert max_diff < 1.0e-5 - - -def test_gaussian_with_grid(tmpdir): - """ - Test do_driz gaussian kernel with grid - """ - output = str(tmpdir.join('output_gaussian_grid.fits')) - output_difference = str(tmpdir.join('difference_gaussian_grid.txt')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') - output_template = os.path.join(DATA_DIR, 'reference_gaussian_grid.fits') - - insci = read_image(input_file) - inwcs = read_wcs(input_file) - insci = make_grid_image(insci, 64, 100.0) - inwht = np.ones(insci.shape,dtype=insci.dtype) - output_wcs = read_wcs(output_template) - - driz = drizzle.Drizzle(outwcs=output_wcs, wt_scl="", kernel='gaussian') - driz.add_image(insci, inwcs, inwht=inwht) - - if ok: - driz.write(output_template) - else: - driz.write(output) - template_data = read_image(output_template) - - min_diff, med_diff, max_diff = centroid_statistics("gaussian with grid", - output_difference, - driz.outsci, - template_data, 20.0, 8) - - assert med_diff < 1.0e-6 - assert max_diff < 2.0e-5 - - -def test_lanczos_with_grid(tmpdir): - """ - Test do_driz lanczos kernel with grid - """ - output = str(tmpdir.join('output_lanczos_grid.fits')) - output_difference = str(tmpdir.join('difference_lanczos_grid.txt')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') - output_template = os.path.join(DATA_DIR, 'reference_lanczos_grid.fits') - - insci = read_image(input_file) - inwcs = read_wcs(input_file) - insci = make_grid_image(insci, 64, 100.0) - inwht = np.ones(insci.shape,dtype=insci.dtype) - output_wcs = read_wcs(output_template) - - driz = drizzle.Drizzle(outwcs=output_wcs, wt_scl="", kernel='lanczos3') - driz.add_image(insci, inwcs, inwht=inwht) - - if ok: - driz.write(output_template) - else: - driz.write(output) - template_data = read_image(output_template) - - min_diff, med_diff, max_diff = centroid_statistics("lanczos with grid", - output_difference, - driz.outsci, - template_data, 20.0, 8) - assert med_diff < 1.0e-6 - assert max_diff < 1.0e-5 - - -def test_tophat_with_grid(tmpdir): - """ - Test do_driz tophat kernel with grid - """ - output = str(tmpdir.join('output_tophat_grid.fits')) - output_difference = str(tmpdir.join('difference_tophat_grid.txt')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') - output_template = os.path.join(DATA_DIR, 'reference_tophat_grid.fits') - - insci = read_image(input_file) - inwcs = read_wcs(input_file) - insci = make_grid_image(insci, 64, 100.0) - inwht = np.ones(insci.shape,dtype=insci.dtype) - output_wcs = read_wcs(output_template) - - driz = drizzle.Drizzle(outwcs=output_wcs, wt_scl="", kernel='tophat') - driz.add_image(insci, inwcs, inwht=inwht) - - if ok: - driz.write(output_template) - else: - driz.write(output) - template_data = read_image(output_template) - - min_diff, med_diff, max_diff = centroid_statistics("tophat with grid", - output_difference, - driz.outsci, - template_data, 20.0, 8) - assert med_diff < 1.0e-6 - assert max_diff < 1.0e-5 - - -def test_point_with_grid(tmpdir): - """ - Test do_driz point kernel with grid - """ - output = str(tmpdir.join('output_point_grid.fits')) - output_difference = str(tmpdir.join('difference_point_grid.txt')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') - output_template = os.path.join(DATA_DIR, 'reference_point_grid.fits') - - insci = read_image(input_file) - inwcs = read_wcs(input_file) - insci = make_grid_image(insci, 64, 100.0) - inwht = np.ones(insci.shape,dtype=insci.dtype) - output_wcs = read_wcs(output_template) - - driz = drizzle.Drizzle(outwcs=output_wcs, wt_scl="", kernel='point') - driz.add_image(insci, inwcs, inwht=inwht) - - if ok: - driz.write(output_template) - else: - driz.write(output) - template_data = read_image(output_template) - - min_diff, med_diff, max_diff = centroid_statistics("point with grid", - output_difference, - driz.outsci, - template_data, 20.0, 8) - assert med_diff < 1.0e-6 - assert max_diff < 1.0e-5 - - -def test_blot_with_point(tmpdir): - """ - Test do_blot with point image - """ - output = str(tmpdir.join('output_blot_point.fits')) - output_difference = str(tmpdir.join('difference_blot_point.txt')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') - output_template = os.path.join(DATA_DIR, 'reference_blot_point.fits') - - outsci = read_image(input_file) - outwcs = read_wcs(input_file) - outsci = make_point_image(outsci, (500, 200), 40.0) - inwcs = read_wcs(output_template) - - driz = drizzle.Drizzle(outwcs=outwcs) - driz.outsci = outsci - - driz.blot_image(inwcs) - - if ok: - driz.write(output_template) - else: - driz.write(output) - template_data = read_image(output_template) - - min_diff, med_diff, max_diff = centroid_statistics("blot with point", - output_difference, - driz.outsci, - template_data, 20.0, 16) - assert med_diff < 1.0e-6 - assert max_diff < 1.0e-5 - - -def test_blot_with_default(tmpdir): - """ - Test do_blot with default grid image - """ - output = str(tmpdir.join('output_blot_default.fits')) - output_difference = str(tmpdir.join('difference_blot_default.txt')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') - output_template = os.path.join(DATA_DIR, 'reference_blot_default.fits') - - outsci = read_image(input_file) - outsci = make_grid_image(outsci, 64, 100.0) - outwcs = read_wcs(input_file) - inwcs = read_wcs(output_template) - - driz = drizzle.Drizzle(outwcs=outwcs) - driz.outsci = outsci - - driz.blot_image(inwcs) - - if ok: - driz.write(output_template) - else: - driz.write(output) - template_data = read_image(output_template) - - min_diff, med_diff, max_diff = centroid_statistics("blot with defaults", - output_difference, - driz.outsci, - template_data, 20.0, 16) - - assert med_diff < 1.0e-6 - assert max_diff < 1.0e-5 - - -def test_blot_with_lan3(tmpdir): - """ - Test do_blot with lan3 grid image - """ - output = str(tmpdir.join('output_blot_lan3.fits')) - output_difference = str(tmpdir.join('difference_blot_lan3.txt')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') - output_template = os.path.join(DATA_DIR, 'reference_blot_lan3.fits') - - outsci = read_image(input_file) - outsci = make_grid_image(outsci, 64, 100.0) - outwcs = read_wcs(input_file) - inwcs = read_wcs(output_template) - - driz = drizzle.Drizzle(outwcs=outwcs) - driz.outsci = outsci - - driz.blot_image(inwcs, interp="lan3") - - if ok: - driz.write(output_template) - else: - driz.write(output) - template_data = read_image(output_template) - - min_diff, med_diff, max_diff = centroid_statistics("blot with lan3", - output_difference, - driz.outsci, - template_data, 20.0, 16) - assert med_diff < 1.0e-6 - assert max_diff < 1.0e-5 - - -def test_blot_with_lan5(tmpdir): - """ - Test do_blot with lan5 grid image - """ - output = str(tmpdir.join('output_blot_lan5.fits')) - output_difference = str(tmpdir.join('difference_blot_lan5.txt')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') - output_template = os.path.join(DATA_DIR, 'reference_blot_lan5.fits') - - outsci = read_image(input_file) - outsci = make_grid_image(outsci, 64, 100.0) - outwcs = read_wcs(input_file) - inwcs = read_wcs(output_template) - - driz = drizzle.Drizzle(outwcs=outwcs) - driz.outsci = outsci - - driz.blot_image(inwcs, interp="lan5") - - if ok: - driz.write(output_template) - else: - driz.write(output) - template_data = read_image(output_template) - - min_diff, med_diff, max_diff = centroid_statistics("blot with lan5", - output_difference, - driz.outsci, - template_data, 20.0, 16) - assert med_diff < 1.0e-6 - assert max_diff < 1.0e-5 - - -def test_context_planes(): - """Reproduce error seen in issue #50""" - shape = [10, 10] - outwcs = wcs.WCS() - outwcs.pixel_shape = shape - driz = drizzle.Drizzle(outwcs=outwcs) - - image = np.ones(shape) - inwcs = wcs.WCS() - inwcs.pixel_shape = shape - - for i in range(32): - driz.add_image(image, inwcs) - assert driz.outcon.shape == (1, 10, 10) - - driz.add_image(image, inwcs) - assert driz.outcon.shape == (2, 10, 10) - - -@pytest.mark.parametrize( - 'kernel', - [ - 'square', - 'point', - 'turbo', - pytest.param( - 'lanczos2', - marks=pytest.mark.xfail(reason='Not a flux-conserving kernel'), - ), - pytest.param( - 'lanczos3', - marks=pytest.mark.xfail(reason='Not a flux-conserving kernel'), - ), - pytest.param( - 'gaussian', - marks=pytest.mark.xfail(reason='Not a flux-conserving kernel'), - ), - ], -) -def test_flux_conservation_nondistorted(kernel): - n = 200 - in_shape = (n, n) - - # input coordinate grid: - y, x = np.indices(in_shape, dtype=np.float64) - - # simulate a gaussian "star": - fwhm = 2.9 - x0 = 50.0 - y0 = 68.0 - sig = fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0 * fwhm))) - sig2 = sig * sig - star = np.exp(-0.5 / sig2 * ((x.astype(np.float32) - x0)**2 + (y.astype(np.float32) - y0)**2)) - in_sci = (star / np.sum(star)).astype(np.float32) # normalize to 1 - in_wht = np.ones(in_shape, dtype=np.float32) - - # linear shift: - xp = x + 0.5 - yp = y + 0.2 - - pixmap = np.dstack([xp, yp]) - - out_shape = (int(yp.max()) + 1, int(xp.max()) + 1) - # make sure distorion is not moving flux out of the image towards negative - # coordinates (just because of the simple way of how we account for output - # image size) - assert np.min(xp) > -0.5 and np.min(yp) > -0.5 - - out_sci = np.zeros(out_shape, dtype=np.float32) - out_ctx = np.zeros(out_shape, dtype=np.int32) - out_wht = np.zeros(out_shape, dtype=np.float32) - - cdrizzle.tdriz( - in_sci, - in_wht, - pixmap, - out_sci, - out_wht, - out_ctx, - pixfrac=1.0, - scale=1.0, - kernel=kernel, - in_units="cps", - expscale=1.0, - wtscale=1.0, - ) - - assert np.allclose( - np.sum(out_sci * out_wht), - np.sum(in_sci), - atol=0.0, - rtol=0.0001, - ) - -@pytest.mark.parametrize( - 'kernel', - [ - 'square', - 'point', - 'turbo', - pytest.param( - 'lanczos2', - marks=pytest.mark.xfail(reason='Not a flux-conserving kernel'), - ), - pytest.param( - 'lanczos3', - marks=pytest.mark.xfail(reason='Not a flux-conserving kernel'), - ), - pytest.param( - 'gaussian', - marks=pytest.mark.xfail(reason='Not a flux-conserving kernel'), - ), - ], -) -def test_flux_conservation_distorted(kernel): - n = 200 - in_shape = (n, n) - - # input coordinate grid: - y, x = np.indices(in_shape, dtype=np.float64) - - # simulate a gaussian "star": - fwhm = 2.9 - x0 = 50.0 - y0 = 68.0 - sig = fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0 * fwhm))) - sig2 = sig * sig - star = np.exp(-0.5 / sig2 * ((x.astype(np.float32) - x0)**2 + (y.astype(np.float32) - y0)**2)) - in_sci = (star / np.sum(star)).astype(np.float32) # normalize to 1 - in_wht = np.ones(in_shape, dtype=np.float32) - - # linear shift: - xp = x + 0.5 - yp = y + 0.2 - # add distortion: - xp += 1e-4 * x**2 + 1e-5 * x * y - yp += 1e-3 * y**2 - 2e-5 * x * y - - pixmap = np.dstack([xp, yp]) - - out_shape = (int(yp.max()) + 1, int(xp.max()) + 1) - # make sure distorion is not moving (pixels with) flux out of the image - # towards negative coordinates (just because of the simple way of how we - # account for output image size): - assert np.min(xp) > -0.5 and np.min(yp) > -0.5 - - out_sci = np.zeros(out_shape, dtype=np.float32) - out_ctx = np.zeros(out_shape, dtype=np.int32) - out_wht = np.zeros(out_shape, dtype=np.float32) - - cdrizzle.tdriz( - in_sci, - in_wht, - pixmap, - out_sci, - out_wht, - out_ctx, - pixfrac=1.0, - scale=1.0, - kernel=kernel, - in_units="cps", - expscale=1.0, - wtscale=1.0, - ) - - assert np.allclose( - np.sum(out_sci * out_wht), - np.sum(in_sci), - atol=0.0, - rtol=0.0001, - ) - - -def test_no_context_array(): - """ - Test that providing None for "outctx" is supported. - """ - # initialize input: - insci = np.ones((200, 400), dtype=np.float32) - inwht = np.ones((200, 400), dtype=np.float32) - inwht[:, 150:155] = 0 - - # initialize output: - outsci = np.zeros((210, 410), dtype=np.float32) - outwht = np.zeros((210, 410), dtype=np.float32) - outctx = None - - # define coordinate mapping: - pixmap = np.moveaxis(np.mgrid[1:201, 1:401][::-1], 0, -1) - - # resample: - cdrizzle.tdriz( - insci, inwht, pixmap, - outsci, outwht, outctx, - uniqid=1, - xmin=0, xmax=400, - ymin=0, ymax=200, - pixfrac=1, - kernel='square', - in_units='cps', - expscale=1, - wtscale=1, - fillstr='INDEF', - ) - - # check that no pixel with 0 weight has any counts: - assert np.sum(np.abs(outsci[(outwht == 0)])) == 0.0 diff --git a/drizzle/tests/test_file_io.py b/drizzle/tests/test_file_io.py deleted file mode 100644 index b44008b4..00000000 --- a/drizzle/tests/test_file_io.py +++ /dev/null @@ -1,173 +0,0 @@ -import os - -import pytest -import numpy as np -from astropy import wcs -from astropy.io import fits - -from drizzle import drizzle -from drizzle import util - -TEST_DIR = os.path.abspath(os.path.dirname(__file__)) -DATA_DIR = os.path.join(TEST_DIR, 'data') - - -def read_header(filename): - """ - Read the primary header from a fits file - """ - fileroot, extn = util.parse_filename(os.path.join(DATA_DIR, filename)) - with fits.open(fileroot) as hdulist: - header = hdulist[0].header - return header - - -def read_image(filename): - """ - Read the image from a fits file - """ - fileroot, extn = util.parse_filename(os.path.join(DATA_DIR, filename)) - with fits.open(fileroot, memmap=False) as hdulist: - data = hdulist[1].data.copy() - return data - - -def read_wcs(filename): - """ - Read the wcs of a fits file - """ - fileroot, extn = util.parse_filename(os.path.join(DATA_DIR, filename)) - with fits.open(fileroot) as hdulist: - the_wcs = wcs.WCS(hdulist[1].header) - return the_wcs - - -@pytest.fixture -def run_drizzle_reference_square_points(tmpdir): - """Create an empty drizzle image""" - output_file = str(tmpdir.join('output_null_run.fits')) - output_template = os.path.join(DATA_DIR, 'reference_square_point.fits') - - output_wcs = read_wcs(output_template) - driz = drizzle.Drizzle(outwcs=output_wcs, wt_scl="expsq", pixfrac=0.5, - kernel="turbo", fillval="NaN") - driz.write(output_file) - - return output_file - - -def test_null_run(run_drizzle_reference_square_points): - output_file = run_drizzle_reference_square_points - with fits.open(output_file) as hdulist: - - assert hdulist.index_of("SCI") == 1 - assert hdulist.index_of("WHT") == 2 - assert hdulist.index_of("CTX") == 3 - - pheader = hdulist["PRIMARY"].header - - assert pheader['DRIZOUDA'] == 'SCI' - assert pheader['DRIZOUWE'] == 'WHT' - assert pheader['DRIZOUCO'] == 'CTX' - assert pheader['DRIZWTSC'] == 'expsq' - assert pheader['DRIZKERN'] == 'turbo' - assert pheader['DRIZPIXF'] == 0.5 - assert pheader['DRIZFVAL'] == 'NaN' - assert pheader['DRIZOUUN'] == 'cps' - assert pheader['EXPTIME'] == 0.0 - assert pheader['DRIZEXPT'] == 1.0 - - -def test_file_init(run_drizzle_reference_square_points): - """ - Initialize drizzle object from a file - """ - input_file = run_drizzle_reference_square_points - - driz = drizzle.Drizzle(infile=input_file) - - assert driz.outexptime == 1.0 - assert driz.wt_scl == 'expsq' - assert driz.kernel == 'turbo' - assert driz.pixfrac == 0.5 - assert driz.fillval == 'NaN' - - -@pytest.fixture -def add_header(tmpdir): - """Add extra keywords read from the header""" - output_file = str(tmpdir.join('output_add_header.fits')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') - output_template = os.path.join(DATA_DIR, 'reference_square_point.fits') - - driz = drizzle.Drizzle(infile=output_template) - image = read_image(input_file) - inwcs = read_wcs(input_file) - driz.add_image(image, inwcs) - - header = fits.header.Header() - header['ONEVAL'] = (1.0, 'test value') - header['TWOVAL'] = (2.0, 'test value') - - driz.write(output_file, outheader=header) - - return output_file - - -def test_add_header(add_header): - output_file = add_header - header = read_header(output_file) - assert header['ONEVAL'] == 1.0 - assert header['TWOVAL'] == 2.0 - assert header['DRIZKERN'] == 'square' - - -def test_add_file(add_header, tmpdir): - """ - Add an image read from a file - """ - test_file = add_header - output_file = str(tmpdir.join('output_add_file.fits')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits[1]') - output_template = os.path.join(DATA_DIR, 'reference_square_point.fits') - - driz = drizzle.Drizzle(infile=output_template) - driz.add_fits_file(input_file) - driz.write(output_file) - - output_image = read_image(output_file) - test_image = read_image(test_file) - diff_image = np.absolute(output_image - test_image) - - assert np.amax(diff_image) == 0.0 - - -def test_blot_file(tmpdir): - """ - Blot an image read from a file - """ - output_file = str(tmpdir.join('output_blot_file.fits')) - test_file = str(tmpdir.join('output_blot_image.fits')) - - input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits[1]') - output_template = os.path.join(DATA_DIR, 'reference_blot_image.fits') - - blotwcs = read_wcs(input_file) - - driz = drizzle.Drizzle(infile=output_template) - driz.add_fits_file(input_file) - driz.blot_image(blotwcs) - driz.write(test_file) - - driz = drizzle.Drizzle(infile=output_template) - driz.add_fits_file(input_file) - driz.blot_fits_file(input_file) - driz.write(output_file) - - output_image = read_image(output_file) - test_image = read_image(test_file) - diff_image = np.absolute(output_image - test_image) - - assert np.amax(diff_image) == 0.0 diff --git a/drizzle/tests/test_overlap_calc.py b/drizzle/tests/test_overlap_calc.py index 9f644bf3..26cf168a 100644 --- a/drizzle/tests/test_overlap_calc.py +++ b/drizzle/tests/test_overlap_calc.py @@ -1,8 +1,8 @@ -import pytest -from math import sqrt from itertools import product +from math import sqrt import numpy as np +import pytest from drizzle.cdrizzle import clip_polygon, invert_pixmap diff --git a/drizzle/tests/test_pixmap.py b/drizzle/tests/test_pixmap.py deleted file mode 100644 index c6215bf6..00000000 --- a/drizzle/tests/test_pixmap.py +++ /dev/null @@ -1,76 +0,0 @@ -import os - -import numpy as np -from numpy.testing import assert_equal, assert_almost_equal -from astropy import wcs -from astropy.io import fits - -from drizzle import calc_pixmap - -TEST_DIR = os.path.abspath(os.path.dirname(__file__)) -DATA_DIR = os.path.join(TEST_DIR, 'data') - - -def test_map_rectangular(): - """ - Make sure the initial index array has correct values - """ - naxis1 = 1000 - naxis2 = 10 - - pixmap = np.indices((naxis1, naxis2), dtype='float32') - pixmap = pixmap.transpose() - - assert_equal(pixmap[5,500], (500,5)) - - -def test_map_to_self(): - """ - Map a pixel array to itself. Should return the same array. - """ - input_file = os.path.join(DATA_DIR, 'input1.fits') - input_hdu = fits.open(input_file) - - input_wcs = wcs.WCS(input_hdu[1].header) - naxis1, naxis2 = input_wcs.pixel_shape - input_hdu.close() - - ok_pixmap = np.indices((naxis1, naxis2), dtype='float32') - ok_pixmap = ok_pixmap.transpose() - - pixmap = calc_pixmap.calc_pixmap(input_wcs, input_wcs) - - # Got x-y transpose right - assert_equal(pixmap.shape, ok_pixmap.shape) - # Mapping an array to itself - assert_almost_equal(pixmap, ok_pixmap, decimal=5) - - -def test_translated_map(): - """ - Map a pixel array to at translated array. - """ - first_file = os.path.join(DATA_DIR, 'input1.fits') - first_hdu = fits.open(first_file) - first_header = first_hdu[1].header - - first_wcs = wcs.WCS(first_header) - naxis1, naxis2 = first_wcs.pixel_shape - first_hdu.close() - - second_file = os.path.join(DATA_DIR, 'input3.fits') - second_hdu = fits.open(second_file) - second_header = second_hdu[1].header - - second_wcs = wcs.WCS(second_header) - second_hdu.close() - - ok_pixmap = np.indices((naxis1, naxis2), dtype='float32') - 2.0 - ok_pixmap = ok_pixmap.transpose() - - pixmap = calc_pixmap.calc_pixmap(first_wcs, second_wcs) - - # Got x-y transpose right - assert_equal(pixmap.shape, ok_pixmap.shape) - # Mapping an array to a translated array - assert_almost_equal(pixmap, ok_pixmap, decimal=5) diff --git a/drizzle/tests/test_resample.py b/drizzle/tests/test_resample.py new file mode 100644 index 00000000..3c02cf0c --- /dev/null +++ b/drizzle/tests/test_resample.py @@ -0,0 +1,1437 @@ +import math +import os + +import numpy as np +import pytest + +from astropy import wcs +from astropy.io import fits +from drizzle import cdrizzle, resample, utils + +TEST_DIR = os.path.abspath(os.path.dirname(__file__)) +DATA_DIR = os.path.join(TEST_DIR, 'data') + + +def bound_image(image): + """ + Compute region where image is non-zero + """ + coords = np.nonzero(image) + ymin = coords[0].min() + ymax = coords[0].max() + xmin = coords[1].min() + xmax = coords[1].max() + return (ymin, ymax, xmin, xmax) + + +def centroid(image, size, center): + """ + Compute the centroid of a rectangular area + """ + ylo = int(center[0] - size / 2) + yhi = min(ylo + size, image.shape[0]) + xlo = int(center[1] - size / 2) + xhi = min(xlo + size, image.shape[1]) + + yx1 = np.mgrid[ylo:yhi, xlo:xhi, 1:2] + center = (yx1[..., 0] * image[ylo:yhi, xlo:xhi]).sum( + axis=(1, 2), + dtype=np.float64, + ) + + if center[2] == 0.0: + return None + + center[0] /= center[2] + center[1] /= center[2] + return center + + +def centroid_close(list_of_centroids, size, point): + """ + Find if any centroid is close to a point + """ + for i in range(len(list_of_centroids) - 1, -1, -1): + if (abs(list_of_centroids[i][0] - point[0]) < int(size / 2) and + abs(list_of_centroids[i][1] - point[1]) < int(size / 2)): + return 1 + + return 0 + + +def centroid_compare(centroid): + return centroid[1] + + +def centroid_distances(image1, image2, amp, size): + """ + Compute a list of centroids and the distances between them in two images + """ + distances = [] + list_of_centroids = centroid_list(image2, amp, size) + for center2 in list_of_centroids: + center1 = centroid(image1, size, center2) + if center1 is None: + continue + + disty = center2[0] - center1[0] + distx = center2[1] - center1[1] + dist = math.sqrt(disty * disty + distx * distx) + dflux = abs(center2[2] - center1[2]) + distances.append([dist, dflux, center1, center2]) + + distances.sort(key=centroid_compare) + return distances + + +def centroid_list(image, amp, size): + """ + Find the next centroid + """ + list_of_centroids = [] + points = np.transpose(np.nonzero(image > amp)) + for point in points: + if not centroid_close(list_of_centroids, size, point): + center = centroid(image, size, point) + list_of_centroids.append(center) + + return list_of_centroids + + +def centroid_statistics(title, fname, image1, image2, amp, size): + """ + write centroid statistics to compare differences btw two images + """ + stats = ("minimum", "median", "maximum") + images = (None, None, image1, image2) + im_type = ("", "", "test", "reference") + + diff = [] + distances = centroid_distances(image1, image2, amp, size) + indexes = (0, int(len(distances) / 2), len(distances) - 1) + fd = open(fname, 'w') + fd.write(f"*** {title:s} ***\n") + + if len(distances) == 0: + diff = [0.0, 0.0, 0.0] + fd.write("No matches!!\n") + + elif len(distances) == 1: + diff = [distances[0][0], distances[0][0], distances[0][0]] + + fd.write("1 match\n") + fd.write( + f"distance = {distances[0][0]:f} " + f"flux difference = {distances[0][1]:f}\n" + ) + + for j in range(2, 4): + ylo = int(distances[0][j][0]) - 1 + yhi = int(distances[0][j][0]) + 2 + xlo = int(distances[0][j][1]) - 1 + xhi = int(distances[0][j][1]) + 2 + subimage = images[j][ylo:yhi, xlo:xhi] + fd.write( + f"\n{im_type[j]} image centroid = " + f"({distances[0][j][0]:f}, {distances[0][j][1]:f}) " + f"image flux = {distances[0][j][2]:f}\n" + ) + fd.write(str(subimage) + "\n") + + else: + fd.write(f"{len(distances)} matches\n") + + for k in range(3): + i = indexes[k] + diff.append(distances[i][0]) + fd.write( + f"\n{stats[k]} distance = {distances[i][0]:f} " + f"flux difference = {distances[i][1]:f}\n" + ) + + for j in range(2, 4): + ylo = int(distances[i][j][0]) - 1 + yhi = int(distances[i][j][0]) + 2 + xlo = int(distances[i][j][1]) - 1 + xhi = int(distances[i][j][1]) + 2 + subimage = images[j][ylo:yhi, xlo:xhi] + fd.write( + f"\n{stats[k]} {im_type[j]} image centroid = " + f"({distances[i][j][0]:f}, {distances[i][j][1]:f}) " + f"image flux = {distances[i][j][2]:f}\n" + ) + fd.write(str(subimage) + "\n") + + fd.close() + return tuple(diff) + + +def make_point_image(input_image, point, value): + """ + Create an image with a single point set + """ + output_image = np.zeros(input_image.shape, dtype=input_image.dtype) + output_image[point] = value + return output_image + + +def make_grid_image(input_image, spacing, value): + """ + Create an image with points on a grid set + """ + output_image = np.zeros(input_image.shape, dtype=input_image.dtype) + + shape = output_image.shape + half_space = int(spacing / 2) + for y in range(half_space, shape[0], spacing): + for x in range(half_space, shape[1], spacing): + output_image[y, x] = value + + return output_image + + +def read_image(filename): + """ + Read the image from a fits file + """ + path = os.path.join(DATA_DIR, filename) + hdu = fits.open(path) + + image = hdu[1].data + hdu.close() + return image + + +def read_wcs(filename): + """ + Read the wcs of a fits file + """ + path = os.path.join(DATA_DIR, filename) + hdu = fits.open(path) + the_wcs = wcs.WCS(hdu[1].header) + hdu.close() + return the_wcs + + +def test_drizzle_defaults(): + n = 200 + in_shape = (n, n) + + # input coordinate grid: + y, x = np.indices(in_shape, dtype=np.float64) + + # simulate data: + in_sci = np.ones(in_shape, dtype=np.float32) + in_wht = np.ones(in_shape, dtype=np.float32) + + # create a Drizzle object using all default parameters (except for 'kernel') + driz = resample.Drizzle( + kernel='square', + ) + + assert driz.out_img is None + assert driz.out_wht is None + assert driz.out_ctx is None + assert driz.total_exptime == 0.0 + + driz.add_image( + in_sci, + exptime=1.0, + pixmap=np.dstack([x, y]), + weight_map=in_wht, + ) + + pixmap = np.dstack([x + 1, y + 2]) + driz.add_image( + 3 * in_sci, + exptime=1.0, + pixmap=pixmap, + weight_map=in_wht, + ) + + assert driz.out_img[0, 0] == 1 + assert driz.out_img[1, 0] == 1 + assert driz.out_img[2, 0] == 1 + assert driz.out_img[1, 1] == 1 + assert driz.out_img[1, 2] == 1 + assert (driz.out_img[2, 1] - 2.0) < 1.0e-14 + + +def test_square_with_point(tmpdir): + """ + Test do_driz square kernel with point + """ + output_difference = str(tmpdir.join('difference_square_point.txt')) + + input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') + output_template = os.path.join(DATA_DIR, 'reference_square_point.fits') + + insci = read_image(input_file) + inwcs = read_wcs(input_file) + insci = make_point_image(insci, (500, 200), 100.0) + inwht = np.ones(insci.shape, dtype=insci.dtype) + output_wcs = read_wcs(output_template) + + pixmap = utils.calc_pixmap( + inwcs, + output_wcs, + ) + + # ignore previous pscale and compute it the old way (only to make + # tests work with old truth files and thus to show that new API gives + # same results when equal definitions of the pixel scale is used): + pscale = np.sqrt( + np.sum(output_wcs.wcs.pc**2, axis=0)[0] / + np.sum(inwcs.wcs.cd**2, axis=0)[0] + ) + + driz = resample.Drizzle( + kernel='square', + out_shape=output_wcs.array_shape, + fillval=0.0, + ) + driz.add_image( + insci, + exptime=1.0, + pixmap=pixmap, + weight_map=inwht, + scale=pscale, + ) + + template_data = read_image(output_template) + + _, med_diff, max_diff = centroid_statistics( + "square with point", + output_difference, + driz.out_img, + template_data, + 20.0, + 8, + ) + + assert med_diff < 1.0e-6 + assert max_diff < 1.0e-5 + + +@pytest.mark.parametrize( + 'kernel,fc', + [ + ('square', True), + ('point', True), + ('turbo', True), + ('lanczos2', False), + ('lanczos3', False), + ('gaussian', False), + ], +) +def test_zero_input_weight(kernel, fc): + """ + Test do_driz square kernel with grid + """ + # initialize input: + insci = np.ones((200, 400), dtype=np.float32) + inwht = np.ones((200, 400), dtype=np.float32) + inwht[:, 150:155] = 0 + + # initialize output: + outsci = np.zeros((210, 410), dtype=np.float32) + outwht = np.zeros((210, 410), dtype=np.float32) + outctx = np.zeros((210, 410), dtype=np.int32) + + # define coordinate mapping: + pixmap = np.moveaxis(np.mgrid[1:201, 1:401][::-1], 0, -1) + + # resample: + if fc: + cdrizzle.tdriz( + insci, + inwht, + pixmap, + outsci, + outwht, + outctx, + uniqid=1, + xmin=0, + xmax=400, + ymin=0, + ymax=200, + pixfrac=1, + kernel=kernel, + in_units='cps', + expscale=1, + wtscale=1, + fillstr='INDEF', + ) + else: + with pytest.warns(Warning): + cdrizzle.tdriz( + insci, + inwht, + pixmap, + outsci, + outwht, + outctx, + uniqid=1, + xmin=0, + xmax=400, + ymin=0, + ymax=200, + pixfrac=1, + kernel=kernel, + in_units='cps', + expscale=1, + wtscale=1, + fillstr='INDEF', + ) + # pytest.xfail("Not a flux-conserving kernel") + + # check that no pixel with 0 weight has any counts: + assert np.sum(np.abs(outsci[(outwht == 0)])) == 0.0 + + +def test_square_with_grid(tmpdir): + """ + Test do_driz square kernel with grid + """ + output_difference = str(tmpdir.join('difference_square_grid.txt')) + + input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') + output_template = os.path.join(DATA_DIR, 'reference_square_grid.fits') + + insci = read_image(input_file) + inwcs = read_wcs(input_file) + insci = make_grid_image(insci, 64, 100.0) + + inwht = np.ones(insci.shape, dtype=insci.dtype) + output_wcs = read_wcs(output_template) + + pixmap = utils.calc_pixmap( + inwcs, + output_wcs, + ) + pscale = utils.estimate_pixel_scale_ratio( + inwcs, + output_wcs, + refpix_from=inwcs.wcs.crpix, + refpix_to=output_wcs.wcs.crpix, + ) + # ignore previous pscale and compute it the old way (only to make + # tests work with old truth files and thus to show that new API gives + # same results when equal definitions of the pixel scale is used): + pscale = np.sqrt( + np.sum(output_wcs.wcs.pc**2, axis=0)[0] / + np.sum(inwcs.wcs.cd**2, axis=0)[0] + ) + + driz = resample.Drizzle( + kernel='square', + out_shape=output_wcs.array_shape, + fillval=0.0, + ) + driz.add_image( + insci, + exptime=1.0, + pixmap=pixmap, + weight_map=inwht, + scale=pscale, + ) + template_data = read_image(output_template) + + _, med_diff, max_diff = centroid_statistics( + "square with grid", + output_difference, + driz.out_img, + template_data, + 20.0, + 8, + ) + assert med_diff < 1.0e-6 + assert max_diff < 1.0e-5 + + +def test_turbo_with_grid(tmpdir): + """ + Test do_driz turbo kernel with grid + """ + output_difference = str(tmpdir.join('difference_turbo_grid.txt')) + + input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') + output_template = os.path.join(DATA_DIR, 'reference_turbo_grid.fits') + + insci = read_image(input_file) + inwcs = read_wcs(input_file) + insci = make_grid_image(insci, 64, 100.0) + inwht = np.ones(insci.shape, dtype=insci.dtype) + output_wcs = read_wcs(output_template) + + pixmap = utils.calc_pixmap( + inwcs, + output_wcs, + ) + pscale = utils.estimate_pixel_scale_ratio( + inwcs, + output_wcs, + refpix_from=inwcs.wcs.crpix, + refpix_to=output_wcs.wcs.crpix, + ) + + # ignore previous pscale and compute it the old way (only to make + # tests work with old truth files and thus to show that new API gives + # same results when equal definitions of the pixel scale is used): + pscale = np.sqrt( + np.sum(output_wcs.wcs.pc**2, axis=0)[0] / + np.sum(inwcs.wcs.cd**2, axis=0)[0] + ) + + driz = resample.Drizzle( + kernel='turbo', + out_shape=output_wcs.array_shape, + fillval=0.0, + ) + driz.add_image( + insci, + exptime=1.0, + pixmap=pixmap, + weight_map=inwht, + scale=pscale, + ) + + template_data = read_image(output_template) + + _, med_diff, max_diff = centroid_statistics( + "turbo with grid", + output_difference, + driz.out_img, + template_data, + 20.0, + 8, + ) + + assert med_diff < 1.0e-6 + assert max_diff < 1.0e-5 + + +def test_gaussian_with_grid(tmpdir): + """ + Test do_driz gaussian kernel with grid + """ + output_difference = str(tmpdir.join('difference_gaussian_grid.txt')) + + input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') + output_template = os.path.join(DATA_DIR, 'reference_gaussian_grid.fits') + + insci = read_image(input_file) + inwcs = read_wcs(input_file) + insci = make_grid_image(insci, 64, 100.0) + inwht = np.ones(insci.shape, dtype=insci.dtype) + output_wcs = read_wcs(output_template) + + pixmap = utils.calc_pixmap( + inwcs, + output_wcs, + ) + pscale = utils.estimate_pixel_scale_ratio( + inwcs, + output_wcs, + refpix_from=inwcs.wcs.crpix, + refpix_to=output_wcs.wcs.crpix, + ) + + # ignore previous pscale and compute it the old way (only to make + # tests work with old truth files and thus to show that new API gives + # same results when equal definitions of the pixel scale is used): + pscale = np.sqrt( + np.sum(output_wcs.wcs.pc**2, axis=0)[0] / + np.sum(inwcs.wcs.cd**2, axis=0)[0] + ) + + driz = resample.Drizzle( + kernel='gaussian', + out_shape=output_wcs.array_shape, + fillval=0.0, + ) + with pytest.warns(Warning): + driz.add_image( + insci, + exptime=1.0, + pixmap=pixmap, + weight_map=inwht, + scale=pscale, + ) + + template_data = read_image(output_template) + + _, med_diff, max_diff = centroid_statistics( + "gaussian with grid", + output_difference, + driz.out_img, + template_data, + 20.0, + 8, + ) + + assert med_diff < 1.0e-6 + assert max_diff < 2.0e-5 + + +def test_lanczos_with_grid(tmpdir): + """ + Test do_driz lanczos kernel with grid + """ + output_difference = str(tmpdir.join('difference_lanczos_grid.txt')) + + input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') + output_template = os.path.join(DATA_DIR, 'reference_lanczos_grid.fits') + + insci = read_image(input_file) + inwcs = read_wcs(input_file) + insci = make_grid_image(insci, 64, 100.0) + inwht = np.ones(insci.shape, dtype=insci.dtype) + output_wcs = read_wcs(output_template) + + pixmap = utils.calc_pixmap( + inwcs, + output_wcs, + ) + pscale = utils.estimate_pixel_scale_ratio( + inwcs, + output_wcs, + refpix_from=inwcs.wcs.crpix, + refpix_to=output_wcs.wcs.crpix, + ) + + # ignore previous pscale and compute it the old way (only to make + # tests work with old truth files and thus to show that new API gives + # same results when equal definitions of the pixel scale is used): + pscale = np.sqrt( + np.sum(output_wcs.wcs.pc**2, axis=0)[0] / + np.sum(inwcs.wcs.cd**2, axis=0)[0] + ) + + driz = resample.Drizzle( + kernel='lanczos3', + out_shape=output_wcs.array_shape, + fillval=0.0, + ) + with pytest.warns(Warning): + driz.add_image( + insci, + exptime=1.0, + pixmap=pixmap, + weight_map=inwht, + scale=pscale, + ) + + template_data = read_image(output_template) + + _, med_diff, max_diff = centroid_statistics( + "lanczos with grid", + output_difference, + driz.out_img, + template_data, + 20.0, + 8, + ) + assert med_diff < 1.0e-6 + assert max_diff < 1.0e-5 + + +def test_point_with_grid(tmpdir): + """ + Test do_driz point kernel with grid + """ + output_difference = str(tmpdir.join('difference_point_grid.txt')) + + input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') + output_template = os.path.join(DATA_DIR, 'reference_point_grid.fits') + + insci = read_image(input_file) + inwcs = read_wcs(input_file) + insci = make_grid_image(insci, 64, 100.0) + inwht = np.ones(insci.shape, dtype=insci.dtype) + output_wcs = read_wcs(output_template) + + pixmap = utils.calc_pixmap(inwcs, output_wcs) + + driz = resample.Drizzle(kernel='point', out_shape=output_wcs.array_shape, fillval=0.0) + driz.add_image(insci, exptime=1.0, pixmap=pixmap, weight_map=inwht) + + template_data = read_image(output_template) + + _, med_diff, max_diff = centroid_statistics( + "point with grid", + output_difference, + driz.out_img, + template_data, + 20.0, + 8, + ) + assert med_diff < 1.0e-6 + assert max_diff < 1.0e-5 + + +def test_blot_with_point(tmpdir): + """ + Test do_blot with point image + """ + output_difference = str(tmpdir.join('difference_blot_point.txt')) + + input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') + output_template = os.path.join(DATA_DIR, 'reference_blot_point.fits') + + outsci = read_image(input_file) + outwcs = read_wcs(input_file) + outsci = make_point_image(outsci, (500, 200), 40.0) + inwcs = read_wcs(output_template) + + pixmap = utils.calc_pixmap(inwcs, outwcs) + + # compute pscale the old way (only to make + # tests work with old truth files and thus to show that new API gives + # same results when equal definitions of the pixel scale is used): + pscale = np.sqrt( + np.sum(inwcs.wcs.pc**2, axis=0)[0] / + np.sum(outwcs.wcs.cd**2, axis=0)[0] + ) + + blotted_image = resample.blot_image( + outsci, + pixmap=pixmap, + pix_ratio=pscale, + exptime=1.0, + output_pixel_shape=inwcs.pixel_shape, + ) + + template_data = read_image(output_template) + + _, med_diff, max_diff = centroid_statistics( + "blot with point", + output_difference, + blotted_image, + template_data, + 20.0, + 16, + ) + assert med_diff < 1.0e-6 + assert max_diff < 1.0e-5 + + +def test_blot_with_default(tmpdir): + """ + Test do_blot with default grid image + """ + output_difference = str(tmpdir.join('difference_blot_default.txt')) + + input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') + output_template = os.path.join(DATA_DIR, 'reference_blot_default.fits') + + outsci = read_image(input_file) + outsci = make_grid_image(outsci, 64, 100.0) + outwcs = read_wcs(input_file) + inwcs = read_wcs(output_template) + + pixmap = utils.calc_pixmap(inwcs, outwcs) + + # compute pscale the old way (only to make + # tests work with old truth files and thus to show that new API gives + # same results when equal definitions of the pixel scale is used): + pscale = np.sqrt( + np.sum(inwcs.wcs.pc**2, axis=0)[0] / + np.sum(outwcs.wcs.cd**2, axis=0)[0] + ) + + blotted_image = resample.blot_image( + outsci, + pixmap=pixmap, + pix_ratio=pscale, + exptime=1.0, + output_pixel_shape=inwcs.pixel_shape, + ) + + template_data = read_image(output_template) + + _, med_diff, max_diff = centroid_statistics( + "blot with defaults", + output_difference, + blotted_image, + template_data, + 20.0, + 16, + ) + + assert med_diff < 1.0e-6 + assert max_diff < 1.0e-5 + + +def test_blot_with_lan3(tmpdir): + """ + Test do_blot with lan3 grid image + """ + output_difference = str(tmpdir.join('difference_blot_lan3.txt')) + + input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') + output_template = os.path.join(DATA_DIR, 'reference_blot_lan3.fits') + + outsci = read_image(input_file) + outsci = make_grid_image(outsci, 64, 100.0) + outwcs = read_wcs(input_file) + inwcs = read_wcs(output_template) + + pixmap = utils.calc_pixmap(inwcs, outwcs) + + # compute pscale the old way (only to make + # tests work with old truth files and thus to show that new API gives + # same results when equal definitions of the pixel scale is used): + pscale = np.sqrt( + np.sum(inwcs.wcs.pc**2, axis=0)[0] / + np.sum(outwcs.wcs.cd**2, axis=0)[0] + ) + + blotted_image = resample.blot_image( + outsci, + pixmap=pixmap, + pix_ratio=pscale, + exptime=1.0, + output_pixel_shape=inwcs.pixel_shape, + interp="lan3", + ) + + template_data = read_image(output_template) + + _, med_diff, max_diff = centroid_statistics( + "blot with lan3", + output_difference, + blotted_image, + template_data, + 20.0, + 16, + ) + + assert med_diff < 1.0e-6 + assert max_diff < 1.0e-5 + + +def test_blot_with_lan5(tmpdir): + """ + Test do_blot with lan5 grid image + """ + output_difference = str(tmpdir.join('difference_blot_lan5.txt')) + + input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') + output_template = os.path.join(DATA_DIR, 'reference_blot_lan5.fits') + + outsci = read_image(input_file) + outsci = make_grid_image(outsci, 64, 100.0) + outwcs = read_wcs(input_file) + inwcs = read_wcs(output_template) + + pixmap = utils.calc_pixmap(inwcs, outwcs) + + # compute pscale the old way (only to make + # tests work with old truth files and thus to show that new API gives + # same results when equal definitions of the pixel scale is used): + pscale = np.sqrt( + np.sum(inwcs.wcs.pc**2, axis=0)[0] / + np.sum(outwcs.wcs.cd**2, axis=0)[0] + ) + + blotted_image = resample.blot_image( + outsci, + pixmap=pixmap, + pix_ratio=pscale, + exptime=1.0, + output_pixel_shape=inwcs.pixel_shape, + interp="lan5", + ) + + template_data = read_image(output_template) + + _, med_diff, max_diff = centroid_statistics( + "blot with lan5", + output_difference, + blotted_image, + template_data, + 20.0, + 16, + ) + + assert med_diff < 1.0e-6 + assert max_diff < 1.0e-5 + + +def test_context_planes(): + """Reproduce error seen in issue #50""" + shape = (10, 10) + output_wcs = wcs.WCS() + output_wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] + output_wcs.wcs.pc = [[1, 0], [0, 1]] + output_wcs.pixel_shape = shape + driz = resample.Drizzle(out_shape=tuple(shape)) + + image = np.ones(shape) + inwcs = wcs.WCS() + inwcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] + inwcs.wcs.cd = [[1, 0], [0, 1]] + inwcs.pixel_shape = shape + + pixmap = utils.calc_pixmap(inwcs, output_wcs) + + # context image must be 2D or 3D: + with pytest.raises(ValueError) as err_info: + resample.Drizzle( + kernel='point', + exptime=0.0, + out_shape=shape, + out_ctx=[0, 0, 0], + ) + assert str(err_info.value).startswith( + "'out_ctx' must be either a 2D or 3D array." + ) + + driz = resample.Drizzle( + kernel='square', + out_shape=output_wcs.array_shape, + fillval=0.0, + ) + + for i in range(32): + assert driz.ctx_id == i + driz.add_image(image, exptime=1.0, pixmap=pixmap) + assert driz.out_ctx.shape == (1, 10, 10) + + driz.add_image(image, exptime=1.0, pixmap=pixmap) + assert driz.out_ctx.shape == (2, 10, 10) + + +def test_no_context_image(): + """Reproduce error seen in issue #50""" + shape = (10, 10) + output_wcs = wcs.WCS() + output_wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] + output_wcs.wcs.pc = [[1, 0], [0, 1]] + output_wcs.pixel_shape = shape + driz = resample.Drizzle( + out_shape=tuple(shape), + begin_ctx_id=-1, + disable_ctx=True + ) + assert driz.out_ctx is None + assert driz.ctx_id is None + + image = np.ones(shape) + inwcs = wcs.WCS() + inwcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] + inwcs.wcs.cd = [[1, 0], [0, 1]] + inwcs.pixel_shape = shape + + pixmap = utils.calc_pixmap(inwcs, output_wcs) + + for i in range(33): + driz.add_image(image, exptime=1.0, pixmap=pixmap) + assert driz.out_ctx is None + assert driz.ctx_id is None + + +def test_init_ctx_id(): + # starting context ID must be positive + with pytest.raises(ValueError) as err_info: + resample.Drizzle( + kernel='square', + exptime=0.0, + begin_ctx_id=-1, + out_shape=(10, 10), + ) + assert str(err_info.value).startswith( + "Invalid context image ID" + ) + + with pytest.raises(ValueError) as err_info: + resample.Drizzle( + kernel='square', + exptime=0.0, + out_shape=(10, 10), + begin_ctx_id=1, + max_ctx_id=0, + ) + assert str(err_info.value).startswith( + "'max_ctx_id' cannot be smaller than 'begin_ctx_id'." + ) + + +def test_context_agrees_with_weight(): + n = 200 + out_shape = (n, n) + + # allocate output arrays: + out_img = np.zeros(out_shape, dtype=np.float32) + out_ctx = np.zeros(out_shape, dtype=np.int32) + out_wht = np.zeros(out_shape, dtype=np.float32) + + # previous data in weight and context must agree: + with pytest.raises(ValueError) as err_info: + out_ctx[0, 0] = 1 + out_ctx[0, 1] = 1 + out_wht[0, 0] = 0.1 + resample.Drizzle( + kernel='square', + out_shape=out_shape, + out_img=out_img, + out_ctx=out_ctx, + out_wht=out_wht, + exptime=1.0, + ) + assert str(err_info.value).startswith( + "Inconsistent values of supplied 'out_wht' and 'out_ctx' " + ) + + +@pytest.mark.parametrize( + 'kernel,fc', + [ + ('square', True), + ('point', True), + ('turbo', True), + ('lanczos2', False), + ('lanczos3', False), + ('gaussian', False), + ], +) +def test_flux_conservation_nondistorted(kernel, fc): + n = 200 + in_shape = (n, n) + + # input coordinate grid: + y, x = np.indices(in_shape, dtype=np.float64) + + # simulate a gaussian "star": + fwhm = 2.9 + x0 = 50.0 + y0 = 68.0 + sig = fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0 * fwhm))) + sig2 = sig * sig + star = np.exp(-0.5 / sig2 * ((x.astype(np.float32) - x0)**2 + (y.astype(np.float32) - y0)**2)) + in_sci = (star / np.sum(star)).astype(np.float32) # normalize to 1 + in_wht = np.ones(in_shape, dtype=np.float32) + + # linear shift: + xp = x + 0.5 + yp = y + 0.2 + + pixmap = np.dstack([xp, yp]) + + out_shape = (int(yp.max()) + 1, int(xp.max()) + 1) + # make sure distorion is not moving flux out of the image towards negative + # coordinates (just because of the simple way of how we account for output + # image size) + assert np.min(xp) > -0.5 and np.min(yp) > -0.5 + + out_img = np.zeros(out_shape, dtype=np.float32) + out_ctx = np.zeros(out_shape, dtype=np.int32) + out_wht = np.zeros(out_shape, dtype=np.float32) + + if fc: + cdrizzle.tdriz( + in_sci, + in_wht, + pixmap, + out_img, + out_wht, + out_ctx, + pixfrac=1.0, + scale=1.0, + kernel=kernel, + in_units="cps", + expscale=1.0, + wtscale=1.0, + ) + else: + with pytest.warns(Warning): + cdrizzle.tdriz( + in_sci, + in_wht, + pixmap, + out_img, + out_wht, + out_ctx, + pixfrac=1.0, + scale=1.0, + kernel=kernel, + in_units="cps", + expscale=1.0, + wtscale=1.0, + ) + pytest.xfail("Not a flux-conserving kernel") + + assert np.allclose( + np.sum(out_img * out_wht), + np.sum(in_sci), + atol=0.0, + rtol=0.0001, + ) + + +@pytest.mark.parametrize( + 'kernel,fc', + [ + ('square', True), + ('point', True), + ('turbo', True), + ('lanczos2', False), + ('lanczos3', False), + ('gaussian', False), + ], +) +def test_flux_conservation_distorted(kernel, fc): + n = 200 + in_shape = (n, n) + + # input coordinate grid: + y, x = np.indices(in_shape, dtype=np.float64) + + # simulate a gaussian "star": + fwhm = 2.9 + x0 = 50.0 + y0 = 68.0 + sig = fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0 * fwhm))) + sig2 = sig * sig + star = np.exp(-0.5 / sig2 * ((x.astype(np.float32) - x0)**2 + (y.astype(np.float32) - y0)**2)) + in_sci = (star / np.sum(star)).astype(np.float32) # normalize to 1 + in_wht = np.ones(in_shape, dtype=np.float32) + + # linear shift: + xp = x + 0.5 + yp = y + 0.2 + # add distortion: + xp += 1e-4 * x**2 + 1e-5 * x * y + yp += 1e-3 * y**2 - 2e-5 * x * y + + pixmap = np.dstack([xp, yp]) + + out_shape = (int(yp.max()) + 1, int(xp.max()) + 1) + # make sure distorion is not moving (pixels with) flux out of the image + # towards negative coordinates (just because of the simple way of how we + # account for output image size): + assert np.min(xp) > -0.5 and np.min(yp) > -0.5 + + out_img = np.zeros(out_shape, dtype=np.float32) + out_ctx = np.zeros(out_shape, dtype=np.int32) + out_wht = np.zeros(out_shape, dtype=np.float32) + + if fc: + cdrizzle.tdriz( + in_sci, + in_wht, + pixmap, + out_img, + out_wht, + out_ctx, + pixfrac=1.0, + scale=1.0, + kernel=kernel, + in_units="cps", + expscale=1.0, + wtscale=1.0, + ) + else: + with pytest.warns(Warning): + cdrizzle.tdriz( + in_sci, + in_wht, + pixmap, + out_img, + out_wht, + out_ctx, + pixfrac=1.0, + scale=1.0, + kernel=kernel, + in_units="cps", + expscale=1.0, + wtscale=1.0, + ) + pytest.xfail("Not a flux-conserving kernel") + + assert np.allclose( + np.sum(out_img * out_wht), + np.sum(in_sci), + atol=0.0, + rtol=0.0001, + ) + + +def test_drizzle_exptime(): + n = 200 + in_shape = (n, n) + + # input coordinate grid: + y, x = np.indices(in_shape, dtype=np.float64) + + # simulate data: + in_sci = np.ones(in_shape, dtype=np.float32) + in_wht = np.ones(in_shape, dtype=np.float32) + + pixmap = np.dstack([x, y]) + + # allocate output arrays: + out_shape = (int(y.max()) + 1, int(x.max()) + 1) + out_img = np.zeros(out_shape, dtype=np.float32) + out_ctx = np.zeros(out_shape, dtype=np.int32) + out_wht = np.zeros(out_shape, dtype=np.float32) + + # starting exposure time must be non-negative: + with pytest.raises(ValueError) as err_info: + driz = resample.Drizzle( + kernel='square', + out_shape=out_shape, + fillval="indef", + exptime=-1.0, + ) + assert str(err_info.value) == "Exposure time must be non-negative." + + driz = resample.Drizzle( + kernel='turbo', + out_shape=out_shape, + fillval="", + out_img=out_img, + out_ctx=out_ctx, + out_wht=out_wht, + exptime=1.0, + ) + assert driz.kernel == 'turbo' + + driz.add_image(in_sci, weight_map=in_wht, exptime=1.03456, pixmap=pixmap) + assert np.allclose(driz.total_exptime, 2.03456, rtol=0, atol=1.0e-14) + + driz.add_image(in_sci, weight_map=in_wht, exptime=3.1415926, pixmap=pixmap) + assert np.allclose(driz.total_exptime, 5.1761526, rtol=0, atol=1.0e-14) + + with pytest.raises(ValueError) as err_info: + driz.add_image(in_sci, weight_map=in_wht, exptime=-1, pixmap=pixmap) + assert str(err_info.value) == "'exptime' *must* be a strictly positive number." + + # exptime cannot be 0 when output data has data: + with pytest.raises(ValueError) as err_info: + out_ctx[0, 0] = 1 + driz = resample.Drizzle( + kernel='square', + out_shape=out_shape, + fillval="indef", + out_img=out_img, + out_ctx=out_ctx, + out_wht=out_wht, + exptime=0.0, + ) + assert str(err_info.value).startswith( + "Inconsistent exposure time and context and/or weight images:" + ) + + # exptime must be 0 when output arrays are not provided: + with pytest.raises(ValueError) as err_info: + driz = resample.Drizzle( + kernel='square', + out_shape=out_shape, + exptime=1.0, + ) + assert str(err_info.value).startswith( + "Exposure time must be 0.0 for the first resampling" + ) + + +def test_drizzle_unsupported_kernel(): + with pytest.raises(ValueError) as err_info: + resample.Drizzle( + kernel='magic_image_improver', + out_shape=(10, 10), + exptime=0.0, + ) + assert str(err_info.value) == "Kernel 'magic_image_improver' is not supported." + + +def test_pixmap_shape_matches_image(): + n = 200 + in_shape = (n, n) + + # input coordinate grid: + y, x = np.indices((n + 1, n), dtype=np.float64) + + # simulate data: + in_sci = np.ones(in_shape, dtype=np.float32) + in_wht = np.ones(in_shape, dtype=np.float32) + + pixmap = np.dstack([x, y]) + + driz = resample.Drizzle( + kernel='square', + fillval=0.0, + exptime=0.0, + ) + + # last two sizes of the pixelmap must match those of input images: + with pytest.raises(ValueError) as err_info: + driz.add_image( + in_sci, + exptime=1.0, + pixmap=pixmap, + weight_map=in_wht, + ) + assert str(err_info.value) == "'pixmap' shape is not consistent with 'data' shape." + + +def test_drizzle_fillval(): + n = 200 + in_shape = (n, n) + + # input coordinate grid: + y, x = np.indices(in_shape, dtype=np.float64) + + # simulate a gaussian "star": + fwhm = 2.9 + x0 = 50.0 + y0 = 68.0 + sig = fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0 * fwhm))) + sig2 = sig * sig + star = np.exp(-0.5 / sig2 * ((x.astype(np.float32) - x0)**2 + (y.astype(np.float32) - y0)**2)) + in_sci = (star / np.sum(star)).astype(np.float32) # normalize to 1 + in_wht = np.zeros(in_shape, dtype=np.float32) + mask = np.where((x.astype(np.float32) - x0)**2 + (y.astype(np.float32) - y0)**2 <= 10) + in_wht[mask] = 1.0 + + # linear shift: + xp = x + 50 + yp = y + 50 + + pixmap = np.dstack([xp, yp]) + + out_shape = (int(yp.max()) + 1, int(xp.max()) + 1) + # make sure distorion is not moving flux out of the image towards negative + # coordinates (just because of the simple way of how we account for output + # image size) + assert np.min(xp) > -0.5 and np.min(yp) > -0.5 + + out_img = np.zeros(out_shape, dtype=np.float32) - 1.11 + out_ctx = np.zeros((1, ) + out_shape, dtype=np.int32) + out_wht = np.zeros(out_shape, dtype=np.float32) + + driz = resample.Drizzle( + kernel='square', + out_shape=out_shape, + fillval="indef", + exptime=0.0, + ) + + driz.add_image(in_sci, weight_map=in_wht, exptime=1.0, pixmap=pixmap) + assert np.isnan(driz.out_img[0, 0]) + assert driz.out_img[int(y0) + 50, int(x0) + 50] > 0.0 + + driz = resample.Drizzle( + kernel='square', + out_shape=out_shape, + fillval="-1.11", + out_img=out_img.copy(), + out_ctx=out_ctx.copy(), + out_wht=out_wht.copy(), + exptime=0.0, + ) + driz.add_image(in_sci, weight_map=in_wht, exptime=1.0, pixmap=pixmap) + assert np.allclose(driz.out_img[0, 0], -1.11, rtol=0.0, atol=1.0e-7) + assert driz.out_img[int(y0) + 50, int(x0) + 50] > 0.0 + assert set(driz.out_ctx.ravel().tolist()) == {0, 1} + + # test same with numeric fillval: + driz = resample.Drizzle( + kernel='square', + out_shape=out_shape, + fillval=-1.11, + out_img=out_img.copy(), + out_ctx=out_ctx.copy(), + out_wht=out_wht.copy(), + exptime=0.0, + ) + driz.add_image(in_sci, weight_map=in_wht, exptime=1.0, pixmap=pixmap) + assert np.allclose(driz.out_img[0, 0], -1.11, rtol=0.0, atol=1.0e-7) + assert np.allclose(float(driz.fillval), -1.11, rtol=0.0, atol=np.finfo(float).eps) + + # make sure code raises exception for unsuported fillval: + with pytest.raises(ValueError) as err_info: + resample.Drizzle( + kernel='square', + out_shape=out_shape, + fillval="fillval", + exptime=0.0, + ) + assert str(err_info.value) == "could not convert string to float: 'fillval'" + + +def test_resample_get_shape_from_pixmap(): + n = 200 + in_shape = (n, n) + + # input coordinate grid: + y, x = np.indices(in_shape, dtype=np.float64) + + # simulate constant data: + in_sci = np.ones(in_shape, dtype=np.float32) + in_wht = np.ones(in_shape, dtype=np.float32) + + pixmap = np.dstack([x, y]) + + driz = resample.Drizzle( + kernel='point', + exptime=0.0, + ) + + driz.add_image(in_sci, weight_map=in_wht, exptime=1.0, pixmap=pixmap) + assert driz.out_img.shape == in_shape + + +def test_resample_counts_units(): + n = 200 + in_shape = (n, n) + + # input coordinate grid: + y, x = np.indices(in_shape, dtype=np.float64) + pixmap = np.dstack([x, y]) + + # simulate constant data: + in_sci = np.ones(in_shape, dtype=np.float32) + in_wht = np.ones(in_shape, dtype=np.float32) + + driz = resample.Drizzle() + driz.add_image(in_sci, weight_map=in_wht, exptime=1.0, pixmap=pixmap, in_units='cps') + cps_max_val = driz.out_img.max() + + driz = resample.Drizzle() + driz.add_image(in_sci, weight_map=in_wht, exptime=2.0, pixmap=pixmap, in_units='counts') + counts_max_val = driz.out_img.max() + + assert abs(counts_max_val - cps_max_val / 2.0) < 1.0e-14 + + +def test_resample_inconsistent_output(): + n = 200 + out_shape = (n, n) + + # different shapes: + out_img = np.zeros((n, n), dtype=np.float32) + out_ctx = np.zeros((1, n, n + 1), dtype=np.int32) + out_wht = np.zeros((n + 1, n + 1), dtype=np.float32) + + # shape from out_img: + driz = resample.Drizzle( + kernel='point', + exptime=0.0, + out_img=out_img, + ) + assert driz.out_img.shape == out_shape + + # inconsistent shapes: + out_shape = (n + 1, n) + with pytest.raises(ValueError) as err_info: + resample.Drizzle( + kernel='point', + exptime=0.0, + out_shape=out_shape, + out_img=out_img, + out_ctx=out_ctx, + out_wht=out_wht, + ) + assert str(err_info.value).startswith("Inconsistent data shapes specified") diff --git a/drizzle/tests/test_utils.py b/drizzle/tests/test_utils.py new file mode 100644 index 00000000..68dac2bb --- /dev/null +++ b/drizzle/tests/test_utils.py @@ -0,0 +1,193 @@ +import os + +import numpy as np +import pytest +from numpy.testing import assert_almost_equal, assert_equal + +from astropy import wcs +from astropy.io import fits +from drizzle.utils import ( + _estimate_pixel_scale, + calc_pixmap, + decode_context, + estimate_pixel_scale_ratio, +) + +TEST_DIR = os.path.abspath(os.path.dirname(__file__)) +DATA_DIR = os.path.join(TEST_DIR, 'data') + + +def test_map_rectangular(): + """ + Make sure the initial index array has correct values + """ + naxis1 = 1000 + naxis2 = 10 + + pixmap = np.indices((naxis1, naxis2), dtype='float32') + pixmap = pixmap.transpose() + + assert_equal(pixmap[5, 500], (500, 5)) + + +def test_map_to_self(): + """ + Map a pixel array to itself. Should return the same array. + """ + input_file = os.path.join(DATA_DIR, 'input1.fits') + input_hdu = fits.open(input_file) + + input_wcs = wcs.WCS(input_hdu[1].header) + naxis1, naxis2 = input_wcs.pixel_shape + input_hdu.close() + + ok_pixmap = np.indices((naxis1, naxis2), dtype='float32') + ok_pixmap = ok_pixmap.transpose() + + pixmap = calc_pixmap(input_wcs, input_wcs) + + # Got x-y transpose right + assert_equal(pixmap.shape, ok_pixmap.shape) + + # Mapping an array to itself + assert_almost_equal(pixmap, ok_pixmap, decimal=5) + + # user-provided shape + pixmap = calc_pixmap(input_wcs, input_wcs, (12, 34)) + assert_equal(pixmap.shape, (12, 34, 2)) + + # Check that an exception is raised for WCS without pixel_shape or + # bounding_box: + input_wcs.pixel_shape = None + with pytest.raises(ValueError): + calc_pixmap(input_wcs, input_wcs) + + # user-provided shape when array_shape is not set: + pixmap = calc_pixmap(input_wcs, input_wcs, (12, 34)) + assert_equal(pixmap.shape, (12, 34, 2)) + + # from bounding box: + input_wcs.bounding_box = ((5.3, 33.5), (2.8, 11.5)) + pixmap = calc_pixmap(input_wcs, input_wcs) + assert_equal(pixmap.shape, (12, 34, 2)) + + # from bounding box and pixel_shape (the later takes precedence): + input_wcs.pixel_shape = (naxis1, naxis2) + pixmap = calc_pixmap(input_wcs, input_wcs) + assert_equal(pixmap.shape, ok_pixmap.shape) + + +def test_translated_map(): + """ + Map a pixel array to at translated array. + """ + first_file = os.path.join(DATA_DIR, 'input1.fits') + first_hdu = fits.open(first_file) + first_header = first_hdu[1].header + + first_wcs = wcs.WCS(first_header) + naxis1, naxis2 = first_wcs.pixel_shape + first_hdu.close() + + second_file = os.path.join(DATA_DIR, 'input3.fits') + second_hdu = fits.open(second_file) + second_header = second_hdu[1].header + + second_wcs = wcs.WCS(second_header) + second_hdu.close() + + ok_pixmap = np.indices((naxis1, naxis2), dtype='float32') - 2.0 + ok_pixmap = ok_pixmap.transpose() + + pixmap = calc_pixmap(first_wcs, second_wcs) + + # Got x-y transpose right + assert_equal(pixmap.shape, ok_pixmap.shape) + # Mapping an array to a translated array + assert_almost_equal(pixmap, ok_pixmap, decimal=5) + + +def test_estimate_pixel_scale_ratio(): + input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') + + with fits.open(input_file) as h: + w = wcs.WCS(h[1].header) + + pscale = estimate_pixel_scale_ratio(w, w, w.wcs.crpix, (0, 0)) + + assert abs(pscale - 0.9999999916964737) < 1.0e-9 + + +def test_estimate_pixel_scale_no_refpix(): + # create a WCS without higher order (polynomial) distortions: + fits_file = os.path.join(DATA_DIR, 'input1.fits') + with fits.open(fits_file) as h: + w = wcs.WCS(h[1].header, h) + w.sip = None + w.det2im1 = None + w.det2im2 = None + w.cpdis1 = None + w.cpdis2 = None + pixel_shape = w.pixel_shape[:] + + ref_pscale = _estimate_pixel_scale(w, w.wcs.crpix) + + if hasattr(w, 'bounding_box'): + del w.bounding_box + pscale1 = _estimate_pixel_scale(w, None) + assert np.allclose(ref_pscale, pscale1, atol=0.0, rtol=1.0e-8) + + w.bounding_box = None + w.pixel_shape = None + pscale2 = _estimate_pixel_scale(w, None) + assert np.allclose(pscale1, pscale2, atol=0.0, rtol=1.0e-8) + + w.pixel_shape = pixel_shape + pscale3 = _estimate_pixel_scale(w, None) + assert np.allclose(pscale1, pscale3, atol=0.0, rtol=1.0e-14) + + w.bounding_box = ((-0.5, pixel_shape[0] - 0.5), (-0.5, pixel_shape[1] - 0.5)) + pscale4 = _estimate_pixel_scale(w, None) + assert np.allclose(pscale3, pscale4, atol=0.0, rtol=1.0e-8) + + +def test_decode_context(): + ctx = np.array( + [[[0, 0, 0, 0, 0, 0], + [0, 0, 0, 36196864, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 537920000, 0, 0, 0]], + [[0, 0, 0, 0, 0, 0,], + [0, 0, 0, 67125536, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 163856, 0, 0, 0]], + [[0, 0, 0, 0, 0, 0], + [0, 0, 0, 8203, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 32865, 0, 0, 0]]], + dtype=np.int32 + ) + + idx1, idx2 = decode_context(ctx, [3, 2], [1, 4]) + + assert sorted(idx1) == [9, 12, 14, 19, 21, 25, 37, 40, 46, 58, 64, 65, 67, 77] + assert sorted(idx2) == [9, 20, 29, 36, 47, 49, 64, 69, 70, 79] + + # context array must be 3D: + with pytest.raises(ValueError): + decode_context(ctx[0], [3, 2], [1, 4]) + + # pixel coordinates must be integer: + with pytest.raises(ValueError): + decode_context(ctx, [3.0, 2], [1, 4]) + + # coordinate lists must be equal in length: + with pytest.raises(ValueError): + decode_context(ctx, [3, 2], [1, 4, 5]) + + # coordinate lists must be 1D: + with pytest.raises(ValueError): + decode_context(ctx, [[3, 2]], [[1, 4]]) diff --git a/drizzle/util.py b/drizzle/util.py index 090de2ba..e9c05d11 100644 --- a/drizzle/util.py +++ b/drizzle/util.py @@ -1,138 +1,15 @@ -import numpy as np +""" +Module ``util`` has been deprecated. +""" +import warnings - -def find_keyword_extn(fimg, keyword, value=None): - """ - This function will return the index of the extension in a multi-extension - FITS file which contains the desired keyword with the given value. - - Parameters - ---------- - - fimg : hdulist - A list of header data units - - keyword : str - The keyword to search for - - value : str or number, optional - If set, the value the keyword must have to match - - Returns - ------- - - The index of the extension - """ - - i = 0 - extnum = -1 - # Search through all the extensions in the FITS object - for chip in fimg: - hdr = chip.header - # Check to make sure the extension has the given keyword - if keyword in hdr: - if value is not None: - # If it does, then does the value match the desired value - # MUST use 'str.strip' to match against any input string! - if hdr[keyword].strip() == value: - extnum = i - break - else: - extnum = i - break - i += 1 - # Return the index of the extension which contained the - # desired EXTNAME value. - return extnum - - -def get_extn(fimg, extn=''): - """ - Returns the FITS extension corresponding to extension specified in - filename. Defaults to returning the first extension with data or the - primary extension, if none have data. - - Parameters - ---------- - - fimg : hdulist - A list of header data units - - extn : str - The extension name and version to match - - Returns - ------- - - The matching header data unit - """ - - if extn: - try: - _extn = parse_extn(extn) - if _extn[0]: - _e = fimg.index_of(_extn) - else: - _e = _extn[1] - - except KeyError: - _e = None - - if _e is None: - _extn = None - else: - _extn = fimg[_e] - - else: - # If no extension is provided, search for first extension - # in FITS file with data associated with it. - - # Set up default to point to PRIMARY extension. - _extn = fimg[0] - # then look for first extension with data. - for _e in fimg: - if _e.data is not None: - _extn = _e - break - - return _extn - - -def get_keyword(fimg, keyword, default=None): - """ - Return a keyword value from the header of an image, - or the default if the keyword is not found. - - Parameters - ---------- - - fimg : hdulist - A list of header data units - - keyword : hdulist - The keyword value to search for - - default : str or number, optional - The default value if not found - - Returns - ------- - - The value if found or default if not - """ - - value = None - if keyword: - _nextn = find_keyword_extn(fimg, keyword) - try: - value = fimg[_nextn].header[keyword] - except KeyError: - value = None - - if value is None and default is not None: - value = default - - return value +warnings.warn( + "Module 'drizzle.util' has been deprecated since version 2.0.0 " + "and it will be removed in a future release. " + "Please replace calls to 'util.is_blank()' with alternative " + "implementation.", + DeprecationWarning +) def is_blank(value): @@ -141,116 +18,17 @@ def is_blank(value): Parameters ---------- - value : str The value to check Returns ------- - True or False """ + warnings.warn( + "'is_blank()' has been deprecated since version 2.0.0 " + "and it will be removed in a future release. " + "Please replace calls to 'is_blank()' with alternative implementation.", + DeprecationWarning + ) return value.strip() == "" - - -def parse_extn(extn=''): - """ - Parse a string representing a qualified fits extension name as in the - output of parse_filename and return a tuple (str(extname), int(extver)), - which can be passed to pyfits functions using the 'ext' kw. - Default return is the first extension in a fits file. - - Examples - -------- - >>> parse_extn('sci,2') - ('sci', 2) - >>> parse_extn('2') - ('', 2) - >>> parse_extn('sci') - ('sci', 1) - - Parameters - ---------- - extn : str - The extension name - - Returns - ------- - A tuple of the extension name and value - """ - if not extn: - return ('', 0) - - try: - lext = extn.split(',') - except Exception: - return ('', 1) - - if len(lext) == 1 and lext[0].isdigit(): - return ("", int(lext[0])) - elif len(lext) == 2: - return (lext[0], int(lext[1])) - else: - return (lext[0], 1) - - -def parse_filename(filename): - """ - Parse out filename from any specified extensions. - Returns rootname and string version of extension name. - - Parameters - ---------- - - filename : str - The filename to be parsed - - Returns - ------- - - A tuple with the filename root and extension - """ - # Parse out any extension specified in filename - _indx = filename.find('[') - if _indx > 0: - # Read extension name provided - _fname = filename[:_indx] - _extn = filename[_indx + 1:-1] - else: - _fname = filename - _extn = '' - - return _fname, _extn - - -def set_pscale(the_wcs): - """ - Calculates the plate scale from cdelt and the pc matrix and adds - it to the WCS. Plate scale is not part of the WCS standard, but is - required by the drizzle code - - Parameters - ---------- - - the_wcs : wcs - A WCS object - """ - try: - cdelt = the_wcs.wcs.get_cdelt() - pc = the_wcs.wcs.get_pc() - - except Exception: - try: - # for non-celestial axes, get_cdelt doesnt work - cdelt = the_wcs.wcs.cd * the_wcs.wcs.cdelt - except AttributeError: - cdelt = the_wcs.wcs.cdelt - - try: - pc = the_wcs.wcs.pc - except AttributeError: - pc = 1 - - pccd = np.array(cdelt * pc) - scales = np.sqrt((pccd ** 2).sum(axis=0, dtype=float)) - the_wcs.pscale = scales[0] diff --git a/drizzle/utils.py b/drizzle/utils.py new file mode 100644 index 00000000..a3362af3 --- /dev/null +++ b/drizzle/utils.py @@ -0,0 +1,239 @@ +import math + +import numpy as np + +__all__ = ["calc_pixmap", "decode_context", "estimate_pixel_scale_ratio"] + +_DEG2RAD = math.pi / 180.0 + + +def calc_pixmap(wcs_from, wcs_to, shape=None): + """ + Calculate a discretized on a grid mapping between the pixels of two images + using provided WCS of the original ("from") image and the destination ("to") + image. + + .. note:: + This function assumes that output frames of ``wcs_from`` and ``wcs_to`` + WCS have the same units. + + Parameters + ---------- + wcs_from : wcs + A WCS object representing the coordinate system you are + converting from. This object's ``array_shape`` (or ``pixel_shape``) + property will be used to define the shape of the pixel map array. + If ``shape`` parameter is provided, it will take precedence + over this object's ``array_shape`` value. + + wcs_to : wcs + A WCS object representing the coordinate system you are + converting to. + + shape : tuple, None, optional + A tuple of integers indicating the shape of the output array in the + ``numpy.ndarray`` order. When provided, it takes precedence over the + ``wcs_from.array_shape`` property. + + Returns + ------- + pixmap : numpy.ndarray + A three dimensional array representing the transformation between + the two. The last dimension is of length two and contains the x and + y coordinates of a pixel center, repectively. The other two coordinates + correspond to the two coordinates of the image the first WCS is from. + + Raises + ------ + ValueError + A `ValueError` is raised when output pixel map shape cannot be + determined from provided inputs. + + Notes + ----- + When ``shape`` is not provided and ``wcs_from.array_shape`` is not set + (i.e., it is `None`), `calc_pixmap` will attempt to determine pixel map + shape from the ``bounding_box`` property of the input ``wcs_from`` object. + If ``bounding_box`` is not available, a `ValueError` will be raised. + + """ + if shape is None: + shape = wcs_from.array_shape + if shape is None: + if (bbox := getattr(wcs_from, "bounding_box", None)) is not None: + if (nd := np.ndim(bbox)) == 1: + bbox = (bbox, ) + if nd > 1: + shape = tuple( + int(math.ceil(lim[1] + 0.5)) for lim in bbox[::-1] + ) + + if shape is None: + raise ValueError( + 'The "from" WCS must have pixel_shape property set.' + ) + + y, x = np.indices(shape, dtype=np.float64) + x, y = wcs_to.world_to_pixel_values(*wcs_from.pixel_to_world_values(x, y)) + pixmap = np.dstack([x, y]) + return pixmap + + +def estimate_pixel_scale_ratio(wcs_from, wcs_to, refpix_from=None, refpix_to=None): + """ + Compute the ratio of the pixel scale of the "to" WCS at the ``refpix_to`` + position to the pixel scale of the "from" WCS at the ``refpix_from`` + position. Pixel scale ratio, + when requested, is computed near the centers of the bounding box + (a property of the WCS object) or near ``refpix_*`` coordinates + if supplied. + + Pixel scale is estimated as the square root of pixel's area, i.e., + pixels are assumed to have a square shape at the reference + pixel position. If input reference pixel position for a WCS is `None`, + it will be taken as the center of the bounding box + if ``wcs_*`` has a bounding box defined, or as the center of the box + defined by the ``pixel_shape`` attribute of the input WCS if + ``pixel_shape`` is defined (not `None`), or at pixel coordinates + ``(0, 0)``. + + Parameters + ---------- + wcs_from : wcs + A WCS object representing the coordinate system you are + converting from. This object *must* have ``pixel_shape`` property + defined. + + wcs_to : wcs + A WCS object representing the coordinate system you are + converting to. + + refpix_from : numpy.ndarray, tuple, list + Image coordinates of the reference pixel near which pixel scale should + be computed in the "from" image. In FITS WCS this could be, for example, + the value of CRPIX of the ``wcs_from`` WCS. + + refpix_to : numpy.ndarray, tuple, list + Image coordinates of the reference pixel near which pixel scale should + be computed in the "to" image. In FITS WCS this could be, for example, + the value of CRPIX of the ``wcs_to`` WCS. + + Returns + ------- + pixel_scale_ratio : float + Estimate the ratio of "to" to "from" WCS pixel scales. This value is + returned only when ``estimate_pixel_scale_ratio`` is `True`. + + """ + pscale_ratio = (_estimate_pixel_scale(wcs_to, refpix_to) / + _estimate_pixel_scale(wcs_from, refpix_from)) + return pscale_ratio + + +def _estimate_pixel_scale(wcs, refpix): + # estimate pixel scale (in rad) using approximate algorithm + # from https://trs.jpl.nasa.gov/handle/2014/40409 + if refpix is None: + if hasattr(wcs, 'bounding_box') and wcs.bounding_box is not None: + refpix = np.mean(wcs.bounding_box, axis=-1) + else: + if wcs.pixel_shape: + refpix = np.array([(i - 1) // 2 for i in wcs.pixel_shape]) + else: + refpix = np.zeros(wcs.pixel_n_dim) + + else: + refpix = np.asarray(refpix) + + l1, phi1 = wcs.pixel_to_world_values(*(refpix - 0.5)) + l2, phi2 = wcs.pixel_to_world_values(*(refpix + [-0.5, 0.5])) + l3, phi3 = wcs.pixel_to_world_values(*(refpix + 0.5)) + l4, phi4 = wcs.pixel_to_world_values(*(refpix + [0.5, -0.5])) + area = _DEG2RAD * abs( + 0.5 * ( + (l4 - l2) * (math.sin(_DEG2RAD * phi1) - math.sin(_DEG2RAD * phi3)) + + (l1 - l3) * (math.sin(_DEG2RAD * phi2) - math.sin(_DEG2RAD * phi4)) + ) + ) + return math.sqrt(area) + + +def decode_context(context, x, y): + """Get 0-based indices of input images that contributed to (resampled) + output pixel with coordinates ``x`` and ``y``. + + Parameters + ---------- + context: numpy.ndarray + A 3D `~numpy.ndarray` of integral data type. + + x: int, list of integers, numpy.ndarray of integers + X-coordinate of pixels to decode (3rd index into the ``context`` array) + + y: int, list of integers, numpy.ndarray of integers + Y-coordinate of pixels to decode (2nd index into the ``context`` array) + + Returns + ------- + A list of `numpy.ndarray` objects each containing indices of input images + that have contributed to an output pixel with coordinates ``x`` and ``y``. + The length of returned list is equal to the number of input coordinate + arrays ``x`` and ``y``. + + Examples + -------- + An example context array for an output image of array shape ``(5, 6)`` + obtained by resampling 80 input images. + + >>> import numpy as np + >>> from drizzle.utils import decode_context + >>> ctx = np.array( + ... [[[0, 0, 0, 0, 0, 0], + ... [0, 0, 0, 36196864, 0, 0], + ... [0, 0, 0, 0, 0, 0], + ... [0, 0, 0, 0, 0, 0], + ... [0, 0, 537920000, 0, 0, 0]], + ... [[0, 0, 0, 0, 0, 0,], + ... [0, 0, 0, 67125536, 0, 0], + ... [0, 0, 0, 0, 0, 0], + ... [0, 0, 0, 0, 0, 0], + ... [0, 0, 163856, 0, 0, 0]], + ... [[0, 0, 0, 0, 0, 0], + ... [0, 0, 0, 8203, 0, 0], + ... [0, 0, 0, 0, 0, 0], + ... [0, 0, 0, 0, 0, 0], + ... [0, 0, 32865, 0, 0, 0]]], + ... dtype=np.int32 + ... ) + >>> decode_context(ctx, [3, 2], [1, 4]) + [array([ 9, 12, 14, 19, 21, 25, 37, 40, 46, 58, 64, 65, 67, 77]), + array([ 9, 20, 29, 36, 47, 49, 64, 69, 70, 79])] + + """ + if context.ndim != 3: + raise ValueError("'context' must be a 3D array.") + + x = np.atleast_1d(x) + y = np.atleast_1d(y) + + if x.size != y.size: + raise ValueError("Coordinate arrays must have equal length.") + + if x.ndim != 1: + raise ValueError("Coordinates must be scalars or 1D arrays.") + + if not (np.issubdtype(x.dtype, np.integer) and + np.issubdtype(y.dtype, np.integer)): + raise ValueError('Pixel coordinates must be integer values') + + nbits = 8 * context.dtype.itemsize + one = np.array(1, context.dtype) + flags = np.array([one << i for i in range(nbits)]) + + idx = [] + for xi, yi in zip(x, y): + idx.append( + np.flatnonzero(np.bitwise_and.outer(context[:, yi, xi], flags)) + ) + + return idx diff --git a/pyproject.toml b/pyproject.toml index e3bae4bf..73595aac 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,6 +30,8 @@ Documentation = "http://spacetelescope.github.io/drizzle/" [project.optional-dependencies] test = [ "pytest", + "pytest-cov", + "pytest-doctestplus", ] docs = [ "tomli; python_version<'3.11'", @@ -37,6 +39,7 @@ docs = [ "sphinx-automodapi", "sphinx-rtd-theme", "matplotlib", + "pytest-doctestplus", ] [build-system] @@ -64,10 +67,9 @@ norecursedirs = [ [tool.setuptools_scm] [tool.ruff] +lint.select = ["ALL"] line-length = 100 - -[tool.ruff.lint] -select = [ +lint.extend-select = [ "F", "W", "E", @@ -76,6 +78,127 @@ select = [ ] exclude = [ "docs", - ".eggs", ".tox", + ".eggs", + "build", +] + +lint.ignore = [ + # flake8-builtins (A) : shadowing a Python built-in. + # New ones should be avoided and is up to maintainers to enforce. + "A00", + + # flake8-annotations (ANN) + "ANN101", # No annotation for `self`. + "ANN102", # No annotation for `cls`. + + # flake8-bugbear (B) + "B008", # FunctionCallArgumentDefault + + # flake8-commas (COM) + "COM812", # TrailingCommaMissing + "COM819", # TrailingCommaProhibited + + # pydocstyle (D) + # Missing Docstrings + "D102", # Missing docstring in public method. Don't check b/c docstring inheritance. + "D105", # Missing docstring in magic method. Don't check b/c class docstring. + # Whitespace Issues + "D200", # FitsOnOneLine + # Docstring Content Issues + "D410", # BlankLineAfterSection. Using D412 instead. + "D400", # EndsInPeriod. NOTE: might want to revisit this. + + # pycodestyle (E, W) + "E711", # NoneComparison (see unfixable) + "E741", # AmbiguousVariableName. Physics variables are often poor code variables + + # flake8-fixme (FIX) + "FIX002", # Line contains TODO | notes for improvements are OK iff the code works + + # pep8-naming (N) + "N803", # invalid-argument-name. Physics variables are often poor code variables + "N806", # non-lowercase-variable-in-function. Physics variables are often poor code variables + + # pandas-vet (PD) + "PD", + + # flake8-self (SLF) + "SLF001", # private member access + + # flake8-todos (TD) + "TD002", # Missing author in TODO + + # Ruff-specific rules (RUF) + "RUF005", # unpack-instead-of-concatenating-to-collection-literal -- it's not clearly faster. + + "Q000", # Single quotes found but double quotes preferred +] + +[tool.ruff.lint.extend-per-file-ignores] +"setup.py" = [ + "INP001", # Part of configuration, not a package. + "ICN001", + "UP018" +] +"test_*.py" = [ + "B018", # UselessExpression + "D", # pydocstyle + "PGH001", # No builtin eval() allowed + "S101", # Use of assert detected +] +".pyinstaller/*.py" = ["INP001"] # Not a package. +"conftest.py" = ["INP001"] # Part of configuration, not a package. +"docs/*.py" = [ + "INP001", # implicit-namespace-package. The examples are not a package. +] +"examples/*.py" = [ + "E402", # Imports are done as needed. + "INP001", # implicit-namespace-package. The examples are not a package. + "T203" # pprint found +] + +[tool.ruff.lint.flake8-annotations] +ignore-fully-untyped = true +mypy-init-return = true + +[tool.ruff.lint.flake8-comprehensions] +allow-dict-calls-with-keyword-arguments = true + +[tool.ruff.lint.flake8-type-checking] +exempt-modules = [] + +[tool.ruff.lint.isort] +known-first-party = ["astropy", "extension_helpers"] + +[tool.ruff.lint.pydocstyle] +convention = "numpy" + +[tool.coverage.run] +omit = [ + "drizzle/__init__*", + "drizzle/**/conftest.py", + "drizzle/**/setup*", + "drizzle/**/tests/*", + "drizzle/version*", + "drizzle/_version.py", + "*/drizzle/__init__*", + "*/drizzle/**/conftest.py", + "*/drizzle/**/setup*", + "*/drizzle/**/tests/*", +] + +[tool.coverage.report] +exclude_lines = [ + # Have to re-enable the standard pragma + "pragma: no cover", + # Don't complain about packages we have installed + "except ImportError", + # Don't complain if tests don't hit defensive assertion code: + "raise AssertionError", + "raise NotImplementedError", + # Don't complain about script hooks + "'def main(.*):'", + # Ignore branches that don't pertain to this version of Python + "pragma: py{ignore_python_version}", ] diff --git a/setup.py b/setup.py index 11580984..73cd4cae 100755 --- a/setup.py +++ b/setup.py @@ -3,9 +3,9 @@ import os import sys -from setuptools import setup, Extension import numpy +from setuptools import Extension, setup def get_extensions(): diff --git a/src/cdrizzleapi.c b/src/cdrizzleapi.c index b16e7fbe..0fcacb02 100644 --- a/src/cdrizzleapi.c +++ b/src/cdrizzleapi.c @@ -84,7 +84,7 @@ tdriz(PyObject *obj UNUSED_PARAM, PyObject *args, PyObject *keywords) struct driz_error_t error; struct driz_param_t p; integer_t isize[2], psize[2], wsize[2]; - char warn_msg[96]; + char warn_msg[128]; driz_log_handle = driz_log_init(driz_log_handle); driz_log_message("starting tdriz"); @@ -130,7 +130,9 @@ tdriz(PyObject *obj UNUSED_PARAM, PyObject *args, PyObject *keywords) goto _exit; } - if (ocon != Py_None) { + if (ocon == Py_None) { + con = NULL; + } else { con = (PyArrayObject *)PyArray_ContiguousFromAny(ocon, NPY_INT32, 2, 2); if (!con) { driz_error_set_message(&error, "Invalid context array"); @@ -138,7 +140,7 @@ tdriz(PyObject *obj UNUSED_PARAM, PyObject *args, PyObject *keywords) } }; - /* Convert t`he fill value string */ + /* Convert the fill value string */ if (fillstr == NULL || *fillstr == 0 || @@ -185,15 +187,8 @@ tdriz(PyObject *obj UNUSED_PARAM, PyObject *args, PyObject *keywords) goto _exit; } - if (kernel == kernel_tophat) { - if (sprintf(warn_msg, - "Kernel '%s' has been deprecated and it will be removed in a future release.", - kernel_str) < 1) { - strcpy(warn_msg, "Selected kernel has been deprecated and it will be removed in a future release."); - } - PyErr_WarnEx(PyExc_DeprecationWarning, warn_msg, 1); - } else if (kernel == kernel_gaussian || kernel == kernel_lanczos2 || kernel == kernel_lanczos3) { - if (sprintf(warn_msg, + if (kernel == kernel_gaussian || kernel == kernel_lanczos2 || kernel == kernel_lanczos3) { + if (snprintf(warn_msg, 128, "Kernel '%s' is not a flux-conserving kernel.", kernel_str) < 1) { strcpy(warn_msg, "Selected kernel is not a flux-conserving kernel."); @@ -208,7 +203,7 @@ tdriz(PyObject *obj UNUSED_PARAM, PyObject *args, PyObject *keywords) /* If the input image is not in CPS we need to divide by the exposure */ if (inun != unit_cps) { - inv_exposure_time = 1.0f / p.exposure_time; + inv_exposure_time = 1.0f / expin; scale_image(img, inv_exposure_time); } @@ -245,14 +240,24 @@ tdriz(PyObject *obj UNUSED_PARAM, PyObject *args, PyObject *keywords) get_dimensions(p.pixmap, psize); if (psize[0] != isize[0] || psize[1] != isize[1]) { - driz_error_set_message(&error, "Pixel map dimensions != input dimensions"); + if (snprintf(warn_msg, 128, + "Pixel map dimensions (%d, %d) != input dimensions (%d, %d).", + psize[0], psize[1], isize[0], isize[1]) < 1) { + strcpy(warn_msg, "Pixel map dimensions != input dimensions."); + } + driz_error_set_message(&error, warn_msg); goto _exit; } if (p.weights) { get_dimensions(p.weights, wsize); if (wsize[0] != isize[0] || wsize[1] != isize[1]) { - driz_error_set_message(&error, "Weights array dimensions != input dimensions"); + if (snprintf(warn_msg, 128, + "Weights array dimensions (%d, %d) != input dimensions (%d, %d).", + wsize[0], wsize[1], isize[0], isize[1]) < 1) { + strcpy(warn_msg, "Weights array dimensions != input dimensions."); + } + driz_error_set_message(&error, warn_msg); goto _exit; } } @@ -315,6 +320,7 @@ tblot(PyObject *obj, PyObject *args, PyObject *keywords) struct driz_error_t error; struct driz_param_t p; integer_t psize[2], osize[2]; + char warn_msg[128]; driz_log_handle = driz_log_init(driz_log_handle); driz_log_message("starting tblot"); @@ -355,7 +361,12 @@ tblot(PyObject *obj, PyObject *args, PyObject *keywords) get_dimensions(out, osize); if (psize[0] != osize[0] || psize[1] != osize[1]) { - driz_error_set_message(&error, "Pixel map dimensions != output dimensions"); + if (snprintf(warn_msg, 128, + "Pixel map dimensions (%d, %d) != output dimensions (%d, %d).", + psize[0], psize[1], osize[0], osize[1]) < 1) { + strcpy(warn_msg, "Pixel map dimensions != output dimensions."); + } + driz_error_set_message(&error, warn_msg); goto _exit; } @@ -470,6 +481,7 @@ invert_pixmap_wrap(PyObject *self, PyObject *args) struct driz_param_t par; double *xy, *xyin; npy_intp *ndim, xyin_dim = 2; + const double half = 0.5 - DBL_EPSILON; xyin = (double *) malloc(2 * sizeof(double)); @@ -491,19 +503,19 @@ invert_pixmap_wrap(PyObject *self, PyObject *args) ndim = PyArray_DIMS(pixmap_arr); if (bbox == Py_None) { - par.xmin = -0.5; - par.xmax = ndim[1] - 0.5; - par.ymin = -0.5; - par.ymax = ndim[0] - 0.5; + par.xmin = 0; + par.xmax = ndim[1] - 1; + par.ymin = 0; + par.ymax = ndim[0] - 1; } else { bbox_arr = (PyArrayObject *)PyArray_ContiguousFromAny(bbox, NPY_DOUBLE, 2, 2); if (!bbox_arr) { return PyErr_Format(gl_Error, "Invalid input bounding box."); } - par.xmin = *(double*) PyArray_GETPTR2(bbox_arr, 0, 0); - par.xmax = *(double*) PyArray_GETPTR2(bbox_arr, 0, 1); - par.ymin = *(double*) PyArray_GETPTR2(bbox_arr, 1, 0); - par.ymax = *(double*) PyArray_GETPTR2(bbox_arr, 1, 1); + par.xmin = (integer_t)(*(double*) PyArray_GETPTR2(bbox_arr, 0, 0) - half); + par.xmax = (integer_t)(*(double*) PyArray_GETPTR2(bbox_arr, 0, 1) + half); + par.ymin = (integer_t)(*(double*) PyArray_GETPTR2(bbox_arr, 1, 0) - half); + par.ymax = (integer_t)(*(double*) PyArray_GETPTR2(bbox_arr, 1, 1) + half); } xy = (double *)PyArray_DATA(xyout_arr); @@ -577,9 +589,9 @@ clip_polygon_wrap(PyObject *self, PyObject *args) static struct PyMethodDef cdrizzle_methods[] = { {"tdriz", (PyCFunction)tdriz, METH_VARARGS|METH_KEYWORDS, - "tdriz(image, weight, output, outweight, context, uniqid, xmin, ymin, scale, pfract, kernel, inun, expin, wtscl, fill, nmiss, nskip, pixmap)"}, + "tdriz(image, weights, pixmap, output, counts, context, uniqid, xmin, xmax, ymin, ymax, scale, pixfrac, kernel, in_units, expscale, wtscale, fillstr)"}, {"tblot", (PyCFunction)tblot, METH_VARARGS|METH_KEYWORDS, - "tblot(image, output, xmin, xmax, ymin, ymax, scale, kscale, interp, ef, misval, sinscl, pixmap)"}, + "tblot(image, pixmap, output, xmin, xmax, ymin, ymax, scale, kscale, interp, exptime, misval, sinscl)"}, {"test_cdrizzle", test_cdrizzle, METH_VARARGS, "test_cdrizzle(data, weights, pixmap, output_data, output_counts)"}, {"invert_pixmap", invert_pixmap_wrap, METH_VARARGS, "invert_pixmap(pixmap, xyout, bbox)"}, diff --git a/src/cdrizzlebox.c b/src/cdrizzlebox.c index a77b5468..f6b401d1 100644 --- a/src/cdrizzlebox.c +++ b/src/cdrizzlebox.c @@ -12,9 +12,6 @@ #include -static char buf[1024]; - - /** -------------------------------------------------------------------------------------------------- * Update the flux and counts in the output image using a weighted average * @@ -251,7 +248,7 @@ static int do_kernel_point(struct driz_param_t* p) { struct scanner s; integer_t i, j, ii, jj; - integer_t ybounds[2], osize[2]; + integer_t osize[2]; float scale2, vc, d, dow; integer_t bv; int xmin, xmax, ymin, ymax, n; @@ -330,123 +327,6 @@ do_kernel_point(struct driz_param_t* p) { return 0; } -/** -------------------------------------------------------------------------------------------------- - * This kernel assumes flux is distrubuted evenly across a circle around the center of a pixel - * - * p: structure containing options, input, and output - */ - -static int -do_kernel_tophat(struct driz_param_t* p) { - struct scanner s; - integer_t bv, i, j, ii, jj, nhit, nxi, nxa, nyi, nya; - integer_t ybounds[2], osize[2]; - float scale2, pfo, pfo2, vc, d, dow; - double xxi, xxa, yyi, yya, ddx, ddy, r2; - int xmin, xmax, ymin, ymax, n; - - scale2 = p->scale * p->scale; - pfo = p->pixel_fraction / p->scale / 2.0; - pfo2 = pfo * pfo; - bv = compute_bit_value(p->uuid); - - if (init_image_scanner(p, &s, &ymin, &ymax)) return 1; - - p->nskip = (p->ymax - p->ymin) - (ymax - ymin); - p->nmiss = p->nskip * (p->xmax - p->xmin); - - /* This is the outer loop over all the lines in the input image */ - - get_dimensions(p->output_data, osize); - for (j = ymin; j <= ymax; ++j) { - /* Check the overlap with the output */ - n = get_scanline_limits(&s, j, &xmin, &xmax); - if (n == 1) { - // scan ended (y reached the top vertex/edge) - p->nskip += (ymax + 1 - j); - p->nmiss += (ymax + 1 - j) * (p->xmax - p->xmin); - break; - } else if (n == 2 || n == 3) { - // pixel centered on y is outside of scanner's limits or image [0, height - 1] - // OR: limits (x1, x2) are equal (line width is 0) - p->nmiss += (p->xmax - p->xmin); - ++p->nskip; - continue; - } else { - // limits (x1, x2) are equal (line width is 0) - p->nmiss += (p->xmax - p->xmin) - (xmax + 1 - xmin); - } - - for (i = xmin; i <= xmax; ++i) { - double ox, oy; - - if (map_pixel(p->pixmap, i, j, &ox, &oy)) { - nhit = 0; - - } else { - /* Offset within the subset */ - xxi = ox - pfo; - xxa = ox + pfo; - yyi = oy - pfo; - yya = oy + pfo; - - nxi = MAX(fortran_round(xxi), 0); - nxa = MIN(fortran_round(xxa), osize[0]-1); - nyi = MAX(fortran_round(yyi), 0); - nya = MIN(fortran_round(yya), osize[1]-1); - - nhit = 0; - - /* Allow for stretching because of scale change */ - d = get_pixel(p->data, i, j) * scale2; - - /* Scale the weighting mask by the scale factor and inversely by - the Jacobian to ensure conservation of weight in the output */ - if (p->weights) { - dow = get_pixel(p->weights, i, j) * p->weight_scale; - } else { - dow = 1.0; - } - - /* Loop over output pixels which could be affected */ - for (jj = nyi; jj <= nya; ++jj) { - ddy = oy - (double)jj; - - /* Check it is on the output image */ - for (ii = nxi; ii <= nxa; ++ii) { - ddx = ox - (double)ii; - - /* Radial distance */ - r2 = ddx*ddx + ddy*ddy; - - /* Weight is one within the specified radius and zero outside. - Note: weight isn't conserved in this case */ - if (r2 <= pfo2) { - /* Count the hits */ - nhit++; - vc = get_pixel(p->output_counts, ii, jj); - - /* If we are create or modifying the context image, - we do so here. */ - if (p->output_context && dow > 0.0) { - set_bit(p->output_context, ii, jj, bv); - } - - if (update_data(p, ii, jj, d, vc, dow)) { - return 1; - } - } - } - } - } - - /* Count cases where the pixel is off the output image */ - if (nhit == 0) ++ p->nmiss; - } - } - - return 0; -} /** -------------------------------------------------------------------------------------------------- * This kernel assumes the flux is distributed acrass a gaussian around the center of an input pixel @@ -458,7 +338,7 @@ static int do_kernel_gaussian(struct driz_param_t* p) { struct scanner s; integer_t bv, i, j, ii, jj, nxi, nxa, nyi, nya, nhit; - integer_t ybounds[2], osize[2]; + integer_t osize[2]; float vc, d, dow; double gaussian_efac, gaussian_es; double pfo, ac, scale2, xxi, xxa, yyi, yya, w, ddx, ddy, r2, dover; @@ -586,7 +466,7 @@ static int do_kernel_lanczos(struct driz_param_t* p) { struct scanner s; integer_t bv, i, j, ii, jj, nxi, nxa, nyi, nya, nhit, ix, iy; - integer_t ybounds[2], osize[2]; + integer_t osize[2]; float scale2, vc, d, dow; double pfo, xx, yy, xxi, xxa, yyi, yya, w, dx, dy, dover; int kernel_order; @@ -987,7 +867,6 @@ kernel_handler_map[] = { do_kernel_square, do_kernel_gaussian, do_kernel_point, - do_kernel_tophat, do_kernel_turbo, do_kernel_lanczos, do_kernel_lanczos diff --git a/src/cdrizzlemap.c b/src/cdrizzlemap.c index bbc7189f..f3c8e64b 100644 --- a/src/cdrizzlemap.c +++ b/src/cdrizzlemap.c @@ -105,7 +105,6 @@ shrink_image_section(PyArrayObject *pixmap, int *xmin, int *xmax, int *ymin, int interpolate_point(struct driz_param_t *par, double xin, double yin, double *xout, double *yout) { - int ipix, jpix, npix, idim; int i0, j0, nx2, ny2; npy_intp *ndim; double x, y, x1, y1, f00, f01, f10, f11, g00, g01, g10, g11; @@ -207,6 +206,7 @@ map_point(struct driz_param_t *par, double xin, double yin, double *xout, if (i >= par->xmin && i <= par->xmax && j >= par->ymin && j <= par->ymax) { status = map_pixel(par->pixmap, i, j, xout, yout); + return status; } else { return 1; } @@ -492,14 +492,14 @@ clip_polygon_to_window(const struct polygon *p, const struct polygon *wnd, int v1_inside, v2_inside; struct polygon p1, p2, *ppin, *ppout, *tpp; struct vertex *pv, *pv_, *wv, *wv_, dp, dw, vi; - double t, u, d, signed_area, app_, aww_; + double d, app_, aww_; if ((p->npv < 3) || (wnd->npv < 3)) { return 1; } - orient_ccw(p); - orient_ccw(wnd); + orient_ccw((struct polygon *)p); + orient_ccw((struct polygon *)wnd); p1 = *p; @@ -906,9 +906,6 @@ int init_image_scanner(struct driz_param_t *par, struct scanner *s, int *ymin, int *ymax) { struct polygon p, q, pq, inpq; - double xyin[2], xyout[2]; - integer_t isize[2], osize[2]; - int ipoint; int k, n; npy_intp *ndim; diff --git a/src/cdrizzleutil.c b/src/cdrizzleutil.c index 88d0ef56..4ea1f3cc 100644 --- a/src/cdrizzleutil.c +++ b/src/cdrizzleutil.c @@ -150,7 +150,6 @@ static const char* kernel_string_table[] = { "square", "gaussian", "point", - "tophat", "turbo", "lanczos2", "lanczos3", diff --git a/src/cdrizzleutil.h b/src/cdrizzleutil.h index 64a71b76..d230bc5f 100644 --- a/src/cdrizzleutil.h +++ b/src/cdrizzleutil.h @@ -81,7 +81,6 @@ enum e_kernel_t { kernel_square, kernel_gaussian, kernel_point, - kernel_tophat, kernel_turbo, kernel_lanczos2, kernel_lanczos3, diff --git a/src/tests/utest_cdrizzle.c b/src/tests/utest_cdrizzle.c index d4449f19..b32e2fd9 100644 --- a/src/tests/utest_cdrizzle.c +++ b/src/tests/utest_cdrizzle.c @@ -615,7 +615,6 @@ FCT_BGN_FN(utest_cdrizzle) /* Test for half overlap */ - integer_t ybounds[2]; /* start of in-bounds */ struct driz_param_t *p; /* parameter structure */ p = setup_parameters(); @@ -639,7 +638,6 @@ FCT_BGN_FN(utest_cdrizzle) /* Test for negative half overlap */ - integer_t ybounds[2]; /* start of in-bounds */ struct driz_param_t *p; /* parameter structure */ p = setup_parameters();