diff options
| -rw-r--r-- | doc/Makefile | 153 | ||||
| -rw-r--r-- | doc/make.bat | 190 | ||||
| -rw-r--r-- | doc/source/api.rst | 30 | ||||
| -rw-r--r-- | doc/source/conf.py | 246 | ||||
| -rw-r--r-- | doc/source/cutlines.png | bin | 0 -> 54023 bytes | |||
| -rw-r--r-- | doc/source/index.rst | 23 | ||||
| -rw-r--r-- | doc/source/inside.png | bin | 0 -> 86567 bytes | |||
| -rw-r--r-- | doc/source/user.rst | 160 | ||||
| -rw-r--r-- | lib/__init__.py | 0 | ||||
| -rw-r--r-- | lib/graph.py | 792 | ||||
| -rw-r--r-- | lib/great_circle_arc.py | 318 | ||||
| -rw-r--r-- | lib/polygon.py | 679 | ||||
| -rw-r--r-- | lib/skyline.py | 48 | ||||
| -rw-r--r-- | lib/test/__init__.py | 1 | ||||
| -rw-r--r-- | lib/test/benchmarks.py | 39 | ||||
| -rw-r--r-- | lib/test/data/1904-66_TAN.fits | 375 | ||||
| -rw-r--r-- | lib/test/data/2chipA.fits.REMOVED.git-id | 1 | ||||
| -rw-r--r-- | lib/test/data/2chipB.fits.REMOVED.git-id | 1 | ||||
| -rw-r--r-- | lib/test/data/2chipC.fits.REMOVED.git-id | 1 | ||||
| -rw-r--r-- | lib/test/test.py | 260 | ||||
| -rw-r--r-- | lib/test/test_intersection.py | 167 | ||||
| -rw-r--r-- | lib/test/test_union.py | 169 | ||||
| -rw-r--r-- | lib/test/test_util.py | 14 | ||||
| -rw-r--r-- | lib/vector.py | 213 | ||||
| -rwxr-xr-x | setup.py | 45 | ||||
| -rw-r--r-- | src/math_util.c | 374 |
26 files changed, 4299 insertions, 0 deletions
diff --git a/doc/Makefile b/doc/Makefile new file mode 100644 index 0000000..6dd7fc4 --- /dev/null +++ b/doc/Makefile @@ -0,0 +1,153 @@ +# Makefile for Sphinx documentation +# + +# You can set these variables from the command line. +SPHINXOPTS = +SPHINXBUILD = sphinx-build +PAPER = +BUILDDIR = build + +# Internal variables. +PAPEROPT_a4 = -D latex_paper_size=a4 +PAPEROPT_letter = -D latex_paper_size=letter +ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source +# the i18n builder cannot share the environment and doctrees with the others +I18NSPHINXOPTS = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source + +.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest gettext + +help: + @echo "Please use \`make <target>' where <target> is one of" + @echo " html to make standalone HTML files" + @echo " dirhtml to make HTML files named index.html in directories" + @echo " singlehtml to make a single large HTML file" + @echo " pickle to make pickle files" + @echo " json to make JSON files" + @echo " htmlhelp to make HTML files and a HTML help project" + @echo " qthelp to make HTML files and a qthelp project" + @echo " devhelp to make HTML files and a Devhelp project" + @echo " epub to make an epub" + @echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter" + @echo " latexpdf to make LaTeX files and run them through pdflatex" + @echo " text to make text files" + @echo " man to make manual pages" + @echo " texinfo to make Texinfo files" + @echo " info to make Texinfo files and run them through makeinfo" + @echo " gettext to make PO message catalogs" + @echo " changes to make an overview of all changed/added/deprecated items" + @echo " linkcheck to check all external links for integrity" + @echo " doctest to run all doctests embedded in the documentation (if enabled)" + +clean: + -rm -rf $(BUILDDIR)/* + +html: + $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html + @echo + @echo "Build finished. The HTML pages are in $(BUILDDIR)/html." + +dirhtml: + $(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml + @echo + @echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml." + +singlehtml: + $(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml + @echo + @echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml." + +pickle: + $(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle + @echo + @echo "Build finished; now you can process the pickle files." + +json: + $(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json + @echo + @echo "Build finished; now you can process the JSON files." + +htmlhelp: + $(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp + @echo + @echo "Build finished; now you can run HTML Help Workshop with the" \ + ".hhp project file in $(BUILDDIR)/htmlhelp." + +qthelp: + $(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp + @echo + @echo "Build finished; now you can run "qcollectiongenerator" with the" \ + ".qhcp project file in $(BUILDDIR)/qthelp, like this:" + @echo "# qcollectiongenerator $(BUILDDIR)/qthelp/SphericalGeometryToolkit.qhcp" + @echo "To view the help file:" + @echo "# assistant -collectionFile $(BUILDDIR)/qthelp/SphericalGeometryToolkit.qhc" + +devhelp: + $(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp + @echo + @echo "Build finished." + @echo "To view the help file:" + @echo "# mkdir -p $$HOME/.local/share/devhelp/SphericalGeometryToolkit" + @echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/SphericalGeometryToolkit" + @echo "# devhelp" + +epub: + $(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub + @echo + @echo "Build finished. The epub file is in $(BUILDDIR)/epub." + +latex: + $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + @echo + @echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex." + @echo "Run \`make' in that directory to run these through (pdf)latex" \ + "(use \`make latexpdf' here to do that automatically)." + +latexpdf: + $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + @echo "Running LaTeX files through pdflatex..." + make -C $(BUILDDIR)/latex all-pdf + @echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex." + +text: + $(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text + @echo + @echo "Build finished. The text files are in $(BUILDDIR)/text." + +man: + $(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man + @echo + @echo "Build finished. The manual pages are in $(BUILDDIR)/man." + +texinfo: + $(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo + @echo + @echo "Build finished. The Texinfo files are in $(BUILDDIR)/texinfo." + @echo "Run \`make' in that directory to run these through makeinfo" \ + "(use \`make info' here to do that automatically)." + +info: + $(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo + @echo "Running Texinfo files through makeinfo..." + make -C $(BUILDDIR)/texinfo info + @echo "makeinfo finished; the Info files are in $(BUILDDIR)/texinfo." + +gettext: + $(SPHINXBUILD) -b gettext $(I18NSPHINXOPTS) $(BUILDDIR)/locale + @echo + @echo "Build finished. The message catalogs are in $(BUILDDIR)/locale." + +changes: + $(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes + @echo + @echo "The overview file is in $(BUILDDIR)/changes." + +linkcheck: + $(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck + @echo + @echo "Link check complete; look for any errors in the above output " \ + "or in $(BUILDDIR)/linkcheck/output.txt." + +doctest: + $(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest + @echo "Testing of doctests in the sources finished, look at the " \ + "results in $(BUILDDIR)/doctest/output.txt." diff --git a/doc/make.bat b/doc/make.bat new file mode 100644 index 0000000..d9ed58d --- /dev/null +++ b/doc/make.bat @@ -0,0 +1,190 @@ +@ECHO OFF + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set BUILDDIR=build +set ALLSPHINXOPTS=-d %BUILDDIR%/doctrees %SPHINXOPTS% source +set I18NSPHINXOPTS=%SPHINXOPTS% source +if NOT "%PAPER%" == "" ( + set ALLSPHINXOPTS=-D latex_paper_size=%PAPER% %ALLSPHINXOPTS% + set I18NSPHINXOPTS=-D latex_paper_size=%PAPER% %I18NSPHINXOPTS% +) + +if "%1" == "" goto help + +if "%1" == "help" ( + :help + echo.Please use `make ^<target^>` where ^<target^> is one of + echo. html to make standalone HTML files + echo. dirhtml to make HTML files named index.html in directories + echo. singlehtml to make a single large HTML file + echo. pickle to make pickle files + echo. json to make JSON files + echo. htmlhelp to make HTML files and a HTML help project + echo. qthelp to make HTML files and a qthelp project + echo. devhelp to make HTML files and a Devhelp project + echo. epub to make an epub + echo. latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter + echo. text to make text files + echo. man to make manual pages + echo. texinfo to make Texinfo files + echo. gettext to make PO message catalogs + echo. changes to make an overview over all changed/added/deprecated items + echo. linkcheck to check all external links for integrity + echo. doctest to run all doctests embedded in the documentation if enabled + goto end +) + +if "%1" == "clean" ( + for /d %%i in (%BUILDDIR%\*) do rmdir /q /s %%i + del /q /s %BUILDDIR%\* + goto end +) + +if "%1" == "html" ( + %SPHINXBUILD% -b html %ALLSPHINXOPTS% %BUILDDIR%/html + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The HTML pages are in %BUILDDIR%/html. + goto end +) + +if "%1" == "dirhtml" ( + %SPHINXBUILD% -b dirhtml %ALLSPHINXOPTS% %BUILDDIR%/dirhtml + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The HTML pages are in %BUILDDIR%/dirhtml. + goto end +) + +if "%1" == "singlehtml" ( + %SPHINXBUILD% -b singlehtml %ALLSPHINXOPTS% %BUILDDIR%/singlehtml + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The HTML pages are in %BUILDDIR%/singlehtml. + goto end +) + +if "%1" == "pickle" ( + %SPHINXBUILD% -b pickle %ALLSPHINXOPTS% %BUILDDIR%/pickle + if errorlevel 1 exit /b 1 + echo. + echo.Build finished; now you can process the pickle files. + goto end +) + +if "%1" == "json" ( + %SPHINXBUILD% -b json %ALLSPHINXOPTS% %BUILDDIR%/json + if errorlevel 1 exit /b 1 + echo. + echo.Build finished; now you can process the JSON files. + goto end +) + +if "%1" == "htmlhelp" ( + %SPHINXBUILD% -b htmlhelp %ALLSPHINXOPTS% %BUILDDIR%/htmlhelp + if errorlevel 1 exit /b 1 + echo. + echo.Build finished; now you can run HTML Help Workshop with the ^ +.hhp project file in %BUILDDIR%/htmlhelp. + goto end +) + +if "%1" == "qthelp" ( + %SPHINXBUILD% -b qthelp %ALLSPHINXOPTS% %BUILDDIR%/qthelp + if errorlevel 1 exit /b 1 + echo. + echo.Build finished; now you can run "qcollectiongenerator" with the ^ +.qhcp project file in %BUILDDIR%/qthelp, like this: + echo.^> qcollectiongenerator %BUILDDIR%\qthelp\SphericalGeometryToolkit.qhcp + echo.To view the help file: + echo.^> assistant -collectionFile %BUILDDIR%\qthelp\SphericalGeometryToolkit.ghc + goto end +) + +if "%1" == "devhelp" ( + %SPHINXBUILD% -b devhelp %ALLSPHINXOPTS% %BUILDDIR%/devhelp + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. + goto end +) + +if "%1" == "epub" ( + %SPHINXBUILD% -b epub %ALLSPHINXOPTS% %BUILDDIR%/epub + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The epub file is in %BUILDDIR%/epub. + goto end +) + +if "%1" == "latex" ( + %SPHINXBUILD% -b latex %ALLSPHINXOPTS% %BUILDDIR%/latex + if errorlevel 1 exit /b 1 + echo. + echo.Build finished; the LaTeX files are in %BUILDDIR%/latex. + goto end +) + +if "%1" == "text" ( + %SPHINXBUILD% -b text %ALLSPHINXOPTS% %BUILDDIR%/text + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The text files are in %BUILDDIR%/text. + goto end +) + +if "%1" == "man" ( + %SPHINXBUILD% -b man %ALLSPHINXOPTS% %BUILDDIR%/man + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The manual pages are in %BUILDDIR%/man. + goto end +) + +if "%1" == "texinfo" ( + %SPHINXBUILD% -b texinfo %ALLSPHINXOPTS% %BUILDDIR%/texinfo + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The Texinfo files are in %BUILDDIR%/texinfo. + goto end +) + +if "%1" == "gettext" ( + %SPHINXBUILD% -b gettext %I18NSPHINXOPTS% %BUILDDIR%/locale + if errorlevel 1 exit /b 1 + echo. + echo.Build finished. The message catalogs are in %BUILDDIR%/locale. + goto end +) + +if "%1" == "changes" ( + %SPHINXBUILD% -b changes %ALLSPHINXOPTS% %BUILDDIR%/changes + if errorlevel 1 exit /b 1 + echo. + echo.The overview file is in %BUILDDIR%/changes. + goto end +) + +if "%1" == "linkcheck" ( + %SPHINXBUILD% -b linkcheck %ALLSPHINXOPTS% %BUILDDIR%/linkcheck + if errorlevel 1 exit /b 1 + echo. + echo.Link check complete; look for any errors in the above output ^ +or in %BUILDDIR%/linkcheck/output.txt. + goto end +) + +if "%1" == "doctest" ( + %SPHINXBUILD% -b doctest %ALLSPHINXOPTS% %BUILDDIR%/doctest + if errorlevel 1 exit /b 1 + echo. + echo.Testing of doctests in the sources finished, look at the ^ +results in %BUILDDIR%/doctest/output.txt. + goto end +) + +:end diff --git a/doc/source/api.rst b/doc/source/api.rst new file mode 100644 index 0000000..5cce592 --- /dev/null +++ b/doc/source/api.rst @@ -0,0 +1,30 @@ +API Documentation +================= + +Vectors +------- + +.. automodule:: vector + :members: + +Great circle arcs +----------------- + +.. automodule:: great_circle_arc + :members: + +Spherical polygons +------------------ + +.. currentmodule:: polygon + +.. automodule:: polygon + :members: + +Graph operations on polygons +---------------------------- + +.. currentmodule:: graph + +.. automodule:: graph + :members: diff --git a/doc/source/conf.py b/doc/source/conf.py new file mode 100644 index 0000000..5ee1814 --- /dev/null +++ b/doc/source/conf.py @@ -0,0 +1,246 @@ +# -*- coding: utf-8 -*- +# +# Spherical Geometry Toolkit documentation build configuration file, created by +# sphinx-quickstart on Tue Oct 11 09:04:48 2011. +# +# This file is execfile()d with the current directory set to its containing dir. +# +# Note that not all possible configuration values are present in this +# autogenerated file. +# +# All configuration values have a default; values that are commented out +# serve to show the default. + +from stsci.sphinxext.conf import * + +import sys, os + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +#sys.path.insert(0, os.path.abspath('.')) + +# -- General configuration ----------------------------------------------------- + +# If your documentation needs a minimal Sphinx version, state it here. +#needs_sphinx = '1.0' + +# Add any Sphinx extension module names here, as strings. They can be extensions +# coming with Sphinx (named 'sphinx.ext.*') or your custom ones. +#extensions += ['sphinx.ext.autodoc', 'sphinx.ext.intersphinx', 'sphinx.ext.mathjax', 'sphinx.ext.viewcode'] +extensions += ['sphinx.ext.autodoc', 'sphinx.ext.intersphinx', 'sphinx.ext.pngmath', 'sphinx.ext.viewcode'] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# The suffix of source filenames. +source_suffix = '.rst' + +# The encoding of source files. +#source_encoding = 'utf-8-sig' + +# The master toctree document. +master_doc = 'index' + +# General information about the project. +project = u'Spherical Geometry Toolkit' +copyright = u'2011, Michael Droettboom, STScI' + +# The version info for the project you're documenting, acts as replacement for +# |version| and |release|, also used in various other places throughout the +# built documents. +# +# The short X.Y version. +version = '0.1' +# The full version, including alpha/beta/rc tags. +release = '0.1' + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +#language = None + +# There are two options for replacing |today|: either, you set today to some +# non-false value, then it is used: +#today = '' +# Else, today_fmt is used as the format for a strftime call. +#today_fmt = '%B %d, %Y' + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +exclude_patterns = [] + +# The reST default role (used for this markup: `text`) to use for all documents. +#default_role = None + +# If true, '()' will be appended to :func: etc. cross-reference text. +#add_function_parentheses = True + +# If true, the current module name will be prepended to all description +# unit titles (such as .. function::). +#add_module_names = True + +# If true, sectionauthor and moduleauthor directives will be shown in the +# output. They are ignored by default. +#show_authors = False + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = 'sphinx' + +# A list of ignored prefixes for module index sorting. +#modindex_common_prefix = [] + + +# -- Options for HTML output --------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +#html_theme = 'default' + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +#html_theme_options = {} + +# Add any paths that contain custom themes here, relative to this directory. +#html_theme_path = [] + +# The name for this set of Sphinx documents. If None, it defaults to +# "<project> v<release> documentation". +#html_title = None + +# A shorter title for the navigation bar. Default is the same as html_title. +#html_short_title = None + +# The name of an image file (relative to this directory) to place at the top +# of the sidebar. +#html_logo = None + +# The name of an image file (within the static path) to use as favicon of the +# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 +# pixels large. +#html_favicon = None + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] + +# If not '', a 'Last updated on:' timestamp is inserted at every page bottom, +# using the given strftime format. +#html_last_updated_fmt = '%b %d, %Y' + +# If true, SmartyPants will be used to convert quotes and dashes to +# typographically correct entities. +#html_use_smartypants = True + +# Custom sidebar templates, maps document names to template names. +#html_sidebars = {} + +# Additional templates that should be rendered to pages, maps page names to +# template names. +#html_additional_pages = {} + +# If false, no module index is generated. +#html_domain_indices = True + +# If false, no index is generated. +#html_use_index = True + +# If true, the index is split into individual pages for each letter. +#html_split_index = False + +# If true, links to the reST sources are added to the pages. +#html_show_sourcelink = True + +# If true, "Created using Sphinx" is shown in the HTML footer. Default is True. +#html_show_sphinx = True + +# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True. +#html_show_copyright = True + +# If true, an OpenSearch description file will be output, and all pages will +# contain a <link> tag referring to it. The value of this option must be the +# base URL from which the finished HTML is served. +#html_use_opensearch = '' + +# This is the file name suffix for HTML files (e.g. ".xhtml"). +#html_file_suffix = None + +# Output file base name for HTML help builder. +htmlhelp_basename = 'SphericalGeometryToolkitdoc' + + +# -- Options for LaTeX output -------------------------------------------------- + +# The paper size ('letter' or 'a4'). +#latex_paper_size = 'letter' + +# The font size ('10pt', '11pt' or '12pt'). +#latex_font_size = '10pt' + +# Grouping the document tree into LaTeX files. List of tuples +# (source start file, target name, title, author, documentclass [howto/manual]). +latex_documents = [ + ('index', 'SphericalGeometryToolkit.tex', u'Spherical Geometry Toolkit Documentation', + u'Michael Droettboom, STScI', 'manual'), +] + +# The name of an image file (relative to this directory) to place at the top of +# the title page. +#latex_logo = None + +# For "manual" documents, if this is true, then toplevel headings are parts, +# not chapters. +#latex_use_parts = False + +# If true, show page references after internal links. +#latex_show_pagerefs = False + +# If true, show URL addresses after external links. +#latex_show_urls = False + +# Additional stuff for the LaTeX preamble. +#latex_preamble = '' + +# Documents to append as an appendix to all manuals. +#latex_appendices = [] + +# If false, no module index is generated. +#latex_domain_indices = True + + +# -- Options for manual page output -------------------------------------------- + +# One entry per manual page. List of tuples +# (source start file, name, description, authors, manual section). +man_pages = [ + ('index', 'sphericalgeometrytoolkit', u'Spherical Geometry Toolkit Documentation', + [u'Michael Droettboom, STScI'], 1) +] + +# If true, show URL addresses after external links. +#man_show_urls = False + + +# -- Options for Texinfo output ------------------------------------------------ + +# Grouping the document tree into Texinfo files. List of tuples +# (source start file, target name, title, author, +# dir menu entry, description, category) +texinfo_documents = [ + ('index', 'SphericalGeometryToolkit', u'Spherical Geometry Toolkit Documentation', u'Michael Droettboom, STScI', + 'SphericalGeometryToolkit', 'One line description of project.', 'Miscellaneous'), +] + +# Documents to append as an appendix to all manuals. +#texinfo_appendices = [] + +# If false, no module index is generated. +#texinfo_domain_indices = True + +# How to display URL addresses: 'footnote', 'no', or 'inline'. +#texinfo_show_urls = 'footnote' + + +# Example configuration for intersphinx: refer to the Python standard library. +intersphinx_mapping = {'http://docs.python.org/': None} diff --git a/doc/source/cutlines.png b/doc/source/cutlines.png Binary files differnew file mode 100644 index 0000000..6766848 --- /dev/null +++ b/doc/source/cutlines.png diff --git a/doc/source/index.rst b/doc/source/index.rst new file mode 100644 index 0000000..680a41a --- /dev/null +++ b/doc/source/index.rst @@ -0,0 +1,23 @@ +.. Spherical Geometry Toolkit documentation master file, created by + sphinx-quickstart on Tue Oct 11 09:04:48 2011. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to Spherical Geometry Toolkit's documentation! +====================================================== + +Contents: + +.. toctree:: + :maxdepth: 2 + + user + api + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` + diff --git a/doc/source/inside.png b/doc/source/inside.png Binary files differnew file mode 100644 index 0000000..a56246f --- /dev/null +++ b/doc/source/inside.png diff --git a/doc/source/user.rst b/doc/source/user.rst new file mode 100644 index 0000000..e17db6f --- /dev/null +++ b/doc/source/user.rst @@ -0,0 +1,160 @@ +User documentation +================== + +.. currentmodule:: sphere + +The `sphere` library is a pure Python package for handling spherical +polygons that represent arbitrary regions of the sky. + +Requirements +------------ + +- Python 2.7 + +- Numpy 1.4 or later + +Coordinate representation +------------------------- + +Coordinates in world space are traditionally represented by right +ascension and declination (*ra* and *dec*), or longitude and latitude. +While these representations are convenient, they have discontinuities +at the poles, making operations on them trickier at arbitrary +locations on the sky sphere. Therefore, all internal operations of +this library are done in 3D vector space, where coordinates are +represented as (*x*, *y*, *z*) vectors. The `sphere.vector` module +contains functions to convert between (*ra*, *dec*) and (*x*, *y*, +*z*) representations. + +While any (*x*, *y*, *z*) triple represents a vector and therefore a +location on the sky sphere, a distinction must be made between +normalized coordinates fall exactly on the unit sphere, and +unnormalized coordinates which do not. A normalized coordinate is +defined as a vector whose length is 1, i.e.: + +.. math:: + + \sqrt{x^2 + y^2 + z^2} = 1 + +To prevent unnecessary recomputation, many methods in this library +assume that the vectors passed in are already normalized. If this is +not the case, `sphere.vector.normalize_vector` can be used to +normalize an array of vectors. + +The library allows the user to work in either degrees or radians. All +methods that require or return an angular value have a `degrees` +keyword argument. When `degrees` is `True`, these measurements are in +degrees, otherwise they are in radians. + +Spherical polygons +------------------ + +Spherical polygons are arbitrary areas on the sky sphere enclosed by +great circle arcs. They are represented by the +`~sphere.polygon.SphericalPolygon` class. + +Representation +`````````````` + +The points defining the polygon are available from the +`~polygon.SphericalPolygon.points` property. It is a Nx3 array where +each row is an (*x*, *y*, *z*) vector, normalized. The polygon points +are explicitly closed, i.e., the first and last points are the same. + +Where is the inside? +^^^^^^^^^^^^^^^^^^^^ + +The edges of a polygon serve to separate the “inside” from the +“outside” area. On a traditional 2D planar surface, the “inside” is +defined as the finite area and the “outside” is the infinite area. +However, since the surface of a sphere is cyclical, i.e., it wraps +around on itself, the a spherical polygon actually defines two finite +areas. To specify which should be considered the “inside” vs. the +“outside”, the definition of the polygon also has an “inside point” +which is just any point that should be considered inside of the +polygon. + +In the following image, the inside point (marked with the red dot) +declares that the area of the polygon is the green region, and not the +white region. + +.. image:: inside.png + +The inside point of the the polygon can be obtained from the +`~polygon.SphericalPolygon.inside` property. + +Cut lines +^^^^^^^^^ + +If the polygon represents two disjoint areas or the polygon has holes, +those areas will be connected by cut lines. The following image shows +a polygon made from the union of a number of cone areas which has both +a hole and a disjoint region connected by cut lines. + +.. image:: cutlines.png + +Creating spherical polygons +``````````````````````````` + +.. currentmodule:: sphere.polygon + +`SphericalPolygon` objects have 4 different constructors: + + - `SphericalPolygon`: Takes an array of (*x*, *y*, *z*) + points and an inside point. + + - `SphericalPolygon.from_radec`: Takes an array of (*ra*, *dec*) + points and an inside point. + + - `SphericalPolygon.from_cone`: Creates a polygon from a cone on the + sky shere. Takes (*ra*, *dec*, *radius*). + + - `SphericalPolygon.from_wcs`: Creates a polygon from the footprint + of a FITS image using its WCS header keywords. Takes a FITS + filename or a `pyfits.Header` object. + +Operations on Spherical Polygons +```````````````````````````````` + +Once one has a `SphericalPolygon` object, there are a number of +operations available: + + - `~SphericalPolygon.contains_point`: Determines if the given point is inside the polygon. + + - `~SphericalPolygon.intersects_poly`: Determines if one polygon intersects with another. + + - `~SphericalPolygon.area`: Determine the area of a polygon. + + - `~SphericalPolygon.union` and `~SphericalPolygon.multi_union`: + Return a new polygon that is the union of two or more polygons. + + - `~SphericalPolygon.intersection` and + `~SphericalPolygon.multi_intersection`: Return a new polygon that + is the intersection of two or more polygons. + + - `~SphericalPolygon.overlap`: Determine how much a given polygon + overlaps another. + + - `~SphericalPolygon.draw`: Plots the polygon using matplotlib’s + Basemap toolkit. This feature is rather bare and intended + primarily for debugging purposes. + +Great circle arcs +----------------- + +.. currentmodule:: sphere.great_circle_arc + +As seen above, great circle arcs are used to define the edges of the +polygon. The `sphere.great_circle_arc` module contains a number of +functions that are useful for dealing with them. + +- `length`: Returns the angular distance between two points on the sphere. + +- `intersection`: Returns the intersection point between two great + circle arcs. + +- `intersects`: Determines if two great circle arcs intersect. + +- `angle`: Calculate the angle between two great circle arcs. + +- `midpoint`: Calculate the midpoint along a great circle arc. diff --git a/lib/__init__.py b/lib/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/lib/__init__.py diff --git a/lib/graph.py b/lib/graph.py new file mode 100644 index 0000000..e9f03ba --- /dev/null +++ b/lib/graph.py @@ -0,0 +1,792 @@ +# -*- coding: utf-8 -*- + +# Copyright (C) 2011 Association of Universities for Research in +# Astronomy (AURA) +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# +# 1. Redistributions of source code must retain the above +# copyright notice, this list of conditions and the following +# disclaimer. +# +# 2. Redistributions in binary form must reproduce the above +# copyright notice, this list of conditions and the following +# disclaimer in the documentation and/or other materials +# provided with the distribution. +# +# 3. The name of AURA and its representatives may not be used to +# endorse or promote products derived from this software without +# specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR +# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT, +# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +# OF THE POSSIBILITY OF SUCH DAMAGE. + +""" +This contains the code that does the actual unioning of regions. +""" +# TODO: Weak references for memory management problems? + +# STDLIB +import itertools +import weakref + +# THIRD-PARTY +import numpy as np + +# LOCAL +#from . import great_circle_arc +#from . import vector +import great_circle_arc +import vector + +# Set to True to enable some sanity checks +DEBUG = False + + +class Graph: + """ + A graph of nodes connected by edges. The graph is used to build + unions between polygons. + + .. note:: + This class is not meant to be used directly. Instead, use + `sphere.polygon.SphericalPolygon.union` and + `sphere.polygon.SphericalPolygon.intersection`. + """ + + class Node: + """ + A `Node` represents a single point, connected by an arbitrary + number of `Edge` objects to other `Node` objects. + """ + def __init__(self, point, source_polygons=[]): + """ + Parameters + ---------- + point : 3-sequence (*x*, *y*, *z*) coordinate + + source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional + The polygon(s) this node came from. Used for bookkeeping. + """ + self._point = np.asanyarray(point) + self._source_polygons = set(source_polygons) + self._edges = weakref.WeakSet() + + def __repr__(self): + return "Node(%s %d)" % (str(self._point), len(self._edges)) + + def follow(self, edge): + """ + Follows from one edge to another across this node. + + Parameters + ---------- + edge : `Edge` instance + The edge to follow away from. + + Returns + ------- + other : `Edge` instance + The other edge. + """ + edges = list(self._edges) + assert len(edges) == 2 + try: + return edges[not edges.index(edge)] + except IndexError: + raise ValueError("Following from disconnected edge") + + def equals(self, other): + """ + Returns `True` if the location of this and the *other* + `Node` are the same. + """ + return np.array_equal(self._point, other._point) + + + class Edge: + """ + An `Edge` represents a connection between exactly two `Node` + objects. This `Edge` class has no direction. + """ + def __init__(self, A, B, source_polygons=[]): + """ + Parameters + ---------- + A, B : `Node` instances + + source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional + The polygon this edge came from. Used for bookkeeping. + """ + self._nodes = [A, B] + for node in self._nodes: + node._edges.add(self) + self._source_polygons = set(source_polygons) + + def __repr__(self): + nodes = self._nodes + return "Edge(%s -> %s)" % (nodes[0]._point, nodes[1]._point) + + def follow(self, node): + """ + Follow along the edge from the given *node* to the other + node. + + Parameters + ---------- + node : `Node` instance + + Returns + ------- + other : `Node` instance + """ + nodes = self._nodes + try: + return nodes[not nodes.index(node)] + except IndexError: + raise ValueError("Following from disconnected node") + + + def __init__(self, polygons): + """ + Parameters + ---------- + polygons : sequence of `~sphere.polygon.SphericalPolygon` instances + Build a graph from this initial set of polygons. + """ + self._nodes = set() + self._edges = set() + self._source_polygons = set() + + self.add_polygons(polygons) + + def add_polygons(self, polygons): + """ + Add more polygons to the graph. + + .. note:: + Must be called before `union` or `intersection`. + + Parameters + ---------- + polygons : sequence of `~sphere.polygon.SphericalPolygon` instances + Set of polygons to add to the graph + """ + for polygon in polygons: + self.add_polygon(polygon) + + def add_polygon(self, polygon): + """ + Add a single polygon to the graph. + + .. note:: + Must be called before `union` or `intersection`. + + Parameters + ---------- + polygon : `~sphere.polygon.SphericalPolygon` instance + Polygon to add to the graph + """ + points = polygon._points + + if len(points) < 3: + raise ValueError("Too few points in polygon") + + self._source_polygons.add(polygon) + + start_node = nodeA = self.add_node(points[0], [polygon]) + for i in range(1, len(points) - 1): + # Don't create self-pointing edges + if not np.array_equal(points[i], nodeA._point): + nodeB = self.add_node(points[i], [polygon]) + self.add_edge(nodeA, nodeB, [polygon]) + nodeA = nodeB + # Close the polygon + self.add_edge(nodeA, start_node, [polygon]) + + def add_node(self, point, source_polygons=[]): + """ + Add a node to the graph. It will be disconnected until used + in a call to `add_edge`. + + Parameters + ---------- + point : 3-sequence (*x*, *y*, *z*) coordinate + + source_polygon : `~sphere.polygon.SphericalPolygon` instance, optional + The polygon this node came from. Used for bookkeeping. + + Returns + ------- + node : `Node` instance + The new node + """ + node = self.Node(point, source_polygons) + self._nodes.add(node) + return node + + def remove_node(self, node): + """ + Removes a node and all of the edges that touch it. + + .. note:: + It is assumed that *Node* is already a part of the graph. + + Parameters + ---------- + node : `Node` instance + """ + assert node in self._nodes + + for edge in list(node._edges): + nodeB = edge.follow(node) + nodeB._edges.remove(edge) + self._edges.remove(edge) + self._nodes.remove(node) + + def add_edge(self, A, B, source_polygons=[]): + """ + Add an edge between two nodes. + + .. note:: + It is assumed both nodes already belong to the graph. + + Parameters + ---------- + A, B : `Node` instances + + source_polygons : `~sphere.polygon.SphericalPolygon` instance, optional + The polygon(s) this edge came from. Used for bookkeeping. + + Returns + ------- + edge : `Edge` instance + The new edge + """ + assert A in self._nodes + assert B in self._nodes + + edge = self.Edge(A, B, source_polygons) + self._edges.add(edge) + return edge + + def remove_edge(self, edge): + """ + Remove an edge from the graph. The nodes it points to remain intact. + + .. note:: + It is assumed that *edge* is already a part of the graph. + + Parameters + ---------- + edge : `Edge` instance + """ + assert edge in self._edges + + A, B = edge._nodes + A._edges.remove(edge) + if len(A._edges) == 0: + self.remove_node(A) + if A is not B: + B._edges.remove(edge) + if len(B._edges) == 0: + self.remove_node(B) + self._edges.remove(edge) + + def split_edge(self, edge, node): + """ + Splits an `Edge` *edge* at `Node` *node*, removing *edge* and + replacing it with two new `Edge` instances. It is intended + that *E* is along the original edge, but that is not enforced. + + Parameters + ---------- + edge : `Edge` instance + The edge to split + + node : `Node` instance + The node to insert + + Returns + ------- + edgeA, edgeB : `Edge` instances + The two new edges on either side of *node*. + """ + assert edge in self._edges + assert node in self._nodes + + A, B = edge._nodes + edgeA = self.add_edge(A, node, edge._source_polygons) + edgeB = self.add_edge(node, B, edge._source_polygons) + self.remove_edge(edge) + return [edgeA, edgeB] + + def _sanity_check(self, node_is_2=False): + """ + For debugging purposes: assert that edges and nodes are + connected to each other correctly and there are no orphaned + edges or nodes. + """ + if not DEBUG: + return + + try: + unique_edges = set() + for edge in self._edges: + for node in edge._nodes: + assert edge in node._edges + assert node in self._nodes + edge_repr = [tuple(x._point) for x in edge._nodes] + edge_repr.sort() + edge_repr = tuple(edge_repr) + # assert edge_repr not in unique_edges + unique_edges.add(edge_repr) + + for node in self._nodes: + if node_is_2: + assert len(node._edges) == 2 + else: + assert len(node._edges) >= 2 + for edge in node._edges: + assert node in edge._nodes + assert edge in self._edges + except AssertionError as e: + import traceback + traceback.print_exc() + self._dump_graph() + raise + + def _dump_graph(self, title=None, lon_0=0, lat_0=90, + projection='ortho', func=lambda x: len(x._edges)): + from mpl_toolkits.basemap import Basemap + from matplotlib import pyplot as plt + fig = plt.figure() + m = Basemap(projection="ortho", + lon_0=lon_0, + lat_0=lat_0) + m.drawparallels(np.arange(-90., 90., 20.)) + m.drawmeridians(np.arange(0., 420., 20.)) + m.drawmapboundary(fill_color='white') + + counts = {} + for node in self._nodes: + count = func(node) + counts.setdefault(count, []) + counts[count].append(list(node._point)) + + for k, v in counts.items(): + v = np.array(v) + ra, dec = vector.vector_to_radec(v[:, 0], v[:, 1], v[:, 2]) + x, y = m(ra, dec) + m.plot(x, y, 'o', label=str(k)) + + for edge in list(self._edges): + A, B = [x._point for x in edge._nodes] + r0, d0 = vector.vector_to_radec(A[0], A[1], A[2]) + r1, d1 = vector.vector_to_radec(B[0], B[1], B[2]) + m.drawgreatcircle(r0, d0, r1, d1, color='blue') + + if title: + plt.title(title) + plt.legend() + plt.show() + + def union(self): + """ + Once all of the polygons have been added to the graph, + join the polygons together. + + Returns + ------- + points : Nx3 array of (*x*, *y*, *z*) points + This is a list of points outlining the union of the + polygons that were given to the constructor. If the + original polygons are disjunct or contain holes, cut lines + will be included in the output. + """ + self._remove_cut_lines() + self._sanity_check() + self._find_all_intersections() + self._sanity_check() + self._remove_interior_nodes() + self._sanity_check() + self._remove_3ary_edges() + self._sanity_check(True) + return self._trace() + + def intersection(self): + """ + Once all of the polygons have been added to the graph, + calculate the intersection. + + Returns + ------- + points : Nx3 array of (*x*, *y*, *z*) points + This is a list of points outlining the intersection of the + polygons that were given to the constructor. If the + resulting polygons are disjunct or contain holes, cut lines + will be included in the output. + """ + self._remove_cut_lines() + self._sanity_check() + self._find_all_intersections() + self._sanity_check() + self._remove_exterior_edges() + self._sanity_check() + self._remove_3ary_edges(large_first=True) + self._sanity_check(True) + return self._trace() + + def _remove_cut_lines(self): + """ + Removes any cutlines that may already have existed in the + input polygons. This is so any cutlines in the final result + will be optimized to be as short as possible and won't + intersect each other. + + This works by finding coincident edges that are reverse to + each other, and then splicing around them. + """ + # As this proceeds, edges are removed from the graph. It + # iterates over a static list of all edges that exist at the + # start, so each time one is selected, we need to ensure it + # still exists as part of the graph. + + # This transforms the following (where = is the cut line) + # + # \ / + # A' + + B' + # | | + # A +====================+ B + # + # D +====================+ C + # | | + # D' + + C' + # / \ + # + # to this: + # + # \ / + # A' + + B' + # | | + # A + + C + # | | + # D' + + C' + # / \ + # + + edges = list(self._edges) + for i in xrange(len(edges)): + AB = edges[i] + if AB not in self._edges: + continue + A, B = AB._nodes + for j in xrange(i + 1, len(edges)): + CD = edges[j] + if CD not in self._edges: + continue + C, D = CD._nodes + # To be a cutline, the candidate edges need to run in + # the opposite direction, hence A == D and B == C, not + # A == C and B == D. + if (A.equals(D) and B.equals(C)): + # Create new edges A -> D' and C -> B' + self.add_edge( + A, D.follow(CD).follow(D), + AB._source_polygons | CD._source_polygons) + self.add_edge( + C, B.follow(AB).follow(B), + AB._source_polygons | CD._source_polygons) + + # Remove B and D which are identical to C and A + # respectively. We do not need to remove AB and + # CD because this will remove it for us. + self.remove_node(D) + self.remove_node(B) + + break + + def _find_all_intersections(self): + """ + Find all the intersecting edges in the graph. For each + intersecting pair, four new edges are created around the + intersection point. + """ + def get_edge_points(edges): + return (np.array([x._nodes[0]._point for x in edges]), + np.array([x._nodes[1]._point for x in edges])) + + # For speed, we want to vectorize all of the intersection + # calculations. Therefore, there is a list of edges, and two + # arrays containing the end points of those edges. They all + # need to have things added and removed from them at the same + # time to keep them in sync, but of course the interface for + # doing so is different between Python lists and numpy arrays. + + edges = list(self._edges) + starts, ends = get_edge_points(edges) + + while len(edges) > 1: + AB = edges.pop(0) + A = starts[0]; starts = starts[1:] # numpy equiv of "pop(0)" + B = ends[0]; ends = ends[1:] # numpy equiv of "pop(0)" + + # Calculate the intersection points between AB and all + # other remaining edges + with np.errstate(invalid='ignore'): + intersections = great_circle_arc.intersection( + A, B, starts, ends) + # intersects is `True` everywhere intersections has an + # actual intersection + intersects = np.isfinite(intersections[..., 0]) + + intersection_indices = np.nonzero(intersects)[0] + + # Iterate through the candidate intersections, if any -- + # we want to eliminate intersections that only intersect + # at the end points + for j in intersection_indices: + C = starts[j] + D = ends[j] + if (np.array_equal(C, A) or np.array_equal(C, B) or + np.array_equal(D, A) or np.array_equal(D, B)): + continue + + CD = edges[j] + E = intersections[j] + + # This is a bona-fide intersection, and E is the + # point at which the two lines intersect. Make a + # new node for it -- this must belong to the all + # of the source polygons of both of the edges that + # crossed. + + # A + # | + # C--E--D + # | + # B + + E = self.add_node( + E, AB._source_polygons | CD._source_polygons) + newA, newB = self.split_edge(AB, E) + newC, newD = self.split_edge(CD, E) + + new_edges = [newA, newB, newC, newD] + + # Delete CD, and push the new edges to the + # front so they will be tested for intersection + # against all remaining edges. + del edges[j] # CD + edges = new_edges + edges + new_starts, new_ends = get_edge_points(new_edges) + starts = np.vstack( + (new_starts, starts[:j], starts[j+1:])) + ends = np.vstack( + (new_ends, ends[:j], ends[j+1:])) + break + + def _remove_interior_nodes(self): + """ + Removes any nodes that are contained inside other polygons. + What's left is the (possibly disjunct) outline. + """ + polygons = self._source_polygons + + # node._count is the number of polygons that the node is + # inside, other than the polygon that it came from + for node in self._nodes: + node._count = 0 + for polygon in polygons: + if polygon in node._source_polygons: + continue + if polygon.contains_point(node._point): + node._count += 1 + + for node in list(self._nodes): + # Nodes with a count >= 1 are completely contained in the + # outer polygon, so they should be removed. What's left + # are only outer nodes. + if node._count >= 1: + self.remove_node(node) + + def _remove_exterior_edges(self): + """ + Removes any edges that are not contained in all of the source + polygons. What's left is the (possibly disjunct) outline. + """ + polygons = self._source_polygons + + for edge in self._edges: + edge._count = 0 + for polygon in polygons: + if (polygon in edge._source_polygons or + polygon.intersects_arc( + edge._nodes[0]._point, edge._nodes[1]._point)): + edge._count += 1 + + for edge in list(self._edges): + if edge._count < len(polygons): + self.remove_edge(edge) + + def _remove_3ary_edges(self, large_first=False): + """ + Remove edges between pairs of nodes that have more than 3 + edges. This removes triangles that can't be traced. + """ + if large_first: + max_ary = 0 + for node in self._nodes: + max_ary = max(len(node._edges), max_ary) + order = range(max_ary + 1, 2, -1) + else: + order = [3] + + for i in order: + removals = [] + for edge in list(self._edges): + if (len(edge._nodes[0]._edges) >= i and + len(edge._nodes[1]._edges) >= i): + removals.append(edge) + + for edge in removals: + if edge in self._edges: + self.remove_edge(edge) + + def _trace(self): + """ + Given a graph that has had cutlines removed and all + intersections found, traces it to find a resulting single + polygon. + """ + polygons = [] + edges = set(self._edges) # copy + seen_nodes = set() + while len(edges): + polygon = [] + edge = edges.pop() + start_node = node = edge._nodes[0] + while True: + # TODO: Do we need this if clause any more? + if len(polygon): + if not np.array_equal(polygon[-1], node._point): + polygon.append(node._point) + else: + polygon.append(node._point) + edge = node.follow(edge) + edges.discard(edge) + node = edge.follow(node) + if node is start_node: + polygon.append(node._point) + break + + polygons.append(np.asarray(polygon)) + + if len(polygons) == 1: + return polygons[0] + elif len(polygons) == 0: + return [] + else: + return self._join(polygons) + + def _join(self, polygons): + """ + If the graph is disjunct, joins the parts with cutlines. + + The closest nodes between each pair that don't intersect + any other edges are used as cutlines. + + TODO: This is not optimal, because the closest distance + between two polygons may not in fact be between two vertices, + but may be somewhere along an edge. + """ + def do_join(polygons): + all_polygons = polygons[:] + + skipped = 0 + + polyA = polygons.pop(0) + while len(polygons): + polyB = polygons.pop(0) + + # If fewer than 3 edges, it's not a polygon, + # just throw it out + if len(polyB) < 4: + continue + + # Find the closest set of vertices between polyA and + # polyB that don't cross any of the edges in *any* of + # the polygons + closest = np.inf + closest_pair_idx = (None, None) + for a in xrange(len(polyA) - 1): + A = polyA[a] + distances = great_circle_arc.length(A, polyB[:-1]) + b = np.argmin(distances) + distance = distances[b] + if distance < closest: + B = polyB[b] + # Does this candidate line cross other edges? + crosses = False + for poly in all_polygons: + if np.any( + great_circle_arc.intersects( + A, B, poly[:-1], poly[1:])): + crosses = True + break + if not crosses: + closest = distance + closest_pair_idx = (a, b) + + if not np.isfinite(closest): + # We didn't find a pair of points that don't cross + # something else, so we want to try to join another + # polygon. Defer the current polygon until later. + if len(polygons) in (0, skipped): + return None + polygons.append(polyB) + skipped += 1 + else: + # Splice the two polygons together using a cut + # line + a, b = closest_pair_idx + new_poly = np.vstack(( + # polyA up to and including the cut point + polyA[:a+1], + + # polyB starting with the cut point and + # wrapping around back to the cut point. + # Ignore the last point in polyB, because it + # is the same as the first + np.roll(polyB[:-1], -b, 0), + + # The cut point on polyB + polyB[b:b+1], + + # the rest of polyA, starting with the cut + # point + polyA[a:] + )) + + skipped = 0 + polyA = new_poly + + return polyA + + for permutation in itertools.permutations(polygons): + poly = do_join(list(permutation)) + if poly is not None: + return poly + + raise RuntimeError("Could not find cut points") diff --git a/lib/great_circle_arc.py b/lib/great_circle_arc.py new file mode 100644 index 0000000..0fff174 --- /dev/null +++ b/lib/great_circle_arc.py @@ -0,0 +1,318 @@ +# -*- coding: utf-8 -*- + +# Copyright (C) 2011 Association of Universities for Research in +# Astronomy (AURA) +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# +# 1. Redistributions of source code must retain the above +# copyright notice, this list of conditions and the following +# disclaimer. +# +# 2. Redistributions in binary form must reproduce the above +# copyright notice, this list of conditions and the following +# disclaimer in the documentation and/or other materials +# provided with the distribution. +# +# 3. The name of AURA and its representatives may not be used to +# endorse or promote products derived from this software without +# specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR +# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT, +# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +# OF THE POSSIBILITY OF SUCH DAMAGE. + +""" +The `sphere.great_circle_arc` module contains functions for computing +the length, intersection, angle and midpoint of great circle arcs. + +Great circles are circles on the unit sphere whose center is +coincident with the center of the sphere. Great circle arcs are the +section of those circles between two points on the unit sphere. +""" + +from __future__ import with_statement, division, absolute_import + +# THIRD-PARTY +import numpy as np + + +try: + #from . import math_util + import math_util + HAS_C_UFUNCS = True +except ImportError: + HAS_C_UFUNCS = False + + +__all__ = ['angle', 'intersection', 'intersects', 'length', 'midpoint'] + + +if not HAS_C_UFUNCS: + from numpy.core.umath_tests import inner1d + + def _fast_cross(a, b): + """ + This is a reimplementation of `numpy.cross` that only does 3D x + 3D, and is therefore faster since it doesn't need any + conditionals. + """ + cp = np.empty(np.broadcast(a, b).shape) + aT = a.T + bT = b.T + cpT = cp.T + + cpT[0] = aT[1]*bT[2] - aT[2]*bT[1] + cpT[1] = aT[2]*bT[0] - aT[0]*bT[2] + cpT[2] = aT[0]*bT[1] - aT[1]*bT[0] + + return cp + + def _cross_and_normalize(A, B): + T = _fast_cross(A, B) + # Normalization + l = np.sqrt(np.sum(T ** 2, axis=-1)) + l = np.expand_dims(l, 2) + # Might get some divide-by-zeros, but we don't care + with np.errstate(invalid='ignore'): + TN = T / l + return TN + + def intersection(A, B, C, D): + ur""" + Returns the point of intersection between two great circle arcs. + The arcs are defined between the points *AB* and *CD*. Either *A* + and *B* or *C* and *D* may be arrays of points, but not both. + + Parameters + ---------- + A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples + Endpoints of the first great circle arc. + + C, D : (*x*, *y*, *z*) triples or Nx3 arrays of triples + Endpoints of the second great circle arc. + + Returns + ------- + T : (*x*, *y*, *z*) triples or Nx3 arrays of triples + If the given arcs intersect, the intersection is returned. If + the arcs do not intersect, the triple is set to all NaNs. + + Notes + ----- + The basic intersection is computed using linear algebra as follows + [1]_: + + .. math:: + + T = \lVert(A × B) × (C × D)\rVert + + To determine the correct sign (i.e. hemisphere) of the + intersection, the following four values are computed: + + .. math:: + + s_1 = ((A × B) × A) · T + + s_2 = (B × (A × B)) · T + + s_3 = ((C × D) × C) · T + + s_4 = (D × (C × D)) · T + + For :math:`s_n`, if all positive :math:`T` is returned as-is. If + all negative, :math:`T` is multiplied by :math:`-1`. Otherwise + the intersection does not exist and is undefined. + + References + ---------- + + .. [1] Method explained in an `e-mail + <http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271>`_ + by Roger Stafford. + + http://www.mathworks.com/matlabcentral/newsreader/view_thread/276271 + """ + A = np.asanyarray(A) + B = np.asanyarray(B) + C = np.asanyarray(C) + D = np.asanyarray(D) + + A, B = np.broadcast_arrays(A, B) + C, D = np.broadcast_arrays(C, D) + + ABX = _fast_cross(A, B) + CDX = _fast_cross(C, D) + T = _cross_and_normalize(ABX, CDX) + T_ndim = len(T.shape) + + if T_ndim > 1: + s = np.zeros(T.shape[0]) + else: + s = np.zeros(1) + s += np.sign(inner1d(_fast_cross(ABX, A), T)) + s += np.sign(inner1d(_fast_cross(B, ABX), T)) + s += np.sign(inner1d(_fast_cross(CDX, C), T)) + s += np.sign(inner1d(_fast_cross(D, CDX), T)) + if T_ndim > 1: + s = np.expand_dims(s, 2) + + cross = np.where(s == -4, -T, np.where(s == 4, T, np.nan)) + + # If they share a common point, it's not an intersection. This + # gets around some rounding-error/numerical problems with the + # above. + equals = (np.all(A == C, axis = -1) | + np.all(A == D, axis = -1) | + np.all(B == C, axis = -1) | + np.all(B == D, axis = -1)) + + equals = np.expand_dims(equals, 2) + + return np.where(equals, np.nan, cross) +else: + inner1d = math_util.inner1d + _fast_cross = math_util.cross + _cross_and_normalize = math_util.cross_and_norm + intersection = math_util.intersection + + +def length(A, B, degrees=True): + ur""" + Returns the angular distance between two points (in vector space) + on the unit sphere. + + Parameters + ---------- + A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples + The endpoints of the great circle arc, in vector space. + + degrees : bool, optional + If `True` (default) the result is returned in decimal degrees, + otherwise radians. + + Returns + ------- + length : scalar or array of scalars + The angular length of the great circle arc. + + Notes + ----- + The length is computed using the following: + + .. math:: + + \Delta = \arccos(A · B) + """ + A = np.asanyarray(A) + B = np.asanyarray(B) + + dot = inner1d(A, B) + with np.errstate(invalid='ignore'): + result = np.arccos(dot) + + if degrees: + return np.rad2deg(result) + else: + return result + + +def intersects(A, B, C, D): + """ + Returns `True` if the great circle arcs between *AB* and *CD* + intersect. Either *A* and *B* or *C* and *D* may be arrays of + points, but not both. + + Parameters + ---------- + A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples + Endpoints of the first great circle arc. + + C, D : (*x*, *y*, *z*) triples or Nx3 arrays of triples + Endpoints of the second great circle arc. + + Returns + ------- + intersects : bool or array of bool + If the given arcs intersect, the intersection is returned as + `True`. + """ + with np.errstate(invalid='ignore'): + intersections = intersection(A, B, C, D) + return np.isfinite(intersections[..., 0]) + + +def angle(A, B, C, degrees=True): + """ + Returns the angle at *B* between *AB* and *BC*. + + This always returns the shortest angle < π. + + Parameters + ---------- + A, B, C : (*x*, *y*, *z*) triples or Nx3 arrays of triples. + + degrees : bool, optional + If `True` (default) the result is returned in decimal degrees, + otherwise radians. + + Returns + ------- + angle : float or array of floats + The angle at *B* between *AB* and *BC*. + + References + ---------- + + .. [1] Miller, Robert D. Computing the area of a spherical + polygon. Graphics Gems IV. 1994. Academic Press. + """ + A = np.asanyarray(A) + B = np.asanyarray(B) + C = np.asanyarray(C) + + A, B, C = np.broadcast_arrays(A, B, C) + + ABX = _fast_cross(A, B) + ABX = _cross_and_normalize(B, ABX) + BCX = _fast_cross(C, B) + BCX = _cross_and_normalize(B, BCX) + with np.errstate(invalid='ignore'): + angle = np.arccos(inner1d(ABX, BCX)) + + if degrees: + angle = np.rad2deg(angle) + return angle + + +def midpoint(A, B): + """ + Returns the midpoint on the great circle arc between *A* and *B*. + + Parameters + ---------- + A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples + The endpoints of the great circle arc. It is assumed that + these points are already normalized. + + Returns + ------- + midpoint : (*x*, *y*, *z*) triple or Nx3 arrays of triples + The midpoint between *A* and *B*, normalized on the unit + sphere. + """ + P = (A + B) / 2.0 + # Now normalize... + l = np.sqrt(np.sum(P * P, axis=-1)) + l = np.expand_dims(l, 2) + return P / l diff --git a/lib/polygon.py b/lib/polygon.py new file mode 100644 index 0000000..b9b5d8c --- /dev/null +++ b/lib/polygon.py @@ -0,0 +1,679 @@ +# -*- coding: utf-8 -*- + +# Copyright (C) 2011 Association of Universities for Research in +# Astronomy (AURA) +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# +# 1. Redistributions of source code must retain the above +# copyright notice, this list of conditions and the following +# disclaimer. +# +# 2. Redistributions in binary form must reproduce the above +# copyright notice, this list of conditions and the following +# disclaimer in the documentation and/or other materials +# provided with the distribution. +# +# 3. The name of AURA and its representatives may not be used to +# endorse or promote products derived from this software without +# specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR +# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT, +# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +# OF THE POSSIBILITY OF SUCH DAMAGE. + +""" +The `polygon` module defines the `SphericalPolygon` class for managing +polygons on the unit sphere. +""" +from __future__ import division, print_function + +# STDLIB +from copy import copy + +# THIRD-PARTY +import numpy as np + +# LOCAL +#from . import great_circle_arc +#from . import vector +import great_circle_arc +import vector + +__all__ = ['SphericalPolygon'] + + +class SphericalPolygon(object): + ur""" + Polygons are represented by both a set of points (in + Cartesian (*x*, *y*, *z*) normalized on the unit sphere), + and an inside point. The inside point is necessary, because + both the inside and outside of the polygon are finite areas + on the great sphere, and therefore we need a way of + specifying which is which. + """ + + def __init__(self, points, inside): + ur""" + Parameters + ---------- + points : An Nx3 array of (*x*, *y*, *z*) triples in vector space + These points define the boundary of the polygon. It must + be "closed", i.e., the last point is the same as the first. + + It may contain zero points, in which it defines the null + polygon. It may not contain one, two or three points. + Four points are needed to define a triangle, since the + polygon must be closed. + + inside : An (*x*, *y*, *z*) triple + This point must be inside the polygon. + """ + if len(points) == 0: + pass + elif len(points) < 3: + raise ValueError("Polygon made of too few points") + else: + assert np.array_equal(points[0], points[-1]), 'Polygon is not closed' + self._points = np.asanyarray(points) + self._inside = np.asanyarray(inside) + + # TODO: Detect self-intersection and fix + + def __copy__(self): + return self.__class__(np.copy(self._points), np.copy(self._inside)) + + @property + def points(self): + """ + The points defining the polygon. It is an Nx3 array of + (*x*, *y*, *z*) vectors. The polygon will be explicitly + closed, i.e., the first and last points are the same. + """ + return self._points + + @property + def inside(self): + """ + Get the inside point of the polygon. + """ + return self._inside + + @classmethod + def from_radec(cls, ra, dec, center=None, degrees=True): + ur""" + Create a new `SphericalPolygon` from a list of (*ra*, *dec*) + points. + + Parameters + ---------- + ra, dec : 1-D arrays of the same length + The vertices of the polygon in right-ascension and + declination. It must be \"closed\", i.e., that is, the + last point is the same as the first. + + center : (*ra*, *dec*) pair, optional + A point inside of the polygon to define its inside. If no + *center* point is provided, the mean of the polygon's + points in vector space will be used. That approach may + not work for concave polygons. + + degrees : bool, optional + If `True`, (default) *ra* and *dec* are in decimal degrees, + otherwise in radians. + + Returns + ------- + polygon : `SphericalPolygon` object + """ + # Convert to Cartesian + x, y, z = vector.radec_to_vector(ra, dec, degrees=degrees) + + if center is None: + xc = x.mean() + yc = y.mean() + zc = z.mean() + center = vector.normalize_vector(xc, yc, zc) + else: + center = vector.radec_to_vector(*center, degrees=degrees) + + return cls(np.dstack((x, y, z))[0], center) + + @classmethod + def from_cone(cls, ra, dec, radius, degrees=True, steps=16.0): + ur""" + Create a new `SphericalPolygon` from a cone (otherwise known + as a "small circle") defined using (*ra*, *dec*, *radius*). + + The cone is not represented as an ideal circle on the sphere, + but as a series of great circle arcs. The resolution of this + conversion can be controlled using the *steps* parameter. + + Parameters + ---------- + ra, dec : float scalars + This defines the center of the cone + + radius : float scalar + The radius of the cone + + degrees : bool, optional + If `True`, (default) *ra*, *dec* and *radius* are in + decimal degrees, otherwise in radians. + + steps : int, optional + The number of steps to use when converting the small + circle to a polygon. + + Returns + ------- + polygon : `SphericalPolygon` object + """ + u, v, w = vector.radec_to_vector(ra, dec, degrees=degrees) + if degrees: + radius = np.deg2rad(radius) + + # Get an arbitrary perpendicular vector. This be be obtained + # by crossing (u, v, w) with any unit vector that is not itself. + which_min = np.argmin([u, v, w]) + if which_min == 0: + perp = np.cross([u, v, w], [1., 0., 0.]) + elif which_min == 1: + perp = np.cross([u, v, w], [0., 1., 0.]) + else: + perp = np.cross([u, v, w], [0., 0., 1.]) + + # Rotate by radius around the perpendicular vector to get the + # "pen" + x, y, z = vector.rotate_around( + u, v, w, perp[0], perp[1], perp[2], radius, degrees=False) + + # Then rotate the pen around the center point all 360 degrees + C = np.linspace(0, np.pi * 2.0, steps) + # Ensure that the first and last elements are exactly the + # same. 2π should equal 0, but with rounding error that isn't + # always the case. + C[-1] = 0 + X, Y, Z = vector.rotate_around(x, y, z, u, v, w, C, degrees=False) + + return cls(np.dstack((X, Y, Z))[0], (u, v, w)) + + @classmethod + def from_wcs(cls, fitspath, steps=1, crval=None): + ur""" + Create a new `SphericalPolygon` from the footprint of a FITS + WCS specification. + + This method requires having `pywcs` and `pyfits` installed. + + Parameters + ---------- + fitspath : path to a FITS file, `pyfits.Header`, or `pywcs.WCS` + Refers to a FITS header containing a WCS specification. + + steps : int, optional + The number of steps along each edge to convert into + polygon edges. + + Returns + ------- + polygon : `SphericalPolygon` object + """ + import pywcs + import pyfits + + if isinstance(fitspath, pyfits.Header): + header = fitspath + wcs = pywcs.WCS(header) + elif isinstance(fitspath, pywcs.WCS): + wcs = fitspath + else: + wcs = pywcs.WCS(fitspath) + if crval is not None: + wcs.wcs.crval = crval + xa, ya = [wcs.naxis1, wcs.naxis2] + + length = steps * 4 + 1 + X = np.empty(length) + Y = np.empty(length) + + # Now define each of the 4 edges of the quadrilateral + X[0 :steps ] = np.linspace(1, xa, steps, False) + Y[0 :steps ] = 1 + X[steps :steps*2] = xa + Y[steps :steps*2] = np.linspace(1, ya, steps, False) + X[steps*2:steps*3] = np.linspace(xa, 1, steps, False) + Y[steps*2:steps*3] = ya + X[steps*3:steps*4] = 1 + Y[steps*3:steps*4] = np.linspace(ya, 1, steps, False) + X[-1] = 1 + Y[-1] = 1 + + # Use wcslib to convert to (ra, dec) + ra, dec = wcs.all_pix2sky(X, Y, 1) + + # Convert to Cartesian + x, y, z = vector.radec_to_vector(ra, dec) + + # Calculate an inside point + ra, dec = wcs.all_pix2sky(xa / 2.0, ya / 2.0, 1) + xc, yc, zc = vector.radec_to_vector(ra, dec) + + return cls(np.dstack((x, y, z))[0], (xc[0], yc[0], zc[0])) + + def contains_point(self, point): + ur""" + Determines if this `SphericalPolygon` contains a given point. + + Parameters + ---------- + point : an (*x*, *y*, *z*) triple + The point to test. + + Returns + ------- + contains : bool + Returns `True` if the polygon contains the given *point*. + """ + P = self._points + r = self._inside + point = np.asanyarray(point) + + intersects = great_circle_arc.intersects(P[:-1], P[1:], r, point) + crossings = np.sum(intersects) + + return (crossings % 2) == 0 + + def intersects_poly(self, other): + ur""" + Determines if this `SphericalPolygon` intersects another + `SphericalPolygon`. + + This method is much faster than actually computing the + intersection region between two polygons. + + Parameters + ---------- + other : `SphericalPolygon` + + Returns + ------- + intersects : bool + Returns `True` if this polygon intersects the *other* + polygon. + + Notes + ----- + + The algorithm proceeds as follows: + + 1. Determine if any single point of one polygon is contained + within the other. + + 2. Deal with the case where only the edges overlap as in:: + + : o---------o + : o----+---------+----o + : | | | | + : o----+---------+----o + : o---------o + + In this case, an edge from one polygon must cross an + edge from the other polygon. + """ + assert isinstance(other, SphericalPolygon) + + # The easy case is in which a point of one polygon is + # contained in the other polygon. + for point in other._points: + if self.contains_point(point): + return True + for point in self._points: + if other.contains_point(point): + return True + + # The hard case is when only the edges overlap, as in: + # + # o---------o + # o----+---------+----o + # | | | | + # o----+---------+----o + # o---------o + # + for i in range(len(self._points) - 1): + A = self._points[i] + B = self._points[i+1] + if np.any(great_circle_arc.intersects( + A, B, other._points[:-1], other._points[1:])): + return True + return False + + def intersects_arc(self, a, b): + """ + Determines if this `SphericalPolygon` intersects or contains + the given arc. + """ + return self.contains_point(great_circle_arc.midpoint(a, b)) + + if self.contains_point(a) and self.contains_point(b): + return True + + P = self._points + intersects = great_circle_arc.intersects(P[:-1], P[1:], a, b) + if np.any(intersects): + return True + return False + + def area(self): + ur""" + Returns the area of the polygon on the unit sphere. + + The algorithm is not able to compute the area of polygons + that are larger than half of the sphere. Therefore, the + area will always be less than 2π. + + The area is computed based on the sum of the radian angles: + + .. math:: + + S = \sum^n_{i=0} \theta_i - (n - 2) \pi + """ + if len(self._points) < 3: + return np.array(0.0) + + A = self._points[:-1] + B = np.roll(A, 1, 0) + C = np.roll(B, 1, 0) + + angles = great_circle_arc.angle(A, B, C, degrees=False) + + midpoints = great_circle_arc.midpoint(A, C) + interior = [self.contains_point(x) for x in midpoints] + angles = np.where(interior, angles, 2.*np.pi - angles) + + sum = np.sum(angles) + area = sum - float(len(angles) - 2) * np.pi + + return area + + def union(self, other): + """ + Return a new `SphericalPolygon` that is the union of *self* + and *other*. + + If the polygons are disjoint, they result will be connected + using cut lines. For example:: + + : o---------o + : | | + : o---------o=====o----------o + : | | + : o----------o + + Parameters + ---------- + other : `SphericalPolygon` + + Returns + ------- + polygon : `SphericalPolygon` object + + See also + -------- + multi_union + + Notes + ----- + For implementation details, see the :mod:`~sphere.graph` + module. + """ + import graph + g = graph.Graph([self, other]) + + polygon = g.union() + + return self.__class__(polygon, self._inside) + + @classmethod + def multi_union(cls, polygons, method='parallel'): + """ + Return a new `SphericalPolygon` that is the union of all of the + polygons in *polygons*. + + Parameters + ---------- + polygons : sequence of `SphericalPolygon` + + method : 'parallel' or 'serial', optional + Specifies the method that is used to perform the unions: + + - 'parallel' (default): A graph is built using all of + the polygons, and the union operation is computed on + the entire thing globally. + + - 'serial': The polygon is built in steps by adding one + polygon at a time and unioning at each step. + + This option is provided because one may be faster than the + other depending on context, but it primarily exposed for + testing reasons. Both modes should theoretically provide + equivalent results. + + Returns + ------- + polygon : `SphericalPolygon` object + + See also + -------- + union + """ + assert len(polygons) + for polygon in polygons: + assert isinstance(polygon, SphericalPolygon) + + import graph + + if method.lower() == 'parallel': + g = graph.Graph(polygons) + polygon = g.union() + return cls(polygon, polygons[0]._inside) + elif method.lower() == 'serial': + result = copy(polygons[0]) + for polygon in polygons[1:]: + result = result.union(polygon) + return result + else: + raise ValueError("method must be 'parallel' or 'serial'") + + @staticmethod + def _find_new_inside(points, polygons): + """ + Finds an acceptable inside point inside of *points* that is + also inside of *polygons*. Used by the intersection + algorithm, and is really only useful in that context because + it requires existing polygons with known inside points. + """ + if len(points) < 4: + return np.array([0, 0, 0]) + + # Special case for a triangle + if len(points) == 4: + return np.sum(points[:3]) / 3.0 + + for i in range(len(points) - 1): + A = points[i] + # Skip the adjacent point, since it is by definition on + # the edge of the polygon, not potentially running through + # the middle. + for j in range(i + 2, len(points) - 1): + B = points[j] + C = great_circle_arc.midpoint(A, B) + in_all = True + for polygon in polygons: + if not polygon.contains_point(C): + in_all = False + break + if in_all: + return C + + raise RuntimeError("Suitable inside point could not be found") + + def intersection(self, other): + """ + Return a new `SphericalPolygon` that is the intersection of + *self* and *other*. + + If the intersection is empty, a `SphericalPolygon` with zero + points will be returned. + + If the result is disjoint, the pieces will be connected using + cut lines. For example:: + + : o---------o + : | | + : o---------o=====o----------o + : | | + : o----------o + + Parameters + ---------- + other : `SphericalPolygon` + + Returns + ------- + polygon : `SphericalPolygon` object + + Notes + ----- + For implementation details, see the :mod:`~sphere.graph` + module. + """ + # if not self.intersects_poly(other): + # return self.__class__([], [0, 0, 0]) + + import graph + g = graph.Graph([self, other]) + + polygon = g.intersection() + + inside = self._find_new_inside(polygon, [self, other]) + + return self.__class__(polygon, inside) + + @classmethod + def multi_intersection(cls, polygons, method='parallel'): + """ + Return a new `SphericalPolygon` that is the intersection of + all of the polygons in *polygons*. + + Parameters + ---------- + polygons : sequence of `SphericalPolygon` + + method : 'parallel' or 'serial', optional + Specifies the method that is used to perform the + intersections: + + - 'parallel' (default): A graph is built using all of + the polygons, and the intersection operation is computed on + the entire thing globally. + + - 'serial': The polygon is built in steps by adding one + polygon at a time and computing the intersection at + each step. + + This option is provided because one may be faster than the + other depending on context, but it primarily exposed for + testing reasons. Both modes should theoretically provide + equivalent results. + + Returns + ------- + polygon : `SphericalPolygon` object + """ + assert len(polygons) + for polygon in polygons: + assert isinstance(polygon, SphericalPolygon) + + # for i in range(len(polygons)): + # polyA = polygons[i] + # for j in range(i + 1, len(polygons)): + # polyB = polygons[j] + # if not polyA.intersects_poly(polyB): + # return cls([], [0, 0, 0]) + + import graph + + if method.lower() == 'parallel': + g = graph.Graph(polygons) + polygon = g.intersection() + inside = cls._find_new_inside(polygon, polygons) + return cls(polygon, inside) + elif method.lower() == 'serial': + result = copy(polygons[0]) + for polygon in polygons[1:]: + result = result.intersection(polygon) + return result + else: + raise ValueError("method must be 'parallel' or 'serial'") + + def overlap(self, other): + ur""" + Returns the fraction of *self* that is overlapped by *other*. + + Let *self* be *a* and *other* be *b*, then the overlap is + defined as: + + .. math:: + + \frac{S_a}{S_{a \cap b}} + + Parameters + ---------- + other : `SphericalPolygon` + + Returns + ------- + frac : float + The fraction of *self* that is overlapped by *other*. + """ + s1 = self.area() + intersection = self.intersection(other) + s2 = intersection.area() + return s2 / s1 + + def draw(self, m, **plot_args): + """ + Draws the polygon in a matplotlib.Basemap axes. + + Parameters + ---------- + m : Basemap axes object + + **plot_args : Any plot arguments to pass to basemap + """ + if not len(self._points): + return + if not len(plot_args): + plot_args = {color: 'blue'} + points = self._points + ra, dec = vector.vector_to_radec( + points[:, 0], points[:, 1], points[:, 2], + degrees=True) + for r0, d0, r1, d1 in zip(ra[0:-1], dec[0:-1], ra[1:], dec[1:]): + m.drawgreatcircle(r0, d0, r1, d1, **plot_args) + ra, dec = vector.vector_to_radec( + *self._inside, degrees=True) + x, y = m(ra, dec) + m.scatter(x, y, 1, **plot_args) + diff --git a/lib/skyline.py b/lib/skyline.py new file mode 100644 index 0000000..05fd279 --- /dev/null +++ b/lib/skyline.py @@ -0,0 +1,48 @@ +""" +skyline -- Manage outlines on the sky + +This module provides support for working with footprints on the sky. +Primary use case would use the following generalized steps:: + +1. Initialize SkyLine objects for each input image. This object would be the + union of all the input image's individual chips WCS footprints. +2. Determine overlap between all images. The determination would employ a + recursive operation to return the extended list of all overlap values + computed as [img1 vs [img2,img3,...,imgN],img2 vs [img3,...,imgN],...] +3. Select the pair with the largest overlap, or the pair which produces the + largest overlap with the first input image. This defines the initial + reference SkyLine object. +4. Perform some operation on the 2 images: for example, match sky in intersecting + regions, or aligning second image with the first (reference) image. +5. Update the second image, either apply the sky value or correct the WCS, then + generate a new SkyLine object for that image. +6. Create a new reference SkyLine object as the union of the initial reference + object and the newly updated SkyLine object. +7. Repeat Steps 2-6 for all remaining input images. + +This process will work reasonably fast as most operations are performed using +the SkyLine objects and WCS information solely, not image data itself. +""" +import pyfits + +import sphere + +class SkyLine(object): + def __init__(self,fname): + pass + + def add_image(self,skyline): + pass + + def compute_overlap(self, skyline): + pass + + def find_intersection(self, skyline): + pass + + def within(self, pos): + pass + + def create_wcs(self): + pass + diff --git a/lib/test/__init__.py b/lib/test/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/lib/test/__init__.py @@ -0,0 +1 @@ + diff --git a/lib/test/benchmarks.py b/lib/test/benchmarks.py new file mode 100644 index 0000000..9711d7b --- /dev/null +++ b/lib/test/benchmarks.py @@ -0,0 +1,39 @@ +import os +import sys +import time + +import numpy as np +from sphere import * +from test_util import * + +def point_in_poly_lots(): + poly1 = SphericalPolygon.from_wcs( + os.path.join(ROOT_DIR, '1904-66_TAN.fits'), 64, crval=[0, 87]) + poly2 = SphericalPolygon.from_wcs( + os.path.join(ROOT_DIR, '1904-66_TAN.fits'), 64, crval=[20, 89]) + poly3 = SphericalPolygon.from_wcs( + os.path.join(ROOT_DIR, '1904-66_TAN.fits'), 64, crval=[180, 89]) + + points = get_point_set(density=25) + + count = 0 + for point in points: + if poly1.contains_point(point) or poly2.contains_point(point) or \ + poly3.contains_point(point): + count += 1 + + assert count == 5 + assert poly1.intersects_poly(poly2) + assert not poly1.intersects_poly(poly3) + assert not poly2.intersects_poly(poly3) + +if __name__ == '__main__': + for benchmark in [point_in_poly_lots]: + t = time.time() + sys.stdout.write(benchmark.__name__) + sys.stdout.write('...') + sys.stdout.flush() + benchmark() + sys.stdout.write(' %.03fs\n' % (time.time() - t)) + sys.stdout.flush() + diff --git a/lib/test/data/1904-66_TAN.fits b/lib/test/data/1904-66_TAN.fits new file mode 100644 index 0000000..09c3929 --- /dev/null +++ b/lib/test/data/1904-66_TAN.fits @@ -0,0 +1,375 @@ +SIMPLE = T BITPIX = -32 / IEEE (big-endian) 32-bit floating point data NAXIS = 2 NAXIS1 = 192 NAXIS2 = 192 BUNIT = 'JY/BEAM ' CTYPE1 = 'RA---TAN' CRPIX1 = -2.680658087122E+02 CDELT1 = -6.666666666667E-02 CRVAL1 = 0.000000000000E+00 CTYPE2 = 'DEC--TAN' CRPIX2 = -5.630437201085E-01 CDELT2 = 6.666666666667E-02 CRVAL2 = -9.000000000000E+01 LONPOLE = 1.800000000000E+02 / Native longitude of celestial pole LATPOLE = -9.000000000000E+01 / Native latitude of celestial pole EQUINOX = 2.000000000000E+03 / Equinox of equatorial coordinates BMAJ = 2.399999936422E-01 / Beam major axis in degrees BMIN = 2.399999936422E-01 / Beam minor axis in degrees BPA = 0.000000000000E+00 / Beam position angle in degrees RESTFRQ = 1.420405750000E+09 / Line rest frequency, Hz HISTORY Parkes Multibeam continuum map HISTORY Formed on Mon 2004/02/09 01:23:37 GMT by "pksgridzilla" which was HISTORY compiled on Feb 9 2004 12:08:02 (local time) within HISTORY AIPS++ version 19.405.00 dated . HISTORY Polarization mode: A and B aggregated HISTORY Gridding parameters: HISTORY Method: WGTMED HISTORY Clip fraction: 0.000 HISTORY Tsys weighting: applied HISTORY Beam weight order: 1 HISTORY Beam FWHM: 14.4 arcmin HISTORY Beam normalization: applied HISTORY Smoothing kernel type: TOP-HAT HISTORY Kernel FWHM: 12.0 arcmin HISTORY Cutoff radius: 6.0 arcmin HISTORY Beam RSS cutoff: 0.0 HISTORY Input data sets: HISTORY 97-10-09_0356_193558-66_206a.sdfits HISTORY 97-10-12_0142_182123-66_193a.sdfits HISTORY 97-10-12_0151_182707-66_194a.sdfits HISTORY 97-10-12_0200_183252-66_195a.sdfits HISTORY 97-11-07_0510_183836-66_196a.sdfits HISTORY 97-11-07_0519_184420-66_197a.sdfits HISTORY 97-11-07_0528_185004-66_198a.sdfits HISTORY 97-11-07_0537_185548-66_199a.sdfits HISTORY 97-11-07_0546_190132-66_200a.sdfits HISTORY 97-11-07_0556_190717-66_201a.sdfits HISTORY 97-11-07_0645_191301-66_202a.sdfits HISTORY 97-11-07_0654_191845-66_203a.sdfits HISTORY 97-11-07_0703_192429-66_204a.sdfits HISTORY 97-11-07_0712_193013-66_205a.sdfits HISTORY 97-11-07_0724_194142-66_207a.sdfits HISTORY 97-11-18_0256_193815-66_206c.sdfits HISTORY 97-11-18_0306_194359-66_207c.sdfits HISTORY 97-11-19_0447_182341-66_193c.sdfits HISTORY 97-11-19_0456_182925-66_194c.sdfits HISTORY 97-11-19_0507_190350-66_200c.sdfits HISTORY 97-11-19_0516_190934-66_201c.sdfits HISTORY 97-11-19_0525_191519-66_202c.sdfits HISTORY 97-11-19_0534_192103-66_203c.sdfits HISTORY 97-11-19_0544_192647-66_204c.sdfits HISTORY 97-11-19_0553_193231-66_205c.sdfits HISTORY 97-11-19_0602_183509-66_195c.sdfits HISTORY 97-11-19_0612_184053-66_196c.sdfits HISTORY 97-11-19_0622_184638-66_197c.sdfits HISTORY 97-11-19_0631_185222-66_198c.sdfits HISTORY 97-11-19_0640_185806-66_199c.sdfits HISTORY 98-03-24_2107_193706-66_206b.sdfits HISTORY 98-03-24_2116_194251-66_207b.sdfits HISTORY 98-03-25_2020_190826-66_201b.sdfits HISTORY 98-03-25_2029_191410-66_202b.sdfits HISTORY 98-03-25_2038_191954-66_203b.sdfits HISTORY 98-03-25_2047_192538-66_204b.sdfits HISTORY 98-03-25_2056_193122-66_205b.sdfits HISTORY 98-03-26_2048_190459-66_200d.sdfits HISTORY 98-03-27_2034_191627-66_202d.sdfits HISTORY 98-03-27_2043_192212-66_203d.sdfits HISTORY 98-03-27_2052_192756-66_204d.sdfits HISTORY 98-03-27_2102_193340-66_205d.sdfits HISTORY 98-03-27_2111_193924-66_206d.sdfits HISTORY 98-03-27_2120_194508-66_207d.sdfits HISTORY 98-03-27_2130_191043-66_201d.sdfits HISTORY 98-05-10_2123_182232-66_193b.sdfits HISTORY 98-05-10_2133_182816-66_194b.sdfits HISTORY 98-05-10_2142_183400-66_195b.sdfits HISTORY 98-05-10_2151_183945-66_196b.sdfits HISTORY 98-05-10_2200_184529-66_197b.sdfits HISTORY 98-05-10_2209_185113-66_198b.sdfits HISTORY 98-05-10_2219_185657-66_199b.sdfits HISTORY 98-05-10_2228_190241-66_200b.sdfits HISTORY 98-05-13_2132_182450-66_193d.sdfits HISTORY 98-05-13_2151_183034-66_194d.sdfits HISTORY 98-05-13_2200_183618-66_195d.sdfits HISTORY 98-05-13_2210_184202-66_196d.sdfits HISTORY 98-05-13_2219_184746-66_197d.sdfits HISTORY 98-05-13_2228_185331-66_198d.sdfits HISTORY 98-05-13_2237_185915-66_199d.sdfits HISTORY 98-05-25_1711_182559-66_193e.sdfits HISTORY 98-05-25_1720_183143-66_194e.sdfits HISTORY 98-05-25_1729_183727-66_195e.sdfits HISTORY 98-05-25_1738_184311-66_196e.sdfits HISTORY 98-05-25_1747_184855-66_197e.sdfits HISTORY 98-05-25_1756_185439-66_198e.sdfits HISTORY 98-05-25_1806_190024-66_199e.sdfits HISTORY 98-05-25_1815_190608-66_200e.sdfits HISTORY 98-05-25_1824_191152-66_201e.sdfits HISTORY 98-05-25_1833_191736-66_202e.sdfits HISTORY 98-05-25_1842_192320-66_203e.sdfits HISTORY 98-05-25_1851_192905-66_204e.sdfits HISTORY 98-05-25_1901_193449-66_205e.sdfits HISTORY 98-05-25_1910_194033-66_206e.sdfits HISTORY 98-05-25_1919_194617-66_207e.sdfits HISTORY Original FITS filename "1904-66_TAN.continuum.fits". HISTORY Noise level of continuum map: 59 mJy (RMS) END s Pwľ'O$ `zģz@ +żF*p,E־ F/[ k>V詤Q"UXҨHp^K+F5;*nҽD`:HtCʍmaO<߽X<P=tY=8;B=yH=ѯf=>k>68½@ѽj1wXLϽ
-`n˔Yۋм߽g<wH=й=-<c ڽ5 a<<H ={ Ƚ`K _WojAC;r+p?;4sn<k<L=0L*=>μj<С< <#= +̂;5===u>/N>H>?y?.z)?,,>_>>W[$T>?1 +CfTv¼=[=.={<(;<
=Uf'=~=cMq3 O + l, c{7yt<2 o=AJ= =
>8=1N=y< <$=]==T> +5>5>GuC>W#?S?v?K>>,=Ӽw2rT|
<\uFGk=!5=zO=Jͺ=Vq=HK{|k+4~7:=`d= Q=.P= <=><=t;el>&t=@<ޠ=M=%>y> >o0>==g3==&=Gi<>j>l=Ħ=17}I<cIͽ*. |V z\sЍK=0=X`>= =g; *J=f~=m>}6>>9>p>??S?=> m=\t=' n䱽`A;^=S: =&i=/C<ܽo7m #c oL@9&ɽOV<ټdh=kb=zoi0{y6"=OJ= 9= =&&-<=<;r= =<Һm
<=؈=>Cv>X%>c
s>j> > x= =>!s>ɝ>r>Zko=-<O,<y ќ>& on75*$%y=N=n>=l<Z[ϲASQ- =e>=>90`>>>û-f=qS&P!ƽ kT弨ۚ< &=t='=S<eZ$£-)mviأ܈=VS=;1u=KR8(h: z Wqgս7G<Q&K_c ==@"=T ;=&R==P<T=Xh=p=4=QF=Q=;/==l>-΅>H>>b->:<= >>;>0)T>2>,j=V5<Ԭ`<'G1e1Y罖7ԽRAfv l=%XE>(>P=m=K $Ž{Qo佩ýo8+ =W_ <:AϽq{ F%ɽ?= 11]=!=>Λ>,=<H,G= AO!nԾOC<KJs-7{<Žc˽CuM}p⽚L< ¼]Añq,;Iyp_ k Cjkaw= =;>x=
= = 9=p6=
J +6l<s==ae=)=a==
|=מ<=>>>> =&H=v +=> w>E=<l &; L; =PbJl}V˽늽9=7><䗣XHą`=7=='A=> k==B酐𫽟"xoỼ3 : ʰTｘhҫ/ֽ ?J ս+:y=3=sV>p>T +D>,v> A=R9w;Jl<=]='i<?;7ˌuzD>%zv +ȽV9pɽQAU~4rUY>R~e}:9gzO϶D;"bwμ.=
W=ކ= =z=+=.<D̽Fi|ݽ = +jHO<y\ ̽GG弧 e(R8J2m=^q>Z>Je> m>n==Lt==X='k= +K +ٻ/==xϽUF7ae]>:====D?( F=^ =^<> +M^>I=<ϼ&R\Fw/C
ܾZƾ(IHF$WF`[6crG@kyAἢGvA7==h=ʙ= =y=lɻɻJ;r'۾ JQ
=}<=oɼASitzn0[pؽ>HS<n<=[Ѽ*Z@=>g>UN%>u>zA>a>Z?>b=> '2>&L>2^==M$>3+k%Z]<=7===놡=Y>! +t4ۼh= ="C;};⋯9wG߫,?z-cd=c<=!뼚\=> +[=){=>l=o:<_<-< <Za< ;.el:LhӖ^Y@$m<<Ӡ=D=9>Ic>vx>vc>3=V2=c_<,<&h{F<[<wP;悽P]p4P"ƽ\Tl<B-!<&7%#K!w>|Qn<u(Y5^>U&p: +/A=Y=3r=4 <<Nl1< =:*=r==sc=YG5;Z\1W .QͲjG ffHS],CNGfn'þDA!iE9 +~.m˜<=.= H=y=8;0<=v=9=@=5Y= =ϛz>x֯> > >>b->T>Q>N>A6 === +N=Gi=Y=2=hD=2<C=HJ=\u=<w< =dR3<R=U =#=y=\=CU> :==&3Rp>i?#z>h>U>;g>"ql=П=0= +=Ze9=>M?>Yn>y.>s2>Z>"]v=؊(<}弇:_ۮH$u _;Fcͼe\`F R{p<ѭ'=y=]> |>e>I>=ƌ; ׀Dz / AyƋ +ڽn]rS9Ge[ѽ\~Erټ61弒*T ~:H# Cv ƾT"ݼ7Ǣ- +$=L7=<66<W<A
==)1===!<6;kL:%YaF- +тx㽽e +Gݽ#b'ΑP꽃u7뢽N{ƽ9ý#!,rǽ(, .ν'@f'"Q| +(ҽjKl)XKz8vD⽋%3u^pg;WXusykKU;=O=B> 9=M==Fy<g<`=1=)v=&,0:Zo 콡[9y }BHuo+6<FdOO澽#` >*9 z7k@˽Umýq%ݽ/:(#
T; rz;N<rS<=(===~\=79qE3n҆c-:Ǘͽ䧽
쎽G$pYtڽjT:(;~<= ^=w=$n<iS8O +")\#õ 9Q=!>&ڸ>,>t>>?+? 0>> S>k>>>9`>)>
%E>0 +w
ė½ޓݽ@~[G0wd:FgmCU8ǽtW0{:a<x=5= =4*p=ݷ=o%==o=V=7>Js żȼħ}5`ܽ|$.^4ǽJxk}d>(*)㟽 QC7f7*) "e 1
<+;:`<z<=$Z=h{<#<ѕx Q$pW>H_EϽٽ;[1 y4ܽ:-kvlN*z <6<V]$<s=/ލ=`X<;/V;Z!?@JoM'֯jڻ3<<=3`>>@,>0>" >To@>c\>'=1 =^= +>%I>>cC>c> ;=͐ +=n;=ཟ=4=#== K;EK=#'{<H2{^27H"ֽg [ GG˪}N<CRɾ}|a=uX>3>Ll\=<ϻ>,t>6V>2[>1v>1=@f=aA<(=;><>v}>@0>A +놽Ŗ6˅ƎhɃн +M=$>;>_G>r>? 91?Vl>_>[>y== +ż Է &<VOWDx9)3ҽ%5M; <] <üaV"4>)]> 4>>aO>.=/Pq=>N>-t><={>==C<j=n?==}$=qI=)= ټ(@(iŽ/Hʽ#.ְK9T3eklu9oMiWf@ɽж2XP7.}:-~Tǻ<GQ<0=j<tGux5u.OT; ;]*4`S,O/ڽ+ ]{J(D] tq# 2L'35,g mcTw 3fL<M2= =&%=*< <-h/ge|5:Ƚ/㽜NkλC:<3:ǻu&VLYkܯxbh8ChjٽD̞
<;`ٺ<!#91N<j:$Oe<'-=җ=n>3nj><>)? Y@ 5@x @@U)@/?a>ٰ>\=!;>D0> G=y<m<d9=/<SoG<;&.;<Nw4<yV=N<w= W;*k=YȻ_Ύ<I<= f=,<u_l=}m=G|=>=+_=b=jE=1=>b>2`}>2b2> O=/>= d=Q=m=\
=^==;F5<=P&
<=E~=u==);OգmYƽ܇~!?aýS4l +?p6ս֑ }
jdEYyӱ:xͼ'S8}=0*Lj-<I<oaU;!{x Sd2 ʜL ;k; `<0O;H(:ߵ<!<Ѧw=*~=d=pD<>ƻXB! /#Xs hVϽ
gX̼{ȼ +ͻ =4:= %=E=vO<}(?Hd0˽x0^5<@>z;B<= #*==e;dϽ扽!cڼ n< Tƽ_ {h<9mO;қ<<TM<L <g<1=,=Ż>=$>k>?~@?r@^,@Ͱ@_@' @ +@6W?/>>&
>Hi>J>19=Q=|AW= &=z/<;yۼUmi2gټdx;Yo_<ym +H=(O=\r=ki>9T=c=;m>"|~<=<*<R;;e< #һU ؼWO廓#<$<<x;:B_R\yBQ?,/kuH?ǘʽg~>
X!VxJ. y'鵽p@mּ⽊V7JNnOһＭ5/xd<485`%;RdB<1=n +>^R>?T@ Y(@!A +t5ht"V:DR:)8 &7S `> <>O=R=$R=P]=9;z==S=K=JE<F;: P! +GLlAV,dPB／Ni¼#5(:$KAF +pe"<#=2Y=(Q=5^=<==@=>{>m?x?T?1?\*>x=jM4C˺ٽVYϽ}ltos,꽞2Dۦ<Ƽ*q\ dm\ofyK:;<g=o=[+= <&=T=->n9>>Ņ>f=bw<ٽLGEK轀NF,Iz>g8;rX=4 =Ķ>"">j>(?yVz@'@@͞AfA-q@~W@ز?`>f>b>, =ya>
s>*> t>/ܳ>*=]:nĽr૽ھ u߾ 2BX=>p/'eJP "پc,^>L CJؽ[`Rtm ^<n=#=Z2=?!=*=jD=SEsi+ 9. %Bݽ5#L 5#Fw. :G=^e>I<NL<(N<@<.<f}*ul 7$QG +;(meWD;ܼt_mh ZRS??oF 3+_ eW0ڲ༼<>=9=́=8= =l=ņ> >?p??? p?q?"K`>W9]˼A+ƽˍ̝XzI T䢽 Aʽ+@Ҡp}>6qNAPK1l#rażϠ<<)<O<'=o>UO>֙>I>>>I?>;fB E+q)"3;,+̽L'b?7t;I<[%=K>m>?.?d@<@ckM@@Q@X@<@E?|>=P="ʼ;g-=[Y>5>C>.>N> *< 彅UѧF'.<[KF+9p6 uƆr_Aa<|+y6v^dn3~{zRg=P:=i=)=3=ȹ=
p(;Y
eRyyry竽@QX <a=EN=OV=2===A|꼮X[l!l1<=]=̬=WZM<[20ĽK*潔65ڷ|ޠ4>;Z369[1DCgw96lf~ 6v<B=z= =ev=3=[ ==>>Uh?n$?3@N@?X?a;>."wLSG G W +<ɒ6!PqٽĽL +>
>+?J?X?Q?)}>?Q>?p=1< .+1Ent9 ý\7!h[q29a=^&>ׄ?d ??R@'@/A@O}@@@Q?>T=Zl=18!9ý(\1<<=sT> 4N>:&> =`O &.ăwdlǠF9I3 7NqоR1v9vm +W"h? -*+Mǽ(WJU4>o>0> o<dഽL?5 ~;9zFͼ<=L=LR=6=O1==T=]I='d}R;&u=>)> ӌ=Ja< 7̽7>Fo;P<{G<(ҽJ;GVF +SϼA?wyir;i\Y%Ii< =F7=;A< <O=a=>s?[d? ??̟?? dY>qû]佾 #})k DʽҨ3>:. Ui 9vmI ὧzu-=h>X%>?&?ch?Zy?U>>#s=Ѿ=<*11ؽ⽡}彅k@G +Dc; }>~fO?4?ǒ(?E@@,@ Z?Z?ːr???Z>}' K +i`.w\ۧ A4 8y +Rl]={6>=?)?V?+?B>}>6.><rTlWνH@Yײ'VDC`ڽD뽔 |e^R:Oψ I̽nھоt"Ӿ S W`Zf|<T =l>R +Fj
S0͟ = ݫ}oֽ6%-,s}e
sQ<4J#=>]>!? t?5? +۹>>>>N=h=GSսhǽg5,SjL<kB>#?.V?)~@܁@@
??U?X> >'=8VTԼ|,ܗ½.o'Mft ;d<ח</&+4 +d<i<Ļߒ0rļx0r н~4&ܾ }F]ʼ8C'6-Z^p瀼k[>8T{==(V=k=3=->u=s>&=^= .< =9Ι= ɺN_/N¼s%B#<G4<t*< <3=' =#%<6μ{ш)ޛ:^<[===S= +1ѽϼ5<O'<
= +x;6żFпs?֟%߽?|1F\{0I7^>|˾Ӿ +,)aUE; =l%>>>^F>>Z>6L><ۼHм{nq GG +JC&}3]4ʽ >½S?cf1 )k}fjNԼf{!}<{=?>l'>@>&>Z>a->|>>>Lp>9| ߺr2~ؽHf"$ Ƃ`< >mz? l)?GG?e ??4%>=`=.x +q)սmG]{gH+\i@/ :z;=,<7<No8;?\=[Tgcr
dyx_,*& 2iAYs ۼ6Y>C8b> >">-O==>=> $>tE>XR>SH>#>s==^,E=qK=n=;C:@;ť<Zhb="Lz=X=Pk=N ={=b;ҽ \e-%w<]V<¼n<>x;7 S3x۠B!pNɮs1\2guTսlRJRk؈}t8F@ &ƅ52=@+>>e>>?>[9>%*=f;gm^ԽF|;b;LD( " t*: +B~d gZp|h+:;;.r;I;;+ͼjF?8<=N$>&{>Y>}>2L>>>U>r>Yz>5H=x<m=7 +Xٽ +41o'}sͽkS]QvT=>,>$?
>3>~> =W^)?3k*<ˏ~ھ&;߽Ű۽4*Ľ;/% XؼC <jT;B<H8Hfh5hKp /y +><GP# a' <FK<<fO ⼪nFF 8vu_3.l⼘k"<1==`=oT=@s8=<<m=& =>'>, >+>#<>>3R>,ݸ=R4=ʭ<Q
Fzf
좽Cy[~_nY潿(F
no'=s_>#>`= B&zlͽHȯJ; +`= =S4=j<u ;t/gls:FUخƼoOZ}JKEu;MT-kGbc Wz>nk%n[% Խp +i
<R =9N==@:rļ dp̽R^1I +.g½Eb= J$:g:?},?[?<?EA>4>E>>}F>7a>F>"~==Q =<}5<<w=B<A=;;ف=6E=8=kGH==*+U<I<k4J=<:\` +<;>;]7>f~hhW/ZI
R 伿gQ $Ӽ0!˟ +ȼB(dt턼26.lV/Q?YֽE<|ֽ!RB]H+<D =1>=Ԙ=:;KTȼ~ݽ C)^Vn*{-g=7t +=d>(>0> =`=j=-"=mo=9=Jwe==Ax=>1>8>"=q>t=hq=#z<\
=MZ<u;1AwU'+˼{5dk/&;_T;%Z<y;;Sν D^~b{cRнc_GVk6{O5ȁ < =4 <m<LZug$2żQl~=g= > O=;*~ EǶ뽥!LϾ1 @u|P)l!<=p?y ?@W> 9>"f>ط>p%>%=C=+==o=m2<<U<4<<k=4>=x +=XS_=`=A3=?= +6=¼>#♼=7\;:t];<%E<{cw" _\ɹϭ_'c;[EOP.8>a.̷MCHEˉ~X +f +λ2<]<4<|;s_a<7Z½zνm[DApR{~⼢)]νg@r+7z< +$Z=@ +< .;; ;= <J <B6 `FIG3!4 dAQϼ' y1 +O㲽WA""216 n٭k5=M<->}> >};)>y>U.>B"><?=7={:=4s<E4'⚼~ey'&;<=8< 86=Bu>uJZSUܼ/W!|<R%=@P=:<D;{7Pp=gZ/9.JhHdzؽ\ǴƻI~4ػ,;?byzn>6)ǽG)gM/ ּ|DihӼq;=7<<;jv f;]4佑11sS>м:=W>>uX>f>*J>#>>V +>2֥>!w>8>~>/>
> >?>` =Ѹ=;;][<z:̼TV}2U(
K@9bcY+ս.vTK4ᢻx;u;I /y"ƽ Ƚ>($$ʽ 8) ؽ.`9/; zbzrCyO:@{T.5]1D;2^ܼQ:sm;K;%:|8;<=NNo=ANY<tťO⽴' <{ĽZӷOXo=t<A=>@w>:XF=,>MrP>CF->x>2d=tn=T=C@="HNA 85/M{QoDj*onV *.<K`<١=)<t7RF 1Xi4QO_-8mO<=R(D=J==V;ּB +ܼGph_fX!pPa埐]땼oɻyuT¼ ޭHDԩx)t}EMۼqFIQ`_ռL`;<uemӽBitW+GloV=5=q[>4\>w>j>>1>{X>D>1>YR>U> + %t;tEk<ɲ<#=U<$< =EO q+q)uQ<[ +<kl(<]<"=`ޫ==[>Ŝ>ف==fSݸ +==]̼1 +[TA*3y7*-N<H9<=5>ej> y=<z¼ĵpi n=ƛ2{սo6Bń<=T=($=&8<9Ƽ>x: 'XhYCy ü"-2n
/Q̪(KԛRpS=(oiݿiDr& ; 9ۼǃUZʽU9 3°q;ͼyM;Ƹ<= =Q>!n>X>DG>>u >MS>"h>R8>) +>Q>??>~n>]=t0<}Fz49Dշ#w;@b㽼FyQ彤Ľ n_"Yս ms
㨽y^;!D< #<];WJ4m7輧μz11;s;]-t$6R.Pڽ7=v ( ս <N5.IѽOp)G/3;J֚<- +>tG> H>L>>U>;=z{;J/ D2Ahü3½IUݽ}Ž@?/ +vm<La=7=nF*=Evt<=*Q<ۡu_#xq"s)[8T @֗])q}Ijuj\n;7ԻiK<P<\<Y&ʻ%y
<AF<֠<oFJR1ܩCʽ:m*]}c揽ݽVm%ֽ7bY=
>Ć??? +>Y=C~νh&D6 +` +==Q<;a81n +MAXG0'-;d==0b>^*>bdy>c=I=q`κ ڼ<#JG=17=ͫ>8|>v>Q>2o=y====.==yf<qW8۽ !\Y0.]sp?ݼ t_SNc:WY.hKg2ud3нý@ D !;D[<l<l<P<=ڑ;Ԗp;el<D3<^d<ghM;<8=e=8x===?6;|ϼpW-I\}kE7҂1 6R ZpyQ;ze=F-мGgr=ds= ;σ;K:2z= +n^̐FCYC ;=0y=s(=4f=JQN<>;K:Ó I⼞'Ͻ.ݽ0,naoW=>z>>h?B>B>A+=m.G;Rt- +R;m=>>B=<w5><3;я:=.=y³==z<w<Op< oe1uLS +R8r!nooPI⽑%;ǽB^ BDH ಼ 45<>h<<D<;ld#<F <;7@<<r=$J< +(<:#=B'=t`=Q>=C=3<Y H;lнfS=Buk<g*R em]mFק}6Hu+y(5 +M뷿;<|!=#o?=M=&<⹙+<0<B=)ө<io<;sɻ~`q E<=y=< <:wglAdDJt뽏
½6>(~? ;?v??-P?>>T =u`T·ʽlr#2j: t 2 `;=y=K =RV=a=p(<GȽH&_ܳ +{j::;T<?= D=+=9S=eD=Sb=c=?=s=|<^;ٽ?c;ZQ;Ov<<FɼuU:B~mJٰ{
F$V[ؽ$Ԍ8<">N$?Q%?G?ֶ?ћ?v>?= 6 ukнٽٽϽp:W+}`=<>=> >%>>=[R'GY |)< =5n=qC=jԸT;;go9۹ʽ蔧tr(EY4r{5］4~f<%<i~u<k;_;H;n<pOv= |=;T=l)=1AF=<Cf=O=PH=3;3@k@t' 1,\`gqqj8ܭ 틽S߈OO87µ&$=<^<='-=&=.=Z=2f=P5'Uoj=<8~,0-)$ 8ܾ | Hž/9&},"6ts(w7սOm +^ M=Ͻ|ӽPySmB=" ¼*; G; ; < g<g=8O='=U=5====E h=<G</< <n<J;}zU #ƼK9A8)V<9eS=eb<+ڼQ8ȼglQ +#8]T3J)KY6BS2ozsb8, )C0sjڼEź[)j<{=6e=b=?֣ }<TJsٽc$CqIWgG +;=|I=o=B==Ҫ=GI=!= +<Uo<-ltR<ROB,b##'"k;CF<<UW<)=K=Ah= +='=K=,V:Ps˻̻`h1'JR. +Rν"$lX>I\^q嶼5_$t<<̜d<L/=2; =U==,=^,<Qe *Iٽ]9iAPC<<B=<j黟hv>+<֡i
KsC +c>"$>ua>> +J>U`>\P=?=H<|8<=="r<(x{l"ׯ#;y;w!<Im=
<c+Vؼ!i
bB; uuJ6G/{ڃӽ6-4YDBneo}Ԗ>u=Ӻ=/>->a/>_9>"=$v{<h=Ė>(>$> +BidD5;(<0*<Q;ͼ_ƼѽBd$#PK ˼۔f~Hb<FѺC۹@m;d-t`A +(͟C:QܽSͽ %I{< =;\hg)'=Dp |~[7`<:3<<s<ypƼ箼SՏ#>#T>r=ʯ=<8:M|яSu
ӽF[1ν 牽FȔr"ҽKHh @ɬP-=z%=>>C>*>
έ>Q=u=W=8a=%=hA=y=sv==^=5(= =w=H%=9===)>
=>c>M5>]1>|>4>kV >1g=j~ <BJ<.G<J<<\ǺLϼ+伄ü i:.si;A==d<eK#L̼{و?PӽWve܆=-EMBCѻbGoi+сcxǽI$<u> }> >>>y>0=H;4< MH< =B_> +\._dY4155<= I=U8=b==۩<'g +(<7,#Lhz,0__=5ڽ~hNxM5%*{;!$-Bg|k2:_=d= >s> ?>rD>va=7=4==n=|y>,>vJb>j>G=+<݄d2 +h>@%=(=s=A}C>a</ةLs<`$QJ,yhFs'C +ݽWp1/Y_ӽM*R<'h<=n +9R98H^J2$Ӽ>Z.5Oͽsܽq\˽,|rp!@:PFh:P;<}[b +R>\ >Vׂ>- +=E<~(gLZlnO2]&sv5s{E8Q齆 ;H ͼ7@P芼k^SIGT\d9UJmD!Q8wݼ ƽ0cw̲(";:;
5:ț;0<=8\= 1<Nh<=4y=W=>,g>h&8>zZ>Sl> 4!=1V=<7<6y˼d<rNǜv۽8b=7L=jbѪq< <*= ==v(<YVD/溿Zh<=4=G=Q5=*<= +Z@=O=r>P>d>{S?> >>j= ===Β==
=E= +%<ɎU= ^=2h=e.=X=E= <;tQQA:;<E<Ӓ_<;ǻ]}N!9j }V(S Y2x ki~bFr k9-KF0}@|9<VH ;\T<N<< <%1^)I驱mй =m>; >
>>>Ey= =@R< <r(=!F= n>e>S==<raݶ@tv! ;"MfKRSu̽aCHi꽚Z +_l v/9½9E%[S7kSM ;kkSCʗ;U䭝*mBP_;i<<X<a<hKQ< <=b=Sb> >b>u>b)e> >!Kp=>=SՇ NցƼ3ýuvt̽4mqCӽ @ڼ<<*<[<<a< +<I̻7 jzڼ&M< +=)=W}=3*=Q6:ļ < =q==3=a:T<g<E3}VθYFv|OY>0EݽAc< +R;#:
-6.lg3:zJ[X N;L;3Y|= eC=+W=> >ET3>?8> >9==+g얽6 ?0ɓ5Rὔ罖D%!½HO{A9h= +=b=HI<ͪ/<nD9<.+$
dE @N^pt <b +'<=v=Wn=g\>\??N?Wg?: + D¼ >7=Ub= > +L> +>#=|<<\*P;7VF橼:39l8!27VD~L] @~PNfY",ڽ|'#kVؽϽ?[P;M@@1Vݽ^~O?ڽI/f'Cn{;ڽ +9 +|q(Co +T=\#TS^6HxM^;= >2%>=^>=5>9R >
|t=:*==~=~E<Y ?#$?Lb95M 0żHe_w:qd;{$V<}b<H#.{FvIUw}Ժ=i +۽ ے;_< 9==[=y4=p>,->h{>uk>U> Z=fv:;e)U:\KQѦM(|wUk8l,e<<<9T; ;VZHqL֪;UctȔJp +.UYҽ ==>==<=\w;琽5w
n
c >&sD^i<oTPo2];~P<YU< <ե<橫=ס=`2=!.<g=X(=u= D<*˼ *:ޥ]h̼＃>ּܔ*racEFzTJȞڽZ彷ƥ +<!C<Ȗ'</z;6Xq35ټ:N
!!`^ +, =Mؤi<P +mLW`qUԷ> G= =:=>Ў=L <j= = ===m:=R=|Q=/<<}@ E-O8B1<A,;Jq<*E<67<o<|=l= =6=!G5<6)< ;ݼ`( r~һ蓻c_$u8o<~L=>4Z>gگ>n!>#= r~K#a!!M6)G,K[kT>5
S<:<z<㯦<翤<z<></{ +=^="=> a=.= =ݩ=S<jF<M/SǿQ齂zBQ b
wFF%0zNt0k~u=7n=w=5= +V== +O%s;B$=5{= +S,yܔcH;R<=gK=&= q@F鼊¼r=u}<4̡<!<el< :!lӒj P<(<U==F=T~=X==v=K> + >#}">,G>==<!٪;fwJ Z9o:E憼 I. gLѽ(@e{w)&
![Wʽ">g۽Ru=A~<f=WΈ> ̏>SD>E<>94>=s>> >v>p===&[<^<~<N=,;^{Cx[ԽP34+`;<Mh <f{u<<M<3E
WƼL'. 4ݠ &i"3*a2]= =ԩ=x=j= ==9==黠==P==]==v=<E:;WVmX <L= += +C= =T==ż==B=>y<Z59̻ ;<T1=i=q:=r=C =%; $<(1<E =
G`==='=~;WG ͽj2NAj`zw +;T=9s=R=NH=S=2;P4p2FO
%Q"巼*yB[*U;%h ;9Vvgk꼳Pj ~4Gs <%= L=>g5<CAy;Hn<+ۥ<1=s=) =0[=a==m=J=fn=c=+=a=
mc<"|;p}&Rmω|`)<U/xX }i$6ȼ6H +D3ŽLEO3s;p5\۽s;a=\P> &>O>[.>Oi>>`>v=ܳ=(^=a=Z:8^Zؼ \6' -`vY7ⲽzVӼ$q;Y=$"=UM="==bʒ<y<<O#(;ʍ<F ;;h㢒p=$==T=1>t=ذ==J=<J>ϧ>V>= D=<&<Z:<(L== =Q==g(7=W=G=t{=U=I=_<`+<I'='=_>ߓ>E|>1~>=K6<3&>=< =C!j=]=8? m4ڵqB +JYL( ֽ&⽇AԼ%<= +Q<mлi}y } +6h 1bÕ3;Vl2=P=9=i8=6>7=§=_=3L==G=Z=J]= E=kx fٶU6n#<==_= +<)Լ
9Ƽ-x ZZҼ& ?缣=NqKu%fý gܽ=E$Fͼ.80Tz; +;4 <Fz y;%[<fU_<,<z5<AM;O<Z <c͠pPh2v߹աjp#XD :pK}#ͽqf|Jf+bkj>y $U =Ak u½ iT ;̻0z<=h=R=c=ƶ=J==h=(><`^j24B
?M\5+J\ʼBkݽ0% a3%伊ٺ +s=r>> \> b=.&=O|==== Z>!jy>*`齛&<S<G= <=d y== ===>qo><> +"=K%:\Xc 50Z<u2Ĩ( +k銽k"Dq e3FlȽ$OMgiL&=<====8=;=@Ϝ=; ;;)<T<5<Ssl;"㻢 +0hμ@xѩfϼ'v0 )9Cp Q;;i6`ټEJƼKΙ<Z^<{ <;E;4]"~] ؝+j?Ua7MuOhh`.ɽwRoo4ܽUA?P"=QC̽C鞼g８; RҼ ü5＃Iu䕼绉<5 <B =3$8= +\;gvW&L?Ɔ͙< X A.r%arӼt+FW5-XicX<<\=J1=#0>@> *=O&=l^=h=Jy=J}=7= [>P=
<e1f<<<=8 =k*/=@=b=Z=ni= +%=5s=Q\=e +=P==e;寧V);g<k<D<^<<7~<_4<<V<u<+_C;N;IX;X < >(>Ma2>~=>>P>
u=;YIC26ƽ(= +O\<ʟ̄ƻsmQt*n< +9콳b\ǔd$1;ɽPI4½罶t큽rXϋڽný<j=A>R>,>>i>>n֥> <=Ń=xٹ= =Ys=w= e="=XK<;Tȼzr69&_軁!E `MO½% ռj.:i8;*;ϼ;!;< =Cn=^W=OK<gA<Կ;ڃFgz{Y5{) 뙽eؽ`9Ǽ?>_] +<b<G;Ǽ@Вpdm48 ʼ + %T|"E;
jS&I4q\u_ulSI I8=R> > +<X<b<8*N<_=`==<2= +97<Ky%w( B/N
]o<<<ງ @ +s96 +7>c>/ =^;==/t`[8krh`-'мKD*o9mx)`<<d< Q<h<j<j:< <=JM==zI>M*=F=(<< ڛ<!b<:A Gdo5;QIetR켆*<:ǹ*2:<<H=3=>>>;=26<Kj༡*< <<<¼|tHN}KH:ùo: o;N;A;W< +$.5\GuC=dU弰i#'üy`aAټWII%dnex +_HAk=qLCh=qYV/QۘLQl{i2s G1̔uʛY!罠ѻY=>@>#>>>0>=ﲕ<)>1̑ ¼b1F)>ha⽈;sc&\ϼ _:C;1<<@<<s;Y.9:jcvī:Źv<d<e=# =A=w7==,=;;̧7;;Ý;Q%8k,&:_^9?g;զ;<%j:ܛ<:y;;լiz<=y>[>9w=e<(PfT|Q) & ^A42R4 ]Վ r?Lļ1&[KiJ'`(rFY3/FJ_zܽ\ jT<lHS~_3݈)4-,V@2vQ;C[ݼ #ŻEĭ<<<\c^A㖼ү;߽X F1'_¼μ,ۢ;$3:<XG< </9.]Zm<ҽqi ?l8*%ƽ
Q,@"Q;|`\:Ė
aď ֩
q]#sͽqнI۽P@ Hi
mn5)
A>] [5ҽhսh> +>;>AEf>5=Bm<>ERcC٢--: ~V2( 팽N +g;N%.q* Q8~J'3(ڼW:_I5s :Q=;^fO1 yyý r-rDqkL|K<.Ǡ<Z<<;,r_hm`)",~Zx: +L 57'Eؽš||fU~>⽡9ϽQ2l'н3Z F +=zzr=<Aϼ^;-;*j Ž@! %C<@Ľd=cS#>fz9>*><ԕ>]=6=Hr=(=<=L<,<o =';UμP<SH<w=HN={<փm }~Җ^J{zgd}>I "v G#7࠻O;;H,^H8
:lٻ<q<=ʽ<yw:#\j:쐾M/v !J;
= =>&>h>O>}@= =b=h==x.XԼt+q)~XkO;0heI a + DԽXps1 +[幷 *_?c\Dt==>1HF>>>g>R}>W<KD ļʻ_<=$ >:p>TQ>l>j>e >=U =0eꩡCK*»oǻ T6Aqj0xeh sC +>Zu> >Z=G==_={ =B<<< lzk ټ⛽1u１u<N z=%"齴X^%o i 0 iG齂UaQ +=3<=k= ==<=,v=N<<A3 o!
н odǼF@]D
<pe?="_=(>'٪>Iō>Z>;ۿ>g=R=SS;]*4.A&"wmhO (@r UPq=MRj~+^&9:92;1<*=w> g>=>g >{UH>1G =n<;Z$Z7ҽMBf;i-:@#ߴ +<M;=ò>>2? F?Z?.>b>X= C1<< +Lik W,T;(;<{j'<$Q༟9μm6ؼż xȼgҼdkZ(+}v{" <<Ȳ=n +=ʞX>J>`>s~>>v6><01==t=V0=1<Ƶ:*[j:!l^~Ԣ;Z;Rb<>'d>o?<}?h[?"??:>=W;_|3
<N=t>2>OZ>a>,>ik>3j>
=:==5<><=7$2==T^U<N-#ؽZHsڽ`YXd>d=ʢ)=~==*=#[<<!;8 t@髛g"'fN%_i:-,x@`GNZ@zؽ\
?;+=\>N>Wa>r4>s>: > >N>*
>r=LU=k@=Bˏ<s<Cn9b`kJp驼,Ѽ %E<5S=S=y>44>==$'= ;#O5]L:o<<C<|8LIh [Q 7<x,;':q[<</=>`?$!?P?`?}?2>e=Zn= <P<;< +,+<N]<=!<J< żzwz{D{?pVxX0dhB$ I"L[;'<V=4=D=;>Rn>=և>x>Ю>B,>=GZ=.1=og3=9R#<<<@:e%<,W<l <<?PF< vL<k>}>y?!?u? ?ږ?:,>(>D ;6lT;&^ļo<E.=>/2>8>=V>====M@\;˭*@<0,=)=w={=E}<X^'5yl+WN=ݿ=$=*!<<3-v
Y>A>]ҽQ䨻
ϻ<:@=/ ==Ͳn=°; 2fP2h
#^B=-%>9C>{>>~>>>K>2>A +>N4]>YӜ>7==AJ+<<<fx:79i/uuX:PCۼM< y <ğ==T;F== =g< ]r]OyR;<x<<; b,Ի55 ; < <=U;nu/<LpW<]K=i>>U?o3\?b?X?]?c>[>7&=Q<Ƕ<> *ѻ;>S;P< +c ==>\MO>.?)?;K?1 +>?>j=:=阽H,xὧ>Cyif<Is= =Z+=Ee=_;=O}=i<<WR<2;KϼZ<r]=Y=z^>mE==6N<w'Br2 >p
=!|<a;2;7DcW4Dpnq >?di@PC r՟3<8=f>%$>,>=|C J`Y}x6# =4>XP>0H>>T>X>f>Z>V&1>]Ŝ>P>![==7<|^ +> >!``<\<-)<<`jd;Ph<D<U +;~J;(1(jZݐ} r
k+Źbӽڗy<%]<$<<;< <[=R=A7= ==|\=ڼ=R}=o==~Z=6*<<@<% +< w=2x=%F<]E< {=.&> |>[>8>G> f=_ ;[rN۽~$v}&4yI<
<Lv<^|J;vɼ1o:҂;x;==R!=>3yk>BO>'$=|x<&; (J)~ۅǽ]=xkP=_<Ϲ<Mb<9C轑Ľg59}˽Z͕3;=" K=>S">B > 2U= +! glR6>;9=Ts> D>{>\>z>>9>zG>G >7AE>5=0=a=!z;s_ ȦT7eWCD|5+4Bc
4QJ-HQ)災}q-n\<Y`<`5<a<U";;^:ٹ<G==:L=Ic=g=M<!<g<g|< ; Z<W /<C?w<<e:=D>H>_??!IY? >>=y< E<;4<5\0❼' :(~ļ0 W} 轖Iet v#59 &I[<G=21=5<<,';<Nd֭<)<@$=%=W-=J= W=ι=p==x= = ,<= =oCV=w =L\=@]ؽ)KH;= =%1==d=G<Cɺ_ *ogY5.vs7bsIm >[>#}6<<S<= =@C=9>ax>j>op>I +ʾv8M|=Wh=>R> ;>) )<!;q Neв2̻0sZ]n4Ǎ3,+i<߫=x>
\B>=פ=Z). ؼ=ʼ?XN=S=>3=>oFO>5>:>z7>N><ј=R=5<;%<ũ<'9;:`hM^ 8mt9[% OY$k9`o$<=M=]=^L=N=v 6=~M=f= A?=5 ++=X[=\Cf<%:ʼn1 +^XĽ1 +kAteI;<f<=UN=Nl= N==p=0ȴb. 7j»gFYx3̐ E=Vug&cTT8~UN~gu=Ul4)=/;=~=\K=o;DOⰼ,.W_ZW/:<n<3</@=MJ=8|=%=u=;=i=3= +3<^ +;< <`;J! <<#6;5{䘄:;52B;XcT>!;AdҽgDF +F))~ʃ輟T8:iS幝Ch6<x=O_'=,ϼR ~G$lݮ
7[Yu}j;';ұ<=)= =Q=;=_K='<<n[R<<A~I伪 += k =KQ#=$=V=]-<hqWߐU7oCq~_I&mr;J< =! ={p=ɀr=R =݆=v=15<}y?* Dq<=
;=>=!{<N]BG?]t]7hٽt ̽LTo&rJ>ӽ8ξ\<8=.|2=w=[==L.=l<;;YKѼ}wC+=(+ԅ;λ# 9; ;<<$; qr:|;n<9ݷ<zXM<3f +]ż_j;gT9@~ < 6%;T:$ڍãN2hNZ:?fM-/4R' gϽ9é)ΛŘ #C158 4+<!=#=K=xУ=~C=L<<V8_< xrےwka}۽PyZAd +|NV& G:s]џ: )govp 5*\μ睬`nvq:w$J4<+)<A<<X< <~84笽4x v4==â=~YF=Hٕ= ?]<==^=8뻖r y̳W 4Q(R*?l +~wi; +I<k=A=?oZ={=J= <9"dȻT+鼑:˽P˛DH8v<i%<+;K +q`. | Tcg;h<ﻉ?R :ArMF;eV;g<<<+<N<I;T9B @h%ݶ bPP'ǽ.н Qbbܼ
<fO<;t5;U:[ViD6q +pkսa1_p̽B~y=AKw ԼBۼ樽Լm4B74<]o<Xo=F=t#=u=Z<<1<<A<v i$'X iz,,/T;v-g%?r)p kp"7kۑ FBH/"ۘ½ AʡН(:;<<%< <<]<<Uy퉽ʽĻ`=P=@=y=ĥ;d=I|===9unǐ(<< *8]J Yoz: 2=%]=4 =/<B T<u!<&< ;l; O6G&<GɽR}P ) B%F6 JsEJճl|;G= +*ʽrT)Jh-s{Ӈ0弼vWPlr߽ c0K/] +g +8 5饼}MHi]զ)ϼGpD;R<=E =V3=kC<<o.;mt ( _Խ=7v;<=\s" Kg3p=KqҢ C#*9<O<ϼ +̄qz+fZvb;< @Zzɜ$G9; _:SXjoF>?6`r +7qP3!Q +%=Me=:=8^>X=i=e<C D +ڼϸ؋ l? ༳8jSƼTfM˼ 53F `3M%뼎UyzV9\<E=5==o<:1ꩽǽC & =a20=!=1<L0^wc(/rR i 8{aI:~K;-6@~iKKӟxw;*%o)% +MkJDJ k "Qc]`뻭3 ;P_(O颽<: w su]n;/<K<+S<<=j=@>f(>0> rc=g= +繽༸ռb;;6b<ш%<<=i5=j=e=a=<}y$crKÆ\$h8亽RDj$ +`ՅL/U!;<}#ۑj:|23<p<z_='=-s< `= <%]<4<_?9ԩ:W]o=!^ 99N;z5[=*Υ=U=>J=ˬr= 佼S9p&＂#+C 3nH<l o<h =L=kR==dZ=]=&='=A%<:}Ͻ Ecp HA2;7ꈼ E.;&⼱V}Q3/ t +k==j>8> ?!3v?D>ŚI>sPb=1=< 8ͅZ߂vtJu #;#W8ܼ8*鏮tT:н$ ݽPEoh,Qh̽wT>XRb7y<!e<2<}<e{m7
ཉZݽ\< 9< :<=&C{>۽Oƽ w{#W<
+=:f=N=T@r +c<l< ;Y<+;J<b;[I ͐d9`Ү9 %:#S0=#> >&2>(G!> @=AAeżn_&Ѽ#8;1 <e<CqN:hg< =Y=,==h= =H=cm=5h<z;njP|dBk4$f +,p2<G^,݂Ѽ!BqHu=N=Uȼv7?'D˼n Ľ3fS<=]W==c@=',<:͇˼>GS,ORݼ5l1軆嬼RgN5 +Uُ+|ؼ;Hw HPF~pMהļV&d ּ5=og;֊<~c<$<ߛ<<a9A;0<:p;c;: :` ֆ,ƪ&y<9=-> x>D>C+>g=l +=$Vk9BF+7k5P|Lk +3>==Cu: +HڗK(!DT<L<c9<]::<RU<=vw=q=1=N=r<;N`Z ʗ:Ç3:-s\|<<lo<Ɵ?< u)U Ö0r{-CH>*mϿ*IOüN2<=^4>'>>>>]? L>&8>]>=G'\<+ZlShxʼmݼ91 h0;
f<Ir<vs<G; +L OhR);v<Ր=CS=, +/iMc$<= s<z=j<=J<ν=XD<<;;~C='QjOx~;Dt<.;Ɋ<=V<=5==>F>j>@>>_D:> +7;RD<ڝ<6:Y;sF<Q<l+<&4<~(= +<<<2c)` ;Q;;;'(3_P(7'PU2ͽN/@;XݽkUd\g}< +>iE>hu>o>,>>MA>3==h=v=(s+eRXb;/K$O"ev:S;ovjr$Xӽ&J!mXl(9K;;L;
qBet ba+(L^%c:OL3I_:9x;a|L[!=6Ď +Gj{ 8m::O;є;!l;)v$3%/#z+*B8LPaQz) RYI：k½2伅?:]/<C< ]ǻfP6vʎD5e(Gw1`9Y /M ȼugB<\<P<><*<j<s<c<k<d<~"<m=̺<S`;i_U%ޮ& 8}$ֻ-s +ử;qWwWNu^̽cE署Od,ȧV<!\<NH< a<<ɸ=%mi=Ic=El<<zC;`;=ڠYrջ,d 0Ŭ d%ѽ2GrѼ¯A:q +>z~t߽_ڽ
`ƒ ;t=-> C>y&s>a'a>o>K=,Y>{>t>0p>*e>#M == [=4iTռX zϼnV +ݽ! +4<M<{H< <o<<w<ĻM%mz
8 +N(\r.I 9μ+J՞.:-;A<F2=g=#<<>M<Y<Rc<=B< 0ܼ̂ a༴sDNX~*`YԼWPad.<'ߔܯꀫ볽[`ogk~|ܽ\GCGl/-#y2A6(QE.5IB%rJ_/ǽILm; LJ4n @M b$֢Z*z7m; +Z#.k;k[;Q<,=%=(<x<i<H7<><B#;0K<<}k0=T<Jt<=}ܻTae'?rre63L0HD ͽ%I\ *$>&?* 쇽w 3➼"r;8:i YRM=@)@)]M<x<T<`F=kv=u=#
=CG=^<-=-
=R=_<P )->&ͼvJ ƽ0(:f 6:Q(̱ ; +e|CT|[Pc%#Pv <}=.==kz=4bz<=>==뛐> '>"2=)M=9:ݽ ~H IŽi;IyBǽ&ډŦFbfSP/rս&Ԝ˼7{/ +:=G< ;;ek;;JAh"f wJ)I=)+L琽^P|2Kt.:z-6~R]2쵽pD0b(6>rǼr>} zr
7Z * -EBI)jE4necO+15(%C*WV o)O8E7&t|=)K s[܌v=ļ̽7q蘼*"?nHX;K%<^0= ?~=&n(= #2=6=1=<<Mg;j4<<n%<4><<Mv ?Y[~"EXw"}BtF:DiW!?龽 +4H +Q^'aY:E +l;0<aD<yS<Ǻo VPԻ<ro* ֽ +~DEZ7WPVOȽG ?m h:J9'G?M#dۻ2'X G &ҽ +;q<Is,;elN:z伒2)T$Q +M ⲽl;"-Lˈ<:=ds<Õ<`<oY< >3F>~sX=> ==r=X2==o&=#=@`= .S Gquh*f80z7H<-uX:t:ek ZqzBo}콆ͽD%Sλ5<5< +<h<sNyƽfܼuv +<U +1>?>j>;ai>e=ۋD=R:wz]q::襻^PLIcmQ0]˽WTż}Sn
* f~bcּ)<r<y;B<ފ<Wz==M}>/K7>[6>]q>>=ŦE= =l> =ڙ={l=<VJw~zi؊s#iŽ*TU4g<5`=_)=i==<QksTuؽ2"l˽;<s=`*H=A8=%K= ;Z < ;w#l/ABR ҽOFwüDֻgX#wZ^a1i磽k^1q1,vwxR1%ҽ*? Ja-y;d"<=WFQ>8 >>1>>yX>?@=w=:<<;&!:`F+/$◼?꼑uax <l:=QD=_=<"<k;a;;><<մ<;qBFOc&T<<P<<;p;;Nq켔8;S6Udܦ +waf3k𧽄$rEӽR\e(7|K2t`q\kw/़LuF"x:;6fUjＭH)ϼT*<N=_ +=, >==>D5>o=;==C|du'ǯ;Jw;yL6<NZ8c ż[ż aQfż: + +& K8m `{*ŽH:<U=*ڄ>:n>:|> }q="=y; WKXʼ;V_"躅c<Ƣ7=4= ==K=9'<?f;;3dRt(( ࣽ>#[L 0t~)̽`ji9ýa_ㇽet) + NȽYq;˽u2|'aO Sfa~/p꼝 Ru<DW=K>T?eC?'? >Y>W>>=g<;=Mɼ4O<#o=8on}K<8l=@=v?=Զ=="</<V<L7<#< ȼ _9u%:<T='<(w<&;R;z;:o <b;ok/ ~c[qb 2h'tzX ;ݾ0;-H)owqY<_T +^{3;O;K;"bԼ]TL )S? ;4<y =7=C==c=1==Gb;ƭ q+F.<er=, =b^=y<9}BA**ʼ.1F@k=gO=#C<I;%; ;9<?<=T7=Q>}>*=ۈ>AO>>^)>0I>"
=H=<F==J=" <uZwaT oŤĠ< ?k=:>(>T>>4>W8==S e<y;́λ+
|:=>=_=o>=3=;ӻe*<>St<<*2AeUJnؽ jŽ +|O(ag%7A*7O'4ͽm.wu6 <j=˝=f$=*,="X=95-= +6igHDBTP@^Bfy>3C98~&;͘<%=O= &==m=;==h==X=h=y{==l<<sEԼ3 ^VEj%=='<;c
E;O?;<K;a< +=B= +>81&>D==+=>5Β==l=V= <<|;<n<mՑаQ"(>f3-=#=>I>>T>>.> $=bv<fa<;u<6z<G~=+V=\l=:h=i=r=7u<7;n;/(< 4<< s1ٽ.ԽgܽPս;ܽ#ƽN80yeVx^fsD9>WY ?˧=<R<=$>)>>6>>KA>>< =e=?S= 9\伹漗M nK2&Uh6Ƚ|Xak;BB^=*3==䝌=!9=y=v<
=
< <D<=|Gӫ{9W<T2<I==Jd=J<<)<fC<<C<*z=G= !=&<)eT <=U=5===h<:%G^( +"y +p̼ v6( (:Qլмх\i};ͷ+&;$ ˤ;+1<]=]=+>p> |.=ޥ`=&=AN=pW={=j= = + < ; 3<N̼ޘ꼨μ<B<ݥ0=+:ս&};"Q=^<|L<^<,s=> +<8=
Yd<K<<,=&;"[ؼY켯k1<!}=3>>1h>=#<\wM~;Յ<<:2BJbc3ROt;n;TZ;|36A;;P<2<=C=s>!R>XA>N9>)d==
N=%%= +u<1=&=4$<V<>`uF.jMѼ>cͺnN^{;ZH<4;ʈ=+̽<<<1=v>,Vۼƻo;NMbVo].Z6GمTf^dzڻ=
=H=>"D>m>>> +< I<GUּkpֽ괽6 YRrʼ2w[Y a½wSUY[%8C([5]x}1Q<_;=Xǻ>mIط7ٻ"Q<Y7*Nh
9;3:K;e<F=^ +bcס@K۽0ν@i69B3 Ƚ2"_غ R<=i=_B=l -= => x>F4>=Im<B=
D0<O ݹ|| +\:=ϼ20u7Sw c{ |h1WX悔 {헫=z07+ʽZƂ;ἙU;!+< H < 9t:G=b=>=J<hL<$ <<<x=n=>P=?:Y X0'R;c +
ȼYkK<x"=6h=$Q=<x<Y<u=K y= 2
<Ʀ<<;ƻ_ +꼜sL
<@0<: <'<'<ž<<ސ <=!=k< [;ՙ<7 +;U;*{;;S<'ē=5=-fs=*g=%<λ4i;[x")bAy;|;jqb;[5;8;ɻ= K:<a͈=w +> s>D==n<z<5 <8< <[ϐR2/=2͖T<9<Z;@ʻۻ3;7<f<:Z<=)7=dn>?p>c>>>;>]@L>q5= +ҏ<0;$< +<:͵<q;cs<':j7{мGR=[û +=<B<Q<%=,-=$=
=3i=l<R<pUL]~,(x.۽^f.%x^߆o㸮<.=^=l((=`1@=;:e<<A[=;<O<k<0ٸ_W/4)ql 2~{nȂʽ
?ZxUҽr Tֹ_C3Q4Lxf=l=!Z=~<zw WY-j<|=Ѕ=J=D<=Jn= Q5< =<M===^O?<ɸ;W <;v4 M9|+<c<ņ<%<B;l:\<o<B<o<\<};ΓR$hhN<9<h<=\/=<C\<< :;=I<7<9(<<|;Y <q)<Rm;R;ҩ;s;<0<tf;E7;U.eWxG6X aMw'! :#;Șc: u߮
džӼ)<q=}= ="=L<S<e<! ;k;}:zٻ* uMbJˉ&o +>>j>N>M>$=tM;WbJ92ځ t O$y3=;&_ +yU;u<<=?<ٿ=Q=-?<l<=/=!z*c pD]l]#=R8ݽ߽π";r<oa=K 3=2:<Ŵ/8t; <O<53QM\C + +E\VB. !,PϽ쌽Nl^:.0t[x=^=zx=<Cɽ);m%<QL=yq===Q=:[=n=O=>== 1=c<2<<,cL<]<?-<n#;pF;;!: *.c9 ;8';\< I7 y1Wzֻ!Z+< = 9==3=a< ;PKiqm^G< i6 Eeۃ?}n;R +9 Ļ7A))ؼb;a)<= f=[T=F<G<֮<qHd;%4M+F0<nP< 0Ɨ|s>ǽ8 +jR ApR v C 2=>J>3>}E>V[>>{9 +nO:<7:u~/= =ȼ +Yzsiɼ6)7 [:¼a?PTIŻU#; $<o<6ϼIY0A ,y' H0C 9}K&U 5H0= ;ѻռY!7="=6Xw<ٍ<T8j=PxjwZ==V==j<<=,K=3=x=i=4<ݫX<<&<i<~R< ;u4ἎBλ;KC8;co)ݼ^.B;x;;@=+M==>|== :`lP kҊ$ Ƽ<։.hq; #<I'<]4\<9;o>'Tg+mۻr0u֩;ά:<><=h==?= #S<c; mD;;Ӥ8Hթ; +:z;<-x= 1=&2O<n{<Aia<*Qo/ =n<P,/ۣW_:C8ֽ, 5m̼Gn !2nA=)>>@>|7>w>H$V> M`==ǐ6 -暽 ټy
o&;I7L E@[ +1)~(FΫ;} +l=2=L=7>; 4=1>k> *=6=>=1S<=D=#= =v= aB$OR-HAiq-=^Խؽaʽ6+4F~</= H<[<ƃl?ڵ"ͽ[;+3;C~<7M%<2N9# ʼH=J&_%v<:!h;9yoL ;?SS"rȼ9m8#C dg$ȼDjcܻ6FX=& =;,p<[=Wj=W=;=zl<P <I<,V;;L<;;Dċdy+/-m
:İ 1o
rT;;=ٳ;5 <33K<=s=n>(%==A=2.:x9jzA V+1+ +s<A;h:|;Sq;^K<1\<'Hu<.562 ++òѼVf?=+6:ԽIGG^*`TE ݼTF!<0J=K>n4> u>J%>*#;> ={=< 3\<
<5029\wyNOepż +Z}3K(x/,<-<< <::9{# r";!g + 4;<<J_=U=<;]"46N-AB" eX?sѼ>JU` +:ӼD@ۅIA ػMӻμF'V_<aW<v<~ = ٹ=S==R=|%=R=I<`= D='Ej=:H;wcܼ@-FqoD};<!<w<?= =
a<%w<EB<2==C> 6>0>'S==@~;XTI( Fͽ:h<l{0<I;!uQ69<s_;A<E<<i<Z<d"<YF<<ȶ<z<s=D,==np=^]= "<bI):ݼ+zf;3G;_ARZ +yD<;{;O`Իм&WE DyԔX;h4Nk^mf%0nin=w=b=(==f=u=7==<;=I;<;<9oknjJ:Z:O.<j&;[f;O;BL<;)<<=$>=]= /=x8=</<v +O(=R==@<fϼfv
x- o(
Suj'|{;2V};;;ڔ<tZؼD ` =(8F<X=*=)?=h&$Pٽ#Z⽉ ܳ{\oC8z~X" +k=z==e= ==P<e5
< +={=b;0˼/ +Ye"\;J<~<<E;S8h;<>d=Zc=z=w=uA=/7=
<w
`<h;< <Sc<Lj=W='w= i<%NI 㼸8`a_dw!TI2⼾1f{?\| 6kټ*X09c۽*;$Aٽs>`m~s}
ONڻ(=Wff=*= +r=q==A=r<< <fz;@<X<cm:o-;_`y=7=*V% ;;=;9B< +=#=<QgN=I<x/yU]+itؽֽ퍱o|ғʽn9 +=3=&=Y=rIT<8
So!=ʠO'u <[H<ܙ|=<= +=, ,>
4"$ DXw60]39=(zA˖w欏奻0q;
k<N ;[ˉ]YlvD;<Y= J=8<<H<K<= =\>$*>Y +E>K7>=sA=s~<Xz;g;|;Bؼ.'3^< {<V<5<AJ):J< O=TO=>1> >5W>1>y={<S#F#c(ٽQ!C<< <Bu~: =? ===z,==5y<X;xq ;0u3;w;7;Ԗ< <{ 5,V` l Zٽ_ټż5Y{˽^{(|I|ѽv+ؽ +w\(p@uTm=bp=<0<ܮ===^b=j<n+ؼIH;6$5<!<g,7=̻qi ⛽* +~ȼ#q NzG*6+&/==U)<&)
ޤx[J̎Ҽ\M `oΙʀH˶ᅟ 罿 ƽQ;߶>0Y>Kg>)>ӳ>PH=Y=;k*<&3c<F +<ˣ\=.>2>t1>zf> == =$&==04<Քcduۼ':d4o"_<=VH=> ~>=§=<;(;Žd\ ˽ ⼐'&g;zg;tA*4; + + % spɻq&ݼuXսPr;xl<ٜ<'<<0.<=~=Ҽ=p=CI^<CJ<~c<,<3#PP;SX+Dҽsｰp.5X` wVPȽ\*}=!= + *(l7p̼:⻥Ӽ-#¼0Oڨ| 4R +0< <Y< +<ػ[\0- Ty7:2{"c,(F0dNB@_}]s4u^h( ~`)+;<l-8:d=ƻ 弢ct` &<vaٺ)Ey0= 5=Ll=<'<~=@@=zQ==/=RН=?<R<<Yp<＿;<%>(8t6$מo3z +ｵb=]j<:J +ͻ<;o< <tr8= +FL>! +Pټ/ռ ;W<Q< $<{] f +;齜:̯<D=K=a# =>>Z>}>2>K(= =n;I;)b<@?=<J<t6bͺYؼծu{R *s-''-;1< +;J )V7k$=>@p>1=l=/;M`ں:e<<2ڼPuO
ý +X+=<>=fp=UL=<;p';+= =4>$f(>P_>pj>~>[>}><A=m= +mPh 1>>Z
=c<;<<;<u<(<,]y& p,?x]i0b +r缇ɼ*Z=G18=r-=٫==q>'=s=Y=t==}>/:>=[!U=%=ݷ=B|==SZ=h=<c</$;q;Ɯ;k +l7<ou=< +=2 =u=],u< ;+/d==N +=);i2<p<;k<I< 0ON
D[VKbNJj Sx kupV<F<L<Q<s<a==b= =!<VջYP̐x@żwo7<Y+1=LA=lW=O<XȻ[c%{]<ý+̬}s\]/:YpD|Ei U>< +2
A㼽7;@< +<_-;A8o2'ցQ L
#KXKȽ/;ѽҳ?7)ơ.ɽcU;1Hv&Լw˼"E;<n=y=ˌ>6= ,= <9v,;=J=>S>i>g >M],>;==כU=<[L<+<F<}=6<t=y<-Ҿȃ)6q;s|ƽ<%$F* tν|ts=uս I/r +鼅~:g8u_+<Z<0Ի/I ~ʽ +U="=Bi<(Xս0BY~&;k;&0[͜Xl]^g Xs R.j~nݽ/L}T@"<s>>_z%>&= =bn,wQȽ +=*=<˸ɽ3D'z)648j؟|(mZY cp;z<w<`3= 5=a|{=>=`$<d<1);~ ;><Z<CѼJj1!#E('8Ƚ9 +Sy\,m;6ԺI9nH9)_֚iɼӷ(DQ'L.~s7WIf;+DI<Un<=ld=V=Fe=L*8=V1;͘CC:Ӽ)Gm+<"6=F=q
==DMm<% +]3Mo:v,% + ځm k;7 )Qr:ڼDӫ<A1<=z =|B=A#=U1<`"P;&iffn<;3{> ?(1h lt
8#,şqt-9૽ cר/dT潛p O
H㼄a{R#Ӽ +MD";h=z> >Bw8>1>e>8 =]&;;Cw˼!='===7=Wl=,=o,N~:e `̽ ڽYk27e=@yqzΰd,JXM:LP녽7.(N4V<j;`&>q>d>=i
=wB=`O# y<OIWX,i*ý9+u<tcq#J JCRQVʛ;i=</)=gr==U=Y=<A;R8x wx{.x_Rb٪4B 5 +ЧD7 b
+Il# }Nl\Ҽ9 &[$p1ͼ<!<yg(=N= +=ÌT=x=8H<y<<Z< <<0<$<F<x,&F +p ۴ V&B R3ⷽVUF y9S@AA^<t =2$= 7s=<<5\BȀy \<=RI>t>Fo>z>z>1$=\;%cHI=J=n=v=ʖ=}Q=J<l;;C)q\~+ нia0٫{{N`|j cO<m~y)i0G=A YӽR쳽{"GE>7>?>q==W>#P=J;[vI84:~( +U6g7`4~,r տ*tC<AW$=y*==_<i<:P<QEnW .k'-]C3>. e)\Sj^ +; +<T0<@ +hz4Ѽ^ SC=@=-D@<C=3=!<q>ZpR1J= !> &>]>,j>]=>hy>!0o T兽IqK:U/Yz"w^auI'"쓽-;mq< <9%i .R* čT?=AN튽<Pqr=&{MjѽcS%)e| +*B'<=˽6~,7 T? j１]ɯ<ܻ<E="HP<:<#<Ui<Il<)<;MJ<<_<R!<l<6`=Z<פ<C:[ܽ"t.%l!%3yx*8W2(нgwwq/ ;=x=N<> +4=C= =*"<ݼt:1s-_9=u4=\ҼD\6h<bC<=JA=n =`<@^:\ӼQ+e]vSn舽p9ǽ8zGa,S8
RüA93=S={=ͣ=.a=g%=C;澽z%wӽ >b0==Ec>,4>>X>^M><e=1p< 8ݸk >?zIǽY_e#ņ}9;9Pd(j F9䝼s.N꾽a{l:zyѽj a&`uzUz{R ^FݽIǽ<:}9样)<=B
>p8><==U=O|H +ǼV;V< ;3Z;G69C~,L"kVEX!x+f2$Հ; +S t +Qνѽכ/`-2VDV:<<g6=E=J=W{< +Ⱥ;)#{; <;8KARX& +g~ebd!ż)s-@<.9!dV, +3l}F 4 ý#)uĽ9hgBrdof$'9<^H=EN=I=9=h8=Kb;_;!:'.:t;
<rp='}o=\=T== =!<J c̡,Hټ
锻 rű<1<<<`n|JhW"(k; <Z=<=+>3>?h>_e{>g==<s8zl)ܽXZl֜NXVI"ټ)WNJ*O;W;]kj nPao}QD xwe:<g!<Je֔;;SA?m<a<66<ܵ<pZ<^:=<F=T=9==U=4= <Y>5K=ú>DvqEoƾ G}̼Žy1";H;/t:[/$A; < <[FmF!= +"ڙ=nf ;]`F{?g$t"7Ts.D~w=D&t>V>N>w>>>=S+溝7MBcM#@("?|k;8<<<<tT\4]ݠ7熼S~T:!;}JmKNSW&!rNջw7ni:M:̓A=Y=a===T=U=o; ;B:Wv<Ou< =-?===<S=mr=I";\N: _P;"<{ϻམ=43=5W=.)=$9n<@"ݻ~8u5P<Ć=\=P>|>>/?>|>wۜ>=o;ѼR'N! ǽ҂
ڇ?ɲj2_LSJ0_. e T;U;h:;UN~9Mm KEӽmսCݼVp듽iQaɼ?qb< <p</7ѻ<$<_<<iZ<ȑ==N=~=.t=r=}=l=e<QIཤ ⡼pbn RYXOja㽒.v`C>m<i ;lO:ʏ;+'ż+c2UUٽQn{=μ@~Q2uc ] +F{/7<*F=?3<Q=V<Fm<C;cȺW3-vC彀 #W|dc6 n9>$= m='A=y=65A=*2=<a<<=W=w<%< &=%=l=v2=qe]=#IGK +p_}ԼkϽ +\1ٽ2(; S</mY<<O7(<#<[<J\Ye.RQz;<-;<q:,-<3k!"Bx :=/I>>>->>L>$> =>%<dcjAl UbxJ<=aC=ƅ >Z=Kn==v<<G=Rp<nj:#x2
u\kp%M6-Zͽ4Zi4= G= /=o#==SZ=H`=='s=M=t=x=};F*<h*= +4<<c<f<fO]/V]$p"ϼPؼhfARv \8ϓ䶴 ˒֨@'ꗻ +=+E=)=R=]x<]
<6> == D<۪<eM<vT"(<Dc==Q6ݽܸn7r <n%Yo4ֽ75z<<= |EӼNֿe=ف"഼ k|y;}:I<= +q(v +)Ѽ.,q:@<`K=AdU=fc= =S=ϫ=J=$N;sl0+?l7
'JEȽkwٶZB,'K;rM<.:= eO= r#<=t=4 +'=n=oX<Ῥ<<:l#:;<<t<c;76:xMg wÈͽN'B 屽<ٽ@q-uJڋ9UͽY н?zH;мҜ2g ٽQcʬ<c=>f`>}>xJ>]C>m>g;=C=8<kPgDov0[=<@=E@i= Z==[<ֲ<< :al=n>+>"b=]==K:k;q=.=P=v>==Ǭ=Fo" \? :(<h n0Ld<y;m-=Jս'f <$=<B<y ==w$===y='a<=M< Ӽ0Ҽ~<: +2<I=K=I=R8=q=?H=7q=o= +e5Nj0Gh.ҽNd^jp\琽>a9';==n$={j=Ȫ= ==D= v==[-=t=°o=f=L3=xOK5c<K=|5={<;f&iTpغ]y==V=Q=B==LM<:=⬐<&g=<Ưo=5s<F<nd=Y<A<`F +9Ὄabɨ, 1)g^&Qy\NiI+exEo%D~x60K& .<'=ȟ>*s>FR>>s=<= M=q<-S磼iW缳ۭ:87"J< += c=^=<; Ѽ*"<֞=Z=F<Ek<>> +>={=?b=ѧ>>MO>N +$A'p3rJk^hT6;,U;ӎ=!:Pwr̙gqL< $*==;>._">Taq>@> v>V== +=P=ӔS=F=h;#L7-&,* ̼@0*wGa< <=s=Fs= =x8==|= j=Z<?V<====:>:j>>&> $>6a=ӻkT _[\\o¼ ݽ5罽#Q̽O\<_=jz>>$>E>K>J>=< a[L~&&c_@<۶=&=_"=3=4;= +;;9q<o= {=+==R=L"=q=-<|m<Eg8p/:yg;<s=f:> >6>mi>|>b =CА?FL^|n rk~M7x.C=9Y>%97>?4?|?oɸ?,>>=e=}n<L Q1<=>A>I-B== = F }ĬI%^Ѽ-59lLnFAKP M/>4+pN$:ھ$4=D^==2=5< \ɔ;n<\f<rq>=pK====Ê=S=<:8< d;<g(< +M<\ <M<Ą<ɘ== < <(/:(P<J</<ă==<<N%<';QˇF`o]@be;le=8=Wo=_i=a@>)>/ w> l=2Y# rнѼTٻq3= t=R>=-=tJ=
=S=V; =s+;Qq`Lv%_ +Հ=gy= =a=Ա=<<;'H7ua;L 7<=K>>$\> =ڏ<<o>ҽZ\D0Q0FlѭD+;t=>9X?OC??"??b̫?<>O==GLɼ; +=0>;">>">v>6u'=6ɐ* zu5K꼽V MCz=w}p }S4!U,M?ɾ*XVyP9t<ż#:T
=t=G=zD;!;6 S=_= S= +=}=-o<a<xo?<<l<J '<gL<< e8n;:<]j<`=pUU=\<=/<ht,<
qrcl[Όr5N?;(=+= =>!>">L>n=*Z ۽)|}/
:
2=UZ=/> >ؾ=s=<R:"wyq^jbEm"IQV' :%5K>KL) {GU+:}E<@ =-=;==\3<;;lȼع7 +dPӼ"5&ֽ+aɽB+$[w?gSϼSa;d"<mq==9*=,;x':Y] aI|V;W<O<<+s4_I ؽ>|2u{B6WnQq B꽃ڽH}ν)뿜=H>{??;<?ŝ?&~?@?r>s = <PлkJ#3;.<
=e>/)># >V>u A>S=Xp<i ffIV1 Lǽy0/[_ ^M +=8Cs95 =J7=|=U`=f7=-<&=Y$=J=P=C=\$L VI7S L;[ <= GI<5Ћ<S<=ߊ=/=O<p:M볼p40[3,9rl:<' =W=[ =@ ='< )<<Y.J<&Pd||I ]a=!A=ߙ> ># >=w==<9T`[ѽo8b+[B-ǽ< + 7t̽XX1o +<`Q<C<|J&<<1b=/"<ھ3;i Q*7!F +eԼ5ytX-zZ
I=QcVL@p?X0ͽI&wk7Mwٯ%W{/+XW]*eż8νOf*5",j 'P[OýcͽsнbJʖ#ۻlV!X=9+>>>az>'K#=s<<a-פM@ed>y/;l<y<<nE;<(Xpo<^=L==<O\=t 02 ͽA$,7˼9V^<Jѽľ'Ym=$===,|= =c=p[Y4*/G <iub`d
K%[=d= +<uD<lt1U=q)=W<]@<=/%=<\<9 ;GucK/tNҕ,1Qn <;<X\<fl<s; +tYlWʽf@Aܽ7 e2 ;j҆Ӻ+Ż{Y>e#-5iٽYjJ
/3F͊B08w $/ W!Y +Ľh奮> =!j==N=<;ν0
kֽ%;<xM<[W<v< ;ݶ=,o6nIF3<VG<<<i<<= <6E;^;j4+H)Y',Iea国V>DS;:<
) OÚ꽑B + =,=T< +{3Y1==X=;"< ]b IнQuo=58=={> == +뼒6n+;< =h< <:;C<} ǁ^G +ټtĘZyN
W)~:=wٔU;fF7y|韼B.'y&b彡j_^dB%rR%=>3Z>+W>%>I >>Bk>>k>>N==s +s;c:<a>=%=˕=C=Ԙ=;1ͽ'YxxMXsҽ|ͽфƮ =D(_uF;^6 <F<x<<zAH:#vr#ZνrDqEo̽JR⯽w^⽍g;3ӽ +kvK DR<f<_$?kfL$E2VD 2}/= NQ`>n# <*0ݽ@z*FD߽)߽SrIyq1 A@a ':5;_]<+<9lPL], kOؼ =ӼW~kQɺ5._;3+<$B;x;
;y((~){ +6rSW І!87B`+$<)>.CH>;b>Wf?I?yc0?^J?A$?
i>>[~>>m=;u<V=,=y=e=΄=L'a;̽>tmMjsy V C1G @֦|ڢ=c^P<pl<=<`=<w &+4v<*ݽ㺽 !58I^aL2gPld~X~2,wҽ<Z[Q9Lϑ}'qRJJ5;#@+¾:⾅9 +~ O`$:N +4bDͶ]K8 ݱt'/Z*
;-<ںT +=hs>?@5@@@@"]@??>c>^3 =<Q$Nؽ@ M< +#,aּa Fp=ɓ<?[̧$uj<]_<;=݅M__5A3cTT}a
ʼҝ [ Ƚw}GЅT=}Gbd<";µ=W==y< %;ڵ< +:,Md~[$=x<=q^e=Qq?e@@ԳA|CA-A"@@.3?/>;> ~=(;./
]ҽ56ĽÑvld
Iy_:[6Ͻ,w8촚@E;=?<a#<1e}E1Wڼs*ksސ ;<R^;F<<Þ<r&hս?Լ̼`= =KV=->a=(=$= 4սVս2" |y^ZD +i=Xm=(6f@¢n2$<:&uf==\= t&,%ݡ<^|= +R=6=#]< ` 9h'漗u9dHۼz* I;;4p<ID<ۯ<@<'F$q( +43*9$W=ڼ=XY<#ED0e; +X== <iX<6<O;}㒼;% ]+5j65]3W;vk;/|'Fg;@< +!L]R^$
!J*CV43z~@n*1`:¹=rjo=> T>Q>?ܜ?ɾ@7?9?>}>[== =O < +˽S;~ٽV\[wG½ĵȽd걼ë"{=
= <)<@U<NY<%= M=O=23-=<yU^/3).Ob:PJ: ;v[;5;A<u<_<J}:kD
<wZ<f=VA=ս>>t>xb>':%MOF,$־!@<ߦ?a]ZֻEQ=sj=&e=_;b0 o==eQ= +Ľ' +ռ|.ɽ νH̽m=ɼûxxlN< <2=s> >]>q +>x>|>v>><>= =T<pWҽyZ *(y"UB[uŽ z%<r<b9V<Yx<O</<~<<+,=Q<`r;獑5Tշ1`
><>j<i<<<Q<<b-;Z;);<$:AUp<) 3:茌=U==V*>|P>N=>-V.A6ݖKV =7>#=Z}!;٬<:D<=U,=V7<=,c= +O= =m`?=:=-y=2=} =eM=w3%=<A<?!!O&pveݿӼ˧¼4܁ἝWThć wq7 86i; O= +=Us> +3 "p= V ΠH8s):ww:<ڹ= .=u<;[WfGR+´@&E߽{+hr|q"}9<cѼ;G<vC= +'<<'ջu-ӼSoS;M3%;<=#{q==r=>\=A=Ϡ=>!%>l>f>=X0غ <=Ct{== +:j;;=_=2h=1
"/^SuW~Pyw=6=x<A;ėWl<FQ* TL;;0peνld*05_&%KC; +Sρϼ,ȼrX6¥qW md`ri$bUϐE#ժm9jG-v b/K <^<=;l;J=+<X<-=== +5D/`p5 :84I]qCrż~gsWG`ۼּmkqVƅ`Riӽwн
e BQ'齟IUfmXi=S+* +ҽ!~` + ֽv!x}n<> +ߺ>B=ؽ$!eN(<);<A=i?=W+l=*=~=a=i=w=7@<5 +=/=d=9r|=oV= =-=kv> +=0dRqn<U1`dMf+N|<J=!:<c<5=3= Q>>>C diff --git a/lib/test/data/2chipA.fits.REMOVED.git-id b/lib/test/data/2chipA.fits.REMOVED.git-id new file mode 100644 index 0000000..50242ca --- /dev/null +++ b/lib/test/data/2chipA.fits.REMOVED.git-id @@ -0,0 +1 @@ +fe2343477504847a56ac369c96fbe2d1b9d835db
\ No newline at end of file diff --git a/lib/test/data/2chipB.fits.REMOVED.git-id b/lib/test/data/2chipB.fits.REMOVED.git-id new file mode 100644 index 0000000..1d70532 --- /dev/null +++ b/lib/test/data/2chipB.fits.REMOVED.git-id @@ -0,0 +1 @@ +8a32e201693484546f631c52e3774dfbf278151b
\ No newline at end of file diff --git a/lib/test/data/2chipC.fits.REMOVED.git-id b/lib/test/data/2chipC.fits.REMOVED.git-id new file mode 100644 index 0000000..7afd241 --- /dev/null +++ b/lib/test/data/2chipC.fits.REMOVED.git-id @@ -0,0 +1 @@ +167e2502832be1609cdde29c52e9dc8c571cc3ae
\ No newline at end of file diff --git a/lib/test/test.py b/lib/test/test.py new file mode 100644 index 0000000..c6c1fcd --- /dev/null +++ b/lib/test/test.py @@ -0,0 +1,260 @@ +import os +import random + +import numpy as np +from numpy.testing import assert_almost_equal, assert_array_less + +from sphere import graph +from sphere import great_circle_arc +from sphere import polygon +from sphere import vector + +from .test_util import * + + +graph.DEBUG = True + + +def test_normalize_vector(): + x, y, z = np.ogrid[-100:100:11,-100:100:11,-100:100:11] + xn, yn, zn = vector.normalize_vector(x.flatten(), y.flatten(), z.flatten()) + l = np.sqrt(xn ** 2 + yn ** 2 + zn ** 2) + assert_almost_equal(l, 1.0) + + +def test_radec_to_vector(): + npx, npy, npz = vector.radec_to_vector(np.arange(-360, 360, 1), 90) + assert_almost_equal(npx, 0.0) + assert_almost_equal(npy, 0.0) + assert_almost_equal(npz, 1.0) + + spx, spy, spz = vector.radec_to_vector(np.arange(-360, 360, 1), -90) + assert_almost_equal(spx, 0.0) + assert_almost_equal(spy, 0.0) + assert_almost_equal(spz, -1.0) + + eqx, eqy, eqz = vector.radec_to_vector(np.arange(-360, 360, 1), 0) + assert_almost_equal(eqz, 0.0) + + +def test_vector_to_radec(): + ra, dec = vector.vector_to_radec(0, 0, 1) + assert_almost_equal(dec, 90) + + ra, dec = vector.vector_to_radec(0, 0, -1) + assert_almost_equal(dec, -90) + + ra, dec = vector.vector_to_radec(1, 1, 0) + assert_almost_equal(ra, 45.0) + assert_almost_equal(dec, 0.0) + + +def test_intersects_poly_simple(): + ra1 = [-10, 10, 10, -10, -10] + dec1 = [30, 30, 0, 0, 30] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = [-5, 15, 15, -5, -5] + dec2 = [20, 20, -10, -10, 20] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert poly1.intersects_poly(poly2) + + # Make sure it isn't order-dependent + ra1 = ra1[::-1] + dec1 = dec1[::-1] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = ra2[::-1] + dec2 = dec2[::-1] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert poly1.intersects_poly(poly2) + + +def test_intersects_poly_fully_contained(): + ra1 = [-10, 10, 10, -10, -10] + dec1 = [30, 30, 0, 0, 30] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = [-5, 5, 5, -5, -5] + dec2 = [20, 20, 10, 10, 20] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert poly1.intersects_poly(poly2) + + # Make sure it isn't order-dependent + ra1 = ra1[::-1] + dec1 = dec1[::-1] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = ra2[::-1] + dec2 = dec2[::-1] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert poly1.intersects_poly(poly2) + + +def test_hard_intersects_poly(): + ra1 = [-10, 10, 10, -10, -10] + dec1 = [30, 30, 0, 0, 30] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = [-20, 20, 20, -20, -20] + dec2 = [20, 20, 10, 10, 20] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert poly1.intersects_poly(poly2) + + # Make sure it isn't order-dependent + ra1 = ra1[::-1] + dec1 = dec1[::-1] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = ra2[::-1] + dec2 = dec2[::-1] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert poly1.intersects_poly(poly2) + + +def test_not_intersects_poly(): + ra1 = [-10, 10, 10, -10, -10] + dec1 = [30, 30, 5, 5, 30] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = [-20, 20, 20, -20, -20] + dec2 = [-20, -20, -10, -10, -20] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert not poly1.intersects_poly(poly2) + + # Make sure it isn't order-dependent + ra1 = ra1[::-1] + dec1 = dec1[::-1] + poly1 = polygon.SphericalPolygon.from_radec(ra1, dec1) + + ra2 = ra2[::-1] + dec2 = dec2[::-1] + poly2 = polygon.SphericalPolygon.from_radec(ra2, dec2) + + assert not poly1.intersects_poly(poly2) + + +def test_point_in_poly(): + point = np.asarray([-0.27475449, 0.47588873, -0.83548781]) + points = np.asarray([[ 0.04821217, -0.29877206, 0.95310589], + [ 0.04451801, -0.47274119, 0.88007608], + [-0.14916503, -0.46369786, 0.87334649], + [-0.16101648, -0.29210164, 0.94273555], + [ 0.04821217, -0.29877206, 0.95310589]]) + inside = np.asarray([-0.03416009, -0.36858623, 0.9289657]) + poly = polygon.SphericalPolygon(points, inside) + assert not poly.contains_point(point) + + +def test_point_in_poly_lots(): + import pyfits + fits = pyfits.open(os.path.join(ROOT_DIR, '1904-66_TAN.fits')) + header = fits[0].header + + poly1 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[0, 87]) + poly2 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[20, 89]) + poly3 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[180, 89]) + + points = get_point_set() + count = 0 + for point in points: + if poly1.contains_point(point) or poly2.contains_point(point) or \ + poly3.contains_point(point): + count += 1 + + assert count == 5 + assert poly1.intersects_poly(poly2) + assert not poly1.intersects_poly(poly3) + assert not poly2.intersects_poly(poly3) + + +def test_great_circle_arc_intersection(): + A = [-10, -10] + B = [10, 10] + + C = [-25, 10] + D = [15, -10] + + E = [-20, 40] + F = [20, 40] + + correct = [0.99912414, -0.02936109, -0.02981403] + + A = vector.radec_to_vector(*A) + B = vector.radec_to_vector(*B) + C = vector.radec_to_vector(*C) + D = vector.radec_to_vector(*D) + E = vector.radec_to_vector(*E) + F = vector.radec_to_vector(*F) + + assert great_circle_arc.intersects(A, B, C, D) + r = great_circle_arc.intersection(A, B, C, D) + assert r.shape == (3,) + assert_almost_equal(r, correct) + + assert np.all(great_circle_arc.intersects([A], [B], [C], [D])) + r = great_circle_arc.intersection([A], [B], [C], [D]) + assert r.shape == (1, 3) + assert_almost_equal(r, [correct]) + + assert np.all(great_circle_arc.intersects([A], [B], C, D)) + r = great_circle_arc.intersection([A], [B], C, D) + assert r.shape == (1, 3) + assert_almost_equal(r, [correct]) + + assert not np.all(great_circle_arc.intersects([A, E], [B, F], [C], [D])) + r = great_circle_arc.intersection([A, E], [B, F], [C], [D]) + assert r.shape == (2, 3) + assert_almost_equal(r[0], correct) + assert np.all(np.isnan(r[1])) + + # Test parallel arcs + r = great_circle_arc.intersection(A, B, A, B) + assert np.all(np.isnan(r)) + + +def test_great_circle_arc_length(): + A = [90, 0] + B = [-90, 0] + A = vector.radec_to_vector(*A) + B = vector.radec_to_vector(*B) + assert great_circle_arc.length(A, B) == 180.0 + + A = [135, 0] + B = [-90, 0] + A = vector.radec_to_vector(*A) + B = vector.radec_to_vector(*B) + assert great_circle_arc.length(A, B) == 135.0 + + A = [0, 0] + B = [0, 90] + A = vector.radec_to_vector(*A) + B = vector.radec_to_vector(*B) + assert great_circle_arc.length(A, B) == 90.0 + + +def test_great_circle_arc_angle(): + A = [1, 0, 0] + B = [0, 1, 0] + C = [0, 0, 1] + assert great_circle_arc.angle(A, B, C) == 90.0 + + # TODO: More angle tests + + +def test_cone(): + random.seed(0) + for i in range(50): + ra = random.randrange(-180, 180) + dec = random.randrange(20, 90) + cone = polygon.SphericalPolygon.from_cone(ra, dec, 8, steps=64) diff --git a/lib/test/test_intersection.py b/lib/test/test_intersection.py new file mode 100644 index 0000000..81512d1 --- /dev/null +++ b/lib/test/test_intersection.py @@ -0,0 +1,167 @@ +from __future__ import print_function + +# STDLIB +import functools +import itertools +import math +import os +import random +import sys + +# THIRD-PARTY +import numpy as np +from numpy.testing import assert_array_almost_equal + +# LOCAL +from sphere import polygon + +GRAPH_MODE = False +ROOT_DIR = os.path.join(os.path.dirname(__file__), 'data') + + +class intersection_test: + def __init__(self, lon_0, lat_0, proj='ortho'): + self._lon_0 = lon_0 + self._lat_0 = lat_0 + self._proj = proj + + def __call__(self, func): + @functools.wraps(func) + def run(*args, **kwargs): + polys = func(*args, **kwargs) + + intersections = [] + num_permutations = math.factorial(len(polys)) + step_size = int(max(float(num_permutations) / 20.0, 1.0)) + if GRAPH_MODE: + print("%d permutations" % num_permutations) + for method in ('parallel', 'serial'): + for i, permutation in enumerate( + itertools.islice( + itertools.permutations(polys), + None, None, step_size)): + filename = '%s_%s_intersection_%04d.svg' % ( + func.__name__, method, i) + print(filename) + + intersection = polygon.SphericalPolygon.multi_intersection( + permutation, method=method) + intersections.append(intersection) + areas = [x.area() for x in permutation] + intersection_area = intersection.area() + assert np.all(intersection_area < areas) + + if GRAPH_MODE: + fig = plt.figure() + m = Basemap(projection=self._proj, + lon_0=self._lon_0, + lat_0=self._lat_0) + m.drawparallels(np.arange(-90., 90., 20.)) + m.drawmeridians(np.arange(0., 420., 20.)) + m.drawmapboundary(fill_color='white') + + intersection.draw(m, color='red', linewidth=3) + for poly in permutation: + poly.draw(m, color='blue', alpha=0.5) + plt.savefig(filename) + fig.clear() + + lengths = np.array([len(x._points) for x in intersections]) + assert np.all(lengths == [lengths[0]]) + areas = np.array([x.area() for x in intersections]) + assert_array_almost_equal(areas, areas[0]) + + return run + + +@intersection_test(0, 90) +def test1(): + import pyfits + fits = pyfits.open(os.path.join(ROOT_DIR, '1904-66_TAN.fits')) + header = fits[0].header + + poly1 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[0, 87]) + poly2 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[20, 89]) + + return [poly1, poly2] + + +@intersection_test(0, 90) +def test2(): + poly1 = polygon.SphericalPolygon.from_cone(0, 60, 8, steps=16) + poly2 = polygon.SphericalPolygon.from_cone(0, 68, 8, steps=16) + poly3 = polygon.SphericalPolygon.from_cone(12, 66, 8, steps=16) + return [poly1, poly2, poly3] + + +@intersection_test(0, 90) +def test3(): + import pyfits + fits = pyfits.open(os.path.join(ROOT_DIR, '1904-66_TAN.fits')) + header = fits[0].header + + poly1 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[0, 87]) + poly3 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[175, 89]) + + return [poly1, poly3] + + +def test4(): + import pyfits + import pywcs + + A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz')) + B = pyfits.open(os.path.join(ROOT_DIR, '2chipB.fits.gz')) + + wcs = pywcs.WCS(A[1].header, fobj=A) + chipA1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(A[4].header, fobj=A) + chipA2 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(B[1].header, fobj=B) + chipB1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(B[4].header, fobj=B) + chipB2 = polygon.SphericalPolygon.from_wcs(wcs) + + Apoly = chipA1.union(chipA2) + Bpoly = chipB1.union(chipB2) + + X = Apoly.intersection(Bpoly) + + +def test5(): + import pyfits + import pywcs + + A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz')) + B = pyfits.open(os.path.join(ROOT_DIR, '2chipB.fits.gz')) + + wcs = pywcs.WCS(A[1].header, fobj=A) + chipA1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(A[4].header, fobj=A) + chipA2 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(B[1].header, fobj=B) + chipB1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(B[4].header, fobj=B) + chipB2 = polygon.SphericalPolygon.from_wcs(wcs) + + Apoly = chipA1.union(chipA2) + Bpoly = chipB1.union(chipB2) + + Apoly.overlap(chipB1) + + + +if __name__ == '__main__': + if '--profile' not in sys.argv: + GRAPH_MODE = True + from mpl_toolkits.basemap import Basemap + from matplotlib import pyplot as plt + + functions = [(k, v) for k, v in globals().items() if k.startswith('test')] + functions.sort() + for k, v in functions: + v() diff --git a/lib/test/test_union.py b/lib/test/test_union.py new file mode 100644 index 0000000..66abda8 --- /dev/null +++ b/lib/test/test_union.py @@ -0,0 +1,169 @@ +from __future__ import print_function + +# STDLIB +import functools +import itertools +import math +import os +import random +import sys + +# THIRD-PARTY +import numpy as np +from numpy.testing import assert_array_almost_equal + +# LOCAL +from sphere import polygon + +GRAPH_MODE = False +ROOT_DIR = os.path.join(os.path.dirname(__file__), 'data') + + +class union_test: + def __init__(self, lon_0, lat_0, proj='ortho'): + self._lon_0 = lon_0 + self._lat_0 = lat_0 + self._proj = proj + + def __call__(self, func): + @functools.wraps(func) + def run(*args, **kwargs): + polys = func(*args, **kwargs) + + unions = [] + num_permutations = math.factorial(len(polys)) + step_size = int(max(float(num_permutations) / 20.0, 1.0)) + if GRAPH_MODE: + print("%d permutations" % num_permutations) + for method in ('parallel', 'serial'): + for i, permutation in enumerate( + itertools.islice( + itertools.permutations(polys), + None, None, step_size)): + filename = '%s_%s_union_%04d.svg' % ( + func.__name__, method, i) + print(filename) + + union = polygon.SphericalPolygon.multi_union( + permutation, method=method) + unions.append(union) + areas = [x.area() for x in permutation] + union_area = union.area() + assert np.all(union_area >= areas) + + if GRAPH_MODE: + fig = plt.figure() + m = Basemap(projection=self._proj, + lon_0=self._lon_0, + lat_0=self._lat_0) + m.drawmapboundary(fill_color='white') + m.drawparallels(np.arange(-90., 90., 20.)) + m.drawmeridians(np.arange(0., 420., 20.)) + + union.draw(m, color='red', linewidth=3) + for poly in permutation: + poly.draw(m, color='blue', alpha=0.5) + plt.savefig(filename) + fig.clear() + + lengths = np.array([len(x._points) for x in unions]) + assert np.all(lengths == [lengths[0]]) + areas = np.array([x.area() for x in unions]) + # assert_array_almost_equal(areas, areas[0], 1) + + return run + + +@union_test(0, 90) +def test1(): + import pyfits + fits = pyfits.open(os.path.join(ROOT_DIR, '1904-66_TAN.fits')) + header = fits[0].header + + poly1 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[0, 87]) + poly2 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[20, 89]) + poly3 = polygon.SphericalPolygon.from_wcs( + header, 1, crval=[175, 89]) + poly4 = polygon.SphericalPolygon.from_cone( + 90, 70, 10, steps=50) + + return [poly1, poly2, poly3, poly4] + + +@union_test(0, 90) +def test2(): + poly1 = polygon.SphericalPolygon.from_cone(0, 60, 7, steps=16) + poly2 = polygon.SphericalPolygon.from_cone(0, 72, 7, steps=16) + poly3 = polygon.SphericalPolygon.from_cone(20, 60, 7, steps=16) + poly4 = polygon.SphericalPolygon.from_cone(20, 72, 7, steps=16) + poly5 = polygon.SphericalPolygon.from_cone(35, 55, 7, steps=16) + poly6 = polygon.SphericalPolygon.from_cone(60, 60, 3, steps=16) + return [poly1, poly2, poly3, poly4, poly5, poly6] + + +@union_test(0, 90) +def test3(): + random.seed(0) + polys = [] + for i in range(10): + polys.append(polygon.SphericalPolygon.from_cone( + random.randrange(-180, 180), + random.randrange(20, 90), + random.randrange(5, 16), + steps=16)) + return polys + + +@union_test(0, 15) +def test4(): + random.seed(64) + polys = [] + for i in range(10): + polys.append(polygon.SphericalPolygon.from_cone( + random.randrange(-30, 30), + random.randrange(-15, 60), + random.randrange(5, 16), + steps=16)) + return polys + + +def test5(): + import pyfits + import pywcs + + A = pyfits.open(os.path.join(ROOT_DIR, '2chipA.fits.gz')) + + wcs = pywcs.WCS(A[1].header, fobj=A) + chipA1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(A[4].header, fobj=A) + chipA2 = polygon.SphericalPolygon.from_wcs(wcs) + + null_union = chipA1.union(chipA2) + + +def test6(): + import pyfits + import pywcs + + A = pyfits.open(os.path.join(ROOT_DIR, '2chipC.fits.gz')) + + wcs = pywcs.WCS(A[1].header, fobj=A) + chipA1 = polygon.SphericalPolygon.from_wcs(wcs) + wcs = pywcs.WCS(A[4].header, fobj=A) + chipA2 = polygon.SphericalPolygon.from_wcs(wcs) + + null_union = chipA1.union(chipA2) + + +if __name__ == '__main__': + if '--profile' not in sys.argv: + GRAPH_MODE = True + from mpl_toolkits.basemap import Basemap + from matplotlib import pyplot as plt + + functions = [(k, v) for k, v in globals().items() if k.startswith('test')] + functions.sort() + for k, v in functions: + v() diff --git a/lib/test/test_util.py b/lib/test/test_util.py new file mode 100644 index 0000000..47a3444 --- /dev/null +++ b/lib/test/test_util.py @@ -0,0 +1,14 @@ +import os + +import numpy as np +from sphere import vector + +ROOT_DIR = os.path.join(os.path.dirname(__file__), 'data') + +def get_point_set(density=25): + points = [] + for i in np.linspace(-85, 85, density, True): + for j in np.linspace(-180, 180, np.cos(np.deg2rad(i)) * density): + points.append([j, i]) + points = np.asarray(points) + return np.dstack(vector.radec_to_vector(points[:,0], points[:,1]))[0] diff --git a/lib/vector.py b/lib/vector.py new file mode 100644 index 0000000..7323c92 --- /dev/null +++ b/lib/vector.py @@ -0,0 +1,213 @@ +# -*- coding: utf-8 -*- + +# Copyright (C) 2011 Association of Universities for Research in +# Astronomy (AURA) +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# +# 1. Redistributions of source code must retain the above +# copyright notice, this list of conditions and the following +# disclaimer. +# +# 2. Redistributions in binary form must reproduce the above +# copyright notice, this list of conditions and the following +# disclaimer in the documentation and/or other materials +# provided with the distribution. +# +# 3. The name of AURA and its representatives may not be used to +# endorse or promote products derived from this software without +# specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR +# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT, +# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED +# OF THE POSSIBILITY OF SUCH DAMAGE. + +""" +The `sphere.vector` module contains the basic operations for handling +vectors and converting them to and from other representations. +""" + +# THIRD-PARTY +import numpy as np + + +__all__ = ['radec_to_vector', 'vector_to_radec', 'normalize_vector', + 'rotate_around'] + + +def radec_to_vector(ra, dec, degrees=True): + ur""" + Converts a location on the unit sphere from right-ascension and + declination to an *x*, *y*, *z* vector. + + Parameters + ---------- + ra, dec : scalars or 1-D arrays + + degrees : bool, optional + + If `True`, (default) *ra* and *dec* are in decimal degrees, + otherwise in radians. + + Returns + ------- + x, y, z : tuple of scalars or 1-D arrays of the same length + + Notes + ----- + Where right-ascension is *α* and declination is *δ*: + + .. math:: + x = \cos\alpha \cos\delta + + y = \sin\alpha \cos\delta + + z = \sin\delta + """ + ra = np.asanyarray(ra) + dec = np.asanyarray(dec) + + if degrees: + ra_rad = np.deg2rad(ra) + dec_rad = np.deg2rad(dec) + else: + ra_rad = ra + dec_rad = dec + + cos_dec = np.cos(dec_rad) + + return ( + np.cos(ra_rad) * cos_dec, + np.sin(ra_rad) * cos_dec, + np.sin(dec_rad)) + + +def vector_to_radec(x, y, z, degrees=True): + ur""" + Converts a vector to right-ascension and declination. + + Parameters + ---------- + x, y, z : scalars or 1-D arrays + The input vectors + + degrees : bool, optional + If `True` (default) the result is returned in decimal degrees, + otherwise radians. + + Returns + ------- + ra, dec : tuple of scalars or arrays of the same length + + Notes + ----- + Where right-ascension is *α* and declination is + *δ*: + + .. math:: + \alpha = \arctan2(y, x) + + \delta = \arctan2(z, \sqrt{x^2 + y^2}) + """ + x = np.asanyarray(x) + y = np.asanyarray(y) + z = np.asanyarray(z) + + result = ( + np.arctan2(y, x), + np.arctan2(z, np.sqrt(x ** 2 + y ** 2))) + + if degrees: + return np.rad2deg(result[0]), np.rad2deg(result[1]) + else: + return result + + +def normalize_vector(x, y, z, inplace=False): + ur""" + Normalizes a vector so it falls on the unit sphere. + + *x*, *y*, *z* may be scalars or 1-D arrays + + Parameters + ---------- + x, y, z : scalars or 1-D arrays of the same length + The input vectors + + inplace : bool, optional + When `True`, the original arrays will be normalized in place, + otherwise a normalized copy is returned. + + Returns + ------- + X, Y, Z : scalars of 1-D arrays of the same length + The normalized output vectors + """ + x = np.asanyarray(x) + y = np.asanyarray(y) + z = np.asanyarray(z) + + l = np.sqrt(x ** 2 + y ** 2 + z ** 2) + + if inplace: + x /= l + y /= l + z /= l + return (x, y, z) + else: + return (x / l, y / l, z / l) + + +def rotate_around(x, y, z, u, v, w, theta, degrees=True): + r""" + Rotates the vector (*x*, *y*, *z*) around the arbitrary axis defined by + vector (*u*, *v*, *w*) by *theta*. + + It is assumed that both (*x*, *y*, *z*) and (*u*, *v*, *w*) are + already normalized. + + Parameters + ---------- + x, y, z : doubles + The normalized vector to rotate + + u, v, w : doubles + The normalized vector to rotate around + + theta : double, or array of doubles + The amount to rotate + + degrees : bool, optional + When `True`, *theta* is given in degrees, otherwise radians. + + Returns + ------- + X, Y, Z : doubles + The rotated vector + """ + if degrees: + theta = np.deg2rad(theta) + + u2 = u**2 + v2 = v**2 + w2 = w**2 + costheta = np.cos(theta) + sintheta = np.sin(theta) + icostheta = 1.0 - costheta + + det = (-u*x - v*y - w*z) + X = (-u*det)*icostheta + x*costheta + (-w*y + v*z)*sintheta + Y = (-v*det)*icostheta + y*costheta + ( w*x - u*z)*sintheta + Z = (-w*det)*icostheta + z*costheta + (-v*x + u*y)*sintheta + + return X, Y, Z diff --git a/setup.py b/setup.py new file mode 100755 index 0000000..b810728 --- /dev/null +++ b/setup.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python + +CONTACT = "Michael Droettboom" +EMAIL = "help@stsci.edu" +VERSION = "0.1" + +from distutils.core import setup, Extension +import sys + +try: + import numpy +except ImportError: + print("numpy must be installed to build sphere.") + print("ABORTING.") + raise + +major, minor, rest = numpy.__version__.split(".", 2) +if (int(major), int(minor)) < (1, 4): + print("numpy version 1.4 or later must be installed to build pywcs.") + print("ABORTING.") + raise ImportError + +try: + numpy_include = numpy.get_include() +except AttributeError: + numpy_include = numpy.get_numpy_include() + +extensions = [ + Extension('sphere.math_util', + ['src/math_util.c'], + include_dirs = [numpy_include], + libraries = ['m']) + ] + +setup( + name = 'sphere', + version = VERSION, + description = "Python based tools for spherical geometry", + author = CONTACT, + author_email = EMAIL, + packages = ['sphere', 'sphere.test'], + package_dir = {'sphere': 'lib', 'sphere.test': 'lib/test'}, + package_data = {'sphere.test': ['data/*.fits', 'data/*.fits.gz']}, + ext_modules = extensions + ) diff --git a/src/math_util.c b/src/math_util.c new file mode 100644 index 0000000..f67bea5 --- /dev/null +++ b/src/math_util.c @@ -0,0 +1,374 @@ +#include "Python.h" +#include "numpy/arrayobject.h" +#include "numpy/ufuncobject.h" + +#include "numpy/npy_3kcompat.h" + +/* + ***************************************************************************** + ** BASICS ** + ***************************************************************************** + */ + +typedef npy_intp intp; + +#define INIT_OUTER_LOOP_1 \ + intp dN = *dimensions++; \ + intp N_; \ + intp s0 = *steps++; + +#define INIT_OUTER_LOOP_2 \ + INIT_OUTER_LOOP_1 \ + intp s1 = *steps++; + +#define INIT_OUTER_LOOP_3 \ + INIT_OUTER_LOOP_2 \ + intp s2 = *steps++; + +#define INIT_OUTER_LOOP_4 \ + INIT_OUTER_LOOP_3 \ + intp s3 = *steps++; + +#define INIT_OUTER_LOOP_5 \ + INIT_OUTER_LOOP_4 \ + intp s4 = *steps++; + +#define BEGIN_OUTER_LOOP_3 \ + for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1, args[2] += s2) { + +#define BEGIN_OUTER_LOOP_4 \ + for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1, args[2] += s2, args[3] += s3) { + +#define BEGIN_OUTER_LOOP_5 \ + for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1, args[2] += s2, args[3] += s3, args[4] += s4) { + +#define END_OUTER_LOOP } + +static inline void +load_point(const char *in, const intp s, double* out) { + out[0] = (*(double *)in); + in += s; + out[1] = (*(double *)in); + in += s; + out[2] = (*(double *)in); +} + +static inline void +save_point(const double* in, char *out, const intp s) { + *(double *)out = in[0]; + out += s; + *(double *)out = in[1]; + out += s; + *(double *)out = in[2]; +} + +static inline void +cross(const double *A, const double *B, double *C) { + C[0] = A[1]*B[2] - A[2]*B[1]; + C[1] = A[2]*B[0] - A[0]*B[2]; + C[2] = A[0]*B[1] - A[1]*B[0]; +} + +static inline void +normalize(double *A) { + double l = sqrt(A[0]*A[0] + A[1]*A[1] + A[2]*A[2]); + A[0] /= l; + A[1] /= l; + A[2] /= l; +} + +static inline double +dot(const double *A, const double *B) { + return A[0]*B[0] + A[1]*B[1] + A[2]*B[2]; +} + +static inline double +sign(const double A) { + return (A == 0.0) ? 0.0 : ((A < 0.0) ? -1.0 : 1.0); +} + +static inline int +equals(const double *A, const double *B) { + return A[0] == B[0] && A[1] == B[1] && A[2] == B[2]; +} + +static inline void +mult(double *T, const double f) { + T[0] *= f; + T[1] *= f; + T[2] *= f; +} + +/* + ***************************************************************************** + ** UFUNC LOOPS ** + ***************************************************************************** + */ + +/*/////////////////////////////////////////////////////////////////////////// + inner1d +*/ + +char *inner1d_signature = "(i),(i)->()"; + +static void +DOUBLE_inner1d(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) +{ + INIT_OUTER_LOOP_3 + intp di = dimensions[0]; + intp i; + intp is1=steps[0], is2=steps[1]; + BEGIN_OUTER_LOOP_3 + char *ip1=args[0], *ip2=args[1], *op=args[2]; + double sum = 0; + for (i = 0; i < di; i++) { + sum += (*(double *)ip1) * (*(double *)ip2); + ip1 += is1; + ip2 += is2; + } + *(double *)op = sum; + END_OUTER_LOOP +} + +static PyUFuncGenericFunction inner1d_functions[] = { DOUBLE_inner1d }; +static void * inner1d_data[] = { (void *)NULL }; +static char inner1d_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE }; + +/*/////////////////////////////////////////////////////////////////////////// + cross +*/ + +char *cross_signature = "(i),(i)->(i)"; + +static void +DOUBLE_cross(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) +{ + double A[3]; + double B[3]; + double C[3]; + + INIT_OUTER_LOOP_3 + intp is1=steps[0], is2=steps[1], is3=steps[2]; + BEGIN_OUTER_LOOP_3 + char *ip1=args[0], *ip2=args[1], *op=args[2]; + + load_point(ip1, is1, A); + load_point(ip2, is2, B); + + cross(A, B, C); + + save_point(C, op, is3); + END_OUTER_LOOP +} + +static PyUFuncGenericFunction cross_functions[] = { DOUBLE_cross }; +static void * cross_data[] = { (void *)NULL }; +static char cross_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE }; + +/*/////////////////////////////////////////////////////////////////////////// + cross_and_norm +*/ + +char *cross_and_norm_signature = "(i),(i)->(i)"; + +/* + * This implements the function + * out[n] = sum_i { in1[n, i] * in2[n, i] }. + */ +static void +DOUBLE_cross_and_norm(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) +{ + double A[3]; + double B[3]; + double C[3]; + + INIT_OUTER_LOOP_3 + intp is1=steps[0], is2=steps[1], is3=steps[2]; + BEGIN_OUTER_LOOP_3 + char *ip1=args[0], *ip2=args[1], *op=args[2]; + + load_point(ip1, is1, A); + load_point(ip2, is2, B); + + cross(A, B, C); + normalize(C); + + save_point(C, op, is3); + END_OUTER_LOOP +} + +static PyUFuncGenericFunction cross_and_norm_functions[] = { DOUBLE_cross_and_norm }; +static void * cross_and_norm_data[] = { (void *)NULL }; +static char cross_and_norm_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE }; + +/*/////////////////////////////////////////////////////////////////////////// + intersection +*/ + +char *intersection_signature = "(i),(i),(i),(i)->(i)"; + +/* + * Finds the intersection of 2 great circle arcs AB and CD. + */ +static void +DOUBLE_intersection(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) +{ + double A[3]; + double B[3]; + double C[3]; + double D[3]; + + double ABX[3]; + double CDX[3]; + double T[3]; + double tmp[3]; + + double nans[3]; + + double s; + int match; + + nans[0] = nans[1] = nans[2] = NPY_NAN; + + INIT_OUTER_LOOP_5 + intp is1=steps[0], is2=steps[1], is3=steps[2], is4=steps[3], is5=steps[4]; + BEGIN_OUTER_LOOP_5 + char *ip1=args[0], *ip2=args[1], *ip3=args[2], *ip4=args[3], *op=args[4]; + + load_point(ip1, is1, A); + load_point(ip2, is2, B); + load_point(ip3, is3, C); + load_point(ip4, is4, D); + + match = !(equals(A, C) | equals(A, D) | equals(B, C) | equals(B, D)); + + if (match) { + cross(A, B, ABX); + cross(C, D, CDX); + cross(ABX, CDX, T); + normalize(T); + + match = 0; + cross(ABX, A, tmp); + s = sign(dot(tmp, T)); + cross(B, ABX, tmp); + if (s == sign(dot(tmp, T))) { + cross(CDX, C, tmp); + if (s == sign(dot(tmp, T))) { + cross(D, CDX, tmp); + if (s == sign(dot(tmp, T))) { + match = 1; + } + } + } + } + + if (match) { + mult(T, s); + save_point(T, op, is5); + } else { + save_point(nans, op, is5); + } + END_OUTER_LOOP +} + +static PyUFuncGenericFunction intersection_functions[] = { DOUBLE_intersection }; +static void * intersection_data[] = { (void *)NULL }; +static char intersection_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE }; + +/* + ***************************************************************************** + ** MODULE ** + ***************************************************************************** + */ + +static void +addUfuncs(PyObject *dictionary) { + PyObject *f; + + f = PyUFunc_FromFuncAndDataAndSignature( + inner1d_functions, inner1d_data, inner1d_signatures, 2, 2, 1, + PyUFunc_None, "inner1d", + "inner on the last dimension and broadcast on the rest \n" \ + " \"(i),(i)->()\" \n", + 0, inner1d_signature); + PyDict_SetItemString(dictionary, "inner1d", f); + Py_DECREF(f); + + f = PyUFunc_FromFuncAndDataAndSignature( + cross_functions, cross_data, cross_signatures, 2, 2, 1, + PyUFunc_None, "cross", + "cross product of 3-vectors only \n" \ + " \"(i),(i)->(i)\" \n", + 0, cross_signature); + PyDict_SetItemString(dictionary, "cross", f); + Py_DECREF(f); + + f = PyUFunc_FromFuncAndDataAndSignature( + cross_and_norm_functions, cross_and_norm_data, cross_and_norm_signatures, 2, 2, 1, + PyUFunc_None, "cross_and_norm", + "cross_and_norm product of 3-vectors only \n" \ + " \"(i),(i)->(i)\" \n", + 0, cross_and_norm_signature); + PyDict_SetItemString(dictionary, "cross_and_norm", f); + Py_DECREF(f); + + f = PyUFunc_FromFuncAndDataAndSignature( + intersection_functions, intersection_data, intersection_signatures, 2, 4, 1, + PyUFunc_None, "intersection", + "intersection product of 3-vectors only \n" \ + " \"(i),(i),(i),(i)->(i)\" \n", + 0, intersection_signature); + PyDict_SetItemString(dictionary, "intersection", f); + Py_DECREF(f); +} + +#if defined(NPY_PY3K) +static struct PyModuleDef moduledef = { + PyModuleDef_HEAD_INIT, + "math_util", + NULL, + -1, + NULL, + NULL, + NULL, + NULL, + NULL +}; +#endif + +#if defined(NPY_PY3K) +#define RETVAL m +PyObject *PyInit_math_util(void) +#else +#define RETVAL +PyMODINIT_FUNC +initmath_util(void) +#endif +{ + PyObject *m; + PyObject *d; + +#if defined(NPY_PY3K) + m = PyModule_Create(&moduledef); +#else + m = Py_InitModule("math_util", NULL); +#endif + if (m == NULL) + return RETVAL; + + import_array(); + import_ufunc(); + + d = PyModule_GetDict(m); + + /* Load the ufunc operators into the module's namespace */ + addUfuncs(d); + + if (PyErr_Occurred()) { + PyErr_SetString(PyExc_RuntimeError, + "cannot load umath_tests module."); + } + + return RETVAL; +} |
