Index: /tags/rel-0-1-0/pois/.cvsignore =================================================================== --- /tags/rel-0-1-0/pois/.cvsignore (revision 11633) +++ /tags/rel-0-1-0/pois/.cvsignore (revision 11633) @@ -0,0 +1,7 @@ +Makefile Makefile.in +aclocal.m4 autom4te.cache +config.* configure +depcomp install-sh libtool ltmain.sh missing mkinstalldirs compile stamp-h* +pois.pc +pois-*.tar.bz2 +pois-*.tar.gz Index: /tags/rel-0-1-0/pois/GPL.txt =================================================================== --- /tags/rel-0-1-0/pois/GPL.txt (revision 11633) +++ /tags/rel-0-1-0/pois/GPL.txt (revision 11633) @@ -0,0 +1,345 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. + 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Library General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Library General +Public License instead of this License. Index: /tags/rel-0-1-0/pois/INSTALL =================================================================== --- /tags/rel-0-1-0/pois/INSTALL (revision 11633) +++ /tags/rel-0-1-0/pois/INSTALL (revision 11633) @@ -0,0 +1,62 @@ +$ cvs co jhbuild + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +Setup $HOME/.jhbuildrc to install in (say) $HOME/JHB + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +I found that I needed to say + +$ jhbuild boostrap + +in order to get my autoconf macros recognised + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +$ jhbuild build pslib + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +Install + ac_pkg_swig.m4 ac_python_devel.m4 +into + $HOME/JHB/share/aclocal + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +To build/use the code, you'll have to source some script such as: + +export JHTREE=$HOME/JHB/Linux +export PYTHONPATH=$HOME/Python:$HOME/Pan-STARRS/Pois/python:$HOME/Pan-STARRS/Pois/swig:$JHTREE/lib +export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$JHTREE/lib +export PATH=$JHTREE/bin:$(echo $PATH | perl -pe "s|$JHTREE/bin||") +export PKG_CONFIG_PATH=$JHTREE/lib/pkgconfig:$PKG_CONFIG_PATH +hash -r + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +You need to set XPA_DIR to somewhere that the xpa libraries can be found; +or (better) install the xpa stuff in $JHTREE/{lib,include} + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +To build from scratch: + +./autogen.sh --prefix=$JHTREE + +(note that you set JHTREE above) + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +To be able to display images, ensure that ds9 and xpans are in your path (xpans +is a utility used to communicate with ds9) + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +When all is setup, say: + +$ python +>>> import pois +>>> help(pois.run) +>>> pois.run(fileDir="data", verbose=3, display={'stamps':1, 'reference':1, 'subtracted':1}) Index: /tags/rel-0-1-0/pois/Makefile.am =================================================================== --- /tags/rel-0-1-0/pois/Makefile.am (revision 11633) +++ /tags/rel-0-1-0/pois/Makefile.am (revision 11633) @@ -0,0 +1,5 @@ +SUBDIRS = src #swig + +CLEANFILES = *.pyc *~ core core.* + +EXTRA_DIST = autogen.sh Index: /tags/rel-0-1-0/pois/NOTES =================================================================== --- /tags/rel-0-1-0/pois/NOTES (revision 11633) +++ /tags/rel-0-1-0/pois/NOTES (revision 11633) @@ -0,0 +1,153 @@ +Installed + ac_python_devel.m4 + ac_pkg_swig.m4 +in ~/JHB/Tree/share/aclocal + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +Added + AC_PROG_SWIG(1.3.17) + SWIG_MULTI_MODULE_SUPPORT + SWIG_PYTHON + +... + + AC_CONFIG_FILES(pois.pc Makefile src/Makefile swig/Makefile) + +to configure.ac + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +Need to fix swig autoconf so that the configure chooses the target language[s], +but the same Makefile.am works for both + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +Initially running make fails due to: + Makefile:261: .deps/libswigrun..Po: No such file or directory + make: *** No rule to make target `.deps/libswigrun..Po'. Stop. +Changing the +include ./$(DEPDIR)/libswigrun.$(SO).Po +in Makefile (not Makefile.am) to +-include ./$(DEPDIR)/libswigrun.$(SO).Po +and running make successfully bootstrapped things. + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +cp ~/Pan-STARRS/psLib/cvs/src/parseErrorCodes.pl ~/JHB/Tree/bin + +(N.b. this is the RHL version) + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +autogen.sh --prefix=$HOME/JHB/Tree + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +Need autoconf for xpa + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +Need swig package (1.3.24?) + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +Need to add Tree/lib to PYTHONPATH (and also PERL5LIB?) + +export PYTHONPATH=$PYTHONPATH:$HOME/Pan-STARRS/Price/python:$HOME/Pan-STARRS/Price/swig:$HOME/JHB/Tree/lib + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +How do I make a .so file LOCALLY, without having to do a make install? + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +Why do things suddenly get installed as _xpa.0.so?? + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +If I try to build pois as well as libpois.a, the compiler fails to find pslib +while building pois.c + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + + +psBinaryOp isn't implemented for & and | + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +psBinaryOp &/| can be faked with a cast to long, but this will fail for some +unsigned types. Casting to unsigned will fail on some machines. Sigh. + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +psImageSubset SDRS should allow -ve row0/col0 as well as row1/col1 + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +psBinaryOp fails if e.g. data1 is U8 and the (scalar) data2 is U32. +Shouldn't it handle promotions? + +The SDRS reads: + Operations between data structures with different types (e.g., psS32 + and psF32) are not allowed and must raise an error (it is the + responsibility of calling functions to perform type conversions). + +but I think that integers should all be handled via the usual C rules +(which are only tricky for unsigned). + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +The file src/libpois.a has to exist + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +psArray* psArrayAlloc(psU32 nalloc) + +fails to initialise the array to NULLs (not specified in the SDRS) + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + +"bool psListAdd(psList *list, int location, void *data); + +the first function, psListAdd, adds an entry to the list and returns a +boolean giving the success or failure of the operation." + +(SDRS). It'd be better to define + +psList *psListAdd(psList *list, int location, void *data); +and to make + psListAdd(NULL, location, data); +equivalent to + psListAlloc(data); + +That makes the code to append to a list simpler; instead of + + list = NULL; + while(foo) { + ... + foo = psFooAlloc(); + if(list == NULL) { + list = psListAlloc(foo); + } else { + bool ok = psListAdd(list, PS_LIST_TAIL, foo); + assert(ok); + } + psFree(foo); + } + +I can simply write: + + list = NULL; + while(foo) { + ... + foo = psFooAlloc(); + psListAdd(list, PS_LIST_TAIL, foo); + psFree(foo); + } + +-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + Index: /tags/rel-0-1-0/pois/autogen.sh =================================================================== --- /tags/rel-0-1-0/pois/autogen.sh (revision 11633) +++ /tags/rel-0-1-0/pois/autogen.sh (revision 11633) @@ -0,0 +1,103 @@ +#!/bin/sh +# Run this to generate all the initial makefiles, etc. + +srcdir=`dirname $0` +test -z "$srcdir" && srcdir=. + +ORIGDIR=`pwd` +cd $srcdir + +PROJECT=pois +TEST_TYPE=-f +# change this to be a unique filename in the top level dir +FILE=autogen.sh + +DIE=0 + +LIBTOOLIZE=libtoolize +ACLOCAL="aclocal $ACLOCAL_FLAGS" +AUTOHEADER=autoheader +AUTOMAKE=automake +AUTOCONF=autoconf + +($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $LIBTOOlIZE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/" + DIE=1 +} + +($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $ACLOCAL installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOHEADER installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOMAKE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOCONF installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +if test "$DIE" -eq 1; then + exit 1 +fi + +test $TEST_TYPE $FILE || { + echo "You must run this script in the top-level $PROJECT directory" + exit 1 +} + +if test -z "$*"; then + echo "I am going to run ./configure with no arguments - if you wish " + echo "to pass any to it, please specify them on the $0 command line." +fi + +$LIBTOOLIZE --copy --force || echo "$LIBTOOlIZE failed" +$ACLOCAL -I m4 || echo "$ACLOCAL failed" +$AUTOHEADER || echo "$AUTOHEADER failed" +$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE failed" +$AUTOCONF || echo "$AUTOCONF failed" + +cd $ORIGDIR + +run_configure=true +for arg in $*; do + case $arg in + --no-configure) + run_configure=false + ;; + *) + ;; + esac +done + +if $run_configure; then + $srcdir/configure --enable-maintainer-mode "$@" + echo + echo "Now type 'make' to compile $PROJECT." +else + echo + echo "Now run 'configure' and 'make' to compile $PROJECT." +fi Index: /tags/rel-0-1-0/pois/configure.ac =================================================================== --- /tags/rel-0-1-0/pois/configure.ac (revision 11633) +++ /tags/rel-0-1-0/pois/configure.ac (revision 11633) @@ -0,0 +1,48 @@ +dnl Process this file with autoconf to produce a configure script. +AC_PREREQ(2.59) + +AC_INIT([pois], [0.1.0], [ipp-support@ifa.hawaii.edu]) +AC_CONFIG_SRCDIR([autogen.sh]) + +AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2]) +AM_CONFIG_HEADER([src/config.h]) +AM_MAINTAINER_MODE + +IPP_STDCFLAGS + +AC_LANG(C) +AC_PROG_CC +AC_PROG_INSTALL +AC_PROG_LIBTOOL +AC_SYS_LARGEFILE + +PKG_CHECK_MODULES([PSLIB], [pslib >= 1.0.0]) + +AC_C_BIGENDIAN([PS_BIGENDIAN=1],[PS_BIGENDIAN=0]) + +dnl AC_PROG_SWIG(1.3.24) +dnl SWIG_MULTI_MODULE_SUPPORT +dnl SWIG_PYTHON + +dnl lifted from autoconf 2.59's configure.ac +AC_PATH_PROG([PERL], perl, no) +AC_SUBST([PERL]) +if test "$PERL" = no; then + AC_MSG_ERROR([perl is not found]) +fi +$PERL -e 'require 5.005_03;' || { + AC_MSG_ERROR([Perl 5.005_03 or better is required]) +} + +dnl Set CFLAGS for build +IPP_STDOPTS +CFLAGS="${CFLAGS=} -Wall -Werror -std=c99" + +AC_CONFIG_FILES([ + Makefile + src/Makefile +]) +dnl swig/Makefile +dnl python/Makefile + +AC_OUTPUT Index: /tags/rel-0-1-0/pois/m4/ac_pkg_swig.m4 =================================================================== --- /tags/rel-0-1-0/pois/m4/ac_pkg_swig.m4 (revision 11633) +++ /tags/rel-0-1-0/pois/m4/ac_pkg_swig.m4 (revision 11633) @@ -0,0 +1,160 @@ +dnl @synopsis AC_PROG_SWIG([major.minor.micro]) +dnl +dnl This macro searches for a SWIG installation on your system. If +dnl found you should call SWIG via $(SWIG). You can use the optional +dnl first argument to check if the version of the available SWIG is +dnl greater than or equal to the value of the argument. It should have +dnl the format: N[.N[.N]] (N is a number between 0 and 999. Only the +dnl first N is mandatory.) +dnl +dnl If the version argument is given (e.g. 1.3.17), AC_PROG_SWIG checks +dnl that the swig package is this version number or higher. +dnl +dnl In configure.in, use as: +dnl +dnl AC_PROG_SWIG(1.3.17) +dnl SWIG_ENABLE_CXX +dnl SWIG_MULTI_MODULE_SUPPORT +dnl SWIG_PYTHON +dnl +dnl @category InstalledPackages +dnl @author Sebastian Huber +dnl @author Alan W. Irwin +dnl @author Rafael Laboissiere +dnl @author Andrew Collier +dnl @version 2004-09-20 +dnl @license GPLWithACException + +AC_DEFUN([AC_PROG_SWIG],[ + AC_PATH_PROG([SWIG],[swig]) + if test -z "$SWIG" ; then + AC_MSG_WARN([cannot find 'swig' program. You should look at http://www.swig.org]) + SWIG='echo "Error: SWIG is not installed. You should look at http://www.swig.org" ; false' + elif test -n "$1" ; then + AC_MSG_CHECKING([for SWIG version]) + [swig_version=`$SWIG -version 2>&1 | grep 'SWIG Version' | sed 's/.*\([0-9][0-9]*\.[0-9][0-9]*\.[0-9][0-9]*\).*/\1/g'`] + AC_MSG_RESULT([$swig_version]) + if test -n "$swig_version" ; then + # Calculate the required version number components + [required=$1] + [required_major=`echo $required | sed 's/[^0-9].*//'`] + if test -z "$required_major" ; then + [required_major=0] + fi + [required=`echo $required | sed 's/[0-9]*[^0-9]//'`] + [required_minor=`echo $required | sed 's/[^0-9].*//'`] + if test -z "$required_minor" ; then + [required_minor=0] + fi + [required=`echo $required | sed 's/[0-9]*[^0-9]//'`] + [required_patch=`echo $required | sed 's/[^0-9].*//'`] + if test -z "$required_patch" ; then + [required_patch=0] + fi + # Calculate the available version number components + [available=$swig_version] + [available_major=`echo $available | sed 's/[^0-9].*//'`] + if test -z "$available_major" ; then + [available_major=0] + fi + [available=`echo $available | sed 's/[0-9]*[^0-9]//'`] + [available_minor=`echo $available | sed 's/[^0-9].*//'`] + if test -z "$available_minor" ; then + [available_minor=0] + fi + [available=`echo $available | sed 's/[0-9]*[^0-9]//'`] + [available_patch=`echo $available | sed 's/[^0-9].*//'`] + if test -z "$available_patch" ; then + [available_patch=0] + fi + if test $available_major -ne $required_major \ + -o $available_minor -ne $required_minor \ + -o $available_patch -lt $required_patch ; then + AC_MSG_WARN([SWIG version >= $1 is required. You have $swig_version. You should look at http://www.swig.org]) + SWIG='echo "Error: SWIG version >= $1 is required. You have '"$swig_version"'. You should look at http://www.swig.org" ; false' + else + AC_MSG_NOTICE([SWIG executable is '$SWIG']) + SWIG_LIB=`$SWIG -swiglib` + AC_MSG_NOTICE([SWIG library directory is '$SWIG_LIB']) + fi + else + AC_MSG_WARN([cannot determine SWIG version]) + SWIG='echo "Error: Cannot determine SWIG version. You should look at http://www.swig.org" ; false' + fi + fi + AC_SUBST([SWIG_LIB]) +]) + +# SWIG_ENABLE_CXX() +# +# Enable SWIG C++ support. This affects all invocations of $(SWIG). +AC_DEFUN([SWIG_ENABLE_CXX],[ + AC_REQUIRE([AC_PROG_SWIG]) + AC_REQUIRE([AC_PROG_CXX]) + SWIG="$SWIG -c++" +]) + +# SWIG_MULTI_MODULE_SUPPORT() +# +# Enable support for multiple modules. This effects all invocations +# of $(SWIG). You have to link all generated modules against the +# appropriate SWIG runtime library. If you want to build Python +# modules for example, use the SWIG_PYTHON() macro and link the +# modules against $(SWIG_LANG_LIBS). +# +AC_DEFUN([SWIG_MULTI_MODULE_SUPPORT],[ + AC_REQUIRE([AC_PROG_SWIG]) + #SWIG="$SWIG -noruntime" # Not needed for SWIG >= 1.3.24 +]) + +# SWIG_PYTHON([use-shadow-classes = {no, yes}]) +# +# Checks for Python and provides the $(SWIG_PYTHON_CPPFLAGS), +# and $(SWIG_PYTHON_OPT) output variables. +# +# $(SWIG_PYTHON_OPT) contains all necessary SWIG options to generate +# code for Python. Shadow classes are enabled unless the value of the +# optional first argument is exactly 'no'. If you need multi module +# support (provided by the SWIG_MULTI_MODULE_SUPPORT() macro) use +# $(SWIG_PYTHON_LIBS) to link against the appropriate library. It +# contains the SWIG Python runtime library that is needed by the type +# check system for example. +AC_DEFUN([SWIG_PYTHON],[ + AC_REQUIRE([AC_PROG_SWIG]) + AC_REQUIRE([AC_PYTHON_DEVEL]) + test "x$1" != "xno" || swig_shadow=" -noproxy" + AC_SUBST([SWIG_PYTHON_OPT],[-python$swig_shadow]) + AC_SUBST([SWIG_PYTHON_CPPFLAGS],[$PYTHON_CPPFLAGS]) + AC_SUBST([SWIG_PYTHON_LIBS],[$PYTHON_LIBS]) + + AC_SUBST([SWIG_TARGET_OPT], [$SWIG_PYTHON_OPT]) + AC_SUBST([SWIG_TARGET_CPPFLAGS],[$SWIG_PYTHON_CPPFLAGS]) + AC_SUBST([SWIG_TARGET_LIBS], [$SWIG_PYTHON_LIBS]) +]) + + +dnl @synopsis AC_LIB_WAD +dnl +dnl This macro searches for installed WAD library. +dnl +AC_DEFUN([AC_LIB_WAD], +[ + AC_REQUIRE([AC_PYTHON_DEVEL]) + AC_ARG_ENABLE(wad, + AC_HELP_STRING([--enable-wad], [enable wad module]), + [ + case "${enableval}" in + no) ;; + *) if test "x${enableval}" = xyes; + then + check_wad="yes" + fi ;; + esac + ], []) + + if test -n "$check_wad"; + then + AC_CHECK_LIB(wadpy, _init, [WADPY=-lwadpy], [], $PYTHON_LDFLAGS $PYTHON_EXTRA_LIBS) + AC_SUBST(WADPY) + fi +]) Index: /tags/rel-0-1-0/pois/m4/ac_python_devel.m4 =================================================================== --- /tags/rel-0-1-0/pois/m4/ac_python_devel.m4 (revision 11633) +++ /tags/rel-0-1-0/pois/m4/ac_python_devel.m4 (revision 11633) @@ -0,0 +1,64 @@ +dnl @synopsis AC_PYTHON_DEVEL +dnl +dnl Checks for Python and tries to get the include path to 'Python.h'. +dnl It provides the $(PYTHON_CPPFLAGS) and $(PYTHON_LDFLAGS) output +dnl variable. +dnl +dnl @category InstalledPackages +dnl @author Sebastian Huber +dnl @author Alan W. Irwin +dnl @author Rafael Laboissiere +dnl @author Andrew Collier +dnl @version 2004-07-14 +dnl @license GPLWithACException + +AC_DEFUN([AC_PYTHON_DEVEL],[ + # + # should allow for checking of python version here... + # + AC_REQUIRE([AM_PATH_PYTHON]) + + # Check for Python include path + AC_MSG_CHECKING([for Python include path]) + python_path=`echo $PYTHON | sed "s,/bin.*$,,"` + for i in "$python_path/include/python$PYTHON_VERSION/" "$python_path/include/python/" "$python_path/" ; do + python_path=`find $i -type f -name Python.h -print | sed "1q"` + if test -n "$python_path" ; then + break + fi + done + python_path=`echo $python_path | sed "s,/Python.h$,,"` + AC_MSG_RESULT([$python_path]) + if test -z "$python_path" ; then + AC_MSG_ERROR([cannot find Python include path]) + fi + AC_SUBST([PYTHON_CPPFLAGS],[-I$python_path]) + + # Check for Python library path + AC_MSG_CHECKING([for Python library path]) + python_path=`echo $PYTHON | sed "s,/bin.*$,,"` + for i in "$python_path/lib/python$PYTHON_VERSION/config/" "$python_path/lib/python$PYTHON_VERSION/" "$python_path/lib/python/config/" "$python_path/lib/python/" "$python_path/" ; do + python_path=`find $i -type f -name libpython$PYTHON_VERSION.* -print | sed "1q"` + if test -n "$python_path" ; then + break + fi + done + python_path=`echo $python_path | sed "s,/libpython.*$,,"` + AC_MSG_RESULT([$python_path]) + if test -z "$python_path" ; then + AC_MSG_ERROR([cannot find Python library path]) + fi + AC_SUBST([PYTHON_LDFLAGS],["-L$python_path -lpython$PYTHON_VERSION"]) + # + python_site=`echo $python_path | sed "s/config/site-packages/"` + AC_SUBST([PYTHON_SITE_PKG],[$python_site]) + # + # libraries which must be linked in when embedding + # + AC_MSG_CHECKING(python extra libraries) + PYTHON_EXTRA_LIBS=`$PYTHON -c "import distutils.sysconfig; \ + conf = distutils.sysconfig.get_config_var; \ + print conf('LOCALMODLIBS')+' '+conf('LIBS')" + AC_MSG_RESULT($PYTHON_EXTRA_LIBS)` + AC_SUBST(PYTHON_EXTRA_LIBS) +]) Index: /tags/rel-0-1-0/pois/python/.cvsignore =================================================================== --- /tags/rel-0-1-0/pois/python/.cvsignore (revision 11633) +++ /tags/rel-0-1-0/pois/python/.cvsignore (revision 11633) @@ -0,0 +1,2 @@ +Makefile Makefile.in +*.pyc Index: /tags/rel-0-1-0/pois/python/ds9.py =================================================================== --- /tags/rel-0-1-0/pois/python/ds9.py (revision 11633) +++ /tags/rel-0-1-0/pois/python/ds9.py (revision 11633) @@ -0,0 +1,163 @@ +# +# XPA +# +import os, re, math, sys, time + +try: import xpaSwig as xpa +except: print "Cannot import xpa" + +class Ds9Error(IOError): + """Some problem talking to ds9""" + +def ds9Cmd(cmd): + """Issue a ds9 command, raising errors as appropriate""" + + try: + xpa.XPASet(None, "ds9", cmd, "", "", 0) + except IOError: + raise Ds9Error, "XPA: (%s)" % cmd + +def initDS9(execDs9 = 1): + try: + ds9Cmd("iconify no; raise") + except IOError: + if execDs9: + print "ds9 doesn't appear to be running, I'll exec it for you" + if not re.search('xpa', os.environ['PATH']): + raise Ds9Error, 'You need the xpa binaries in your path to use ds9 with python' + + os.system('ds9 &') + for i in range(0,10): + try: + ds9Cmd("frame 0; scale histequ; scale mode minmax") + return + except IOError: + print "waiting for ds9...\r", + time.sleep(0.5) + else: + break + + print " \r", + sys.stdout.flush(); + + raise Ds9Error + +def mtv(data, meta = None, frame = 0, init = 1, WCS = ""): + """Display an OTA or Cell on a DS9 display""" + + if frame == None: + return + + if init: + for i in range(0,3): + try: + initDS9(i == 0) + except IOError: + print "waiting for ds9...\r", ; sys.stdout.flush(); + time.sleep(0.5) + else: + break + + ds9Cmd("frame %d" % frame) + + if 1: + xpa_cmd = "xpaset ds9 fits" + pfd = os.popen(xpa_cmd, "w") + else: + pfd = file("foo.fits", "w") + + xpa.rhlWriteFits(pfd.fileno(), data, meta, WCS) + + try: + pass + except: + print "Error" + pass + + try: + pfd.close() + except: + pass +# +# Graphics commands +# +def erase(frame = 0): + """Erase the specified DS9 frame""" + if frame == None: + return + + ds9Cmd("frame %d; regions delete all" % frame) + +def dot(symb, c, r, frame = 0, size = 2, ctype = 'green'): + """Draw a symbol onto the specfied DS9 frame at (row,col) = (r,c) [0-based coordinates] +Possible values are: + + Draw a + + x Draw an x + o Draw a circle +Any other value is interpreted as a string to be drawn +""" + if frame == None: + return + + cmd = "frame %d; regions physical; " % frame + r += 1; c += 1; # ds9 uses 1-based coordinates + if (symb == '+'): + cmd += 'regions line %g %g %g %g # color=%s; ' % (c, r+size, c, r-size, ctype) + cmd += 'regions line %g %g %g %g ' % (c-size, r, c+size, r) + elif (symb == 'x'): + size = size/math.sqrt(2) + cmd += 'regions line %g %g %g %g # color=%s; ; ' % (c+size, r+size, c-size, r-size, ctype) + cmd += 'regions line %g %g %g %g ' % (c-size, r+size, c+size, r-size) + elif (symb == 'o'): + cmd += 'regions circle %g %g %g ' % (c, r, size) + else: + cmd += 'regions text %g %g \"%s\"' % (c, r, symb) + + cmd += ' # color=%s' % ctype + + ds9Cmd(cmd) + +def line(points, frame = 0, symbs = False, ctype = 'green'): + """Draw a set of symbols or connect the points, a list of (row,col) +If symbs is True, draw points at the specified points using the desired symbol, +otherwise connect the dots. Ctype is the name of a colour (e.g. 'red')""" + + if frame == None: + return + + if symbs: + for (r, c) in points: + dot(symbs, c, r, frame = frame, size = 0.5, ctype = ctype) + else: + cmd = "frame %d; regions image; regions line " % (frame) + + for (r, c) in points: + r += 1; c += 1; # ds9 uses 1-based coordinates + cmd += '%g %g ' % (c, r) + + cmd += ' # color=%s' % ctype + + ds9Cmd(cmd) +# +# Zoom and Pan +# +def zoom(zoomfac = None, rowc = None, colc = None, frame = 0): + """Zoom frame by specified amount, optionally panning also""" + + if frame == None: + return + + if (rowc and not colc) or (not rowc and colc): + raise Ds9Error, "Please specify row and column center to pan about" + + if zoomfac == None and rowc == None: + zoomfac = 2 + + cmd = "" + if zoomfac != None: + cmd += "zoom to %d; " % zoomfac + + if rowc != None: + cmd += "pan to %g %g physical; " % (colc + 1, rowc + 1) # ds9 is 1-indexed. Grrr + + ds9Cmd(cmd) Index: /tags/rel-0-1-0/pois/python/pois.py =================================================================== --- /tags/rel-0-1-0/pois/python/pois.py (revision 11633) +++ /tags/rel-0-1-0/pois/python/pois.py (revision 11633) @@ -0,0 +1,360 @@ +# +# POIS: Pan-STARRS (or Pricey's) Optimal Image Subtraction +# +import os, re, sys, time +import ds9 +import ps +import psLibSwig as psLib; from psLibSwig import psTraceSetLevel, psTrace +import poisSwig as pois + +def getTime(): + return time.clock() + +def displayImage(image, frame = 0, caption = None, ctype = "green", init = False): + global init_display + if init_display: + init = True + + ds9.mtv(image, init = init, frame = frame) + if caption: + ds9.dot(caption, 10, image.numRows + 5, ctype = ctype, frame = frame) + + init_display = False + +def run(showLeaks = True, *args, **nargs): + #"""Wrapper for run, with a check on memory leaks after out-of-scope + #objects have been garbage collected.""" + # + # Only look for leaks of memory allocated after this point + # + id0 = psLib.psMemGetId() + # + # Do the work + # + value = _run(*args, **nargs) + # + # Look for otherwise-undetected problems + # + ps.raisePsError() + + nleak = psLib.psMemCheckLeaks(id0, None, None, False) + if nleak: + if showLeaks != False: + ps.memCheckLeaks(sys.stderr, id0 = id0, addresses = True) + + raise IOError, "Memory Leak (%d pointers); id0 = %d" % (nleak, id0) + + return value + +def _run(fileDir = ".", refFile = "test1.fits", inFile = "test2.fits", outFile=None, + section = "[100:500,100:500]", + testing = False, verbose = 0, psLib_verbose = 0, + display = {}, setTrace = True): + """POIS: Pan-STARRS (or Pricey's) Optimal Image Subtraction. + +Attempts to fit for a general convolution kernel that will match the PSFs between two images. +Arguments: + fileDir Directory for input and output files (default: .) + refFile Reference image (default: test1.fits) + inFile File to be analysed (default: test2.fits) + outFile Write output to this file (default: none) + section Section of datafiles to process, or None (=> all) (default: [100:500,100:500]) + + testing Write some extra test information/files (default: False) + verbose Verbosity passed to psTrace (default: 0).0, psLib_verbose = 0, + setTrace Should the code fiddle with trace levels? (default: True) + + display Dictionary of things that can be displayed; use display="help" for help (default: {}) + +E.g. + pois.run(fileDir="data", verbose=3, display={'stamps':1, 'reference':1, 'subtracted':1}) +""" + + startTime = getTime(); + + if setTrace: + psTraceSetLevel(".", 0) + psTraceSetLevel("pois", verbose) + psTraceSetLevel("pois.config", verbose) + psTraceSetLevel("pois.time", verbose + 2) + psTraceSetLevel("pois.solution", 0) + psTraceSetLevel("pois.kernelBasisFunctions", verbose) + psTraceSetLevel("pois.calculateEquation", verbose) + psTraceSetLevel("pois.convolveImage", verbose) + psTraceSetLevel("pois.findStamps", verbose) + psTraceSetLevel("pois.parseConfig", verbose) + psTraceSetLevel("pois.makeMask", verbose) + psTraceSetLevel("pois.extractKernel", verbose) + psTraceSetLevel("pois.calculateDeviations", verbose) + + psTraceSetLevel("ps.utils", psLib_verbose) + + if verbose > 3: + print "Trace levels:" + psLib.psTracePrintLevels() + + # Set logging level + psLib.psLogSetDestination("dest:stderr") + + # + # Set missing keys in display + keys = ['all', 'pause'] # permitted keywords in display + subkeys = ['convolved', 'input', 'kernel', 'mask', 'reference', 'stamps', 'subtracted'] + + if display == "help": + print """Usage: + display={'key1':1, 'key2':True, 'key3':0, ...} (1==True, 0==False) +Keys: + all Display everything that can be displayed + pause Pause after every iteration over the stamp list + + convolved Show the kernel-convolved reference image + input Show the input image + kernel Show the convolution kernel + mask Show the image mask + reference Show the reference image + stamps Overlay the reference image with the centroids of stamp stars (circle bad ones in red) + subtracted Show the subtracted image + """ + return + + for k in (keys + subkeys): + if not k in display: + display[k] = False + + # Are illegal keys present in display[]? + illegal = [] + for k in display.keys(): + if not k in keys + subkeys: + illegal.append(k) + + if illegal: + raise ("Unrecognised keywords %s are present in display[]" % (illegal.__str__())) + + # Deal with display values + if display['all']: + for k in subkeys: + display[k] = True + + if display['stamps']: + display['input'] = True + + global init_display; + init_display = True # have we initialised the display? + + # Set desired variables from config file + if not re.search("/$", fileDir): + fileDir += "/" + + config = pois.poisConfigAlloc() # Configuration values + config.inFile = fileDir + inFile + config.refFile = fileDir + refFile + try: + config.outFile = fileDir + outFile + except: + config.outFile = outFile + + # Read inputs + imageRegion = psLib.psRegionAlloc(0, 0, 0, 0) + + fd = psLib.psFitsAlloc(config.refFile) + refImage = psLib.psFitsReadImage(None, fd, imageRegion, 0) + del fd + psTrace("pois", 3, ("Read %s: %dx%d\n" % (config.refFile, refImage.numCols, refImage.numRows))) + + fd = psLib.psFitsAlloc(config.inFile) + inImage = psLib.psFitsReadImage(None, fd, imageRegion, 0) + del fd + psTrace("pois", 3, ("Read %s: %dx%d\n" % (config.inFile, inImage.numCols, inImage.numRows))) + + del imageRegion + + if section: + inImage = psLib.psImageSubsection(inImage, section) + refImage = psLib.psImageSubsection(refImage, section) + + if display['reference']: + displayImage(refImage, 1, "Reference Image", "red") + if display['input']: + displayImage(inImage, 2, "Input Image", "red") + + if (refImage.numCols != inImage.numCols) or (refImage.numRows != inImage.numRows): + raise IOError, \ + ("Reference and input images must have the same dimensions: %dx%d vs %dx%d\n" % + (refImage.numCols, refImage.numRows, inImage.numCols, inImage.numRows)) + + config.xImage = refImage.numCols + config.yImage = refImage.numRows + + psTrace("pois.time", 1, ("Inputs read at %f sec\n" % (getTime() - startTime))) + + # Add a background offset to get the image all above zero + # An image that goes negative (e.g., a background-subtracted image) produces + # horrible results --- the matrix goes NaN because the noise is determined + # as the square root of the value. + if config.background != 0.0: + psBinaryOp(refImage, refImage, "+", psScalarAlloc(config.background, PS_TYPE_F32)) + psBinaryOp(inImage, inImage, "+", psScalarAlloc(config.background, PS_TYPE_F32)) + + # Generate the kernel parameters + kernelBasisFunctions = pois.poisKernelBasisFunctions(config) + + # Generate a mask + mask = pois.poisMakeMask(None, refImage, inImage, config) + + stamps = None # Array of stamps + stampMask = None # Mask for stamps + solution = None # Solution vector + + psLib.psMemCheckCorruption(1) + + # Iterate for a good solution + for iterNum in range(0, config.numIter): + psTrace("pois", 1, ("Iteration %d...\n" % (iterNum))) + + # Find stamps + stamps = pois.poisFindStamps(stamps, refImage, mask, config); + + numStamps = 0; + for s in range(0, stamps.n): + stamp = pois.poisStamp_Cast(stamps.get_data(s)) # Stamp of interest + if stamp: + if stamp.status == pois.POIS_STAMP_BAD: + ctype = "red"; dotType = "x" + else: + ctype = "green"; dotType = "+" + numStamps += 1 + + if display['stamps']: + ds9.dot(dotType, stamp.x, stamp.y, size = 5, ctype = ctype, frame = 1) + + psTrace("pois.time", 1, ("%d stamps found at %f sec\n" % (numStamps, getTime() - startTime))) + if numStamps == 0: + raise IOError, "No stamps found. Check input parameters" + + # Calculate equation + pois.poisCalculateEquation(stamps, refImage, inImage, kernelBasisFunctions, config) + psTrace("pois.time", 1, ("Equation calculated at %f sec\n" % (getTime() - startTime))) + + # Solve the equation + solution = pois.poisSolveEquation(solution, stamps, config) + psTrace("pois.time", 1, ("Matrix equation solved at %f sec\n" % (getTime() - startTime))) + + if testing: + # Print solution + print "Solution:" + for i in range(0, solution.n - 1): + kernel = kernelBasisFunctions.get_data(i) # The i-th kernel basis function + kernel = pois.poisKernelBasis_Cast(kernel) + print "%d: (%d,%d) x^%d y^%d - %f" % \ + (i, kernel.u, kernel.v, kernel.xOrder, kernel.yOrder, solution.get_data_F64()[i]) + + + if testing or display['kernel']: + kernelImage = pois.poisExtractKernel(solution, kernelBasisFunctions, 0.0, 0.0, config) + + if display['kernel']: + displayImage(kernelImage, 5, "Kernel") + + if testing: + # Save the kernel postage stamp + kernelName = "kernel%d.fits" % iterNum + kernelFile = psLib.psFitsAlloc(kernelName) + psLib.psFitsWriteImage(kernelFile, None, kernelImage, 0, None) + psTrace("pois", 2, ("Kernel written to %s\n" %kernelName)) + + # Calculate deviations for the stamps + deviations = pois.poisCalculateDeviations(None, stamps, refImage, inImage, mask, + kernelBasisFunctions, solution, config) + # Have we converged? + rejected = pois.poisRejectStamps(stamps, mask, deviations, config) + + if rejected == None: + break + + if display['stamps']: + iter = psLib.psListIteratorAlloc(rejected, psLib.PS_LIST_HEAD, False) + while True: + stamp = pois.poisStamp_Cast(psLib.psListGetAndIncrement(iter)) + if not stamp: + break + ds9.dot("o", stamp.x, stamp.y, size = 5, ctype = "red", frame = 1) + + if display['pause']: + if raw_input("Iteration %d completed: continue? (q to abort) " % iterNum) == "q": + raise "XXX" + + # If there was rejection on the last round, re-solve the equation using only the good stamps + if iterNum == config.numIter: + psTrace("pois", 2, ("Failed to converge in %d iterations\n" % (config.numIter))) + solution = pois.poisSolveEquation(solution, stamps, config); + + # Clear the mask for bad stamp areas + if display['mask']: + displayImage(mask, 0, "Stamps mask") + + psLib.psBinaryOp(mask, mask, "&", psLib.psScalarAlloc(~pois.POIS_MASK_STAMP, psLib.PS_TYPE_U8)) + + # Convolve the input image + convImage = pois.poisConvolveImage(refImage, mask, solution, kernelBasisFunctions, config) + psTrace("pois.time", 1, ("Convolution completed at %f sec\n" % (getTime() - startTime))) + + if testing and config.outFile: + # Save the convolved image + convName = "%s.conv" % config.outFile + convFile = psLib.psFitsAlloc(convName) + psLib.psFitsWriteImage(convFile, None, convImage, 0, None) + psTrace("pois", 2, ("Convolved image written to %s\n" % (convName))) + + if display['convolved']: + displayImage(convImage, 4, "Convolved image", "red") + + # Do the subtraction + subImage = psLib.psImageAlloc(config.xImage, config.yImage, psLib.PS_TYPE_F32) + psLib.psBinaryOp(subImage, inImage, "-", convImage) + psTrace("pois.time", 1, ("Subtraction completed at %f sec\n" % (getTime() - startTime))) + + if testing: + stats = psLib.psStatsAlloc(psLib.PS_STAT_SAMPLE_MEAN | psLib.PS_STAT_SAMPLE_STDEV) + psLib.psImageStats(stats, subImage, mask, 0) + print "Subtracted image statistics: %f %f\n" % (stats.sampleMean, stats.sampleStdev) + #psFree(stats); + + # Mask out the edges + X0 = config.xKernel; X1 = mask.numCols - config.xKernel # non-masked area is (X0,Y0) -- (X1-1, Y1-1) + Y0 = config.yKernel; Y1 = mask.numRows - config.yKernel + + tmp = psLib.psImageSubset(mask, 0, 0, mask.numCols, Y0) + psLib.psBinaryOp(tmp, tmp, "|", psLib.psScalarAlloc(pois.POIS_MASK_BAD, psLib.PS_TYPE_U8)) + + tmp = psLib.psImageSubset(mask, 0, Y1, mask.numCols, mask.numRows) + psLib.psBinaryOp(tmp, tmp, "|", psLib.psScalarAlloc(pois.POIS_MASK_BAD, psLib.PS_TYPE_U8)) + + tmp = psLib.psImageSubset(mask, 0, Y0, X0, Y1) + psLib.psBinaryOp(tmp, tmp, "|", psLib.psScalarAlloc(pois.POIS_MASK_BAD, psLib.PS_TYPE_U8)) + + tmp = psLib.psImageSubset(mask, X1, Y0, mask.numCols, Y1) + psLib.psBinaryOp(tmp, tmp, "|", psLib.psScalarAlloc(pois.POIS_MASK_BAD, psLib.PS_TYPE_U8)) + del tmp + + pois.poisImageSetValInMask(subImage, subImage, mask, 0, pois.POIS_MASK_BAD) + + if False: + for y in range(config.yKernel, subImage.numRows - config.yKernel): + for x in range(config.xKernel, subImage.numCols - config.xKernel): + if mask.data.U8[y][x]: + psLib.subImage.data.F32[y][x] = 0.0 + + + if display['subtracted']: + displayImage(subImage, 3, "Subtracted image", "red") + + if config.outFile: + subFile = psLib.psFitsAlloc(config.outFile) + psLib.psFitsWriteImage(subFile, None, subImage, 0, None) + psTrace("pois", 2, ("Subtracted image written to %s\n" % (config.outFile))) + + psTrace("pois.time", 1, ("I/O completed at %f sec\n" % (getTime() - startTime))) + +if not run.__doc__: + run.__doc__ = _run.__doc__ Index: /tags/rel-0-1-0/pois/python/ps.py =================================================================== --- /tags/rel-0-1-0/pois/python/ps.py (revision 11633) +++ /tags/rel-0-1-0/pois/python/ps.py (revision 11633) @@ -0,0 +1,151 @@ +import psLibSwig as psLib +import sys + +def logcntl(what = "NM"): + """Control logging; disable if what == \"\"""" + + if what == "": + psLib.psLogSetDestination("none") # turn off logging + else: + psLib.psLogSetDestination("dest:stderr") + psLib.psLogSetFormat(what) + +logcntl("NM") + +class PsError(StandardError): + """An error signalled via a call to psError""" + +def raisePsError(clear = True): + """Raise the most recently signalled psError (the psError stack is cleared if clear is True)""" + + err = psLib.psErrorLast() + if err.code != 0: + if clear: # clear the error? + psLib.psErrorClear() + raise PsError, "psError was called without an error return: (%s) %s" % (err.name, err.msg) + +def errStack(): + """Print the current error stack""" + + i = 0; + while 1: + err = psLib.psErrorGetForSwig(i) + if err == None : + return + + errString = psLib.psErrorCodeString(err.code) + if not errString: + errString = "(unknown code: %d)" % err.code + print "%s\t%-40s %s" %(err.name, errString + ':', err.msg) + i += 1 + +def printMD(meta, match=None, file=sys.stdout, hdr=""): + """Print a set of psMetadata""" + + for i in range(0,meta.list.size): + md = psLib.psMetadataGet(meta, i) + + if md.type == psLib.PS_META_STR: + fmt = "%s" + elif md.type == psLib.PS_META_S32: + fmt = "%d" + elif md.type == psLib.PS_META_F32 or md.type == psLib.PS_META_F64: + fmt = "%g" + psLib.psMetadataItemPrint(file, md.name + " " + fmt + "\n", md) + +def memCheckLeaks(fd = None, id0 = 0, getList = 0, unique = 0, addresses = 0, persistent = False): + """Check for memory leaks. + + If file (an fd or a filename) is provided, write the list of + leaked pointers to it. If id0 is specified, only check for leaked + memory allocated after that memory id. + + If getList is true, return a list of tuples for each leaked block: + (isMagic0, id, file, lineNo, refCounter, isMagic) + + If unique is true, print/return a list of places where memory is leaking, along with + the number of pointers leaked at each point. If addresses is true, include the addresses + of the pointers + """ + + if (unique and addresses): + raise "Please choose unique _or_ addresses" + + if isinstance(fd, str): + fd = file(fd, "w") + + if addresses or getList or unique: + arr = [] + else: + arr = None + + if addresses or unique: + nleak = psLib.psMemCheckLeaks(id0, arr, None, persistent) + else: + nleak = psLib.psMemCheckLeaks(id0, arr, fd, persistent) + + if fd == None and not (addresses or unique): + if nleak == 0: + print "No leaks" + else: + print "%d leaked pointers" % (nleak) + + if addresses or unique: + un = {} + for el in arr: + if addresses: + key = "%s:%d:%d" % (el[2], el[3], el[6]) + un[key] = el[1] + else: + key = "%s:%d:" % (el[2], el[3]) + + if un.get(key): + un[key] += 1 + else: + un[key] = 1 + + keys = un.keys() + + if unique: + def cmpKeys(a,b): # sort keys by file then linenumber + a=a.split(':') + b=b.split(':') + c= cmp(a[0],b[0]) + if c: + return c + else: + return int(a[1]) - int(b[1]) + else: + def cmpKeys(a,b): # sort keys by ID + return un[a] - un[b] + + keys.sort(cmpKeys) + + arr = [] + for el in keys: + arr.append((el, un[el])) + + if fd == None and not getList: + fd = sys.stdout + + if fd != None: + if len(arr) == 0: + fd.write("No leaks\n") + else: + if unique: + fd.write("%-20s Nptr\n" % "file:line") + for el in arr: + fd.write("%-20s %d\n" % el) + else: + fd.write("%-15s %4s %-10s ID\n" % ("file", "line", "ADDR")) + for el in arr: + el = el[0].split(":") + [el[1]] + fd.write("%-15s %4s 0x%08x %d\n" % (el[0], el[1], int(el[2]), el[3])) + + try: + if fd != sys.stderr and fd != sys.stdout: + fd.close(); + except: pass + + if getList: + return arr Index: /tags/rel-0-1-0/pois/src/.cvsignore =================================================================== --- /tags/rel-0-1-0/pois/src/.cvsignore (revision 11633) +++ /tags/rel-0-1-0/pois/src/.cvsignore (revision 11633) @@ -0,0 +1,8 @@ +Makefile Makefile.in +.deps .libs +errorCodes.c poisErrorCodes.h +*.la *.lo +pois +config.h +config.h.in +stamp-h1 Index: /tags/rel-0-1-0/pois/src/Makefile.am =================================================================== --- /tags/rel-0-1-0/pois/src/Makefile.am (revision 11633) +++ /tags/rel-0-1-0/pois/src/Makefile.am (revision 11633) @@ -0,0 +1,58 @@ +AM_CFLAGS = $(PSLIB_CFLAGS) +AM_LDFLAGS = $(PSLIB_LIBS) + +pois_CPPFLAGS = $(PSLIB_CFLAGS) +pois_LDFLAGS = $(PSLIB_LIBS) + +noinst_HEADERS = pois.h poisErrorCodes.h +bin_PROGRAMS = pois + +BUILT_SOURCES = \ + poisErrorCodes.h errorCodes.c + +pois_SOURCES = \ + pois.c \ + poisAddPenalty.c \ + poisCalculateDeviations.c \ + poisCalculateEquation.c \ + poisCheckKernel.c \ + poisConvolveImage.c \ + poisExtractKernel.c \ + poisFindStamps.c \ + poisKernelBasisFunctions.c \ + poisMakeMask.c \ + poisReadStamps.c \ + poisRejectStamps.c \ + poisParseConfig.c \ + poisSolveEquation.c \ + poisStamp.c \ + errorCodes.c +### poisMaskOps.c + +# +# Error codes +# +poisErrorCodes.h : errorCodes.dat errorCodes-skl.h + $(PERL) parseErrorCodes.pl --data=errorCodes.dat --outfile=poisErrorCodes.h \ + errorCodes-skl.h +errorCodes.c : errorCodes.dat errorCodes-skl.c poisErrorCodes.h + $(PERL) parseErrorCodes.pl --data=errorCodes.dat --outfile=errorCodes.c \ + errorCodes-skl.c + +EXTRA_DIST = \ + parseErrorCodes.pl \ + errorCodes.dat \ + errorCodes-skl.h \ + errorCodes-skl.c \ + poisErrorCodes.h + +# +# test program +# +test: pois test1.fits test2.fits + -$(RM) testout.fits stamp*.fits kernel*.fits + ./pois -v test1.fits test2.fits testout.fits + +# Tags for emacs +tags: + etags `find . -name \*.[ch] -print` Index: /tags/rel-0-1-0/pois/src/errorCodes-skl.c =================================================================== --- /tags/rel-0-1-0/pois/src/errorCodes-skl.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/errorCodes-skl.c (revision 11633) @@ -0,0 +1,36 @@ +/* + * The lines between + * //~Start TEMPLATE + * and + * //~End + * will be machine generated from errorCodes.dat + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "pslib.h" +#include "poisErrorCodes.h" + +void poisErrorRegister(void) +{ + static psErrorDescription errors[] = { + { POIS_ERR_BASE, "First value we use; lower values belong to psLib" }, + //~Start {POIS_ERR_$1, "$2"}, + {POIS_ERR_UNKNOWN, "Unknown POIS error code"}, + {POIS_ERR_NOT_IMPLEMENTED, "Desired feature is not yet implemented"}, + {POIS_ERR_FITS, "Problem in FITS I/O"}, + {POIS_ERR_FITS_WCS, "Error interpreting FITS WCS information"}, + //~End + }; + static int nerror = POIS_ERR_NERROR - POIS_ERR_BASE; // number of values in enum + + for (int i = 0; i < nerror; i++) { + psErrorDescription *tmp = psAlloc(sizeof(psErrorDescription)); + *tmp = errors[i]; + psErrorRegister(tmp, 1); + psFree(tmp); /* it's on the internal list */ + } + nerror = 0; // don't register more than once +} Index: /tags/rel-0-1-0/pois/src/errorCodes-skl.h =================================================================== --- /tags/rel-0-1-0/pois/src/errorCodes-skl.h (revision 11633) +++ /tags/rel-0-1-0/pois/src/errorCodes-skl.h (revision 11633) @@ -0,0 +1,23 @@ +#if !defined(POIS_ERROR_CODES_H) +#define POIS_ERROR_CODES_H +/* + * The lines between + * //~Start TEMPLATE + * and + * //~End + * will be machine generated from errorCodes.dat + */ +typedef enum { + POIS_ERR_BASE = 512, + //~Start POIS_ERR_$1, + POIS_ERR_UNKNOWN, + POIS_ERR_NOT_IMPLEMENTED, + POIS_ERR_FITS, + POIS_ERR_FITS_WCS, + //~End + POIS_ERR_NERROR +} poisErrorCode; + +void poisErrorRegister(void); + +#endif Index: /tags/rel-0-1-0/pois/src/errorCodes.dat =================================================================== --- /tags/rel-0-1-0/pois/src/errorCodes.dat (revision 11633) +++ /tags/rel-0-1-0/pois/src/errorCodes.dat (revision 11633) @@ -0,0 +1,11 @@ +# +# This file is used to generate poisErrorClasses.h +# +BASE = 512 First value we use; lower values belong to psLib +UNKNOWN Unknown POIS error code +NOT_IMPLEMENTED Desired feature is not yet implemented +FITS Problem in FITS I/O +FITS_WCS Error interpreting FITS WCS information +SWIG Problem in SWIG bindings +UNSUPPORTED_TYPE The datatype is not supported for this operation + Index: /tags/rel-0-1-0/pois/src/parseErrorCodes.pl =================================================================== --- /tags/rel-0-1-0/pois/src/parseErrorCodes.pl (revision 11633) +++ /tags/rel-0-1-0/pois/src/parseErrorCodes.pl (revision 11633) @@ -0,0 +1,109 @@ +#!/usr/bin/perl + +# Provides functions for handling long command line options +use Getopt::Long; + +my @ErrorCodes = (); +my @ErrorDescriptions = (); + +my $data = "psErrorCodes.dat"; + +# Assign variables based on the presence of command line options to the script +GetOptions( + "data=s" => \$data, + "outfile=s" => \$outputfile, + "verbose" => \$verbose, + "help" => \$help + ) or exit(1); + +if ($help) { + print "Usage: parseErrorCodes ", "[--data=$data] ", "[--outfile=$file] ", "[--help] ", + "[--verbose] [filename [filename ...]]\n\n"; + exit(0); +} + +if(defined($outputfile)) { + if (@ARGV > 1) { + die "You may only specify one input filename if you use --outfile\n"; + } +} + +print "Using data file '$data'\n" if $verbose; +unless ( open( DATAFILE, "<", $data ) ) { + die "Can not open data file $data."; +} + +print "Datafile:\n" if $verbose; +while () { + chop; + if (/^\s*\#/) { + print "C $_\n" if $verbose; + } + else { + if (/^\s*(\w+)\s+([\%\w].*)/) { + my $ErrorCode = $1; + my $ErrorDescription = $2; + print " $ErrorCode: '$ErrorDescription'\n" if $verbose; + push( @ErrorCodes, $ErrorCode ); + push( @ErrorDescriptions, $ErrorDescription ); + } else { + print "I $_\n" if $verbose; + } + } +} +close(DATAFILE); + +my $found = $#ErrorCodes + 1; +print "\nFound $found error codes.\n" if $verbose; + +foreach (@ARGV) { + my $filename = $_; + my @result = (); + + die "Failed to open input file '$filename'" + if !open( INFILE, "<", $filename ); + + print "\nOutput File:\n" if $verbose; + while () { + chop; + push( @result, $_ ); + if (/^\s*\/\/~Start(.*)$/) { + $line = $1; + for ( $n = 0 ; $n < $found ; $n++ ) { + $_ = $line; + s/\$1/$ErrorCodes[$n]/g; + s/\$2/$ErrorDescriptions[$n]/g; + s/\$n/$n/g; + push( @result, $_ ); + print "$_\n" if $verbose; + } + + $break = 0; + while ( ( $break == 0 ) && ( $_ = ) ) { + if (/^\s*\/\/~End/) { + $break = 1; + } + } + chop; + push( @result, $_ ); + } + } + + close(INFILE); + + if (defined($outputfile)) { + $ofile = $outputfile; + } else { + $ofile = $filename; + } + die "Failed to overwrite input file" + if !open( OUTFILE, ">", $ofile ); + + foreach (@result) { + print OUTFILE "$_\n"; + print "$_\n" if $verbose; + } + + close(OUTFILE); + +} Index: /tags/rel-0-1-0/pois/src/pois.c =================================================================== --- /tags/rel-0-1-0/pois/src/pois.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/pois.c (revision 11633) @@ -0,0 +1,354 @@ +// POIS: Pan-STARRS (or Pricey's) Optimal Image Subtraction +// +// Attempts to fit for a general convolution kernel that will match the PSFs between two images. +// +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include "pslib.h" +#include "pois.h" + +#define MAXCHAR 80 + +// Gets the current time. +double getTime(void) +{ + struct timeval tv; + + gettimeofday(&tv, NULL); + return(tv.tv_sec + 1.e-6 * tv.tv_usec); +} + + +int main(int argc, char **argv) +{ + double startTime; + startTime = getTime(); + + // Set trace levels + psTraceSetLevel(".", 0); + + // Parse the command line + poisConfig *config = poisParseConfig(argc, argv); // Configuration values + + if (config->verbose) { + psTraceSetLevel("pois", 10); + psTraceSetLevel("pois.config", 10); + psTraceSetLevel("pois.time", 10); + psTraceSetLevel("pois.solution", 0); + psTraceSetLevel("pois.kernelBasisFunctions", 10); + psTraceSetLevel("pois.calculateEquation", 10); + psTraceSetLevel("pois.convolveImage", 10); + psTraceSetLevel("pois.findStamps", 10); + psTraceSetLevel("pois.parseConfig", 10); + psTraceSetLevel("pois.makeMask", 10); + psTraceSetLevel("pois.extractKernel", 10); + psTraceSetLevel("pois.calculateDeviations", 10); + psTraceSetLevel("pois.checkKernel", 10); + if (config->stampFile) { + psTraceSetLevel("pois.checkStamp", 10); + } else { + psTraceSetLevel("pois.checkStamp", 0); + } + // Set logging level + psLogSetLevel(9); + } + + // Read inputs + psMetadata *header = NULL; // Header for output image + + psRegion imageRegion = {0, 0, 0, 0}; + psFits *reference = psFitsOpen(config->refFile, "r"); + if (config->reverse) { + header = psFitsReadHeader(NULL, reference); + } + psImage *refImage = psFitsReadImage(reference, imageRegion, 0); + if (refImage == NULL) { + psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->refFile); + exit(EXIT_FAILURE); + } + psTrace("pois", 3, "Read %s: %dx%d\n", config->refFile, refImage->numCols, refImage->numRows); + psFitsClose(reference); + // Convert to 32-bit floating point, in necessary + if (refImage->type.type != PS_TYPE_F32) { + psTrace("pois", 3, "Converting %s to floating point in memory....\n", config->refFile); + psImage *temp = psImageCopy(NULL, refImage, PS_TYPE_F32); + psFree(refImage); + refImage = temp; + } + + psFits *input = psFitsOpen(config->inFile, "r"); + if (! config->reverse) { + header = psFitsReadHeader(NULL, input); + } + psImage *inImage = psFitsReadImage(input, imageRegion, 0); + if (inImage == NULL) { + psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->inFile); + exit(EXIT_FAILURE); + } + psTrace("pois", 3, "Read %s: %dx%d\n", config->inFile, inImage->numCols, inImage->numRows); + psFitsClose(input); + // Convert to 32-bit floating point, in necessary + if (inImage->type.type != PS_TYPE_F32) { + psTrace("pois", 3, "Converting %s to floating point in memory....\n", config->inFile); + psImage *temp = psImageCopy(NULL, inImage, PS_TYPE_F32); + psFree(inImage); + inImage = temp; + } + + if ((refImage->numCols != inImage->numCols) || (refImage->numRows != inImage->numRows)) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "Reference and input images must have the same dimensions: %dx%d vs %dx%d\n", + refImage->numCols, refImage->numRows, inImage->numCols, inImage->numRows); + exit(EXIT_FAILURE); + } + config->xImage = refImage->numCols; + config->yImage = refImage->numRows; + + psTrace("pois.time", 1, "Inputs read at %f sec\n", getTime() - startTime); + + // Add a background offset to get the image all above zero + // An image that goes negative (e.g., a background-subtracted image) produces + // horrible results --- the matrix goes NaN because the noise is determined + // as the square root of the value. + if (config->background != 0.0) { + (void)psBinaryOp(refImage, refImage, "+", psScalarAlloc(config->background, PS_TYPE_F32)); + (void)psBinaryOp(inImage, inImage, "+", psScalarAlloc(config->background, PS_TYPE_F32)); + } + + // Generate the kernel parameters + psArray *kernelBasisFunctions = poisKernelBasisFunctions(config); + + // Generate a mask + psImage *mask = poisMakeMask(NULL, refImage, inImage, config); + +#ifdef TESTING + // Write the mask out + psFits *maskFile = psFitsOpen("mask.fits", "w"); + if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write mask: mask.fits\n"); + } + psTrace("pois", 1, "Mask written to mask.fits\n"); + psFitsClose(maskFile); +#endif + + psArray *stamps = NULL; // Array of stamps + psVector *solution = NULL; // Solution vector + + FILE *stampsFP = NULL; // File pointer for stamps + if (config->stampFile) { + psTrace("pois", 5, "Opening stamp file, %s\n", config->stampFile); + stampsFP = fopen(config->stampFile, "r"); + if (! stampsFP) { + psError(PS_ERR_IO, true, "Unable to open stamps file, %s!\n", config->stampFile); + exit(EXIT_FAILURE); + } + } + + // Iterate for a good solution + bool badStamps = true; // Do we have bad stamps, such that we need to continue to iterate? + for (int iterNum = 0; iterNum < config->numIter && badStamps; iterNum++) { + psTrace("pois", 1, "Iteration %d...\n", iterNum); + + // Find stamps + if (config->stampFile) { + stamps = poisReadStamps(stamps, &stampsFP, refImage, mask, config); + } else { + stamps = poisFindStamps(stamps, refImage, mask, config); + } + int numStamps = 0; + for (int s = 0; s < stamps->n; s++) { + poisStamp *stamp = stamps->data[s]; // Stamp of interest + if (stamp->status != POIS_STAMP_BAD) { + numStamps++; + } + } + psTrace("pois.time", 1, "%d stamps found at %f sec\n", numStamps, getTime() - startTime); + if (numStamps == 0) { + fprintf(stderr, "No stamps found. Check input parameters.\n"); + exit(EXIT_FAILURE); + } + + // Calculate equation + (void)poisCalculateEquation(stamps, refImage, inImage, kernelBasisFunctions, config); + psTrace("pois.time", 1, "Equation calculated at %f sec\n", getTime() - startTime); + + // Solve the equation + solution = poisSolveEquation(solution, stamps, config); + psTrace("pois.time", 1, "Matrix equation solved at %f sec\n", getTime() - startTime); + +#ifdef TESTING + // Print solution + psTrace("pois.solution", 5, "Solution:\n"); + for (int i = 0; i < solution->n - 1; i++) { + poisKernelBasis *kernel = kernelBasisFunctions->data[i]; // The i-th kernel basis function + psTrace("pois.solution", 7, "%d: (%d,%d) x^%d y^%d --> %f\n", i, kernel->u, kernel->v, + kernel->xOrder, kernel->yOrder, solution->data.F64[i]); + } +#endif + +#ifdef TESTING + // Save the kernel postage stamp + char kernelName[MAXCHAR]; // File name for kernel + snprintf(kernelName, MAXCHAR, "kernel%d.fits", iterNum); + psImage *kernelImage = poisExtractKernel(solution, kernelBasisFunctions, 0.0, 0.0, config); + psFits *kernelFile = psFitsOpen(kernelName, "w"); + if (!psFitsWriteImage(kernelFile, NULL, kernelImage, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write kernel: %s\n", kernelName); + } + psTrace("pois", 1, "Kernel written to %s\n", kernelName); + psFitsClose(kernelFile); + psFree(kernelImage); +#endif + + // Calculate deviations for the stamps + psVector *deviations = poisCalculateDeviations(NULL, stamps, refImage, inImage, mask, + kernelBasisFunctions, solution, config); + badStamps = poisRejectStamps(stamps, mask, deviations, config); + + psFree(deviations); + } // Iterating for good solution + + // If there was rejection on the last round, re-solve the equation using only the good stamps + if (badStamps) { + solution = poisSolveEquation(solution, stamps, config); + } + + if (stampsFP) { + fclose(stampsFP); + } + + // Undo the mask for bad stamp area + for (int j = 0; j < mask->numRows; j++) { + for (int i = 0; i < mask->numCols; i++) { + if (mask->data.U8[j][i] & POIS_MASK_STAMP) { + mask->data.U8[j][i] &= ~POIS_MASK_STAMP; + } + } + } + + // Check the result + (void)poisCheckKernel(solution, kernelBasisFunctions, config); + + // Convolve the input image + psImage *convImage = poisConvolveImage(refImage, mask, solution, kernelBasisFunctions, config); + psTrace("pois.time", 1, "Convolution completed at %f sec\n", getTime() - startTime); + +#ifdef TESTING + // Save the convolved image + char convName[MAXCHAR]; // Filename for convolved image + snprintf(convName, MAXCHAR, "%s.conv", config->outFile); + psFits *convFile = psFitsOpen(convName, "w"); + if (!psFitsWriteImage(convFile, NULL, convImage, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write convolved image: %s\n", convName); + } + psFitsClose(convFile); + psTrace("pois", 1, "Convolved image written to %s\n", convName); +#endif + + // Do the subtraction + psImage *subImage = psImageAlloc(config->xImage, config->yImage, PS_TYPE_F32); + if (config->reverse) { + (void)psBinaryOp(subImage, convImage, "-", inImage); + } else { + (void)psBinaryOp(subImage, inImage, "-", convImage); + } + psTrace("pois.time", 1, "Subtraction completed at %f sec\n", getTime() - startTime); + +#ifdef TESTING + { + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + (void)psImageStats(stats, subImage, mask, 0); + psTrace("pois", 5, "Subtracted image statistics: %f %f\n", stats->sampleMean, stats->sampleStdev); + psFree(stats); + } +#endif + + // Mask out the edges + for (int y = config->yKernel; y < subImage->numRows - config->yKernel; y++) { + for (int x = config->xKernel; x < subImage->numCols - config->xKernel; x++) { + + for (int x = 0; x < config->xKernel; x++) { + subImage->data.F32[y][x] = 0.0; + } + for (int x = subImage->numCols - config->xKernel; x < subImage->numCols; x++) { + subImage->data.F32[y][x] = 0.0; + } + } + } + for (int y = 0; y < config->yKernel; y++) { + for (int x = 0; x < subImage->numCols; x++) { + subImage->data.F32[y][x] = 0.0; + } + } + for (int y = subImage->numRows - config->yKernel; y < subImage->numRows; y++) { + for (int x = 0; x < subImage->numCols; x++) { + subImage->data.F32[y][x] = 0.0; + } + } + + + // Mask out bad pixels + for (int y = config->yKernel; y < subImage->numRows - config->yKernel; y++) { + for (int x = config->xKernel; x < subImage->numCols - config->xKernel; x++) { + if (mask->data.U8[y][x]) { + subImage->data.F32[y][x] = 0.0; + } + } + } + + if (config->mask) { + char maskName[80]; // Name of mask image + sprintf(maskName, "%s.mask", config->outFile); + psFits *maskFile = psFitsOpen(maskName, "w"); + if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName); + } + psTrace("pois", 1, "Mask image written to %s\n", maskName); + psFitsClose(maskFile); + } + + + int numNaN = psImageClipNaN(subImage, 0.0); + if (numNaN > 0) { + printf("Masked %d NAN pixels to zero.\n", numNaN); + } + + psFits *subFile = psFitsOpen(config->outFile, "w"); + if (!psFitsWriteImage(subFile, header, subImage, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write subtracted image: %s\n", config->outFile); + } + psFitsClose(subFile); + psTrace("pois", 1, "Subtracted image written to %s\n", config->outFile); + psTrace("pois.time", 1, "I/O completed at %f sec\n", getTime() - startTime); + + // Clean up + psFree(header); + psFree(kernelBasisFunctions); + psFree(refImage); + psFree(inImage); + psFree(convImage); + psFree(subImage); + psFree(config); +} + Index: /tags/rel-0-1-0/pois/src/pois.h =================================================================== --- /tags/rel-0-1-0/pois/src/pois.h (revision 11633) +++ /tags/rel-0-1-0/pois/src/pois.h (revision 11633) @@ -0,0 +1,234 @@ +#if !defined(POIS_H) +#define POIS_H + +#include "pslib.h" + +#define abs(x) ((x) >= 0 ? (x) : (-(x))) + +/************************************************************************************************************/ + +/* Configuration options */ +typedef struct { + int verbose; // Verbosity level + char *refFile; // Reference image + char *inFile; // Input image filename + char *outFile; // Output convolved image filename + int xKernel, yKernel; // Half the size of the kernel in x and y + int xImage, yImage; // Size of the images in x and y + float threshold; // Threshold for pixel to include in calculation + int footprint; // Footprint size for stamp + int nsx, nsy; // Number of stamps in x and y + float refSat, inSat; // Saturation values for reference and input images + float refBad, inBad; // Bad values for reference and input images + float background; // Value to add to the background in images to get above zero + int spatialOrder; // Order of spatial variations in the kernel + int numIter; // Number of rejection iterations + float sigmaRej; // Limit (in std dev) for rejection + float penalty; // Penalty value + char *stampFile; // Filename for file with stamps + bool reverse; // Reverse the subtraction? + bool mask; // Output the mask? +} poisConfig; + +/************************************************************************************************************/ + +// A kernel basis function +typedef struct { + int u, v; // Offset + int xOrder, yOrder; // x^xOrder * y^yOrder +} poisKernelBasis; + +/************************************************************************************************************/ + +typedef enum { + POIS_STAMP_USED, // Use this stamp + POIS_STAMP_RESET, // This stamp has been rejected --- reset + POIS_STAMP_RECALC, // Having been reset, this stamp needs to be recalculated + POIS_STAMP_BAD // No stamp in this region +} poisStampStatus; + +// A stamp +typedef struct { + int x, y; // Position + psImage *matrix; // Associated matrix + psVector *vector; // Associated vector + poisStampStatus status; // Status of stamp +} poisStamp; + +/************************************************************************************************************/ + +// poisStamp.c + +// Construct a stamp +poisStamp *poisStampAlloc(void); + +// Destroy a stamp +void poisStampFree(poisStamp *stamp); + +// Return true if the stamp is OK +bool poisCheckStamp(const psImage *image, // Image to check for threshold + const psImage *mask, // Mask to check for bad pixels + int x, int y, // Pixel coordinates + const poisConfig *config // Configuration + ); + +/************************************************************************************************************/ + +typedef enum { + POIS_MASK_OK = 0x00, // Fine + POIS_MASK_BAD = 0x01, // Bad pixels + POIS_MASK_NEAR_BAD = 0x02, // Near bad pixels + POIS_MASK_STAMP = 0x04 // Masking a stamp +} poisMask; + +/************************************************************************************************************/ + +// poisParseConfig.c + +// Print usage information +void help(void); + +// Parse the command line arguments +poisConfig *poisParseConfig(int argc, // Number of command-line arguments + char **argv // Command-line arguments + ); + +/************************************************************************************************************/ + +// poisKernelBasisFunctions.c + +// Generate the kernel basis functions for ease of reference in the construction +psArray *poisKernelBasisFunctions(const poisConfig *config // Configuration + ); + +/************************************************************************************************************/ + +// poisMakeMask.c + +// Generate a mask of bad pixels in the reference and input images --- don't use these pixels for stamps! +psImage *poisMakeMask(psImage *mask, // Mask to update, or NULL + const psImage *refImage, // Reference image + const psImage *inImage, // Input image + const poisConfig *config // Configuration + ); + +/************************************************************************************************************/ + +// poisFindStamps.c + +// Find stamps +psArray *poisFindStamps(psArray *stamps, // Existing list of stamps, or NULL + const psImage *image, // Image for which to find stamps + const psImage *mask, // Mask image + const poisConfig *config // Configuration values + ); + +/************************************************************************************************************/ + +// poisExtractKernel.c + +// Create a kernel image from the solution vector +psImage *poisExtractKernel(const psVector *solution, // Vector containing the kernel elements + const psArray *kernelParams, // The kernel parameters + float x, // The relative (-1 to 1) x position for which to get the kernel + float y, // The relative (-1 to 1) y position for which to get the kernel + const poisConfig *config // Configuration + ); + +/************************************************************************************************************/ + +// poisCalculateEquation.c + +// Calculate the matrix and vector for the matrix equation +int poisCalculateEquation(psArray *stamps, // The stamps + const psImage *refImage, // The reference image + const psImage *inImage, // The input image + const psArray *kernelParams, // Kernel parameters + const poisConfig *config // Configuration + ); + +/************************************************************************************************************/ + +// poisConvolveImage.c + +// Convolve an image by the kernel +psImage *poisConvolveImage(const psImage *in, // Input image, to be convolved + const psImage *mask, // Mask image + const psVector *solution, // The solution vector, containing the coefficients + const psArray *kernelParams, // The kernel parameters + const poisConfig *config // The configuration + ); + +/************************************************************************************************************/ + +// poisCalculateDeviations.c + +// Calculate deviations from the best fit for the stamps +psVector *poisCalculateDeviations(psVector *deviations, // Output array of deviations, or NULL + psArray *stamps, // Array of stamps + psImage *refImage, // Reference image + psImage *inImage, // Input image + psImage *mask, // Mask image + psArray *kernelParams, // Array of kernel parameters + psVector *solution, // Solution vector + poisConfig *config // Configuration + ); + +/************************************************************************************************************/ + +// poisRejectStamps.c + +// Reject stamps, based on the deviations +bool poisRejectStamps(psArray *stamps, // Array of stamps + psImage *mask, // Mask image + const psVector *deviations, // Vector of deviations for the stamps + const poisConfig *config // Configuration + ); + +/************************************************************************************************************/ + +// poisSolveEquation.c + +// Solve the matrix equation +psVector *poisSolveEquation(psVector *solution, // Solution vector, or NULL + const psArray *stamps, // Array of stamps + const poisConfig *config + ); + +/************************************************************************************************************/ + +// poisAddPenalty.c + +// Generate penalty matrix and vector +void poisAddPenalty(psImage *matrix, // The matrix + psVector *vector, // The vector + float penalty, // Penalty value + const psArray *kernels, // Kernel basis functions + const poisConfig *config // Configuration + ); + +/************************************************************************************************************/ + +// poisReadStamps.c + +// Read stamps from a file +psArray *poisReadStamps(psArray *stamps,// Stamps + FILE **stampFP, // File pointer for stamps file + const psImage *image, // Image for which to find stamps + const psImage *mask, // Mask image + const poisConfig *config // Configuration + ); + +/************************************************************************************************************/ + +// poisCheckKernel.c + +// Check kernel for power at large radii +bool poisCheckKernel(const psVector *solution,// Kernel solution + const psArray *kernels, // Kernel basis functions + const poisConfig *config // Configuration + ); + +/************************************************************************************************************/ + +#endif Index: /tags/rel-0-1-0/pois/src/poisAddPenalty.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisAddPenalty.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisAddPenalty.c (revision 11633) @@ -0,0 +1,77 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" +#include "pois.h" + +void poisAddPenalty(psImage *matrix, // The matrix + psVector *vector, // The vector + float penalty, // Penalty value + const psArray *kernels, // Kernel basis functions + const poisConfig *config // Configuration + ) +{ + // Check inputs + assert(matrix); + assert(vector); + assert(kernels); + assert(config); + assert(matrix->numRows == matrix->numCols); + assert(vector->n == matrix->numRows); + assert(vector->n == kernels->n + 1);// Additional space for constant background + assert(matrix->type.type == PS_TYPE_F64); + assert(vector->type.type == PS_TYPE_F64); + + int numKernels = kernels->n; + + for (int j = 0; j < numKernels; j++) { + poisKernelBasis *kernel = kernels->data[j]; + if (kernel->xOrder == 0 && kernel->yOrder == 0) { + for (int k = j + 1; k < numKernels; k++) { + poisKernelBasis *compare = kernels->data[k]; + if (compare->xOrder == 0 && compare->yOrder == 0 && + (compare->u == kernel->u + 1 || compare->u == kernel->u - 1 || + compare->v == kernel->v + 1 || compare->v == kernel->v - 1)) { + matrix->data.F64[j][k] -= penalty; + matrix->data.F64[k][j] -= penalty; + } + } + } + + if (abs(kernel->u) == config->xKernel) { + if (abs(kernel->v) == config->yKernel) { + vector->data.F64[j] -= 3.0 * penalty; + } else { + vector->data.F64[j] -= 5.0 * penalty; + } + } else { + if (abs(kernel->v) == config->yKernel) { + vector->data.F64[j] -= 5.0 * penalty; + } else { + vector->data.F64[j] -= 8.0 * penalty; + } + } + + } + + +} Index: /tags/rel-0-1-0/pois/src/poisCalculateDeviations.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisCalculateDeviations.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisCalculateDeviations.c (revision 11633) @@ -0,0 +1,125 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" +#include "pois.h" + +#define MAXCHAR 80 + +psVector *poisCalculateDeviations(psVector *deviations, // Output array of deviations, or NULL + psArray *stamps, // Array of stamps + psImage *refImage, // Reference image + psImage *inImage, // Input image + psImage *mask, // Mask image + psArray *kernelParams, // Array of kernel parameters + psVector *solution, // Solution vector + poisConfig *config // Configuration + ) +{ + // Check inputs + assert(!deviations || deviations->type.type == PS_TYPE_F32); + assert(!deviations || deviations->n == stamps->n); + assert(stamps); + assert(refImage->type.type == PS_TYPE_F32); + assert(inImage->type.type == PS_TYPE_F32); + assert(refImage->numRows == inImage->numRows); + assert(refImage->numCols == inImage->numCols); + assert(mask->type.type == PS_TYPE_U8); + assert(refImage->numRows == mask->numRows); + assert(refImage->numCols == mask->numCols); + assert(solution->n == kernelParams->n + 1); + assert(solution->type.type == PS_TYPE_F64); + assert(config); + + if (!deviations) { + deviations = psVectorAlloc(stamps->n, PS_TYPE_F32); + } + + int footprint = config->footprint; // Size of stamp footprint + int xSize = footprint + config->xKernel; // (Half) size of subimage in x + int ySize = footprint + config->yKernel; // (Half) size of subimage in y + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics + + // Pre-allocate for optimisation: + psImage *subStamp = psImageAlloc(2 * xSize, 2 * ySize, PS_TYPE_F32); // Subtraction of stamp + + for (int s = 0; s < stamps->n; s++) { + poisStamp *stamp = stamps->data[s]; // The coordinates of the stamp of interest + int x = stamp->x; // Stamp x coord + int y = stamp->y; // Stamp y coord + if (stamp->status == POIS_STAMP_USED) { + psRegion stampRegion = {x - xSize, x + xSize, y - ySize, y + ySize}; + psImage *refStamp = psImageSubset(refImage, stampRegion); + psImage *inStamp = psImageSubset(inImage, stampRegion); + psImage *maskStamp = psImageSubset(mask, stampRegion); + psImage *convRefStamp = poisConvolveImage(refStamp, maskStamp, solution, kernelParams, config); + // Calculate chi^2 + (void)psBinaryOp(subStamp, inStamp, "-", convRefStamp); + (void)psBinaryOp(subStamp, subStamp, "*", subStamp); + (void)psBinaryOp(subStamp, subStamp, "/", inStamp); + + // Trim regions. Need separate regions for the subtraction and the mask, since the subtraction + // isn't a child image (it's the subtraction of two child images, but placed in a pre-allocated + // image which isn't a child image), but the mask *is* a child image. Therefore the regions are + // different, since they refer to the region on the parent image (if it exists). + psRegion stampTrim = { config->xKernel, config->xKernel + 2 * footprint, config->yKernel, + config->yKernel + 2 * footprint }; +#if 1 + psRegion maskTrim = { x - xSize + config->xKernel, x + xSize - config->xKernel, + y - ySize + config->yKernel, y + ySize - config->yKernel }; + psImage *maskStampTrim = psImageSubset(maskStamp, maskTrim); +#else + psImage *maskStampTrim = psImageSubset(maskStamp, stampTrim); +#endif + psImage *subStampTrim = psImageSubset(subStamp, stampTrim); + (void)psImageStats(stats, subStampTrim, maskStampTrim, POIS_MASK_BAD | POIS_MASK_NEAR_BAD); + + deviations->data.F32[s] = sqrtf(stats->sampleMean / 2.0); + psTrace("pois.calculateDeviations", 5, "Deviation of stamp %d (%d,%d) is %f\n", s, x, y, + deviations->data.F32[s]); + +#ifdef TESTING + char stampName[MAXCHAR]; // File name for stamp + snprintf(stampName, MAXCHAR, "stamp%d.fits", s); + psFits *stampFile = psFitsOpen(stampName, "w"); + if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write stamp: %s\n", stampName); + } + psFitsClose(stampFile); +#endif + +#if 0 + psFree(convRefStamp); + psFree(maskStampTrim); + psFree(subStampTrim); + psFree(maskStamp); + psFree(refStamp); + psFree(inStamp); +#endif + } + } + + psFree(stats); + psFree(subStamp); + + return deviations; +} Index: /tags/rel-0-1-0/pois/src/poisCalculateEquation.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisCalculateEquation.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisCalculateEquation.c (revision 11633) @@ -0,0 +1,198 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include "pslib.h" +#include "pois.h" + +/* Calculate the matrix and vector for the matrix equation */ +int poisCalculateEquation(psArray *stamps, // The stamps + const psImage *refImage, // The reference image + const psImage *inImage, // The input image + const psArray *kernelParams, // Kernel parameters + const poisConfig *config // Configuration + ) +{ + // Check inputs + assert(stamps); + assert(refImage); + assert(inImage); + assert(kernelParams); + assert(config); + + // In addition to the kernel parameters, need one extra parameter to fit the (assumed constant) difference + // in the background levels. + int numSolveParams = kernelParams->n + 1; + int bgIndex = kernelParams->n; // Index in matrix for the background (it's the last one, after the + // kernel parameters) + float **reference = refImage->data.F32; // The reference pixels + float **input = inImage->data.F32; // The input pixels + int nPix = 0; // Number of pixels used in calculation (for interest) + float nx = 0.5 * (float)refImage->numCols; // Half-size of input image in x + float ny = 0.5 * (float)refImage->numRows; // Half-size of input image in y + int numKernelParams = kernelParams->n; // Number of kernel parameters + + psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, config->spatialOrder + 1, + config->spatialOrder + 1); // Polynomial for evaluation + for (int j = 0; j < config->spatialOrder + 1; j++) { + for (int i = 0; i < config->spatialOrder + 1; i++) { + poly->coeff[j][i] = 1.0; + poly->mask[j][i] = 1; // Mask all coefficients; unmask to evaluate + } + } + + psImage *polyValues = psImageAlloc(config->spatialOrder + 1, config->spatialOrder + 1, + PS_TYPE_F64); // Evaluations for each of the polynomial orders + + /* Iterate over the stamps */ + for (int s = 0; s < stamps->n; s++) { + poisStamp *stamp = stamps->data[s]; // The stamp + + if (stamp->status == POIS_STAMP_RECALC) { + psTrace("pois.calculateEquation", 5, "Calculating for stamp %d: %d,%d\n", s, stamp->x, stamp->y); + + // Initialise the matrix and vector if required + if (! stamp->matrix) { + psTrace("pois.calculateEquation", 9, "Allocating matrix: %dx%d\n", numSolveParams, + numSolveParams); + stamp->matrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64); + } + if (! stamp->vector) { + psTrace("pois.calculateEquation", 3, "Allocating vector: %d\n", numSolveParams); + stamp->vector = psVectorAlloc(numSolveParams, PS_TYPE_F64); + stamp->vector->n = numSolveParams; + } + + // Dereference, for convenience + psImage *matrix = stamp->matrix; + psVector *vector = stamp->vector; + + // Will assert on the type, since can't do much about it if it's wrong. + assert(matrix->type.type == PS_TYPE_F64); + assert(vector->type.type == PS_TYPE_F64); + + // Initialise the matrix and vector + for (int j = 0; j < numSolveParams; j++) { + for (int i = 0; i < numSolveParams; i++) { + matrix->data.F64[j][i] = 0.0; + } + vector->data.F64[j] = 0.0; + } + + // Evaluate the polynomial for each order + for (int j = 0; j < config->spatialOrder + 1; j++) { + for (int i = 0; i < config->spatialOrder + 1 - j; i++) { + poly->mask[j][i] = 0; + polyValues->data.F64[j][i] = psPolynomial2DEval(poly, ((double)stamp->x - nx)/nx, + ((double)stamp->y - ny)/ny); + poly->mask[j][i] = 1; + } + } + + // Iterate over the pixels + for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) { + for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) { + float invNoise2 = 1.0/reference[y][x]; // The inverse of the noise, squared. + nPix++; + + /* Iterate over the first convolution */ + for (int k1 = 0; k1 < numKernelParams; k1++) { + poisKernelBasis *kernel1 = kernelParams->data[k1]; // Kernel basis #1 + + int u1 = kernel1->u; // Offset in x + int v1 = kernel1->v; // Offset in y + int i1 = kernel1->xOrder; // Polynomial order in x + int j1 = kernel1->yOrder; // Polynomial order in y + // First convolution + float conv1 = polyValues->data.F64[j1][i1] * reference[y - v1][x - u1]; + + // Assuming that the first kernel component is 0 order in x and y, and 0 offset + if (k1 != 0) { + conv1 -= reference[y][x]; + } + + /* Iterate over the second convolution */ + for (int k2 = k1; k2 < numKernelParams; k2++) { + poisKernelBasis *kernel2 = kernelParams->data[k2]; // Kernel basis #2 + + int u2 = kernel2->u; // Offset in x + int v2 = kernel2->v; // Offset in y + int i2 = kernel2->xOrder; // Polynomial order in x + int j2 = kernel2->yOrder; // Polynomial order in y + // Second convolution + float conv2 = polyValues->data.F64[j2][i2] * reference[y - v2][x - u2]; + + // Assuming that the first kernel component is 0 order in x and y, and 0 offset + if (k2 != 0) { + conv2 -= reference[y][x]; + } + + // Add into the matrix element + matrix->data.F64[k1][k2] += conv1 * // First convolution + conv2 * // Second convolution + invNoise2; // Divide by sigma^2, ~ poisson. + } // Second convolution + + // Add into the vector element + vector->data.F64[k1] += input[y][x] * // Input image, no convolution + conv1 * // Convolved reference image + invNoise2; // Divide by sigma^2, ~ poisson. + + /* Background term */ + matrix->data.F64[k1][bgIndex] += conv1 * // Convolved reference image + invNoise2; // Divide by sigma^2, ~ poisson. + + } // First convolution + + /* Background only terms */ + matrix->data.F64[bgIndex][bgIndex] += invNoise2; + vector->data.F64[bgIndex] += input[y][x] * invNoise2; + } + } // Iterating over pixels in stamp + + // Fill in other side of symmetric matrix + for (int k1 = 0; k1 < kernelParams->n; k1++) { + for (int k2 = 0; k2 < k1; k2++) { + matrix->data.F64[k1][k2] = matrix->data.F64[k2][k1]; + } + matrix->data.F64[bgIndex][k1] = matrix->data.F64[k1][bgIndex]; + } + + for (int k = 0; k < kernelParams->n; k++) { + poisKernelBasis *kernel = kernelParams->data[k]; // Kernel basis #2 + int u = kernel->u; // Offset in x + int v = kernel->v; // Offset in y + matrix->data.F64[k][k] += config->penalty * (float)(u*u + v*v); + } + + stamp->status = POIS_STAMP_USED; // We've calculated it now, so don't need to do so again. + + } // Legitimate stamps + + } // Iterating over stamps + + psFree(polyValues); + psFree(poly); + + // Return number of pixels over which matrix has been calculated + return nPix; +} Index: /tags/rel-0-1-0/pois/src/poisCheckKernel.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisCheckKernel.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisCheckKernel.c (revision 11633) @@ -0,0 +1,83 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include "pslib.h" +#include "pois.h" + +bool poisCheckKernel(const psVector *solution,// Kernel solution + const psArray *kernels, // Kernel basis functions + const poisConfig *config // Configuration + ) +{ + psTrace("pois.checkKernel", 1, "Checking kernel....\n"); + + // Get kernel value as a function of radius + psVector *radKernel = psVectorAlloc((2*config->xKernel + 1) * (2*config->yKernel + 1), PS_TYPE_F32); + psVector *radii = psVectorAlloc((2*config->xKernel + 1) * (2*config->yKernel + 1), PS_TYPE_F32); // Radius + for (int i = 0; i < solution->n - 1; i++) { + poisKernelBasis *kernel = kernels->data[i]; // The kernel basis function + if (kernel->xOrder == 0 && kernel->yOrder == 0) { + int u = kernel->u; // x offset + int v = kernel->v; // y offset + float distance = sqrtf((float)(u*u) + (float)(v*v)); // Distance from the centre + radKernel->data.F32[i] = solution->data.F64[i]; + radii->data.F32[i] = distance; + } + } + + psVector *sortIndex = psVectorSortIndex(NULL, radii); // Indices from sort + float distance = 0.0; // Distance from the centre that's being examined + float sum = 0.0; // Sum of kernel for that distance + int num = 0; // Number of pixels at that distance + int centreIndex = sortIndex->data.U32[0]; // Index of the centre pixel + int numSuspect = 0; // Number of radii considered suspect + for (int i = 0; i < (2*config->xKernel + 1) * (2*config->yKernel + 1); i++) { + unsigned int index = sortIndex->data.U32[i]; + if (radii->data.F32[index] != distance) { + psTrace("pois.checkKernel", 7, "Radius: %f\tMean: %f\n", distance, sum/(float)num); + if (fabs(distance*sum/(float)num) > radKernel->data.F32[centreIndex]) { + numSuspect++; + } + sum = 0.0; + num = 0; + distance = radii->data.F32[index]; + } + sum += radKernel->data.F32[index]; + num++; + } + psTrace("pois.checkKernel", 7, "Radius: %f\tMean: %f\n", distance, sum/(float)num); + if (fabs(distance*sum/(float)num) > radKernel->data.F32[centreIndex]) { + numSuspect++; + } + + psFree(radKernel); + psFree(radii); + psFree(sortIndex); + + psTrace("pois.checkKernel", 5, "%d radii considered suspect.\n", numSuspect); + + if (numSuspect > 0) { + psWarning("There is significant power at large radii. The kernel may be bad.\n"); + return false; + } + + return true; +} Index: /tags/rel-0-1-0/pois/src/poisConvolveImage.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisConvolveImage.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisConvolveImage.c (revision 11633) @@ -0,0 +1,122 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" +#include "pois.h" + + +/* Convolve an image by the kernel */ +psImage *poisConvolveImage(const psImage *input, // Input image, to be convolved + const psImage *mask, // Mask image + const psVector *solution, // The solution vector, containing the coefficients + const psArray *kernelParams, // The kernel parameters + const poisConfig *config // The configuration + ) +{ + assert(solution->n == kernelParams->n + 1); // Extra parameter for background + assert(solution->type.type == PS_TYPE_F64); + assert(input->numCols == mask->numCols && input->numRows == mask->numRows); + assert(mask->type.type == PS_TYPE_U8); + assert(input->type.type == PS_TYPE_F32); + + int nx = input->numCols; // The size of the image in x + int ny = input->numRows; // The size of the image in y + float nxHalf = 0.5 * (float)nx; // Half the size of the image + float nyHalf = 0.5 * (float)ny; // Half the size of the image + int xKernel = config->xKernel; // The half-size of the kernel in x + int yKernel = config->yKernel; // The half-size of the kernel in y + float background = solution->data.F64[solution->n-1]; // The difference in background + int numParams = kernelParams->n; // Number of kernel parameters + + psImage *convolved = psImageAlloc(nx, ny, PS_TYPE_F32); // The convolved image, to be returned + + psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, config->spatialOrder + 1, + config->spatialOrder + 1); // Polynomial + for (int i = 0; i < config->spatialOrder + 1; i++) { + for (int j = 0; j < config->spatialOrder + 1; j++) { + poly->coeff[i][j] = 1.0; + poly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate + } + } + + // Convolve only the middle part + for (int y = yKernel; y < ny - yKernel; y++) { + for (int x = xKernel; x < nx - xKernel; x++) { + if (! mask->data.U8[y][x] & (POIS_MASK_BAD | POIS_MASK_NEAR_BAD)) { + double conv = background; // Convolved value + + float imageX = (x - nxHalf) / nxHalf; // Normalised position in x + float imageY = (y - nyHalf) / nyHalf; // Normalised position in y + + // Iterate over the kernel basis functions + for (int k = 0; k < numParams; k++) { + poisKernelBasis *kernel = kernelParams->data[k]; // The kernel basis function + int u = kernel->u; + int v = kernel->v; + int xOrder = kernel->xOrder; + int yOrder = kernel->yOrder; + + double polyVal = 1.0; + if (xOrder != 0 || yOrder != 0) { + // Evaluate the polynomial + poly->mask[xOrder][yOrder] = 0; + polyVal = psPolynomial2DEval(poly, imageX, imageY); + poly->mask[xOrder][yOrder] = 1; + } + + if (k == 0) { + conv += solution->data.F64[k] * input->data.F32[y - v][x - u] * polyVal; + } else { + conv += solution->data.F64[k] * + (input->data.F32[y - v][x - u] * polyVal - input->data.F32[y][x]); + } + } + convolved->data.F32[y][x] = (float)conv; + } + } + + /* Pad the rest with zeros */ + for (int x = 0; x < xKernel; x++) { + convolved->data.F32[y][x] = 0.0; + } + for (int x = nx - xKernel; x < nx; x++) { + convolved->data.F32[y][x] = 0.0; + } + } + + for (int y = 0; y < yKernel; y++) { + for (int x = 0; x < nx; x++) { + convolved->data.F32[y][x] = 0.0; + } + } + for (int y = ny - yKernel; y < ny; y++) { + for (int x = 0; x < nx; x++) { + convolved->data.F32[y][x] = 0.0; + } + } + + psFree(poly); + + return convolved; +} + + Index: /tags/rel-0-1-0/pois/src/poisExtractKernel.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisExtractKernel.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisExtractKernel.c (revision 11633) @@ -0,0 +1,90 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" +#include "pois.h" + +/* Create a kernel image from the solution vector */ +psImage *poisExtractKernel(const psVector *solution, // Vector containing the kernel elements + const psArray *kernelParams, // The kernel parameters + float x, // The relative (-1 to 1) x position for which to get the kernel + float y, // The relative (-1 to 1) y position for which to get the kernel + const poisConfig *config // Configuration + ) +{ + assert(solution->n == kernelParams->n + 1); + assert(solution->type.type == PS_TYPE_F64); + assert(config); + assert(x >= -1.0 && x <= 1.0); + assert(y >= -1.0 && y <= 1.0); + + int xKernel = config->xKernel; // Half-size of kernel in x + int yKernel = config->yKernel; // Half-size of kernel in y + float kernelSum = 0.0; // Sum of the kernel + + /* Allocate the image */ + psImage *image = psImageAlloc(2*xKernel + 1, 2*yKernel + 1, PS_TYPE_F32); // The kernel image + for (int j = 0; j < (2*yKernel + 1); j++) { + for (int i = 0; i < (2*xKernel + 1); i++) { + image->data.F32[j][i] = 0.0; + } + } + + psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, config->spatialOrder + 1, + config->spatialOrder + 1); // Polynomial for evaluation + for (int j = 0; j < config->spatialOrder + 1; j++) { + for (int i = 0; i < config->spatialOrder + 1; i++) { + poly->coeff[j][i] = 1.0; + poly->mask[j][i] = 1; // Mask all coefficients; unmask to evaluate + } + } + + psTrace("pois.extractKernel", 4, "Evaluating kernel at %f,%f\n", x, y); + /* Iterate over the kernel parameters */ + for (int k = 0; k < kernelParams->n; k++) { + poisKernelBasis *kernel = kernelParams->data[k]; // The kernel basis function + + int u = kernel->u; + int v = kernel->v; + int xOrder = kernel->xOrder; + int yOrder = kernel->yOrder; + + // Evaluate polynomial + poly->mask[xOrder][yOrder] = 0; + double evaluation = solution->data.F64[k] * psPolynomial2DEval(poly, x, y); + poly->mask[xOrder][yOrder] = 1; + + image->data.F32[v + yKernel][u + xKernel] += (float)evaluation; + kernelSum += (float)evaluation; + } + + // Subtract from the centre pixel to remove the effect of the flux-conservation treatment + // This assumes that the "reference" is the centre pixel: u = 0, v = 0, xOrder = 0, yOrder = 0 + image->data.F32[yKernel][xKernel] -= kernelSum - (float)solution->data.F64[0]; + + psTrace("pois.extractKernel", 4, "Kernel sum: %f\n", kernelSum); + + psFree(poly); + + return image; +} + Index: /tags/rel-0-1-0/pois/src/poisFindStamps.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisFindStamps.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisFindStamps.c (revision 11633) @@ -0,0 +1,101 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include "pslib.h" +#include "pois.h" + + +psArray *poisFindStamps(psArray *stamps, // Existing list of stamps, or NULL + const psImage *image, // Image for which to find stamps + const psImage *mask, // Mask image + const poisConfig *config // Configuration values + ) +{ + int nx = config->xImage; // Size of image in x + int ny = config->yImage; // Size of image in y + + assert(image->numCols == nx && image->numRows == ny); + assert(mask->numCols == nx && mask->numRows == ny); + assert(image->type.type == PS_TYPE_F32); + assert(mask->type.type == PS_TYPE_U8); + assert(!stamps || stamps->n == config->nsx * config->nsy); + + int footprint = config->footprint; // Footprint for stamps + int borderx = footprint + config->xKernel; // Border around image in x + int bordery = footprint + config->yKernel; // Border around image in y + int nsx = config->nsx; // Number of stamps in x + int nsy = config->nsy; // Number of stamps in y + + psF32 **pixels = image->data.F32; // The image pixels + + if (stamps == NULL) { + stamps = psArrayAlloc(nsx * nsy); + // Initialise + for (int i = 0; i < nsx * nsy; i++) { + stamps->data[i] = poisStampAlloc(); + } + } + + // Iterate over the image sections + int num = 0; + for (int i = 0; i < nsx; i++) { + for (int j = 0; j < nsy; j++) { + poisStamp *stamp = stamps->data[num]; // The stamp + + // Only find a new stamp if we need to + if ((stamp->x == 0 && stamp->y == 0) || stamp->status == POIS_STAMP_RESET) { + // Find maximum non-masked value in the image section, but don't include a footprint around + // the edge + float max = - INFINITY; // Negative infinity + int bestx = 0, besty = 0; // Position of maximum + + for (int y = bordery + j * (ny - 2.0 * bordery) / nsy; + y < bordery + (j + 1) * (ny - 2.0 * bordery) / nsy; y++) { + for (int x = borderx + i * (nx - 2.0 * borderx) / nsx; + x < borderx + (i + 1) * (nx - 2.0 * borderx) / nsx; x++) { + if (pixels[y][x] > max && poisCheckStamp(image, mask, x, y, config)) { + max = pixels[y][x]; + bestx = x; + besty = y; + } + } + } // Done iterating over this image section + + // Store the stamp + if (max != - INFINITY) { + stamp->x = bestx; + stamp->y = besty; + stamp->status = POIS_STAMP_RECALC; + } else { + // No stamp in this section + stamp->status = POIS_STAMP_BAD; + } + } // Finding a new stamp + + psTrace("pois.findStamps", 7, "Stamp %d: %d,%d\n", num, stamp->x, stamp->y); + num++; + } + } // Done iterating over sections + + return stamps; +} Index: /tags/rel-0-1-0/pois/src/poisKernelBasisFunctions.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisKernelBasisFunctions.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisKernelBasisFunctions.c (revision 11633) @@ -0,0 +1,97 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" +#include "pois.h" + +psArray *poisKernelBasisFunctions(const poisConfig *config // Configuration + ) +{ + assert(config); + + psTrace("pois.kernelBasisFunctions", 1, "Setting kernel basis functions....\n"); + + int xKernelHalfSize = config->xKernel; // Half size of kernel in x + int yKernelHalfSize = config->yKernel; // Half size of kernel in y + int spatialOrder = config->spatialOrder; // Order of spatial variations + + // Number of basis functions + int nBF = (2 * xKernelHalfSize + 1) * (2 * yKernelHalfSize + 1) * (spatialOrder + 1) * + (spatialOrder + 2) / 2; + + psArray *array = psArrayAlloc(nBF); // Array to hold the basis functions + + int num = 0; // Kernel parameter number + + // Put the (u,v) = (0,0) component right at the start of the list for convenience + for (int order = 0; order <= spatialOrder; order++) { + for (int i = 0; i <= order; i++) { + int j = order - i; + poisKernelBasis *kernel = psAlloc(sizeof(poisKernelBasis)); + kernel->u = 0; + kernel->v = 0; + kernel->xOrder = i; + kernel->yOrder = j; + + // Stuff into the array + array->data[num] = kernel; + + psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", i, kernel->u, + kernel->v, kernel->xOrder, kernel->yOrder); + num++; + } + } + + // Iterate over (u,v) + for (int u = -xKernelHalfSize; u <= xKernelHalfSize; u++) { + for (int v = -yKernelHalfSize; v <= yKernelHalfSize; v++) { + + // Already did (u,v) = (0,0): it's at the start, so skip it now. + if ((u != 0) || (v != 0)) { + + // Iterate over spatial order + for (int order = 0; order <= spatialOrder; order++) { + + for (int i = 0; i <= order; i++) { + int j = order - i; + poisKernelBasis *kernel = psAlloc(sizeof(poisKernelBasis)); + + kernel->u = u; + kernel->v = v; + kernel->xOrder = i; + kernel->yOrder = j; + + // Stuff into the array + array->data[num] = kernel; + + psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", num, + kernel->u, kernel->v, kernel->xOrder, kernel->yOrder); + num++; + } + } // Iterating over spatial order + } // Except (u,v) = (0,0) + } + } // Iterating over (u,v) + + return array; +} + Index: /tags/rel-0-1-0/pois/src/poisMakeMask.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisMakeMask.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisMakeMask.c (revision 11633) @@ -0,0 +1,83 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" +#include "pois.h" + +#define MAX(a,b) ((a) > (b) ? (a) : (b)) +#define MIN(a,b) ((a) < (b) ? (a) : (b)) + +psImage *poisMakeMask(psImage *mask, // Mask to update, or NULL + const psImage *refImage, // Reference image + const psImage *inImage, // Input image + const poisConfig *config // Configuration + ) +{ + assert(!mask || mask->type.type == PS_TYPE_U8); + assert(!mask || (mask->numCols == refImage->numCols && mask->numRows == refImage->numRows)); + assert(refImage->numCols == inImage->numCols && refImage->numRows == inImage->numRows); + assert(refImage->type.type == PS_TYPE_F32); + assert(inImage->type.type == PS_TYPE_F32); + assert(config); + + int nx = config->xImage; // Size in x + int ny = config->yImage; // Size in y + + psTrace("pois.makeMask", 5, "Creating mask: %dx%d\n", nx, ny); + + // Allocate the mask image, if necessary + if (mask == NULL) { + mask = psImageAlloc(nx, ny, PS_TYPE_U8); + // Initialise to zero + for (int j = 0; j < ny; j++) { + for (int i = 0; i < nx; i++) { + mask->data.U8[j][i] = POIS_MASK_OK; + } + } + } + + // Iterate over the images + for (int y = config->yKernel; y < ny - config->yKernel; y++) { + for (int x = config->xKernel; x < nx - config->xKernel; x++) { + // Mask all around something that's saturated in the reference, since it gets convolved + // But only mask single pixel in the input image, since it doesn't get convolved + if (refImage->data.F32[y][x] >= config->refSat || refImage->data.F32[y][x] <= config->refBad || + isnan(refImage->data.F32[y][x])) { + mask->data.U8[y][x] |= POIS_MASK_BAD; + int xlow = MAX(0, x - config->xKernel); + int xhigh = MIN(nx, x + config->xKernel); + int ylow = MAX(0, y - config->yKernel); + int yhigh = MIN(ny, y + config->yKernel); + for (int v = ylow; v <= yhigh; v++) { + for (int u = xlow; u <= xhigh; u++) { + mask->data.U8[v][u] |= POIS_MASK_NEAR_BAD; + } + } + } else if (inImage->data.F32[y][x] >= config->inSat || inImage->data.F32[y][x] <= config->inBad || + isnan(inImage->data.F32[y][x])) { + mask->data.U8[y][x] |= POIS_MASK_BAD; + } + } + } + + return mask; +} Index: /tags/rel-0-1-0/pois/src/poisMaskOps.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisMaskOps.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisMaskOps.c (revision 11633) @@ -0,0 +1,112 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +/* + * Some utility routines that should be in psLib + */ +#include +#include +#include +#include "pslib.h" +#include "pois.h" + +psImage *poisImageSetVal(psImage *restrict in, // Input image + const psC64 val // set to this value + ) +{ + assert (in != NULL); + + switch (in->type.type) { + case PS_TYPE_U8: + { + const int numCols = in->numCols; + psU8 *row0 = in->data.U8[0]; + for (int c = 0; c < numCols; c++) { + row0[c] = val; + } + for (int r = 1; r < in->numRows; r++) { + memcpy(in->data.U8[r], row0, numCols*sizeof(psU8)); + } + } + break; + case PS_TYPE_F32: + { + const int numCols = in->numCols; + psF32 *row0 = in->data.F32[0]; + for (int c = 0; c < numCols; c++) { + row0[c] = val; + } + for (int r = 1; r < in->numRows; r++) { + memcpy(in->data.F32[r], row0, numCols*sizeof(psF32)); + } + } + break; + default: + psError(POIS_ERR_UNSUPPORTED_TYPE, true, + "Type %d is not supported", in->type.type); + return NULL; + } + + return in; +} + +psImage *poisImageSetValInMask(psImage *out, // Image to update, or NULL + const psImage *restrict in, // Input image + const psImage *mask, // mask for image + const psC64 val, // set to this value + const unsigned long bits // set pixels where mask & bits != 0 + ) +{ + assert(mask != NULL); + assert(mask->type.type == PS_TYPE_U8); + const int numCols = mask->numCols; + const int numRows = mask->numRows; + assert(in != NULL); + assert(in->type.type == PS_TYPE_F32); + assert(in->numCols == numCols && in->numRows == numRows); + + // Allocate the output image, if necessary + if (out == NULL) { + out = psImageCopy(NULL, in, in->type.type); + } + + switch (in->type.type) { + case PS_TYPE_F32: + { + const int numCols = in->numCols; + for (int r = 0; r < in->numRows; r++) { + psF32 *row = out->data.F32[r]; // n.b. a copy of in + psU8 *mrow = mask->data.U8[r]; + for (int c = 0; c < numCols; c++) { + if (mrow[c] & bits) { + row[c] = val; + } + } + } + } + break; + default: + psError(POIS_ERR_UNSUPPORTED_TYPE, true, + "Type %d is not supported", in->type.type); + return NULL; + } + + return out; +} Index: /tags/rel-0-1-0/pois/src/poisParseConfig.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisParseConfig.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisParseConfig.c (revision 11633) @@ -0,0 +1,207 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include "pois.h" + +void help(void) +{ + fprintf (stderr, "POIS: Pan-STARRS (or Pricey's) Optimal Image Subtraction\n" + "Usage: pois [-h] [-v] [-k X Y] [-t T] [-s S1 S2] [-b B1 B2] [-B BG] [-o N] [-n NSTAMPS] [-S STAMPS] [-i ITER] [-r REJ] [-m] REF IN OUT\n" + "where\n" + "\t-h Help (this info)\n" + "\t-v Increase verbosity level\n" + "\t-k X Y Kernel half-size in x and y (5)\n" + "\t-t T Threshold for stamps on reference image (0)\n" + "\t-s S1 S2 Saturation level for reference (S1) and input (S2) (50000)\n" + "\t-b B1 B2 Bad level for both images (0)\n" + "\t-B BG Level to add to background to get above zero (0)\n" + "\t-o ORDER Order of spatial variations in the kernel (0)\n" + "\t-n NSTAMPS Number of stamps in each dimension (5)\n" + "\t-S STAMPS File from which to read stamps\n" + "\t-R Reverse the sense of the subtraction (REF - IN)\n" + "\t-i ITER Number of rejection iterations (10)\n" + "\t-r REJ Rejection level in standard deviations (2.5)\n" + "\t-m Output mask frame (OUT.mask)\n" + "\tREF Reference image\n" + "\tIN Input image (to be matched to REF)\n" + "\tOUT Output image\n" + ); +} + +poisConfig *poisParseConfig(int argc, // Number of command-line arguments + char **argv // Command-line arguments + ) +{ + poisConfig *config; // Configuration values + + /* Variables for getopt */ + int opt; /* Option, from getopt */ + extern char *optarg; /* Argument accompanying switch */ + extern int optind; /* getopt variables */ + + /* Initialise configuration */ + config = (poisConfig *) psAlloc(sizeof(poisConfig)); + config->verbose = 0; + config->refFile = NULL; + config->inFile = NULL; + config->outFile = NULL; + config->xKernel = 5; + config->yKernel = 5; + config->footprint = 15; + config->threshold = 0.0; + config->refSat = 60000.0; + config->inSat = 60000.0; + config->refBad = 0.0; + config->inBad = 0.0; + config->background = 0.0; + config->spatialOrder = 0; + config->nsx = 5; + config->nsy = 5; + config->numIter = 10; + config->sigmaRej = 2.5; + config->penalty = 0.0; + config->stampFile = NULL; + config->reverse = false; + config->mask = false; + + /* Parse command-line arguments using getopt */ + while ((opt = getopt(argc, argv, "hvRmk:f:q:t:s:o:b:B:n:i:r:p:S:")) != -1) { + switch (opt) { + case 'h': + help(); + exit(EXIT_SUCCESS); + case 'v': + config->verbose++; + break; + case 'k': + if ((sscanf(argv[optind-1], "%d", &config->xKernel) != 1) || + (sscanf(argv[optind++], "%d", &config->yKernel) != 1)) { + /* + Note: incrementing optind, so I can read more than + one parameter. + */ + help(); + exit(EXIT_FAILURE); + } + break; + case 'f': + if (sscanf(optarg, "%d", &config->footprint) != 1) { + help(); + exit(EXIT_FAILURE); + } + break; + case 't': + if (sscanf(optarg, "%f", &config->threshold) != 1) { + help(); + exit(EXIT_FAILURE); + } + break; + case 's': + if ((sscanf(argv[optind-1], "%f", &config->refSat) != 1) || + (sscanf(argv[optind++], "%f", &config->inSat) != 1)) { + /* + Note: incrementing optind, so I can read more than + one parameter. + */ + help(); + exit(EXIT_FAILURE); + } + break; + case 'o': + if (sscanf(optarg, "%d", &config->spatialOrder) != 1) { + help(); + exit(EXIT_FAILURE); + } + break; + case 'b': + if ((sscanf(argv[optind-1], "%f", &config->refBad) != 1) || + (sscanf(argv[optind++], "%f", &config->inBad) != 1)) { + /* + Note: incrementing optind, so I can read more than + one parameter. + */ + help(); + exit(EXIT_FAILURE); + } + break; + case 'B': + if (sscanf(optarg, "%f", &config->background) != 1) { + help(); + exit(EXIT_FAILURE); + } + break; + case 'n': + if (sscanf(argv[optind-1], "%d", &config->nsx) != 1) { + help(); + exit(EXIT_FAILURE); + } + config->nsy = config->nsx; + break; + case 'i': + if (sscanf(argv[optind-1], "%d", &config->numIter) != 1) { + help(); + exit(EXIT_FAILURE); + } + break; + case 'r': + if (sscanf(argv[optind-1], "%f", &config->sigmaRej) != 1) { + help(); + exit(EXIT_FAILURE); + } + break; + case 'p': + if (sscanf(argv[optind-1], "%f", &config->penalty) != 1) { + help(); + exit(EXIT_FAILURE); + } + break; + case 'S': + config->stampFile = argv[optind-1]; + config->threshold = -INFINITY; + break; + case 'R': + config->reverse = true; + break; + case 'm': + config->mask = true; + break; + default: + help(); + } + } + + /* Get the remaining, mandatory values */ + argc -= optind; + argv += optind; + if (argc < 3) { + help(); + exit(EXIT_FAILURE); + } + + config->refFile = argv[0]; + config->inFile = argv[1]; + config->outFile = argv[2]; + + return config; +} + Index: /tags/rel-0-1-0/pois/src/poisReadStamps.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisReadStamps.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisReadStamps.c (revision 11633) @@ -0,0 +1,81 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include "pslib.h" +#include "pois.h" + +#define MAXLINE 128 // Maximum length of a line in the stamps file + + +psArray *poisReadStamps(psArray *stamps,// Stamps + FILE **stampFP, // File pointer for stamps file + const psImage *image, // Image for which to find stamps + const psImage *mask, // Mask image + const poisConfig *config // Configuration + ) +{ + FILE *fp = *stampFP; // File pointer + + if (stamps == NULL) { + stamps = psArrayAlloc(config->nsx); + // Initialise + for (int i = 0; i < config->nsx; i++) { + stamps->data[i] = poisStampAlloc(); + } + } + + for (int i = 0; i < stamps->n; i++) { + poisStamp *stamp = stamps->data[i]; + + if (stamp->status == POIS_STAMP_RESET) { + if (fp) { + float x, y; // Position of stamp + int xPixel, yPixel; // Integer position + int result = 0; // Result of fscanf + char line[MAXLINE]; // The line from the stamps file + char *lineBack; // The line back from fgets + do { + lineBack = fgets(line, MAXLINE, fp); + result = sscanf(line, "%f %f", &x, &y); + + // Add 0.5 for rounding, and subtract 1.0 for unit-offset to zero-offset + xPixel = (int)(x - 0.5); + yPixel = (int)(y - 0.5); + } while (result == 2 && lineBack && ! poisCheckStamp(image, mask, xPixel, yPixel, config)); + if (result == 2 && lineBack) { + stamp->x = xPixel; + stamp->y = yPixel; + stamp->status = POIS_STAMP_RECALC; + psTrace("pois.readStamps", 7, "Stamp %d: %d,%d\n", i, stamp->x, stamp->y); + } else { + // We've run out of stamps in the file --- probably got EOF + stamp->status = POIS_STAMP_BAD; + } + } else { + // Stamp file isn't open! + stamp->status = POIS_STAMP_BAD; + } + } + } + + return stamps; +} + Index: /tags/rel-0-1-0/pois/src/poisRejectStamps.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisRejectStamps.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisRejectStamps.c (revision 11633) @@ -0,0 +1,86 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" +#include "pois.h" + +#define SQUARE(x) ((x)*(x)) + +bool poisRejectStamps(psArray *stamps, // Array of stamps + psImage *mask, // Mask image + const psVector *deviations, // Vector of deviations for the stamps + const poisConfig *config // Configuration + ) +{ + assert(stamps); + assert(mask); + assert(deviations); + assert(config); + assert(stamps->n == deviations->n); + assert(mask->numCols == config->xImage && mask->numRows == config->yImage); + assert(mask->type.type == PS_TYPE_U8); + assert(deviations->type.type == PS_TYPE_F32); + + bool masked = false; // Have we masked any stamps? + + // Don't want the standard deviation, but rather the deviation from zero + double sum = 0.0; + int numDev = 0; + for (int i = 0; i < deviations->n; i++) { + poisStamp *stamp = stamps->data[i]; // The stamp + // Only interested in the stamps that we're actually using + if (stamp->status == POIS_STAMP_USED) { + sum += SQUARE(deviations->data.F32[i]); + numDev++; + } + } + float meanDev = sqrtf(sum / (float)numDev); + float limit = meanDev * config->sigmaRej; + psTrace("pois.rejectStamps", 2, "Mean RMS deviation: %f\n", meanDev); + psTrace("pois.rejectStamps", 2, "Rejection limit: %f\n", limit); + + // Reject stamps + for (int s = 0; s < deviations->n; s++) { + poisStamp *stamp = stamps->data[s]; + if (stamp->status == POIS_STAMP_USED && + (fabsf(deviations->data.F32[s]) > limit || ! isfinite(deviations->data.F32[s]))) { + masked = true; + psTrace("pois.rejectStamps", 1, "Rejecting stamp %d (%d,%d): %f\n", s, stamp->x, stamp->y, + deviations->data.F32[s]); + + // Mask out the stamp in the image so you don't find it again + for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) { + for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) { + mask->data.U8[y][x] |= POIS_MASK_STAMP; + } + } + + // Set stamp for replacement + stamp->x = 0; + stamp->y = 0; + stamp->status = POIS_STAMP_RESET; + + } // Bad stamps + } // Iterating over stamps + + return masked; +} Index: /tags/rel-0-1-0/pois/src/poisSolveEquation.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisSolveEquation.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisSolveEquation.c (revision 11633) @@ -0,0 +1,89 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" +#include "pois.h" + +psVector *poisSolveEquation(psVector *solution, // Solution vector, or NULL + const psArray *stamps, // Array of stamps + const poisConfig *config + ) +{ + assert(stamps); + assert(config); + assert(!solution || solution->type.type == PS_TYPE_F64); + + int size = (2 * config->xKernel + 1) * (2 * config->yKernel + 1) * (config->spatialOrder + 1) * + (config->spatialOrder + 2) / 2 + 1; // Size of matrix and vector + + assert(!solution || solution->n == size); + + psImage *matrix = psImageAlloc(size, size, PS_TYPE_F64); + psVector *vector = psVectorAlloc(size, PS_TYPE_F64); + for (int j = 0; j < size; j++) { + for (int i = 0; i < size; i++) { + matrix->data.F64[j][i] = 0.0; + } + vector->data.F64[j] = 0.0; + } + + for (int s = 0; s < stamps->n; s++) { + poisStamp *stamp = stamps->data[s]; // Stamp of interest + psImage *stampMatrix = stamp->matrix; // Corresponding matrix + psVector *stampVector = stamp->vector; // Corresponding vector + + // Check inputs + assert(matrix->type.type == PS_TYPE_F64); + assert(vector->type.type == PS_TYPE_F64); + assert(matrix->numCols == size && matrix->numRows == size); + assert(vector->n == size); + + if (stamp->status == POIS_STAMP_USED) { + (void)psBinaryOp(matrix, matrix, "+", stampMatrix); + (void)psBinaryOp(vector, vector, "+", stampVector); + } + } + + // Solution via LU Decomposition + psVector *permutation = NULL; // Permutation vector for LU decomposition + psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix + solution = psMatrixLUSolve(solution, luMatrix, vector, permutation); // Solution in x + psTrace("pois.solveEquation", 1, "Background difference is %f\n", solution->data.F64[solution->n - 1]); + +#if 0 + printf("Matrix:\n"); + for (int i = 0; i < matrix->numCols; i++) { + for (int j = 0; j < matrix->numRows; j++) { + printf("%f ", matrix->data.F64[i][j]); + } + printf("\n"); + } +#endif + + // Clean up + psFree(luMatrix); + psFree(permutation); + psFree(matrix); + psFree(vector); + + return solution; +} Index: /tags/rel-0-1-0/pois/src/poisStamp.c =================================================================== --- /tags/rel-0-1-0/pois/src/poisStamp.c (revision 11633) +++ /tags/rel-0-1-0/pois/src/poisStamp.c (revision 11633) @@ -0,0 +1,81 @@ +// Copyright (C) 2006 Paul A. Price (price@ifa.hawaii.edu) +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include "pslib.h" +#include "pois.h" + +poisStamp *poisStampAlloc(void) +{ + poisStamp *stamp = (poisStamp*)psAlloc(sizeof(poisStamp)); + stamp->x = 0; + stamp->y = 0; + stamp->status = POIS_STAMP_RESET; + stamp->matrix = NULL; + stamp->vector = NULL; + psMemSetDeallocator(stamp, (psFreeFunc)poisStampFree); + + return stamp; +} + +void poisStampFree(poisStamp *stamp) +{ + if (stamp->matrix) { + psFree(stamp->matrix); + } + if (stamp->vector) { + psFree(stamp->vector); + } + psFree(stamp); +} + + +// Return true if the stamp is OK +bool poisCheckStamp(const psImage *image, // Image to check for threshold + const psImage *mask, // Mask to check for bad pixels + int x, int y, // Pixel coordinates + const poisConfig *config // Configuration + ) +{ + if (x < config->footprint + config->xKernel || + y < config->footprint + config->yKernel || + x > config->xImage - config->footprint - config->xKernel || + y > config->yImage - config->footprint - config->yKernel || + image->data.F32[y][x] < config->threshold) { + psTrace("pois.checkStamp", 9, "Rejecting stamp at %d,%d (border/threshold)\n", x, y); + return false; + } + + // Check the footprint + bool ok = true; // Is the footprint OK? + int footprint = config->footprint; // Footprint size + for (int v = y - footprint; v <= y + footprint && ok; v++) { + for (int u = x - footprint; u <= x + footprint && ok; u++) { + if (mask->data.U8[v][u] & + (POIS_MASK_BAD | POIS_MASK_NEAR_BAD | POIS_MASK_STAMP)) { + psTrace("pois.checkStamp", 9, "Rejecting stamp at %d,%d (bad footprint)\n", x, y); + ok = false; + } + } + } + + return ok; +} + Index: /tags/rel-0-1-0/pois/swig/.cvsignore =================================================================== --- /tags/rel-0-1-0/pois/swig/.cvsignore (revision 11633) +++ /tags/rel-0-1-0/pois/swig/.cvsignore (revision 11633) @@ -0,0 +1,9 @@ +Makefile Makefile.in +*.pm *.pyc +*_wrap.c +poisSwig.py psLibSwig.py xpaSwig.py +swigrun.py +xpa.py +.Makefile.depend +*.la *.lo +.deps .libs Index: /tags/rel-0-1-0/pois/swig/Makefile.am =================================================================== --- /tags/rel-0-1-0/pois/swig/Makefile.am (revision 11633) +++ /tags/rel-0-1-0/pois/swig/Makefile.am (revision 11633) @@ -0,0 +1,82 @@ +AM_CFLAGS = @pois_CFLAGS@ + +AM_CFLAGS += $(PSLIB_CFLAGS) +AM_CFLAGS += $(SWIG_TARGET_CPPFLAGS) +AM_CFLAGS += -DPS_ALLOW_MALLOC +AM_CFLAGS += -I../src + +LIBS = $(AM_LIBS) +LIBS += $(SWIG_TARGET_LIBS) + +SWIG += -w451 +SWIG += $(SWIG_TARGET_OPT) +SWIG += $(PSLIB_CFLAGS) +SWIG += -Drestrict= -Dinline= + +AM_CFLAGS += -I../src +# +# Libraries to install +# +SWIG_LIB_NAMES = poisSwig psLibSwig xpaSwig +lib_LTLIBRARIES = libpoisSwig.la libpsLibSwig.la libxpaSwig.la + +BUILT_SOURCES = \ + psLibSwig.py psLibSwig_wrap.c \ + xpaSwig.py xpaSwig_wrap.c +# +# PSLib wrappers +# +libpsLibSwig_la_LDFLAGS = -module +libpsLibSwig_la_LIBADD = $(PSLIB_LIBS) + +libpsLibSwig_la_SOURCES = \ + psLibSwig_wrap.c \ + psLibForSwig.c psLibForSwig.h + +psLibSwig.py psLibSwig_wrap.c : psLibSwig.i p_psSwig.i + $(SWIG) $(PS_INCLUDES) psLibSwig.i +# +# Now for XPA +# +libxpaSwig_la_LIBADD = $(PSLIB_LIBS) +libxpaSwig_la_LIBADD += -lxpa + +libxpaSwig_la_LDFLAGS = -module +libxpaSwig_la_SOURCES = \ + xpaSwig_wrap.c \ + simpleFits.c + +xpaSwig.py xpaSwig_wrap.c : xpaSwig.i p_psSwig.i + $(SWIG) $(PS_INCLUDES) xpaSwig.i + @case "$(SWIG)" in \ + *-perl*) \ + perl -pi -e 's/([a-zA-Z0-9_]+Array)/_xpa_\1/g' xpaSwig_wrap.c;; \ + esac +# +# And the pois modules +# +libpoisSwig_la_LDFLAGS = -module +libpoisSwig_la_SOURCES = \ + poisSwig_wrap.c +libpoisSwig_la_LIBADD = ../src/libpois.la + +poisSwig.py poisSwig_wrap.c : poisSwig.i ../src/pois.h p_psSwig.i + $(SWIG) $(PS_INCLUDES) -I../src poisSwig.i + @case "$(SWIG)" in \ + *-perl*) \ + perl -pi -e 's/([a-zA-Z0-9_]+Array)/_pois_\1/g' poisLib_wrap.c;; \ + esac +# +# Patch up name of .so file for the sake of python +# +install-exec-hook: + cd $(prefix)/lib && \ + for f in $(SWIG_LIB_NAMES); do \ + $(LN_S) -f `perl -ne "if(/dlname='([^']+)'/) { print \\$$1; exit }" lib$$f.la` _$$f.so; \ + done + +uninstall-local: + for f in $(SWIG_LIB_NAMES); do $(RM) $(prefix)/lib/_$$f.so; done + +CLEANFILES = Makefile Makefile.in \ + *_wrap.c *.pyc poisSwig.py psLibSwig.py xpaSwig.py *~ core core.* Index: /tags/rel-0-1-0/pois/swig/p_psSwig.i =================================================================== --- /tags/rel-0-1-0/pois/swig/p_psSwig.i (revision 11633) +++ /tags/rel-0-1-0/pois/swig/p_psSwig.i (revision 11633) @@ -0,0 +1,267 @@ +// -*- panstarrs-c -*- +%nodefault; + +%include "exception.i" +#if !defined(SWIGXML) +%include "typemaps.i" +#endif +%include "carrays.i" + +%array_class(int, intArray); +%array_class(float, floatArray); +%array_class(double, doubleArray); + +/************************************************************************************************************/ + +#if defined(SWIGPYTHON) +%pythoncode { +class PsUtilsError(StandardError): + """An error generated by psUtils""" +} +#endif + +/************************************************************************************************************/ + +%{ +# include "pslib.h" + +static void *exception(int exception_type, const char *fmt, ...) +{ + char message[256]; + va_list ap; + + va_start(ap, fmt); + vsnprintf(message, 256, fmt, ap); + va_end(ap); + + SWIG_exception(exception_type, message) + + fail: // used by SWIG_exception + return NULL; +} + +static void *ps_exception(void) +{ + psErr *err = psErrorLast(); + char message[256]; + snprintf(message, 256, "(%s) %s", err->name, err->msg); + psFree(err); + + SWIG_exception(SWIG_RuntimeError, message) + + fail: // used by SWIG_exception + return NULL; +} + +#if 0 +static PyObject *set_thisown(PyObject *resultobj, int val) +{ + PyObject *n = PyInt_FromLong(val); + PyObject_SetAttrString(resultobj, (char *)"thisown", n); + Py_DECREF(n); + + return resultobj; +} +#endif + +%} + +/************************************************************************************************************/ +/* + * Typemaps + */ +%typemap(python, out) psStatus { + if (result == 0) { + Py_INCREF(Py_None); resultobj = Py_None; + } else { + ps_exception(); + return NULL; + } +} + +%typemap(python,in) FILE * { + if ($input == Py_None) { + $1 = NULL; + } else if (!PyFile_Check($input)) { + PyErr_SetString(PyExc_TypeError, "Need a file!"); + goto fail; + } else { + $1 = PyFile_AsFile($input); + } +} + +%typemap(python,in) char * { + if ($input == Py_None) { + $1 = NULL; + } else if (!PyString_Check($input)) { + PyErr_SetString(PyExc_TypeError, "Need a string!"); + goto fail; + } else { + $1 = PyString_AsString($input); + } +} + +/************************************************************************************************************/ +/* + * Tell our extension language that this function returns a new object (that can't be NULL) + */ +%define NEWOBJECT(func) + %newobject func; + NOTNULL(func); +%enddef + +/************************************************************************************************************/ +/* + * Tell our extension language that it should call psFree to delete an object + */ +%define EXPORT_TYPE(type) + %extend type { + ~type() { + if (psMemGetRefCounter(self) == 1) { + const psMemBlock *mb = ((psMemBlock *)self) - 1; + +#if !defined(PS_NO_TRACE) + p_psTrace("ps.utils.memory", 6, "Freeing %p (%ld: %s:%d)\n", self, (long)mb->id, + mb->file, mb->lineno); +#endif + } + psFree(self); + } + } + + %{ + type *type ## _Cast(void *ptr) { + return (type *)ptr; + } + %} + type *type ## _Cast(void *ptr); + + NEWOBJECT(type ## Alloc); +%enddef + +/************************************************************************************************************/ +/* + * Throw an exception if func returns NULL + */ +%define NOTNULL(func) + %exception func { + $action; + if (result == NULL) { + $cleanup; +#if defined(SWIGPERL) + ps_exception(); +#else + return ps_exception(); +#endif + } + } +%enddef + +/* + * Throw an exception if func returns a negative value + */ +%define NOTNEGATIVE(func) + %exception func { + $action; + if (result < 0) { + $cleanup; +#if defined(SWIGPERL) + ps_exception(); +#else + return ps_exception(); +#endif + } + } +%enddef + +/* + * Throw an exception if func returns non-zero + */ +%define ISZERO(func) + %exception func { + $action; + if (result != 0) { + $cleanup; +#if defined(SWIGPERL) + ps_exception(); +#else + return ps_exception(); +#endif + } + } +%enddef + +/************************************************************************************************************/ +/* + * If an argument is passed in, and also returned by the function, we need to declare + * the function %newobject, and also ensure that the input function has its refCounter + * incremented (as perl/python etc. will think that it's gone out of scope) + * + * This macro defined a typemap to handle the refcounter. + * + * Do this as an argout map as we don't get another reference to the image if there's an error + */ +%define TYPE_INOUT(type, argname) + %typemap(argout) type argname { + if ($1 != NULL) { + const psMemBlock *mb = ((psMemBlock *)$1) - 1; + p_psTrace("ps.utils.memory", 4, "Incrementing reference counter for %p (%ld %s:%d)\n", + $1, (long)mb->id, mb->file, mb->lineno); + psMemIncrRefCounter($1); + } + } +%enddef + +%define CLEAR_TYPE_INOUT(type, argname) + %typemap(argout) type argname; +%enddef + +/************************************************************************************************************/ + +%exception; + +/************************************************************************************************************/ +/* + * Accept a python list of floats if the routine wants an array + */ +%typemap(in) (const float *floatArr, int nel) { + /* Check if is a list */ + if (PyList_Check($input)) { + $2 = PyList_Size($input); + if ($2 == 0) { + $1 = NULL; + } else { + $1 = psAlloc($2*sizeof(float)); + + for (int i = 0; i < $2; i++) { + PyObject *o = PyList_GetItem($input,i); + if (PyFloat_Check(o)) + $1[i] = PyFloat_AsDouble(PyList_GetItem($input, i)); + else { + PyErr_SetString(PyExc_TypeError, "list must contain floats"); + psFree($1); + return NULL; + } + } + } + } else { + PyErr_SetString(PyExc_TypeError,"not a list"); + return NULL; + } +} + +%typemap(freearg) (const float *floatArr, int nel) { + psFree($1); +} + +/************************************************************************************************************/ + +%apply double *OUTPUT { double *irow }; +%apply double *OUTPUT { double *icol }; +%apply double *OUTPUT { double *xPos }; +%apply double *OUTPUT { double *yPos }; +%apply double *OUTPUT { double *alpha }; +%apply double *OUTPUT { double *delta }; +%apply double *OUTPUT { double *alt }; +%apply double *OUTPUT { double *az }; +%apply double *OUTPUT { double *lst }; Index: /tags/rel-0-1-0/pois/swig/poisSwig.i =================================================================== --- /tags/rel-0-1-0/pois/swig/poisSwig.i (revision 11633) +++ /tags/rel-0-1-0/pois/swig/poisSwig.i (revision 11633) @@ -0,0 +1,33 @@ +// -*- C -*- +/* + */ +%module poisSwig + +%{ +#include "pois.h" +%} + +%include "p_psSwig.i" +%import "psLibSwig.i" + +EXPORT_TYPE(poisConfig); +EXPORT_TYPE(poisStamp); +EXPORT_TYPE(poisKernelBasis); +NEWOBJECT(poisKernelBasisFunctions); +NEWOBJECT(poisMakeMask); +NEWOBJECT(poisFindStamps); +NEWOBJECT(poisExtractKernel); +NEWOBJECT(poisConvolveImage); +NEWOBJECT(poisCalculateDeviations); +NEWOBJECT(poisSolveEquation); +%newobject poisRejectStamps; +NEWOBJECT(poisImageSetVal); +NEWOBJECT(poisImageSetValInMask); + +TYPE_INOUT(psArray *, stampsIO); +TYPE_INOUT(psImage *, IO); +TYPE_INOUT(psVector *, solutionIO) +%include "pois.h" +CLEAR_TYPE_INOUT(psArray *, stampsIO); +CLEAR_TYPE_INOUT(psImage *, IO); +CLEAR_TYPE_INOUT(psVector *, solutionIO) Index: /tags/rel-0-1-0/pois/swig/psLibForSwig.c =================================================================== --- /tags/rel-0-1-0/pois/swig/psLibForSwig.c (revision 11633) +++ /tags/rel-0-1-0/pois/swig/psLibForSwig.c (revision 11633) @@ -0,0 +1,41 @@ +/* + * Wrappers for psLib + */ +#include +#include "pslib.h" + +void +psErrorStackPrintForSwig(void) +{ + psErrorStackPrint(stdout, ""); +} + +const psErr *psErrorGetForSwig(int i) +{ + const psErr *err = psErrorGet(i); + + return (err->code == PS_ERR_NONE ? NULL : err); +} + +/************************************************************************************************************/ + +#if defined(psTrace) +# undef psTrace +#endif + +void psTrace(const char *comp, // component being traced + int level, // desired trace level + char *str) // output string +{ + p_psTrace(comp, level, str); +} + +/************************************************************************************************************/ +/* + * SWIG doesn't like varargs functions; so provide non-varargs versions of psMetadataAdd + */ +bool psMetadataAddBool(psMetadata* restrict md, int where, const char *name, const char *comment, int val) +{ + return psMetadataAdd(md, where, name, PS_TYPE_BOOL, comment, val); +} + Index: /tags/rel-0-1-0/pois/swig/psLibForSwig.h =================================================================== --- /tags/rel-0-1-0/pois/swig/psLibForSwig.h (revision 11633) +++ /tags/rel-0-1-0/pois/swig/psLibForSwig.h (revision 11633) @@ -0,0 +1,14 @@ +#if !defined(PS_LIB_FOR_SWIG_H) +#define PS_LIB_FOR_SWIG_H + +void psErrorStackPrintForSwig(void); +const psErr *psErrorGetForSwig(int i); + +#if defined(psTrace) +# undef psTrace +#endif +void psTrace(const char *comp, int level, char *str); + +bool psMetadataAddBool(psMetadata* restrict md, int where, const char *name, const char *comment, int val); + +#endif Index: /tags/rel-0-1-0/pois/swig/psLibSwig.i =================================================================== --- /tags/rel-0-1-0/pois/swig/psLibSwig.i (revision 11633) +++ /tags/rel-0-1-0/pois/swig/psLibSwig.i (revision 11633) @@ -0,0 +1,352 @@ +// -*- C -*- +%module psLibSwig + +%include "p_psSwig.i" + +%{ +# include "pslib.h" +#if 0 +# include "Private/p_psMemory.h" /* not available from prototype psLib include files (in src dir) */ +#else +# define P_PS_MEMMAGIC (void *)0xdeadbeef // Magic number in psMemBlock header +#endif +# include "psLibForSwig.h" +%} + +/************************************************************************************************************/ +/* + * Typemaps + */ +%typemap(in) FILE * { +#if defined(SWIGPERL) + if (!SvOK($input)) { + $1 = NULL; + } else { + $1 = PerlIO_findFILE(IoIFP(sv_2io($input))); + } +#elif defined(SWIGPYTHON) + if ($input == Py_None) { + $1 = NULL; + } else if (!PyFile_Check($input)) { + PyErr_SetString(PyExc_TypeError, "Need a file!"); + goto fail; + } else { + $1 = PyFile_AsFile($input); + } +#endif +} + +/************************************************************************************************************/ + +%typemap(in) psMemBlock ***arr (psMemBlock **arr) { +#if defined(SWIGPERL) + if (!SvOK($input)) { + $1 = NULL; + } else { + arr = NULL; + $1 = &arr; + } + +#elif defined(SWIGPYTHON) + if ($input == Py_None) { + $1 = NULL; + } else { + if (!PyList_Check($input)) { + SWIG_exception(SWIG_RuntimeError, "\"$1_name\" isn't a list"); + } + + arr = NULL; + $1 = &arr; + } +#endif +} + +%typemap(in) psMetadata ** (psMetadata *) { + if +#if defined(SWIGPYTHON) + ($input == Py_None) +#else + ($input == NULL) +#endif + { + $1 = NULL; + } else { + static $*1_type arr[1]; + if ((SWIG_ConvertPtr($input, (void **) &$1, $*1_descriptor, 1)) == -1) return NULL; + arr[0] = ($*1_type)$1; + $1 = arr; + } +} + +/* + * Return a (psMemBlock ***) as a list of tuples; they are _appended_ to the pre-existing list. E.g. + * + * arr = []; IPP.psMemCheckLeaks(0, arr, None, False) + * for el in arr: + * print el + */ +%typemap(argout) psMemBlock ***arr { + if ($1 != NULL) { + psMemBlock **arr = *$1; + + for (int i = 0; i < result; i++) { + const psMemBlock *m = arr[i]; +#if defined(SWIGPERL) + SV *sv = newSVpvf("%d:%ld:%s:%ld:%d:%d:%ld", + (m->startblock == P_PS_MEMMAGIC) ? 1 : 0, (long)m->id, + m->file, (long)m->lineno, (int)m->refCounter, + (m->endblock == P_PS_MEMMAGIC) ? 1 : 0, (long)(m + 1)); + + if (SvROK($input)) { + AV *av = (AV *)SvRV($input); + if (SvTYPE(av) != SVt_PVAV) { + SWIG_croak("Type error in argument $argnum of $symname; Expected an array of psMemBlock"); + } + + av_push(av, sv); + } else { + SWIG_croak("Type error in argument $argnum of $symname; Expected a reference to an array of psMemBlock"); + } +#elif defined(SWIGPYTHON) + const int nel = 7; + PyObject *ty = PyTuple_New(nel); + int j = 0; + PyTuple_SET_ITEM(ty, j++, PyInt_FromLong(m->startblock == P_PS_MEMMAGIC)); + PyTuple_SET_ITEM(ty, j++, PyInt_FromLong(m->id)); + PyTuple_SET_ITEM(ty, j++, PyString_FromString(m->file)); + PyTuple_SET_ITEM(ty, j++, PyInt_FromLong(m->lineno)); + PyTuple_SET_ITEM(ty, j++, PyInt_FromLong(m->refCounter)); + PyTuple_SET_ITEM(ty, j++, PyInt_FromLong(m->endblock == P_PS_MEMMAGIC)); + PyTuple_SET_ITEM(ty, j++, PyInt_FromLong((long)(m + 1))); + assert (j == nel); + + PyList_Append($input, ty); +#endif + } + + psFree(arr); + } +} + +/************************************************************************************************************/ +%{ +typedef struct { float re, im; } swig_psC32; ///< 32-bit complex floating point +typedef struct { double re, im; } swig_psC64; ///< 64-bit complex floating point +%} + +#if 0 +%include // SWIG has trouble with nested typedefs in this file + // at least on Darwin machines +#else +typedef unsigned char uint8_t; +typedef unsigned short uint16_t; +typedef unsigned int uint32_t; +typedef unsigned long uint64_t; +typedef char int8_t; +typedef short int16_t; +typedef int int32_t; +typedef long int64_t; +#endif + +typedef int psBool; +typedef void *psPtr; +typedef uint8_t psU8; ///< 8-bit unsigned int +typedef uint16_t psU16; ///< 16-bit unsigned int +typedef uint32_t psU32; ///< 32-bit unsigned int +typedef uint64_t psU64; ///< 64-bit unsigned int +typedef int8_t psS8; ///< 8-bit signed int +typedef int16_t psS16; ///< 16-bit signed int +typedef int32_t psS32; ///< 32-bit signed int +typedef int64_t psS64; ///< 64-bit signed int +typedef float psF32; ///< 32-bit floating point +typedef double psF64; ///< 64-bit floating point +typedef struct { float re, im; } swig_psC32; ///< 32-bit complex floating point +typedef struct { double re, im; } swig_psC64; ///< 64-bit complex floating point + +%typemap(in) psC32 = float; +%typemap(in) psC64 = double; + +typedef enum { + PS_TYPE_S8 = 0x0101, ///< Character. + PS_TYPE_S16 = 0x0102, ///< Short integer. + PS_TYPE_S32 = 0x0104, ///< Integer. + PS_TYPE_S64 = 0x0108, ///< Long integer. + PS_TYPE_U8 = 0x0301, ///< Unsigned character. + PS_TYPE_U16 = 0x0302, ///< Unsigned short integer. + PS_TYPE_U32 = 0x0304, ///< Unsigned integer. + PS_TYPE_U64 = 0x0308, ///< Unsigned long integer. + PS_TYPE_F32 = 0x0404, ///< Single-precision Floating point. + PS_TYPE_F64 = 0x0408, ///< Double-precision floating point. + PS_TYPE_C32 = 0x0808, ///< Complex numbers consisting of single-precision floating point. + PS_TYPE_C64 = 0x0810 ///< Complex numbers consisting of double-precision floating point. +} psElemType; + +typedef enum { + PS_DIMEN_SCALAR, ///< Scalar. + PS_DIMEN_VECTOR, ///< Vector. + PS_DIMEN_TRANSV, ///< Transposed vector. + PS_DIMEN_IMAGE, ///< Image. + PS_DIMEN_OTHER ///< Something else that's not supported for arithmetic. +} psDimen; + +typedef struct { + psElemType type; ///< Primitive type. + psDimen dimen; ///< Dimensionality. +} +psType; + +%include "psErrorCodes.h" +%include "psLogMsg.h" +%include "psMemory.h" +%ignore psTraceFree; +%include "psTrace.h" + +EXPORT_TYPE(psList); +%include "psList.h" + +EXPORT_TYPE(psImage); +EXPORT_TYPE(psRegion); +%typemap(memberin) psImage * { + psFree($1); + $1 = $input; +#if defined(SWIGPERL) + p_psMemIncrRefCounter($1, "memberin(psImage *)", 0); +#endif +} +NEWOBJECT(psImageReadSection); +NEWOBJECT(psImageSubset); +NEWOBJECT(psImageSubsection); +NEWOBJECT(psImageCopy); +NOTNULL(psImageTrim); +%include "psImage.h" +%include "psImageIO.h" +%include "psImageExtraction.h" + +EXPORT_TYPE(psMetadata) +%typemap(memberin) psMetadata * { + psFree($1); + $1 = $input; +#if defined(SWIGPERL) + p_psMemIncrRefCounter($1, "memberin(psMetadata *)", 0); +#endif +} + +%include "psMetadata.h" + +%apply int *OUTPUT { psU32 *nFail } +TYPE_INOUT(psMetadata *, md); +NEWOBJECT(psMetadataParseConfig); +NOTNULL(psMetadataLookup); +%include "psMetadataIO.h" +CLEAR_TYPE_INOUT(psMetadata *, md); + +%extend psMetadataItem { + const char *get_STR(void) { + if (self->type != PS_META_STR) { + return NULL; + } else { + return self->data.V; + } + } +} + +%include "psLibForSwig.h" + +EXPORT_TYPE(psHash); +%include "psHash.h" + +EXPORT_TYPE(psErr); +NEWOBJECT(psErrorGet); +NEWOBJECT(psErrorLast); +%include "psError.h" + +EXPORT_TYPE(psArray); +%extend psArray { + const void *get_data(int i) { + return self->data[i]; + } + + const psCell *get_cell(int i) { + return (psCell *)(self->data[i]); + } + const psChip *get_chip(int i) { + return (psChip *)(self->data[i]); + } + const psReadout *get_readout(int i) { + return (psReadout *)(self->data[i]); + } +} +%include "psArray.h" + +// N.b. assumes e.g. %array_class(int, intArray); +%define GET_CARRAY(structType, ctype, typeName) + %extend structType { + ctype ## Array *get_data_ ## typeName(void) { + if (self->type.type == PS_TYPE_ ## typeName) { + return self->data.typeName; + } else { + psError(PS_ERR_BAD_PARAMETER_TYPE, 1, "psVector is of type %d, not PS_TYPE_" "typeName", + self->type.type); + return NULL; + } + } + } +%enddef + +EXPORT_TYPE(psVector); +%include "psVector.h" +GET_CARRAY(psVector, float, F32); +GET_CARRAY(psVector, double, F64); + +EXPORT_TYPE(psObservatory); +EXPORT_TYPE(psFPA); +EXPORT_TYPE(psChip); +EXPORT_TYPE(psCell); +EXPORT_TYPE(psReadout); +EXPORT_TYPE(psExposure); +%include "psAstrometry.h" + +%extend psFPA { + void set_chip(int i, psChip *chip) { + self->chips->data[i] = chip; + p_psMemIncrRefCounter(chip, "set_chip", 0); + } +} +%extend psChip { + void set_cell(int i, psCell *cell) { + self->cells->data[i] = cell; + p_psMemIncrRefCounter(cell, "set_cell", 0); + } +} +%extend psCell { + char *extname; + void set_readout(int i, psReadout *readout) { + self->readouts->data[i] = readout; + p_psMemIncrRefCounter(readout, "set_readout", 0); + } +} +%{ + char *psCell_extname_get(psCell *cell) { return NULL; } + void psCell_extname_set(psCell *cell, char *val) { return; } +%} + +NEWOBJECT(psBinaryOp); +TYPE_INOUT(void *, out); +psImage *psBinaryOp(void *out, void *in1, char *op, void *in2); +CLEAR_TYPE_INOUT(void *, out); +%include "psScalar.h" + +EXPORT_TYPE(psTime); +NEWOBJECT(psTimeGetTime); +%include "psTime.h" + +EXPORT_TYPE(psFits); +NEWOBJECT(psFitsReadImage); +TYPE_INOUT(psImage *, out); +%include "psFits.h" +CLEAR_TYPE_INOUT(psImage *, out); + +EXPORT_TYPE(psStats); +EXPORT_TYPE(psHistogram); +%include "psStats.h" +%include "psImageStats.h" Index: /tags/rel-0-1-0/pois/swig/simpleFits.c =================================================================== --- /tags/rel-0-1-0/pois/swig/simpleFits.c (revision 11633) +++ /tags/rel-0-1-0/pois/swig/simpleFits.c (revision 11633) @@ -0,0 +1,678 @@ +/* + * Write a fits image to a file descriptor; useful for talking to DS9 + * + * This version knows about psLib data structures + */ +#include +#include +#include +#include +#include +#include +#include +#include "pslib.h" +#include "poisErrorCodes.h" + +#include "simpleFits.h" + +#define FITS_SIZE 2880 + +typedef struct { + char keyword[9]; + enum { + UNKNOWN_CARD = -1, + LOGICAL_CARD, + STRING_CARD, + INTEGER_CARD, + DOUBLE_CARD, + FLOAT_CARD, + } type; + union { + char *l; + char *s; + int i; + double d; + float f; + } val; + char *commnt; +} CARD; + + +/* + * Write a string value + */ +static int +write_card_s(int fd, + int ncard, + char *record, + const char *keyword, + const char *val, + const char *commnt) +{ + char *card = &record[80*ncard]; + char value[20]; /* blank-padded val, if needed */ + + if(strlen(val) < 8) { /* FITS requires at least 8 chars */ + sprintf(value, "%-8s", val); + val = value; + } + + if(*keyword == '\0' || + !strcmp(keyword,"COMMENT") || !strcmp(keyword,"END") || !strcmp(keyword,"HISTORY")) { + sprintf(card,"%-8.8s%-72s",keyword,val); + } else { + sprintf(card,"%-8.8s= '%s' %c%-*s",keyword,val, + (commnt[0] == '\0' ? ' ' : '/'), + (int)(80 - 14 - strlen(val)), commnt); + } +/* + * Write record if full + */ + if(++ncard == 36) { + if(write(fd, record, FITS_SIZE) != FITS_SIZE) { + psError(POIS_ERR_FITS, true,"Cannot write header record\n"); + } + ncard = 0; + } + + return ncard; +} + +/*****************************************************************************/ +/* + * Write an integer (%d) + */ +static int +write_card_i(int fd, + int ncard, + char *record, + const char *keyword, + const int val, + const char *commnt) +{ + char *card = &record[80*ncard]; + + sprintf(card,"%-8.8s= %20d %c%-48s",keyword,val, + (commnt[0] == '\0' ? ' ' : '/'),commnt); +/* + * Write record if full + */ + if(++ncard == 36) { + if(write(fd, record, FITS_SIZE) != FITS_SIZE) { + psError(POIS_ERR_FITS, true,"Cannot write header record\n"); + } + ncard = 0; + } + + return ncard; +} + +/*****************************************************************************/ +/* + * Write a logical value + */ +static int +write_card_l(int fd, + int ncard, + char *record, + const char *keyword, + const char *val, + const char *commnt) +{ + char *card = &record[80*ncard]; + + if(strcmp(val,"T") != 0 && strcmp(val,"F") != 0) { + psError(POIS_ERR_FITS, true,"Invalid logical %s for keyword %s\n",val,keyword); + val = "?"; + } + sprintf(card,"%-8.8s= %20s %c%-48s",keyword,val, + (commnt[0] == '\0' ? ' ' : '/'),commnt); +/* + * Write record if full + */ + if(++ncard == 36) { + if(write(fd, record, FITS_SIZE) != FITS_SIZE) { + psError(POIS_ERR_FITS, true,"Cannot write header record\n"); + } + ncard = 0; + } + + return ncard; +} + +/*****************************************************************************/ +/* + * Write a double floating value + */ +static int +write_card_d(int fd, + int ncard, + char *record, + const char *keyword, + const double val, + const char *commnt) +{ + char *card = &record[80*ncard]; + + sprintf(card,"%-8.8s= %20g %c%-48s",keyword,val, + (commnt[0] == '\0' ? ' ' : '/'),commnt); +/* + * Write record if full + */ + if(++ncard == 36) { + if(write(fd, record, FITS_SIZE) != FITS_SIZE) { + psError(POIS_ERR_FITS, true,"Cannot write header record\n"); + } + ncard = 0; + } + + return ncard; +} + +/*****************************************************************************/ +/* + * Write a floating value + */ +static int +write_card_f(int fd, + int ncard, + char *record, + const char *keyword, + const float val, + const char *commnt) +{ + char *card = &record[80*ncard]; + + sprintf(card,"%-8.8s= %20g %c%-48s",keyword,val, + (commnt[0] == '\0' ? ' ' : '/'),commnt); +/* + * Write record if full + */ + if(++ncard == 36) { + if(write(fd, record, FITS_SIZE) != FITS_SIZE) { + psError(POIS_ERR_FITS, true,"Cannot write header record\n"); + } + ncard = 0; + } + + return ncard; +} + +/*****************************************************************************/ +/* + * Write a CARD + */ +static int +write_card(int fd, + int ncard, + char *record, + const CARD *card) +{ + switch (card->type) { + case UNKNOWN_CARD: // we should already have complained + return ncard; + case LOGICAL_CARD: + return write_card_l(fd, ncard, record, card->keyword, card->val.l, card->commnt); + case STRING_CARD: + return write_card_s(fd, ncard, record, card->keyword, card->val.s, card->commnt); + case INTEGER_CARD: + return write_card_i(fd, ncard, record, card->keyword, card->val.i, card->commnt); + case DOUBLE_CARD: + return write_card_d(fd, ncard, record, card->keyword, card->val.d, card->commnt); + case FLOAT_CARD: + return write_card_f(fd, ncard, record, card->keyword, card->val.f, card->commnt); + } + printf("Impossible type for card %s, %d\n", card->keyword, card->type); + abort(); + return ncard; /* NOTREACHED */ +} + +/*****************************************************************************/ +/* + * Byte swap ABABAB -> BABABAB in place + */ +static void swap_2(char *arr, // array to swap + int n) // number of bytes +{ + char *end, + t; + + if(n%2 != 0) { + psError(POIS_ERR_FITS, true,"Attempt to byte swap odd number of bytes: %d\n",n); + n = 2*(int)(n/2); + } + + for(end = arr + n;arr < end;arr += 2) { + t = arr[0]; + arr[0] = arr[1]; + arr[1] = t; + } +} + +/* + * Byte swap ABCDABCD -> DCBADCBA in place (e.g. sun <--> vax) + */ +static void swap_4(char *arr, // array to swap + int n) // number of bytes +{ + char *end, + t; + + if(n%4 != 0) { + psError(POIS_ERR_FITS, true,"Attempt to byte swap non-multiple of 4 bytes: %d\n",n); + n = 4*(int)(n/4); + } + + for(end = arr + n;arr < end;arr += 4) { + t = arr[0]; + arr[0] = arr[3]; + arr[3] = t; + t = arr[1]; + arr[1] = arr[2]; + arr[2] = t; + } +} + +/* + * Byte swap ABCDEFGH -> HGFEDCBA in place (e.g. sun <--> vax) + */ +static void swap_8(char *arr, // array to swap + int n) // number of bytes +{ + char *end, + t; + + if(n%8 != 0) { + psError(POIS_ERR_FITS, true,"Attempt to byte swap non-multiple of 8 bytes: %d\n",n); + n = 8*(int)(n/8); + } + + for(end = arr + n;arr < end;arr += 8) { + t = arr[0]; + arr[0] = arr[7]; + arr[7] = t; + t = arr[1]; + arr[1] = arr[6]; + arr[6] = t; + t = arr[2]; + arr[2] = arr[5]; + arr[5] = t; + t = arr[3]; + arr[3] = arr[4]; + arr[4] = t; + } +} + +/*****************************************************************************/ + +static int +write_fits_hdr(int fd, + int bitpix, + int naxis, + int *naxes, + CARD **cards, /* extra header cards */ + const psMetadata *metaData, /* metadata to write; or NULL */ + int primary) /* is this the primary HDU? */ +{ + int i; + int ncard; /* number of cards written */ + char record[FITS_SIZE + 1]; /* write buffer */ + + ncard = 0; + if(primary) { + ncard = write_card_l(fd,ncard,record,"SIMPLE","T",""); + } else { + ncard = write_card_s(fd,ncard,record,"XTENSION","IMAGE",""); + } + ncard = write_card_i(fd,ncard,record,"BITPIX",bitpix,""); + ncard = write_card_i(fd,ncard,record,"NAXIS",naxis,""); + for(i = 0; i < naxis; i++) { + char key[] = "NAXIS."; + sprintf(key, "NAXIS%d", i + 1); + ncard = write_card_i(fd,ncard,record,key,naxes[i],""); + } + if(primary) { + ncard = write_card_l(fd,ncard,record, + "EXTEND","T","There may be extensions"); + } +/* + * Write extra header cards + */ + if(cards != NULL) { + while(*cards != NULL) { + ncard = write_card(fd,ncard,record,*cards); + cards++; + } + } + /* + * Is there metadata to write? + */ + if (metaData != NULL) { + CARD card; // used to write metadata to + card.keyword[8] = '\0'; // we\'ll be using strncpy, so need to guarantee NUL termination + psMetadataItem *meta = NULL; // used to traverse a set of metadata + int new_error = 1; // is this a new error? + + + psListIterator *iter = psListIteratorAlloc(metaData->list, PS_LIST_HEAD, false); + while ((meta = psListGetAndIncrement(iter)) != NULL) { + strncpy(card.keyword, meta->name, 8); + card.commnt = meta->comment; + switch (meta->type) { + case PS_META_F64: + card.type = DOUBLE_CARD; + card.val.d = meta->data.F64; + break; + case PS_META_F32: + card.type = FLOAT_CARD; + card.val.f = meta->data.F32; + break; + case PS_META_S32: + card.type = INTEGER_CARD; + card.val.i = meta->data.S32; + break; + case PS_META_STR: + card.type = STRING_CARD; + card.val.s = meta->data.V; + break; + default: + card.type = UNKNOWN_CARD; + psError(POIS_ERR_FITS, new_error, + "I don't know how to write metadata of type %d", meta->type); + new_error = 0; + break; + } + if (card.type != UNKNOWN_CARD) { + ncard = write_card(fd, ncard, record, &card); + } + } + psFree(iter); + } + + ncard = write_card_s(fd,ncard,record,"END","",""); + while(ncard != 0) { + ncard = write_card_s(fd,ncard,record,"","",""); + } + + return 0; +} + +/* + * Pad out to a FITS record boundary + */ +static void pad_to_fits_record(int fd, // output file descriptor + int npixel, // number of pixels already written to HDU + int bitpix) // bitpix for this datatype +{ + const int bytes_per_pixel = (bitpix > 0 ? bitpix : -bitpix)/8; + int nbyte = npixel*bytes_per_pixel; + + if(nbyte%FITS_SIZE != 0) { + char record[FITS_SIZE + 1]; /* write buffer */ + + nbyte = FITS_SIZE - nbyte%FITS_SIZE; + memset(record, ' ', nbyte); + if(write(fd, record, nbyte) != nbyte) { + psError(POIS_ERR_FITS, true, "error padding file to multiple of fits block size\n"); + } + } +} + +static int write_fits_data(int fd, + int bitpix, + int npix, + char *data, + int pad_record) // Pad data out to a FITS record boundary? +{ + char *buff = NULL; /* I/O buffer */ + const int bytes_per_pixel = (bitpix > 0 ? bitpix : -bitpix)/8; + int nbyte = npix*bytes_per_pixel; + int swap_bytes = 0; /* the default */ + static int warned = 0; /* Did we warn about BZERO/BSCALE? */ + +#if PS_BIGENDIAN == 0 /* we'll need to byte swap FITS */ + if(bytes_per_pixel > 1) { + swap_bytes = 1; + } +#endif + + if(swap_bytes) { + if((buff = psAlloc(FITS_SIZE*bytes_per_pixel)) == NULL) { + psError(POIS_ERR_FITS, true,"Can't allocate storage for buff\n"); + return -1; + } + } + + if(bytes_per_pixel == 2 && !warned) { + warned = 1; + fprintf(stderr,"Worry about BZERO/BSCALE\n"); + } + + while(nbyte > 0) { + int nwrite = (nbyte > FITS_SIZE) ? FITS_SIZE : nbyte; + + if(swap_bytes) { + memcpy(buff, data, nwrite); + if(bytes_per_pixel == 2) { + swap_2((char *)buff, nwrite); + } else if(bytes_per_pixel == 4) { + swap_4((char *)buff, nwrite); + } else if(bytes_per_pixel == 8) { + swap_8((char *)buff, nwrite); + } else { + fprintf(stderr,"You cannot get here\n"); + abort(); + } + } else { + buff = data; + } + + if(write(fd, buff, nwrite) != nwrite) { + perror("Error writing image: "); + break; + } + + nbyte -= nwrite; + data += nwrite; + } + if (nbyte == 0) { /* no write errors */ + if (pad_record) { + pad_to_fits_record(fd, npix, bitpix); + } + } + + if(swap_bytes) { + psFree(buff); + } + + return (nbyte == 0 ? 0 : -1); +} + +int +rhlWriteFits(int fd, /* file descriptor to write to */ + const psImage *data, /* The data to write */ + const psMetadata *meta, /* associated metadata */ + const char *WCS) /* name of WCS for pixel coordinates */ +{ + CARD **cards = NULL; /* extra header info */ + int CRVAL1 = 0, CRVAL2 = 0; /* card indices for keywords */ + int i; + const int naxis = 2; // == NAXIS + int naxes[2]; /* values of NAXIS1 etc */ + int ncard = 12; /* number of extra header cards */ + + psErrorClear(); // clear error stack + /* + * Does the data have a WCS? If so, call our pixel-based WCS "A" + */ + if (WCS != NULL && *WCS == '\0') { + WCS = NULL; + } + if (WCS == NULL && + meta != NULL && psMetadataLookup((psMetadata *)meta, "CTYPE1") != NULL) { + WCS = "A"; + } else { + psErrorClear(); // set by psMetadataLookup XXX + WCS = ""; + } + /* + * What sort if image is it? + */ + int bitpix; // == BITPIX + char **rows = NULL; + switch (data->type.type) { + case PS_TYPE_F64: + bitpix = -64; + rows = (char **)data->data.F64; + break; + case PS_TYPE_F32: + bitpix = -32; + rows = (char **)data->data.F32; + break; + case PS_TYPE_U16: + bitpix = -16; + rows = (char **)data->data.U16; + break; + case PS_TYPE_S8: + bitpix = 8; + rows = (char **)data->data.S8; + break; + case PS_TYPE_U8: + bitpix = 8; + rows = (char **)data->data.U8; + break; + default: + bitpix = 0; // make gcc happy + psAbort(__func__, "Unsupported image data type: %d", data->type.type); + break; // NOTREACHED + } + assert (rows != NULL); + + /* + * Allocate cards for FITS headers + */ + if (WCS != NULL) { + cards = psAlloc((ncard + 1)*sizeof(CARD *)); + assert(cards != NULL); + + cards[0] = psAlloc((ncard + 1)*sizeof(CARD)); + assert(cards[0] != NULL); + + for(i = 0; i < ncard; i++) { + cards[i] = cards[0] + i; + } + cards[i] = NULL; + /* + * Generate cards for WCS, so that pixel (0,0) is correctly labelled + */ + i = 0; + sprintf(cards[i]->keyword, "CRVAL1%s", WCS); CRVAL1 = i; + cards[i]->type = INTEGER_CARD; + cards[i]->val.i = 0; + cards[i]->commnt = "(output) Column pixel of Reference Pixel"; + + i++; + sprintf(cards[i]->keyword, "CRVAL2%s", WCS); CRVAL2 = i; + cards[i]->type = INTEGER_CARD; + cards[i]->val.i = 0; + cards[i]->commnt = "(output) Row pixel of Reference Pixel"; + + i++; + sprintf(cards[i]->keyword, "CRPIX1%s", WCS); CRVAL1 = i; + cards[i]->type = INTEGER_CARD; + cards[i]->val.i = 1; + cards[i]->commnt = "Column Pixel Coordinate of Reference"; + + i++; + sprintf(cards[i]->keyword, "CRPIX2%s", WCS); CRVAL1 = i; + cards[i]->type = INTEGER_CARD; + cards[i]->val.i = 1; + cards[i]->commnt = "Row Pixel Coordinate of Reference"; + + i++; + sprintf(cards[i]->keyword, "CTYPE1%s", WCS); CRVAL1 = i; + cards[i]->type = STRING_CARD; + cards[i]->val.s = "LINEAR"; + cards[i]->commnt = "Type of projection"; + + i++; + sprintf(cards[i]->keyword, "CTYPE1%s", WCS); CRVAL1 = i; + cards[i]->type = STRING_CARD; + cards[i]->val.s = "LINEAR"; + cards[i]->commnt = "Type of projection"; + + i++; + sprintf(cards[i]->keyword, "CUNIT1%s", WCS); CRVAL1 = i; + cards[i]->type = STRING_CARD; + cards[i]->val.s = "PIXEL"; + cards[i]->commnt = "Column unit"; + + i++; + sprintf(cards[i]->keyword, "CUNIT2%s", WCS); CRVAL1 = i; + cards[i]->type = STRING_CARD; + cards[i]->val.s = "PIXEL"; + cards[i]->commnt = "Row unit"; + + cards[++i] = NULL; + assert(i <= ncard); + } + + naxes[0] = data->numCols; + naxes[1] = data->numRows; + + write_fits_hdr(fd, bitpix, naxis, naxes, cards, meta, 1); + for (int r = 0; r < data->numRows; r++) { + if(write_fits_data(fd, bitpix, data->numCols, rows[r], ((r==data->numRows - 1) ? 1 : 0)) < 0) { + psError(POIS_ERR_UNKNOWN, 1, "Error writing data for row %d", r); + break; + } + } + + if (cards != NULL) { + psFree(cards[0]); + psFree(cards); + } + + psErr *err = psErrorLast(); + int code = err->code == PS_ERR_NONE ? 0 : -1; + psFree(err); + + return code; +} + +/************************************************************************************************************/ + +int +rhlWriteFitsFile(const char *filename, /* file to write to (or "| cmd"); NULL => "|xpaset ds9 fits"*/ + const psImage *data, /* The data to write */ + const psMetadata *meta, /* associated metadata */ + const char *WCS) /* name of WCS for pixel coordinates */ +{ + if (filename == NULL) { + filename = "| xpaset ds9 fits"; // i.e. display on DS9 + } + + int fd; + + if (filename[0] == '|') { // a command + const char *cmd = filename + 1; + while (isspace(*cmd)) { + cmd++; + } + + fd = fileno(popen(cmd, "w")); + } else { + fd = creat(filename, 777); + } + + if (fd < 0) { + psError(POIS_ERR_FITS, 1, "Cannot open \"%s\"", filename); + return -1; + } + + int ret = rhlWriteFits(fd, data, meta, WCS); + + (void)close(fd); + + return ret; +} Index: /tags/rel-0-1-0/pois/swig/simpleFits.h =================================================================== --- /tags/rel-0-1-0/pois/swig/simpleFits.h (revision 11633) +++ /tags/rel-0-1-0/pois/swig/simpleFits.h (revision 11633) @@ -0,0 +1,15 @@ +#if !defined(SIMPLE_FITS_H) +#define SIMPLE_FITS_H 1 + +int +rhlWriteFits(int fd, /* file descriptor to write to */ + const psImage *data, /* The data to write */ + const psMetadata *meta, /* associated metadata */ + const char *WCS); /* which WCS to use for pixel coordinates */ +int +rhlWriteFitsFile(const char *filename, /* file to write to (or "| cmd"); NULL => "|xpaset ds9 fits" */ + const psImage *data, /* The data to write */ + const psMetadata *meta, /* associated metadata */ + const char *WCS); /* name of WCS for pixel coordinates */ + +#endif Index: /tags/rel-0-1-0/pois/swig/xpaSwig.i =================================================================== --- /tags/rel-0-1-0/pois/swig/xpaSwig.i (revision 11633) +++ /tags/rel-0-1-0/pois/swig/xpaSwig.i (revision 11633) @@ -0,0 +1,148 @@ +// -*- C -*- +/* + * This would be a very simple module, except that I didn't want to + * deal with (char **) in SWIG; so I wrote C wrappers to pass a single + * (char *) instead. + * + * I also include some utility routines used to write psImages to + * a file descriptor; these are used to display images on ds9 via xpa + */ +%module xpaSwig + +%include "p_psSwig.i" +%import "psLibSwig.i" + +%rename(XPAGet) XPAGet1; +%rename(XPASet) XPASet1; +%rename(XPASetFd) XPASetFd1; + +%{ +#include "xpa.h" + +static char * +xmalloc(long n) +{ + char *ptr = malloc(n); + assert(ptr != NULL); + + return(ptr); +} + +/* + * A binding for XPAGet that talks to only one server, but doesn't have to talk (char **) with SWIG + */ +char * +XPAGet1(XPA xpa, + char *xtemplate, + char *paramlist, + char *mode) +{ + char *buf = NULL; /* desired response */ + int len = 0; /* length of buf; ignored */ + char *error = NULL; /* returned error if any*/ + + int n = XPAGet(xpa, xtemplate, paramlist, mode, + &buf, &len, NULL, &error, 1); + + if(n == 0) { + return(NULL); + } + if(error != NULL) { + char *errStr = xmalloc(strlen(error) + 1); + strcpy(errStr, error); + return(errStr); + } + + return(buf); +} + +/*****************************************************************************/ + +char * +XPASet1(XPA xpa, + char *xtemplate, + char *paramlist, + char *mode, + char *buf, // desired extra data + int len) // length of buf (or -1) +{ + if(len < 0) { + len = strlen(buf); // length of buf + } + char *error = NULL; // returned error if any + + int n = XPASet(xpa, xtemplate, paramlist, mode, + buf, len, NULL, &error, 1); + + if(n == 0) { + return(NULL); + } + if(error != NULL) { + char *errStr = xmalloc(strlen(error) + 1); + strcpy(errStr, error); + return(errStr); + } + + return ""; +} + + +/*****************************************************************************/ + +char * +XPASetFd1(XPA xpa, + char *xtemplate, + char *paramlist, + char *mode, + int fd) /* file descriptor for xpa to read */ +{ + char *error = NULL; /* returned error if any*/ + + int n = XPASetFd(xpa, xtemplate, paramlist, mode, + fd, NULL, &error, 1); + + if(n == 0) { + return(NULL); + } + if(error != NULL) { + char *errStr = xmalloc(strlen(error) + 1); + strcpy(errStr, error); + + return(errStr); + } + + return NULL; +} +%} + +%rename(XPA_in) in; // avoid conflict with python keyword in xpa.h + +%import "prsetup.h" +%import "xpa.h" + +%include "exception.i" + +%exception { + $action + if (result == NULL) { + SWIG_exception(SWIG_IOError, "XPA returned NULL") + } +} + +char *XPAGet1(XPA xpa, char *xtemplate, char *paramlist, char *mode); +char *XPASet1(XPA xpa, char *xtemplate, char *paramlist, char *mode, char *buf, int len); +char *XPASetFd1(XPA xpa, char *xtemplate, char *paramlist, char *mode, int fd); + +/************************************************************************************************************/ +/* + * Talk to ds9 + */ +%{ +# include "pslib.h" +# include "simpleFits.h" +%} + +NOTNEGATIVE(rhlWriteFits); +NOTNEGATIVE(rhlWriteFitsFile); + +%include "simpleFits.h" Index: /tags/rel-0-1-0/ppStac/.cvsignore =================================================================== --- /tags/rel-0-1-0/ppStac/.cvsignore (revision 11633) +++ /tags/rel-0-1-0/ppStac/.cvsignore (revision 11633) @@ -0,0 +1,17 @@ +aclocal.m4 +compile +config.guess +config.log +config.status +config.sub +configure +depcomp +install-sh +libtool +ltmain.sh +Makefile +Makefile.in +missing +ppStac.pc +test +autom4te.cache Index: /tags/rel-0-1-0/ppStac/Makefile.am =================================================================== --- /tags/rel-0-1-0/ppStac/Makefile.am (revision 11633) +++ /tags/rel-0-1-0/ppStac/Makefile.am (revision 11633) @@ -0,0 +1,9 @@ +SUBDIRS = src + +CLEANFILES = *~ core core.* + +pkgconfigdir = $(libdir)/pkgconfig +pkgconfig_DATA= ppStac.pc + +EXTRA_DIST = \ + ppStac.pc.in Index: /tags/rel-0-1-0/ppStac/autogen.sh =================================================================== --- /tags/rel-0-1-0/ppStac/autogen.sh (revision 11633) +++ /tags/rel-0-1-0/ppStac/autogen.sh (revision 11633) @@ -0,0 +1,103 @@ +#!/bin/sh +# Run this to generate all the initial makefiles, etc. + +srcdir=`dirname $0` +test -z "$srcdir" && srcdir=. + +ORIGDIR=`pwd` +cd $srcdir + +PROJECT=ppStac +TEST_TYPE=-f +# change this to be a unique filename in the top level dir +FILE=autogen.sh + +DIE=0 + +LIBTOOLIZE=libtoolize +ACLOCAL="aclocal $ACLOCAL_FLAGS" +AUTOHEADER=autoheader +AUTOMAKE=automake +AUTOCONF=autoconf + +($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $LIBTOOlIZE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/" + DIE=1 +} + +($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $ACLOCAL installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOHEADER installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOMAKE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOCONF installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +if test "$DIE" -eq 1; then + exit 1 +fi + +test $TEST_TYPE $FILE || { + echo "You must run this script in the top-level $PROJECT directory" + exit 1 +} + +if test -z "$*"; then + echo "I am going to run ./configure with no arguments - if you wish " + echo "to pass any to it, please specify them on the $0 command line." +fi + +$LIBTOOLIZE --copy --force || echo "$LIBTOOlIZE failed" +$ACLOCAL || echo "$ACLOCAL failed" +$AUTOHEADER || echo "$AUTOHEADER failed" +$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE failed" +$AUTOCONF || echo "$AUTOCONF failed" + +cd $ORIGDIR + +run_configure=true +for arg in $*; do + case $arg in + --no-configure) + run_configure=false + ;; + *) + ;; + esac +done + +if $run_configure; then + $srcdir/configure --enable-maintainer-mode "$@" + echo + echo "Now type 'make' to compile $PROJECT." +else + echo + echo "Now run 'configure' and 'make' to compile $PROJECT." +fi Index: /tags/rel-0-1-0/ppStac/configure.ac =================================================================== --- /tags/rel-0-1-0/ppStac/configure.ac (revision 11633) +++ /tags/rel-0-1-0/ppStac/configure.ac (revision 11633) @@ -0,0 +1,30 @@ +AC_PREREQ(2.59) + +AC_INIT([ppStac], [0.1.0], [ipp-support@ifa.hawaii.edu]) +AC_CONFIG_SRCDIR([src]) + +AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2]) +AM_CONFIG_HEADER([src/config.h]) +AM_MAINTAINER_MODE + +IPP_STDCFLAGS + +AC_LANG(C) +AC_GNU_SOURCE +AC_PROG_CC +AC_PROG_INSTALL +AC_PROG_LIBTOOL +AC_SYS_LARGEFILE + +PKG_CHECK_MODULES([PSLIB], [pslib >= 1.0.0]) +PKG_CHECK_MODULES([PSMODULE], [psmodules >= 1.0.0]) + +IPP_STDOPTS +CFLAGS="${CFLAGS=} -Wall -Werror -std=c99" + +AC_CONFIG_FILES([ + Makefile + src/Makefile + ppStac.pc +]) +AC_OUTPUT Index: /tags/rel-0-1-0/ppStac/ppStac.pc.in =================================================================== --- /tags/rel-0-1-0/ppStac/ppStac.pc.in (revision 11633) +++ /tags/rel-0-1-0/ppStac/ppStac.pc.in (revision 11633) @@ -0,0 +1,11 @@ +prefix=@prefix@ +exec_prefix=@exec_prefix@ +libdir=@libdir@ +includedir=@includedir@ + +Name: libppStac +Description: Pan-STARRS Pipeline Simultaneous Telescope Array Combination +Version: @VERSION@ +Requires: pslib psmodules +Libs: -L${libdir} +Cflags: -I${includedir} Index: /tags/rel-0-1-0/ppStac/src/.cvsignore =================================================================== --- /tags/rel-0-1-0/ppStac/src/.cvsignore (revision 11633) +++ /tags/rel-0-1-0/ppStac/src/.cvsignore (revision 11633) @@ -0,0 +1,10 @@ +.deps +Makefile +Makefile.in +.libs +*.lo +*.la +ppStac +config.h +config.h.in +stamp-h1 Index: /tags/rel-0-1-0/ppStac/src/Makefile.am =================================================================== --- /tags/rel-0-1-0/ppStac/src/Makefile.am (revision 11633) +++ /tags/rel-0-1-0/ppStac/src/Makefile.am (revision 11633) @@ -0,0 +1,20 @@ +bin_PROGRAMS = ppStac + +ppStac_CFLAGS = $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +ppStac_LDFLAGS = $(PSMODULE_LIBS) $(PSLIB_LIBS) + +ppStac_SOURCES = \ + ppStac.c \ + ppStacOptions.c \ + ppStacCombine.c + +noinst_HEADERS = \ + ppStac.h + +CLEANFILES = *~ + +clean-local: + -rm -f TAGS +# Tags for emacs +tags: + etags `find . -name \*.[ch] -print` Index: /tags/rel-0-1-0/ppStac/src/ppStac.c =================================================================== --- /tags/rel-0-1-0/ppStac/src/ppStac.c (revision 11633) +++ /tags/rel-0-1-0/ppStac/src/ppStac.c (revision 11633) @@ -0,0 +1,43 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#define TIMERNAME "PPSTAC" + +#include "ppStac.h" + +int main(int argc, char *argv[]) +{ + psExit status = PS_EXIT_SUCCESS; + + psLibInit(NULL); + psTimerStart(TIMERNAME); + + // Parse the command-line arguments, open the files + ppStacOptions *options = ppStacOptionsRead(&argc, argv); + if (!options) { + psError(PS_ERR_UNKNOWN, false, "Unable to set program options.\n"); + status = PS_EXIT_CONFIG_ERROR; + goto die; + } + + if (!ppStacCombine(options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to combine images.\n"); + status = PS_EXIT_DATA_ERROR; + goto die; + } + + + // Common code for the death. +die: + psFree(options); + pmConceptsDone(); + pmConfigDone(); + psLibFinalize(); + + return status; +} Index: /tags/rel-0-1-0/ppStac/src/ppStac.h =================================================================== --- /tags/rel-0-1-0/ppStac/src/ppStac.h (revision 11633) +++ /tags/rel-0-1-0/ppStac/src/ppStac.h (revision 11633) @@ -0,0 +1,29 @@ +#ifndef PP_STAC_H +#define PP_STAC_H + +#define PPSTAC_RECIPE "PPSTAC" + +#include + +// Configuration information +typedef struct { + pmConfig *config; // psModules configuration + psArray *inFiles; // Input files + psFits *outFile; // Output file + psStatsOptions combine; // Combination statistic + int iter; // Number of rejection iterations + float rej; // Rejection threshold (standard deviations) + psMaskType maskVal; // Mask value +} ppStacOptions; + + +// Allocator for ppStacOptions +ppStacOptions *ppStacOptionsAlloc(void); + +// Read the options +ppStacOptions *ppStacOptionsRead(int *argc, char *argv[]); + +// Do the combination +bool ppStacCombine(ppStacOptions *options); + +#endif Index: /tags/rel-0-1-0/ppStac/src/ppStacCombine.c =================================================================== --- /tags/rel-0-1-0/ppStac/src/ppStacCombine.c (revision 11633) +++ /tags/rel-0-1-0/ppStac/src/ppStacCombine.c (revision 11633) @@ -0,0 +1,41 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppStac.h" + +bool ppStacCombine(ppStacOptions *options) +{ + PS_ASSERT_PTR_NON_NULL(options, false); + + psArray *inputs = psArrayAlloc(options->inFiles->n); // Input images + psRegion region = psRegionSet(0, 0, 0, 0); // Region to read + for (int i = 0; i < inputs->n; i++) { + psImage *image = psFitsReadImage(options->inFiles->data[i], region, 0); + if (!image) { + psError(PS_ERR_IO, false, "Unable to read image %d\n", i); + return false; + } + inputs->data[i] = image; + } + + psImage *output = pmCombineImages(NULL, NULL, inputs, NULL, NULL, options->maskVal, NULL, + options->iter, options->rej); + + if (!output) { + psError(PS_ERR_UNKNOWN, false, "Unable to combine images.\n"); + return false; + } + + if (!psFitsWriteImage(options->outFile, NULL, output, 0, NULL)) { + psError(PS_ERR_UNKNOWN, false, "Unable to write output image.\n"); + return false; + } + psFree(output); + + return true; +} Index: /tags/rel-0-1-0/ppStac/src/ppStacErrorCodes.c.in =================================================================== --- /tags/rel-0-1-0/ppStac/src/ppStacErrorCodes.c.in (revision 11633) +++ /tags/rel-0-1-0/ppStac/src/ppStacErrorCodes.c.in (revision 11633) @@ -0,0 +1,26 @@ +/* + * The line + { PPSTAC_ERR_$X{ErrorCode}, "$X{ErrorDescription}"}, + * (without the Xs) + * will be replaced by values from errorCodes.dat + */ +#include "pslib.h" +#include "ppStacErrorCodes.h" + +void ppStacErrorRegister(void) +{ + static psErrorDescription errors[] = { + { PPSTAC_ERR_BASE, "First value we use; lower values belong to psLib" }, + { PPSTAC_ERR_${ErrorCode}, "${ErrorDescription}"}, + }; + static int nerror = PPSTAC_ERR_NERROR - PPSTAC_ERR_BASE; // number of values in enum + + for (int i = 0; i < nerror; i++) { + psErrorDescription *tmp = psAlloc(sizeof(psErrorDescription)); + p_psMemSetPersistent(tmp, true); + *tmp = errors[i]; + psErrorRegister(tmp, 1); + psFree(tmp); /* it's on the internal list */ + } + nerror = 0; // don't register more than once +} Index: /tags/rel-0-1-0/ppStac/src/ppStacErrorCodes.dat =================================================================== --- /tags/rel-0-1-0/ppStac/src/ppStacErrorCodes.dat (revision 11633) +++ /tags/rel-0-1-0/ppStac/src/ppStacErrorCodes.dat (revision 11633) @@ -0,0 +1,10 @@ +# +# This file is used to generate ppStacErrorClasses.h +# +BASE = 800 First value we use; lower values belong to psLib +UNKNOWN Unknown PM error code +NOT_IMPLEMENTED Desired feature is not yet implemented +ARGUMENTS Incorrect arguments +CONFIG Problem in configure files +IO Problem in FITS I/O +DATA Problem in data values Index: /tags/rel-0-1-0/ppStac/src/ppStacErrorCodes.h.in =================================================================== --- /tags/rel-0-1-0/ppStac/src/ppStacErrorCodes.h.in (revision 11633) +++ /tags/rel-0-1-0/ppStac/src/ppStacErrorCodes.h.in (revision 11633) @@ -0,0 +1,18 @@ +#if !defined(PPSTAC_ERROR_CODES_H) +#define PPSTAC_ERROR_CODES_H +/* + * The line + * PPSTAC_ERR_$X{ErrorCode}, + * (without the X) + * + * will be replaced by values from errorCodes.dat + */ +typedef enum { + PPSTAC_ERR_BASE = 512, + PPSTAC_ERR_${ErrorCode}, + PPSTAC_ERR_NERROR +} ppStacErrorCode; + +void ppStacErrorRegister(void); + +#endif Index: /tags/rel-0-1-0/ppStac/src/ppStacOptions.c =================================================================== --- /tags/rel-0-1-0/ppStac/src/ppStacOptions.c (revision 11633) +++ /tags/rel-0-1-0/ppStac/src/ppStacOptions.c (revision 11633) @@ -0,0 +1,125 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include +#include + +#include "ppStac.h" + +static void stacOptionsFree(ppStacOptions *options) +{ + psFree(options->config); + psFree(options->inFiles); + if (options->outFile) { + psFitsClose(options->outFile); + } +} + + +ppStacOptions *ppStacOptionsAlloc(void) +{ + ppStacOptions *options = psAlloc(sizeof(ppStacOptions)); + psMemSetDeallocator(options, (psFreeFunc)stacOptionsFree); + + options->config = NULL; + options->inFiles = NULL; + options->outFile = NULL; + options->combine = 0; + options->iter = 1; + options->rej = 3.0; + options->maskVal = 0; + + return options; +} + + + +// Parse a recipe option according to its type +#define OPTION_PARSE(OPTION,MD,NAME,TYPE) \ +{ \ + psMetadataItem *item = psMetadataLookup(MD, NAME); \ + if (item) { \ + OPTION = psMetadataItemParse##TYPE(item); \ + } else { \ + psWarning("Recipe option %s isn't specified; using default.\n", NAME); \ + } \ +} + +ppStacOptions *ppStacOptionsRead(int *argc, char *argv[]) +{ + PS_ASSERT_INT_POSITIVE(*argc, NULL); + PS_ASSERT_PTR_NON_NULL(argv, NULL); + + ppStacOptions *options = ppStacOptionsAlloc(); // Options, to be returned + + // Parse the configuration and arguments + pmConfigReadParamsSet(false); + options->config = pmConfigRead(argc, argv, PPSTAC_RECIPE); + if (!options->config) { + psError(PS_ERR_UNKNOWN, false, "Unable to read configuration.\n"); + psFree(options); + return NULL; + } + + if (*argc <= 1 || psArgumentGet(*argc, argv, "-help") || psArgumentGet(*argc, argv, "-h")) { + fprintf(stderr, + "ppStac: Combine multiple registered images.\n\n" + "Usage: %s OUT.fits IN1.fits IN2.fits ... [-recipe NAME FILE]\n\n", + argv[0]); + psFree(options); + exit(PS_EXIT_SUCCESS); + } + + bool mdok; // Status of MD lookup + psMetadata *recipe = psMetadataLookupMetadata(&mdok, options->config->recipes, + PPSTAC_RECIPE); // Recipe to use + if (!recipe) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find recipe %s.\n", PPSTAC_RECIPE); + psFree(options); + return NULL; + } + + const char *outName = argv[1]; // Output filename + options->outFile = psFitsOpen(outName, "w"); + if (!options->outFile) { + psError(PS_ERR_IO, false, "Unable to open output file %s.\n", outName); + psFree(options); + return NULL; + } + + int numInputs = *options->config->argc - 2; // Number of input files + if (numInputs <= 1) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Not enough input images supplied to combine.\n"); + return NULL; + } + + options->inFiles = psArrayAlloc(numInputs); + for (int i = 0; i < numInputs; i++) { + const char *inName = argv[i + 2]; + psFits *inFile = psFitsOpen(inName, "r"); + if (!inFile) { + psError(PS_ERR_IO, false, "Unable to open input file %s.\n", inName); + psFree(options); + return NULL; + } + options->inFiles->data[i] = inFile; + } + + OPTION_PARSE(options->rej, recipe, "REJ", F32); + OPTION_PARSE(options->iter, recipe, "ITER", S32); + OPTION_PARSE(options->maskVal, recipe, "MASKVAL", U8); + + const char *stat = psMetadataLookupStr(&mdok, recipe, "COMBINE"); // The combination statistic + if (!mdok || !stat || strlen(stat) == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Can't recognise statistic COMBINE.\n"); + psFree(options); + return NULL; + } + options->combine = psStatsOptionFromString(stat); + + return options; +} Index: /tags/rel-0-1-0/pswarp/.cvsignore =================================================================== --- /tags/rel-0-1-0/pswarp/.cvsignore (revision 11633) +++ /tags/rel-0-1-0/pswarp/.cvsignore (revision 11633) @@ -0,0 +1,18 @@ +bin +Makefile +Makefile.in +aclocal.m4 +autom4te.cache +compile +config.log +config.status +configure +depcomp +install-sh +missing +config.guess +libtool +ltmain.sh +stamp-h1 +config.sub +pswarp.pc Index: /tags/rel-0-1-0/pswarp/Makefile.am =================================================================== --- /tags/rel-0-1-0/pswarp/Makefile.am (revision 11633) +++ /tags/rel-0-1-0/pswarp/Makefile.am (revision 11633) @@ -0,0 +1,10 @@ +SUBDIRS = src + +CLEANFILES = *.pyc *~ core core.* + +pkgconfigdir = $(libdir)/pkgconfig +pkgconfig_DATA= pswarp.pc + +EXTRA_DIST = \ + pswarp.pc.in \ + autogen.sh Index: /tags/rel-0-1-0/pswarp/autogen.sh =================================================================== --- /tags/rel-0-1-0/pswarp/autogen.sh (revision 11633) +++ /tags/rel-0-1-0/pswarp/autogen.sh (revision 11633) @@ -0,0 +1,103 @@ +#!/bin/sh +# Run this to generate all the initial makefiles, etc. + +srcdir=`dirname $0` +test -z "$srcdir" && srcdir=. + +ORIGDIR=`pwd` +cd $srcdir + +PROJECT=pswarp +TEST_TYPE=-f +# change this to be a unique filename in the top level dir +FILE=autogen.sh + +DIE=0 + +LIBTOOLIZE=libtoolize +ACLOCAL="aclocal $ACLOCAL_FLAGS" +AUTOHEADER=autoheader +AUTOMAKE=automake +AUTOCONF=autoconf + +($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $LIBTOOlIZE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/" + DIE=1 +} + +($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $ACLOCAL installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOHEADER installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOMAKE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOCONF installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +if test "$DIE" -eq 1; then + exit 1 +fi + +test $TEST_TYPE $FILE || { + echo "You must run this script in the top-level $PROJECT directory" + exit 1 +} + +if test -z "$*"; then + echo "I am going to run ./configure with no arguments - if you wish " + echo "to pass any to it, please specify them on the $0 command line." +fi + +$LIBTOOLIZE --copy --force || echo "$LIBTOOlIZE failed" +$ACLOCAL || echo "$ACLOCAL failed" +$AUTOHEADER || echo "$AUTOHEADER failed" +$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE failed" +$AUTOCONF || echo "$AUTOCONF failed" + +cd $ORIGDIR + +run_configure=true +for arg in $*; do + case $arg in + --no-configure) + run_configure=false + ;; + *) + ;; + esac +done + +if $run_configure; then + $srcdir/configure --enable-maintainer-mode "$@" + echo + echo "Now type 'make' to compile $PROJECT." +else + echo + echo "Now run 'configure' and 'make' to compile $PROJECT." +fi Index: /tags/rel-0-1-0/pswarp/configure.ac =================================================================== --- /tags/rel-0-1-0/pswarp/configure.ac (revision 11633) +++ /tags/rel-0-1-0/pswarp/configure.ac (revision 11633) @@ -0,0 +1,203 @@ +dnl Process this file with autoconf to produce a configure script. +AC_PREREQ(2.59) + +AC_INIT([pswarp], [0.1.0], [ipp-support@ifa.hawaii.edu]) +AC_CONFIG_SRCDIR([src]) + +AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2]) +AM_CONFIG_HEADER([src/config.h]) +AM_MAINTAINER_MODE + +IPP_STDCFLAGS + +AC_LANG(C) +AC_GNU_SOURCE +AC_PROG_CC +AC_PROG_INSTALL +AC_PROG_LIBTOOL +AC_SYS_LARGEFILE + +dnl ------------------------------------------------------------ + +AC_PATH_PROG([ERRORCODES], [psParseErrorCodes], [missing]) +if test "$ERRORCODES" = "missing" ; then + AC_MSG_ERROR([psParseErrorCodes is required]) +fi + +dnl ------------------ kapa,libkapa options ------------------------- +dnl -- libkapa implies the requirement for libpng, libjpeg as well -- + +dnl save LIBS/CFLAGS/LDFLAGS +TMP_LIBS=${LIBS} +TMP_CFLAGS=${CFLAGS} +TMP_LDFLAGS=${LDFLAGS} +TMP_CPPFLAGS=${CPPFLAGS} + +dnl test for command-line options: use ohana-config if not supplied +KAPA_CFLAGS_CONFIG="true" +KAPA_LIBS_CONFIG="true" +AC_ARG_WITH(kapa, +[AS_HELP_STRING(--with-kapa=DIR,Specify location of libkapa)], +[KAPA_CFLAGS="-I$withval/include" KAPA_LIBS="-L$withval/lib" + KAPA_CFLAGS_CONFIG="false" KAPA_LIBS_CONFIG="false"]) +AC_ARG_WITH(kapa-include, +[AS_HELP_STRING(--with-kapa-include=DIR,Specify libkapa include directory.)], +[KAPA_CFLAGS="-I$withval" KAPA_CFLAGS_CONFIG="false"]) +AC_ARG_WITH(kapa-lib, +[AS_HELP_STRING(--with-kapa-lib=DIR,Specify libkapa library directory.)], +[KAPA_LIBS="-L$withval" KAPA_LIBS_CONFIG="false"]) + +echo "KAPA_CFLAGS_CONFIG: $KAPA_CFLAGS_CONFIG" +echo "KAPA_LIBS_CONFIG: $KAPA_LIBS_CONFIG" +echo "KAPA_CFLAGS: $KAPA_CFLAGS" +echo "KAPA_LIBS: $KAPA_LIBS" + +dnl HAVE_KAPA is set to false if any of the tests fail +HAVE_KAPA="true" +AC_MSG_NOTICE([checking for libkapa]) +if test "$KAPA_CFLAGS_CONFIG" = "true" -o "$KAPA_LIBS_CONFIG" = "true"; then + AC_MSG_NOTICE([kapa info supplied by ohana-config]) + KAPA_CONFIG=`which ohana-config` + AC_CHECK_FILE($KAPA_CONFIG,[], + [HAVE_KAPA="false"; AC_MSG_WARN([libkapa is not found: output plots disabled. Obtain libkapa at http://kiawe.ifa.hawaii.edu/Elixir/Ohana or use --with-kapa to specify location])]) + + echo "HAVE_KAPA: $HAVE_KAPA" + echo "KAPA_CFLAGS_CONFIG: $KAPA_CFLAGS_CONFIG" + + if test "$HAVE_KAPA" = "true" -a "$KAPA_CFLAGS_CONFIG" = "true" ; then + AC_MSG_NOTICE([libkapa cflags info supplied by ohana-config]) + AC_MSG_CHECKING([libkapa cflags]) + KAPA_CFLAGS="`${KAPA_CONFIG} --cflags`" + AC_MSG_RESULT([${KAPA_CFLAGS}]) + fi + + if test "$HAVE_KAPA" = "true" -a "$KAPA_LIBS_CONFIG" = "true" ; then + AC_MSG_NOTICE([libkapa ldflags info supplied by ohana-config]) + AC_MSG_CHECKING([libkapa ldflags]) + KAPA_LIBS="`${KAPA_CONFIG} --libs` -lX11" + AC_MSG_RESULT([${KAPA_LIBS}]) + fi +fi + +if test "$HAVE_KAPA" = "true" ; then + AC_MSG_NOTICE([libkapa supplied]) + PSWARP_CFLAGS="${PSWARP_CFLAGS} ${KAPA_CFLAGS}" + PSWARP_LIBS="${PSWARP_LIBS} ${KAPA_LIBS}" +else + AC_MSG_NOTICE([libkapa ignored]) +fi + +dnl restore the CFLAGS/LDFLAGS +LIBS=${TMP_LIBS} +CFLAGS=${TMP_CFLAGS} +LDFLAGS=${TMP_LDFLAGS} +CPPFLAGS=${TMP_CPPFLAGS} + +dnl ------------------ libjpeg options --------------------- + +dnl save LIBS/CFLAGS/LDFLAGS +TMP_LIBS=${LIBS} +TMP_CFLAGS=${CFLAGS} +TMP_LDFLAGS=${LDFLAGS} +TMP_CPPFLAGS=${CPPFLAGS} + +AC_ARG_WITH(jpeg, +[AS_HELP_STRING(--with-jpeg=DIR,Specify location of libjpeg.)], +[JPEG_CFLAGS="-I$withval/include" + JPEG_LDFLAGS="-L$withval/lib"]) +AC_ARG_WITH(jpeg-include, +[AS_HELP_STRING(--with-jpeg-include=DIR,Specify libjpeg include directory.)], +[JPEG_CFLAGS="-I$withval"]) +AC_ARG_WITH(jpeg-lib, +[AS_HELP_STRING(--with-jpeg-lib=DIR,Specify libjpeg library directory.)], +[JPEG_LDFLAGS="-L$withval"]) + +CFLAGS="${CFLAGS} ${JPEG_CFLAGS}" +CPPFLAGS=${CFLAGS} +LDFLAGS="${LDFLAGS} ${JPEG_LDFLAGS}" + +AC_CHECK_HEADERS([jpeglib.h], + [PSWARP_CFLAGS="$PSWARP_CFLAGS $JPEG_CFLAGS" AC_SUBST(JPEG_CFLAGS)], + [HAVE_KAPA=false; AC_MSG_WARN([libjpeg headers not found: output plots disabled. Obtain libjpeg from http://www.ijg.org/ or use --with-jpeg to specify location.])] +) + +AC_CHECK_LIB(jpeg,jpeg_CreateCompress, + [PSWARP_LIBS="$PSWARP_LIBS $JPEG_LDFLAGS -ljpeg"], + [HAVE_KAPA=false; AC_MSG_WARN([libjpeg library not found: output plots disabled. Obtain libjpeg from http://www.ijg.org/ or use --with-jpeg to specify location.])] +) + +dnl restore the CFLAGS/LDFLAGS +LIBS=${TMP_LIBS} +CFLAGS=${TMP_CFLAGS} +LDFLAGS=${TMP_LDFLAGS} +CPPFLAGS=${TMP_CPPFLAGS} + +dnl ------------------ libpng options --------------------- + +dnl save LIBS/CFLAGS/LDFLAGS +TMP_LIBS=${LIBS} +TMP_CFLAGS=${CFLAGS} +TMP_LDFLAGS=${LDFLAGS} +TMP_CPPFLAGS=${CPPFLAGS} + +AC_ARG_WITH(png, +[AS_HELP_STRING(--with-png=DIR,Specify location of libpng.)], +[PNG_CFLAGS="-I$withval/include" + PNG_LDFLAGS="-L$withval/lib"]) +AC_ARG_WITH(png-include, +[AS_HELP_STRING(--with-png-include=DIR,Specify libpng include directory.)], +[PNG_CFLAGS="-I$withval"]) +AC_ARG_WITH(png-lib, +[AS_HELP_STRING(--with-png-lib=DIR,Specify libpng library directory.)], +[PNG_LDFLAGS="-L$withval"]) + +CFLAGS="${CFLAGS} ${PNG_CFLAGS}" +CPPFLAGS=${CFLAGS} +LDFLAGS="${LDFLAGS} ${PNG_LDFLAGS}" + +AC_CHECK_HEADERS([png.h], + [PSWARP_CFLAGS="$PSWARP_CFLAGS $PNG_CFLAGS" AC_SUBST(PNG_CFLAGS)], + [HAVE_KAPA=false; AC_MSG_WARN([libpng headers not found: output plots disabled. Obtain libpng from http://www.ijg.org/ or use --with-png to specify location.])] +) + +AC_CHECK_LIB(png,png_init_io, + [PSWARP_LIBS="$PSWARP_LIBS $PNG_LDFLAGS -lpng"], + [HAVE_KAPA=false; AC_MSG_WARN([libpng library not found: output plots disabled. Obtain libpng from http://www.ijg.org/ or use --with-png to specify location.])] +) + +dnl restore the CFLAGS/LDFLAGS +LIBS=${TMP_LIBS} +CFLAGS=${TMP_CFLAGS} +LDFLAGS=${TMP_LDFLAGS} +CPPFLAGS=${TMP_CPPFLAGS} + +dnl ------------------ use kapa or not? --------------------- + +if test "$HAVE_KAPA" == "true" ; then + AC_MSG_RESULT([including plotting functions]) + AC_DEFINE([HAVE_KAPA],[1],[enable use of libkapa]) +else + AC_MSG_RESULT([skipping plotting functions]) + AC_DEFINE([HAVE_KAPA],[0],[disable use of libkapa]) +fi + +dnl ------------- psLib, psModules --------------- +PKG_CHECK_MODULES(PSLIB, pslib >= 1.0.0) +PKG_CHECK_MODULES(PSMODULE, psmodules >= 1.0.0) + +dnl Set CFLAGS for build +IPP_STDOPTS +CFLAGS="${CFLAGS} -Wall -Werror -std=c99" +echo "PSWARP_CFLAGS: $PSWARP_CFLAGS" +echo "PSWARP_LIBS: $PSWARP_LIBS" + +AC_SUBST([PSWARP_CFLAGS]) +AC_SUBST([PSWARP_LIBS]) + +AC_CONFIG_FILES([ + Makefile + src/Makefile + pswarp.pc +]) + +AC_OUTPUT Index: /tags/rel-0-1-0/pswarp/doc/notes.txt =================================================================== --- /tags/rel-0-1-0/pswarp/doc/notes.txt (revision 11633) +++ /tags/rel-0-1-0/pswarp/doc/notes.txt (revision 11633) @@ -0,0 +1,29 @@ + +- pswarp takes an image (and corresponding mask and weight images) and + warps the pixels to match the specified skycell projection and pixel + grid. the astrometry of the input image is supplied by the header + of the image or by a specified SMP file. + +- the input image, weight, mask are loaded into pmFPA structures with + the layout determined by the camera/format/layout + +- the skycell is loaded into a pmFPA structure. at a minimum, a + single skycell may contain a single chip, cell, and readout, which + in turn correspond to the single set of pixels representing the + image region. + + * is there a need or use for a skycell to have a multilayered + structure equivalent in any way to the fpa->chip->cell->readout + layer? + + - perhaps we will want to use a single fpa for a larger region + with small sky cells representing contiguous pixel regions. + + - perhaps we will want to warp an entire fpa exposure at once and + use a skyfpa to represent all of the possible output skycells + for a single input fpa exposure. + +- the output is a file or set of files corresponding to the skycell + files. the details should be specified by the PSWARP.OUTPUT file + rule. + Index: /tags/rel-0-1-0/pswarp/pswarp.pc.in =================================================================== --- /tags/rel-0-1-0/pswarp/pswarp.pc.in (revision 11633) +++ /tags/rel-0-1-0/pswarp/pswarp.pc.in (revision 11633) @@ -0,0 +1,12 @@ +prefix=@prefix@ +exec_prefix=@exec_prefix@ +libdir=@libdir@ +includedir=@includedir@ + +Name: libpswarp +Description: Pan-STARRS Photometry Library +Version: @VERSION@ +Requires: pslib psmodules +Libs: -L${libdir} -lpswarp +Libs.private: @PSWARP_LIBS@ +Cflags: -I${includedir} @PSWARP_CFLAGS@ Index: /tags/rel-0-1-0/pswarp/src/.cvsignore =================================================================== --- /tags/rel-0-1-0/pswarp/src/.cvsignore (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/.cvsignore (revision 11633) @@ -0,0 +1,14 @@ +*.o +*.lo +.libs +.deps +Makefile +Makefile.in +pswarp +pswarp.loT +pswarpErrorCodes.h +pswarpErrorCodes.c +libpswarp.la +config.h +config.h.in +stamp-h1 Index: /tags/rel-0-1-0/pswarp/src/Makefile.am =================================================================== --- /tags/rel-0-1-0/pswarp/src/Makefile.am (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/Makefile.am (revision 11633) @@ -0,0 +1,47 @@ + +lib_LTLIBRARIES = libpswarp.la +libpswarp_la_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSWARP_CFLAGS) + +bin_PROGRAMS = pswarp +pswarp_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSWARP_CFLAGS) +pswarp_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSWARP_LIBS) +pswarp_LDADD = libpswarp.la + +pswarp_SOURCES = \ + pswarp.c + +libpswarp_la_SOURCES = \ +pswarpArguments.c \ +pswarpCleanup.c \ +pswarpDataLoad.c \ +pswarpDataSave.c \ +pswarpDefine.c \ +pswarpErrorCodes.c \ +pswarpMapGrid.c \ +pswarpMatchRange.c \ +pswarpParseCamera.c \ +pswarpTransformReadout.c \ +pswarpTransformReadout_Opt.c \ +pswarpVersion.c + +include_HEADERS = \ + pswarp.h \ + pswarpErrorCodes.h + +clean-local: + -rm -f TAGS + +# Tags for emacs +tags: + etags `find . -name \*.[ch] -print` + +### Error codes. +BUILT_SOURCES = pswarpErrorCodes.h pswarpErrorCodes.c +CLEANFILES = pswarpErrorCodes.h pswarpErrorCodes.c +EXTRA_DIST = pswarpErrorCodes.dat pswarpErrorCodes.c.in pswarpErrorCodes.h.in + +pswarpErrorCodes.h : pswarpErrorCodes.dat pswarpErrorCodes.h.in + $(ERRORCODES) --data=pswarpErrorCodes.dat --outdir=. pswarpErrorCodes.h + +pswarpErrorCodes.c : pswarpErrorCodes.dat pswarpErrorCodes.c.in pswarpErrorCodes.h + $(ERRORCODES) --data=pswarpErrorCodes.dat --outdir=. pswarpErrorCodes.c Index: /tags/rel-0-1-0/pswarp/src/pswarp.c =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarp.c (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarp.c (revision 11633) @@ -0,0 +1,47 @@ +# include "pswarp.h" + +static void usage (void) { + fprintf (stderr, "USAGE: pswarp [-file image(s)] [-list imagelist] [options] (output) (skycell)\n"); + fprintf (stderr, " options:\n"); + fprintf (stderr, " [-astrom astrom.cmp] : provide an alternative astrometry calibration\n"); + fprintf (stderr, " [-mask mask.fits] : provide a corresponding mask image\n"); + fprintf (stderr, " [-weight weight.fits] : provide a corresponding weight image (pixel varience)\n"); + psErrorStackPrint(stderr, "\n"); + exit (2); +} + +int main (int argc, char **argv) { + + psTimerStart ("complete"); + + pmConfig *config = pswarpArguments (argc, argv); + if (!config) usage (); + + // load identify the data sources + if (!pswarpParseCamera (config)) { + psErrorStackPrint(stderr, "error setting up the camera\n"); + exit (1); + } + + // load the skycell layout information + if (!pswarpDefine (config)) { + psErrorStackPrint(stderr, "error loading output definition\n"); + exit (1); + } + + // load and warp + if (!pswarpDataLoad (config)) { + psErrorStackPrint(stderr, "error reading input data\n"); + exit (1); + } + + // write out result image + if (!pswarpDataSave (config)) { + psErrorStackPrint(stderr, "error writing output data\n"); + exit (1); + } + + psLogMsg ("pswarp", 3, "complete pswarp run: %f sec\n", psTimerMark ("complete")); + pswarpCleanup (config); + exit (EXIT_SUCCESS); +} Index: /tags/rel-0-1-0/pswarp/src/pswarp.h =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarp.h (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarp.h (revision 11633) @@ -0,0 +1,52 @@ +# ifdef HAVE_CONFIG_H +# include +# endif + +# include +# include +# include // for strcasecmp +# include // for unlink +# include +# include + +# include "pswarpErrorCodes.h" +# define PSWARP_RECIPE "PSWARP" // Name of the recipe to use + +// a single pswarpMap converts coordinates from one image to a second image +// the linear model is only valid over a limited range of pixels +typedef struct { + double Xo, Xx, Xy; + double Yo, Yx, Yy; + int xo; + int yo; +} pswarpMap; + +// the pswarpMapGrid carries a collection of pswarpMag structures representing the +// local value of the pswarpMap at different locations in the image. +typedef struct { + pswarpMap ***maps; + int nXpts, nYpts; // number of x,y samples in the grid + int nXpix, nYpix; // x,y spacing in src image pixels of grid samples + int xMin, yMin; // coordinate of first grid sample +} pswarpMapGrid; + +pmConfig *pswarpArguments (int argc, char **argv); +bool pswarpParseCamera (pmConfig *config); +bool pswarpDefine (pmConfig *config); +bool pswarpDataLoad (pmConfig *config); +bool pswarpDataSave (pmConfig *config); +void pswarpCleanup (pmConfig *config); +bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config); +bool pswarpTransformReadout_Opt (pmReadout *output, pmReadout *input, pmConfig *config); + +bool pswarpMatchRange (int *minX, int *minY, int *maxX, int *maxY, pmReadout *dest, pmReadout *src); + +pswarpMap *pswarpMapAlloc (); +pswarpMapGrid *pswarpMapGridAlloc (int Nx, int Ny); + +pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix); +bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY); +bool pswarpMapGridNextGrid (pswarpMapGrid *grid, int gridX, int gridY, int *nextX, int *nextY); +double pswarpMapGridMaxError (pswarpMapGrid *grid); +bool pswarpMapApply (double *outX, double *outY, pswarpMap *map, double inX, double inY); +bool pswarpMapSetLocalModel (pswarpMap *map, pmReadout *dest, pmReadout *src, int ix, int iy); Index: /tags/rel-0-1-0/pswarp/src/pswarpArguments.c =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpArguments.c (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpArguments.c (revision 11633) @@ -0,0 +1,70 @@ +# include "pswarp.h" +# include + +pmConfig *pswarpArguments (int argc, char **argv) { + + bool status; + int N; + + if (argc == 1) { + psError(PSWARP_ERR_ARGUMENTS, true, "No arguments supplied"); + return NULL; + } + + // load config data from default locations + pmConfig *config = pmConfigRead(&argc, argv, PSWARP_RECIPE); + if (config == NULL) { + psError(PSWARP_ERR_CONFIG, false, "Can't read site configuration"); + return NULL; + } + + // save the following additional recipe values based on command-line options + // these options override the PSWARP recipe values loaded from recipe files + psMetadata *options = pmConfigRecipeOptions (config, PSWARP_RECIPE); + + // astrom : used to supply astrometric parameters + if ((N = psArgumentGet (argc, argv, "-astrom"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "ASTROM_SOURCE", PS_META_REPLACE, "", argv[N]); + psArgumentRemove (N, &argc, argv); + } + + // chip selection is used to limit chips to be processed + if ((N = psArgumentGet (argc, argv, "-chip"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_DATA_STRING, "", psStringCopy(argv[N])); + psArgumentRemove (N, &argc, argv); + } + + // drop the local view on the options (saved in config->arguments) + psFree (options); + + status = pmConfigFileSetsMD (config->arguments, config, "INPUT", "-file", "-list"); + if (!status) { + psError(PSWARP_ERR_ARGUMENTS, true, "Missing -file (input) or -list (input)"); + return NULL; + } + + // the mask and weight entries are optional (build from gain?) + pmConfigFileSetsMD (config->arguments, config, "MASK", "-mask", "-masklist"); + pmConfigFileSetsMD (config->arguments, config, "WEIGHT", "-weight", "-weightlist"); + + if (argc != 3) { + psError(PSWARP_ERR_ARGUMENTS, true, "Incorrect arguments supplied"); + return NULL; + } + + psArray *array; + + // output position is fixed + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "", argv[1]); + + // skycell position is fixed + array = psArrayAlloc(1); + array->data[0] = psStringCopy (argv[2]); + status = psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "SKYCELL", PS_DATA_ARRAY, "", array); + psFree (array); + + psTrace("pswarp", 1, "Done with pswarpArguments...\n"); + return (config); +} Index: /tags/rel-0-1-0/pswarp/src/pswarpCleanup.c =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpCleanup.c (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpCleanup.c (revision 11633) @@ -0,0 +1,17 @@ +# include "pswarp.h" + +void pswarpCleanup (pmConfig *config) { + + psFree (config); + + psTimerStop (); + psMemCheckCorruption(stderr, true); + // XXX not used in pswarp pmModelGroupCleanup (); + psTimeFinalize (); + pmConceptsDone (); + pmConfigDone (); + // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "pswarp"); + fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "pswarp"); + + return; +} Index: /tags/rel-0-1-0/pswarp/src/pswarpDataLoad.c =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpDataLoad.c (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpDataLoad.c (revision 11633) @@ -0,0 +1,98 @@ +# include "pswarp.h" + +// this loop loads the data from the input files and selects the +// brighter stars for astrometry +// at the end of this function, the complete stellar data is loaded +// into the correct fpa structure locations (readout.analysis:PSPHOT.SOURCES) + +// all of the different astrometry analysis modes use the same data load loop +bool pswarpDataLoad (pmConfig *config) { + + pmChip *chip; + pmCell *cell; + pmReadout *readout; + pmFPAview *view; + + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSWARP.INPUT"); + if (!input) { + psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n"); + return false; + } + + // select the output readout + view = pmFPAviewAlloc (0); + view->chip = 0; + view->cell = 0; + view->readout = 0; + pmReadout *output = pmFPAfileThisReadout (config->files, view, "PSWARP.OUTPUT"); + if (!output) { + psError(PSWARP_ERR_CONFIG, true, "Can't find output data!\n"); + return false; + } + psFree (view); + + // de-activate PSWARP.SKYCELL and PSWARP.OUTPUT + pmFPAfileActivate (config->files, false, "PSWARP.SKYCELL"); + pmFPAfileActivate (config->files, false, "PSWARP.OUTPUT"); + + view = pmFPAviewAlloc (0); + + // find the FPA phu + // XXX wrap the test below into a function + // XXX need to read only the headers for the skycell + // XXX need to optionally load the astrometry datafile + bool bilevelAstrometry = false; + pmHDU *phu = pmFPAviewThisPHU (view, input->fpa); + if (phu) { + char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); + if (ctype) { + bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); + } + } + if (bilevelAstrometry) { + pmAstromReadBilevelMosaic (input->fpa, phu->header); + } + + // files associated with the science image + pmFPAfileIOChecks (config, view, PM_FPA_BEFORE); + + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + pmFPAfileIOChecks (config, view, PM_FPA_BEFORE); + + // read WCS data from the corresponding header + pmHDU *hdu = pmFPAviewThisHDU (view, input->fpa); + if (bilevelAstrometry) { + pmAstromReadBilevelChip (chip, hdu->header); + } else { + // XXX get pixelScale from recipes. does it matter? + float pixelScale = 13.5; + pmAstromReadWCS (input->fpa, chip, hdu->header, pixelScale); + } + + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + pmFPAfileIOChecks (config, view, PM_FPA_BEFORE); + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + pmFPAfileIOChecks (config, view, PM_FPA_BEFORE); + if (! readout->data_exists) { continue; } + + // XXX _Opt version uses locally-linear map + // pswarpTransformReadout (output, readout, config); + pswarpTransformReadout_Opt (output, readout, config); + + pmFPAfileIOChecks (config, view, PM_FPA_AFTER); + } + pmFPAfileIOChecks (config, view, PM_FPA_AFTER); + } + pmFPAfileIOChecks (config, view, PM_FPA_AFTER); + } + pmFPAfileIOChecks (config, view, PM_FPA_AFTER); + psFree (view); + return true; +} Index: /tags/rel-0-1-0/pswarp/src/pswarpDataSave.c =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpDataSave.c (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpDataSave.c (revision 11633) @@ -0,0 +1,60 @@ +# include "pswarp.h" + +// this loop loads the data from the input files and selects the +// brighter stars for astrometry +// at the end of this function, the complete stellar data is loaded +// into the correct fpa structure locations (readout.analysis:PSPHOT.SOURCES) + +// all of the different astrometry analysis modes use the same data load loop +bool pswarpDataSave (pmConfig *config) { + + pmChip *chip; + pmCell *cell; + pmReadout *readout; + + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSWARP.INPUT"); + if (!input) { + psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n"); + return false; + } + + // de-activate all files except PSWARP.OUTPUT + pmFPAfileActivate (config->files, false, NULL); + pmFPAfileActivate (config->files, true, "PSWARP.OUTPUT"); + pmFPAview *view = pmFPAviewAlloc (0); + + // files associated with the science image + pmFPAfileIOChecks (config, view, PM_FPA_BEFORE); + + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + pmFPAfileIOChecks (config, view, PM_FPA_BEFORE); + + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + pmFPAfileIOChecks (config, view, PM_FPA_BEFORE); + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + pmFPAfileIOChecks (config, view, PM_FPA_BEFORE); + if (! readout->data_exists) { continue; } + + // pswarpConvertReadout (readout, config); + + pmFPAfileIOChecks (config, view, PM_FPA_AFTER); + } + pmFPAfileIOChecks (config, view, PM_FPA_AFTER); + } + pmFPAfileIOChecks (config, view, PM_FPA_AFTER); + } + pmFPAfileIOChecks (config, view, PM_FPA_AFTER); + psFree (view); + + // re-activate all pmFPAfiles + pmFPAfileActivate (config->files, true, NULL); + + return true; +} Index: /tags/rel-0-1-0/pswarp/src/pswarpDefine.c =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpDefine.c (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpDefine.c (revision 11633) @@ -0,0 +1,57 @@ +# include "pswarp.h" +// this function loads the skycell metadata + +bool pswarpDefine (pmConfig *config) { + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE); + if (!recipe) { + psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n"); + return false; + } + + // select the input data sources + pmFPAfile *skycell = psMetadataLookupPtr (NULL, config->files, "PSWARP.SKYCELL"); + if (!skycell) { + psError(PSWARP_ERR_CONFIG, true, "Can't find skycell data!\n"); + return false; + } + + // open the full skycell file; no need to defer different depths + pmFPAview *view = pmFPAviewAlloc (0); + pmFPAfileOpen (skycell, view, config); + pmFPAfileRead (skycell, view, config); + + // XXX this is not a sufficient test + view->chip = 0; + pmChip *chip = pmFPAviewThisChip (view, skycell->fpa); + pmHDU *phu = pmFPAviewThisPHU (view, skycell->fpa); + pmHDU *hdu = pmFPAviewThisHDU (view, skycell->fpa); + bool bilevelAstrometry = false; + if (phu) { + char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); + if (ctype) { + bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); + } + } + if (bilevelAstrometry) { + pmAstromReadBilevelMosaic (skycell->fpa, phu->header); + pmAstromReadBilevelChip (chip, hdu->header); + } else { + // XXX get pixelScale from recipes. does it matter? + float pixelScale = 13.5; + pmAstromReadWCS (skycell->fpa, chip, hdu->header, pixelScale); + } + + // construct the output fpa here + pmFPAfile *output = psMetadataLookupPtr (NULL, config->files, "PSWARP.OUTPUT"); + if (!output) { + psError(PSWARP_ERR_CONFIG, true, "Can't find output data!\n"); + return false; + } + + pmFPAfileCreate(output, view, config); + + psFree (view); + return true; +} Index: /tags/rel-0-1-0/pswarp/src/pswarpErrorCodes.c.in =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpErrorCodes.c.in (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpErrorCodes.c.in (revision 11633) @@ -0,0 +1,25 @@ +/* + * The line + { PSWARP_ERR_$X{ErrorCode}, "$X{ErrorDescription}"}, + * (without the Xs) + * will be replaced by values from errorCodes.dat + */ +#include "pslib.h" +#include "pswarpErrorCodes.h" + +void pswarpErrorRegister(void) +{ + static psErrorDescription errors[] = { + { PSWARP_ERR_BASE, "First value we use; lower values belong to psLib" }, + { PSWARP_ERR_${ErrorCode}, "${ErrorDescription}"}, + }; + static int nerror = PSWARP_ERR_NERROR - PSWARP_ERR_BASE; // number of values in enum + + for (int i = 0; i < nerror; i++) { + psErrorDescription *tmp = psAlloc(sizeof(psErrorDescription)); + *tmp = errors[i]; + psErrorRegister(tmp, 1); + psFree(tmp); /* it's on the internal list */ + } + nerror = 0; // don't register more than once +} Index: /tags/rel-0-1-0/pswarp/src/pswarpErrorCodes.dat =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpErrorCodes.dat (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpErrorCodes.dat (revision 11633) @@ -0,0 +1,10 @@ +# +# This file is used to generate pswarpErrorClasses.h +# +BASE = 800 First value we use; lower values belong to psLib +UNKNOWN Unknown PM error code +NOT_IMPLEMENTED Desired feature is not yet implemented +ARGUMENTS Incorrect arguments +CONFIG Problem in configure files +IO Problem in FITS I/O +DATA Problem in data values Index: /tags/rel-0-1-0/pswarp/src/pswarpErrorCodes.h.in =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpErrorCodes.h.in (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpErrorCodes.h.in (revision 11633) @@ -0,0 +1,18 @@ +#if !defined(PSWARP_ERROR_CODES_H) +#define PSWARP_ERROR_CODES_H +/* + * The line + * PSWARP_ERR_$X{ErrorCode}, + * (without the X) + * + * will be replaced by values from errorCodes.dat + */ +typedef enum { + PSWARP_ERR_BASE = 600, + PSWARP_ERR_${ErrorCode}, + PSWARP_ERR_NERROR +} pswarpErrorCode; + +void pswarpErrorRegister(void); + +#endif Index: /tags/rel-0-1-0/pswarp/src/pswarpMapGrid.c =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpMapGrid.c (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpMapGrid.c (revision 11633) @@ -0,0 +1,204 @@ +# include "pswarp.h" + +// construct a grid with superpixel spacing of nXpix, nYpix +// XXX for the moment, ignore readout->cell->chip offsets +pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix) { + + int i, ni; + int j, nj; + + // split the difference of the remainder + int xMin = 0.5*(src->image->numCols % nXpix); + int yMin = 0.5*(src->image->numRows % nYpix); + + int nXpts = src->image->numCols / nXpix; + if (src->image->numCols % nXpix == 0) nXpts ++; + int nYpts = src->image->numRows / nYpix; + if (src->image->numRows % nYpix == 0) nYpts ++; + + pswarpMapGrid *grid = pswarpMapGridAlloc (nXpts, nYpts); + + for (ni = 0, i = xMin; ni < nXpts; i += nXpix, ni++) { + for (nj = 0, j = yMin; nj < nYpts; j += nYpix, nj++) { + pswarpMapSetLocalModel (grid->maps[ni][nj], dest, src, i, j); + } + } + + grid->nXpix = nXpix; + grid->nYpix = nYpix; + grid->xMin = xMin; + grid->yMin = yMin; + return grid; +} + +bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY) { + + *gridX = (ix - grid->xMin + 0.5*grid->nXpix) / grid->nXpix; + *gridY = (iy - grid->yMin + 0.5*grid->nYpix) / grid->nYpix; + return true; +} + +bool pswarpMapGridNextGrid (pswarpMapGrid *grid, int gridX, int gridY, int *nextX, int *nextY) { + + *nextX = gridX*grid->nXpix + grid->xMin + 0.5*grid->nXpix; + *nextY = gridY*grid->nYpix + grid->yMin + 0.5*grid->nYpix; + return true; +} + +// measure the max error accumulated in appling one grid point to its neighbors +double pswarpMapGridMaxError (pswarpMapGrid *grid) { + + double xRaw, yRaw; + double xRef, yRef; + double maxError = 0; + + for (int i = 0; i < grid->nXpts - 1; i++) { + for (int j = 0; j < grid->nYpts - 1; j++) { + + // measure the output coordinates for the next grid position using the current grid map + // compare with the coordinates measured using the next grid map + pswarpMapApply (&xRaw, &yRaw, grid->maps[i][j], grid->maps[i][j]->xo + grid->nXpix, grid->maps[i][j]->yo); + pswarpMapApply (&xRef, &yRef, grid->maps[i+1][j], grid->maps[i][j]->xo + grid->nXpix, grid->maps[i][j]->yo); + + double posError = hypot (xRaw-xRef, yRaw-yRef); + maxError = PS_MAX (maxError, posError); + } + } + return maxError; +} + +bool pswarpMapApply (double *outX, double *outY, pswarpMap *map, double inX, double inY) { + + *outX = map->Xo + map->Xx*inX + map->Xy*inY; + *outY = map->Yo + map->Yx*inX + map->Yy*inY; + + return true; +} + +// determine the map for the given pixel from src to dest. pixel is in src coords +bool pswarpMapSetLocalModel (pswarpMap *map, pmReadout *dest, pmReadout *src, int ix, int iy) { + + pmCell *cell = NULL; + + cell = src->parent; + pmChip *chipSrc = cell->parent; + pmFPA *fpaSrc = chipSrc->parent; + + cell = dest->parent; + pmChip *chipDest = cell->parent; + pmFPA *fpaDest = chipDest->parent; + + // XXX save these as static for speed? + psPlane *offset = psPlaneAlloc(); + + psPlane *FP = psPlaneAlloc(); + psPlane *TP = psPlaneAlloc(); + psSphere *sky = psSphereAlloc(); + + psPlane *V00 = psPlaneAlloc(); + psPlane *V10 = psPlaneAlloc(); + psPlane *V01 = psPlaneAlloc(); + + // XXX need to include readout->cell->chip offsets + + // V(0,0) position + offset->x = ix; + offset->y = iy; + psPlaneTransformApply(FP, chipSrc->toFPA, offset); + psPlaneTransformApply (TP, fpaSrc->toTPA, FP); + psDeproject (sky, TP, fpaSrc->toSky); + psProject (TP, sky, fpaDest->toSky); + psPlaneTransformApply (FP, fpaDest->fromTPA, TP); + psPlaneTransformApply (V00, chipDest->fromFPA, FP); + + // V(1,0) position + offset->x = ix + 1; + offset->y = iy; + psPlaneTransformApply(FP, chipSrc->toFPA, offset); + psPlaneTransformApply (TP, fpaSrc->toTPA, FP); + psDeproject (sky, TP, fpaSrc->toSky); + psProject (TP, sky, fpaDest->toSky); + psPlaneTransformApply (FP, fpaDest->fromTPA, TP); + psPlaneTransformApply (V10, chipDest->fromFPA, FP); + + // V(0,1) position + offset->x = ix; + offset->y = iy + 1; + psPlaneTransformApply(FP, chipSrc->toFPA, offset); + psPlaneTransformApply (TP, fpaSrc->toTPA, FP); + psDeproject (sky, TP, fpaSrc->toSky); + psProject (TP, sky, fpaDest->toSky); + psPlaneTransformApply (FP, fpaDest->fromTPA, TP); + psPlaneTransformApply (V01, chipDest->fromFPA, FP); + + map->Xx = V10->x - V00->x; + map->Xy = V01->x - V00->x; + map->Xo = V00->x - map->Xx*ix - map->Xy*iy; + + map->Yx = V10->y - V00->y; + map->Yy = V01->y - V00->y; + map->Yo = V00->y - map->Yx*ix - map->Yy*iy; + + map->xo = ix; + map->yo = iy; + + psFree (offset); + psFree (FP); + psFree (TP); + psFree (sky); + + psFree (V00); + psFree (V10); + psFree (V01); + + return true; +} + +static void pswarpMapFree (pswarpMap *map) { + return; +} + +pswarpMap *pswarpMapAlloc () { + + pswarpMap *map = (pswarpMap *) psAlloc (sizeof(pswarpMap)); + psMemSetDeallocator(map, (psFreeFunc) pswarpMapFree); + + return map; +} + +static void pswarpMapGridFree (pswarpMapGrid *grid) { + + if (grid == NULL) return; + if (grid->maps == NULL) return; + + for (int i = 0; i < grid->nXpts; i++) { + for (int j = 0; j < grid->nYpts; j++) { + psFree (grid->maps[i][j]); + } + psFree (grid->maps[i]); + } + psFree (grid->maps); + return; +} + +pswarpMapGrid *pswarpMapGridAlloc (int nXpts, int nYpts) { + + pswarpMapGrid *grid = (pswarpMapGrid *) psAlloc (sizeof(pswarpMapGrid)); + psMemSetDeallocator(grid, (psFreeFunc) pswarpMapGridFree); + + grid->maps = psAlloc (nXpts*sizeof(void **)); + for (int i = 0; i < nXpts; i++) { + grid->maps[i] = psAlloc (nYpts*sizeof(void *)); + for (int j = 0; j < nYpts; j++) { + grid->maps[i][j] = pswarpMapAlloc(); + } + } + grid->nXpts = nXpts; + grid->nYpts = nYpts; + + grid->nXpix = 0; + grid->nYpix = 0; + + return grid; +} + Index: /tags/rel-0-1-0/pswarp/src/pswarpMatchRange.c =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpMatchRange.c (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpMatchRange.c (revision 11633) @@ -0,0 +1,90 @@ +# include "pswarp.h" + +# define TEST_MINMAX \ + psPlaneTransformApply(FP, chipDest->toFPA, destPix); \ + psPlaneTransformApply (TP, fpaDest->toTPA, FP); \ + psDeproject (sky, TP, fpaDest->toSky); \ + psProject (TP, sky, fpaSrc->toSky); \ + psPlaneTransformApply (FP, fpaSrc->fromTPA, TP); \ + psPlaneTransformApply (srcPix, chipSrc->fromFPA, FP); \ + *minX = PS_MIN (*minX, srcPix->x); \ + *minY = PS_MIN (*minY, srcPix->y); \ + *maxX = PS_MAX (*maxX, srcPix->x); \ + *maxY = PS_MAX (*maxY, srcPix->y); + +// we are warping from src to dest. find the max overlapping pixels in the INPUT (src) +// coordinate frame +bool pswarpMatchRange (int *minX, int *minY, int *maxX, int *maxY, pmReadout *dest, pmReadout *src) { + + // transform input corners and edge centers to output coords + // find max overlapping region of output image + + pmCell *cell = NULL; + + cell = src->parent; + pmChip *chipSrc = cell->parent; + pmFPA *fpaSrc = chipSrc->parent; + + cell = dest->parent; + pmChip *chipDest = cell->parent; + pmFPA *fpaDest = chipDest->parent; + + *minX = src->image->numCols; + *minY = src->image->numRows; + *maxX = 0; + *maxY = 0; + + // XXX save these as static for speed? + psPlane *srcPix = psPlaneAlloc(); + psPlane *destPix = psPlaneAlloc(); + + psPlane *FP = psPlaneAlloc(); + psPlane *TP = psPlaneAlloc(); + psSphere *sky = psSphereAlloc(); + + destPix->x = 0; + destPix->y = 0; + TEST_MINMAX; + + destPix->x = dest->image->numCols; + destPix->y = 0; + TEST_MINMAX; + + destPix->x = dest->image->numCols; + destPix->y = dest->image->numRows; + TEST_MINMAX; + + destPix->x = 0; + destPix->y = dest->image->numRows; + TEST_MINMAX; + + destPix->x = 0.5*dest->image->numCols; + destPix->y = 0.0; + TEST_MINMAX; + + destPix->x = 0.0; + destPix->y = 0.5*dest->image->numRows; + TEST_MINMAX; + + destPix->x = 0.5*dest->image->numCols; + destPix->y = dest->image->numRows; + TEST_MINMAX; + + destPix->x = dest->image->numCols; + destPix->y = 0.5*dest->image->numRows; + TEST_MINMAX; + + *minX = PS_MAX (*minX, 0); + *minY = PS_MAX (*minY, 0); + *maxX = PS_MIN (*maxX, src->image->numCols); + *maxY = PS_MIN (*maxY, src->image->numRows); + + psFree (srcPix); + psFree (destPix); + psFree (FP); + psFree (TP); + psFree (sky); + + return true; +} + Index: /tags/rel-0-1-0/pswarp/src/pswarpParseCamera.c =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpParseCamera.c (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpParseCamera.c (revision 11633) @@ -0,0 +1,60 @@ +# include "pswarp.h" + +bool pswarpParseCamera (pmConfig *config) { + + bool status; + + // the input image(s) are required arguments; they define the camera + status = false; + pmFPAfile *input = pmFPAfileDefineFromArgs (&status, config, "PSWARP.INPUT", "INPUT"); + if (!input) { + psError(PSWARP_ERR_CONFIG, false, "Failed to build FPA from PSWARP.INPUT"); + return false; + } + + // the input skycell is a required argument: it defines the output image + status = false; + pmFPAfile *skycell = pmFPAfileDefineFromArgs (&status, config, "PSWARP.SKYCELL", "SKYCELL"); + if (!skycell) { + psError(PSWARP_ERR_CONFIG, false, "Failed to build FPA from PSWARP.SKYCELL"); + return false; + } + + // the mask is not required; the skycell is a required argument: it defines the output image + status = false; + pmFPAfileDefineFromArgs (&status, config, "PSWARP.MASK", "MASK"); + if (!status) { + // build the basic mask concept here + } + + // the weight is not required; the skycell is a required argument: it defines the output image + status = false; + pmFPAfileDefineFromArgs (&status, config, "PSWARP.WEIGHT", "WEIGHT"); + if (!status) { + // build the basic weight concept here + } + + // these calls bind the I/O handle to the specified fpa + if (!pmFPAfileDefineOutput (config, skycell->fpa, "PSWARP.OUTPUT")) { + psError(PSWARP_ERR_CONFIG, false, "Failed to build FPA from PSWARP.OUTPUT"); + return false; + } + + // Chip selection: turn on only the chips specified + char *chipLine = psMetadataLookupStr(NULL, config->arguments, "CHIP_SELECTIONS"); + psArray *chips = psStringSplitArray (chipLine, ",", false); + if (chips->n > 0) { + pmFPASelectChip (input->fpa, -1, true); // deselect all chips + for (int i = 0; i < chips->n; i++) { + int chipNum = atoi(chips->data[i]); + if (! pmFPASelectChip(input->fpa, chipNum, false)) { + psError(PSWARP_ERR_CONFIG, true, "Chip number %d doesn't exist in camera.\n", chipNum); + return false; + } + } + } + psFree (chips); + + psTrace("pswarp", 1, "Done with pswarpParseCamera...\n"); + return true; +} Index: /tags/rel-0-1-0/pswarp/src/pswarpTransformReadout.c =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpTransformReadout.c (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpTransformReadout.c (revision 11633) @@ -0,0 +1,78 @@ +# include "pswarp.h" + +bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config) { + + // XXX this implementation currently ignores the use of the region + psImage *region = NULL; + pmCell *cell = NULL; + + // select the current recipe + // psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE); + + int outNx = output->image->numCols; + int outNy = output->image->numRows; + + psPlane *inPix = psPlaneAlloc(); // Coordinates on the input detector + psPlane *outPix = psPlaneAlloc(); // Coordinates on the output detector + + psPlane *FP = psPlaneAlloc(); // Coordinates on the focal plane + psPlane *TP = psPlaneAlloc(); // Coordinates on the tangent plane + psSphere *sky = psSphereAlloc(); // Coordinates on the sky + + cell = input->parent; + pmChip *chipInput = cell->parent; + pmFPA *fpaInput = chipInput->parent; + + cell = output->parent; + pmChip *chipOutput = cell->parent; + pmFPA *fpaOutput = chipOutput->parent; + + psF32 **outData = output->image->data.F32; + psImage *inImage = input->image; + + // Iterate over the output image pixels + for (int y = 0; y < outNy; y++) { + for (int x = 0; x < outNx; x++) { + // Only transform those pixels requested + if (region && region->data.U8[y][x]) continue; + + // XXX double check this 1/2 pixel offset + outPix->x = (double)x + 0.5; + outPix->y = (double)y + 0.5; + + psPlaneTransformApply(FP, chipOutput->toFPA, outPix); + psPlaneTransformApply (TP, fpaOutput->toTPA, FP); + psDeproject (sky, TP, fpaOutput->toSky); + + psProject (TP, sky, fpaInput->toSky); + psPlaneTransformApply (FP, fpaInput->fromTPA, TP); + psPlaneTransformApply (inPix, chipInput->fromFPA, FP); + + // XXX get interpolation method from the recipe + outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR); + // outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, mask, 1, NAN, PS_INTERPOLATE_BILINEAR); + // modify zero and scale? + } + } + + psFree (inPix); + psFree (outPix); + psFree (FP); + psFree (TP); + psFree (sky); + + return true; +} + +# if (0) + if (error) { + // Error is actually the variance + outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error, + detector->x, + detector->y, + mask, 1, NAN); + } + if (error) { + outError->data.F32[y][x] = outError->data.F32[y][x] / SQUARE(scale); + } +# endif Index: /tags/rel-0-1-0/pswarp/src/pswarpTransformReadout_Opt.c =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpTransformReadout_Opt.c (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpTransformReadout_Opt.c (revision 11633) @@ -0,0 +1,82 @@ +# include "pswarp.h" + +// NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT +bool pswarpTransformReadout_Opt (pmReadout *output, pmReadout *input, pmConfig *config) { + + int minX, minY, maxX, maxY; + int gridX, gridY, nextGridX, nextGridY; + pswarpMap *map = NULL; + + // XXX this implementation currently ignores the use of the region + psImage *region = NULL; + + // select the current recipe + // psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + + psPlane *inPix = psPlaneAlloc(); // Coordinates on the input detector + psImage *inImage = input->image; + psF32 **outData = output->image->data.F32; + + // we might want to do the rectangular regions outside of the selection independently + psImageInit (output->image, NAN); + + // find the output pixel range + pswarpMatchRange (&minX, &minY, &maxX, &maxY, input, output); + + pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, 128, 128); + + // XXX need to modify the grid based on this result and force the maxError < XXX + double maxError = pswarpMapGridMaxError (grid); + fprintf (stderr, "maximum error using this grid sampling: %f\n", maxError); + + pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY); + pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY); + map = grid->maps[gridX][gridY]; + + // Iterate over the output image pixels + for (int y = minY; y < maxY; y++) { + if (y >= nextGridY) { + gridY ++; + pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY); + map = grid->maps[gridX][gridY]; + } + for (int x = minX; x < maxX; x++) { + if (x >= nextGridX) { + gridX ++; + pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY); + map = grid->maps[gridX][gridY]; + } + + // Only transform those pixels requested + if (region && region->data.U8[y][x]) continue; + + // XXX double check this 1/2 pixel offset + // XXX subtract 0.5,0.5 from result? + pswarpMapApply (&inPix->x, &inPix->y, map, x + 0.5, y + 0.5); + + // XXX get interpolation method from the recipe + // XXX include mask + // XXX apply scale and offset? + outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR); + } + } + + psFree (inPix); + psFree (grid); + return true; +} + +# if (0) + if (error) { + // Error is actually the variance + outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error, + detector->x, + detector->y, + mask, 1, NAN); + } + if (error) { + outError->data.F32[y][x] = outError->data.F32[y][x] / SQUARE(scale); + } +# endif + + Index: /tags/rel-0-1-0/pswarp/src/pswarpVersion.c =================================================================== --- /tags/rel-0-1-0/pswarp/src/pswarpVersion.c (revision 11633) +++ /tags/rel-0-1-0/pswarp/src/pswarpVersion.c (revision 11633) @@ -0,0 +1,27 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include "pswarp.h" + +static const char *cvsTag = "$Name: not supported by cvs2svn $";// CVS tag name + +psString pswarpVersion(void) +{ + psString version = NULL; // Version, to return + psStringAppend(&version, "%s-%s",PACKAGE_NAME,PACKAGE_VERSION); + return version; +} + +psString pswarpVersionLong(void) +{ + psString version = pswarpVersion(); // Version, to return + psString tag = psStringStripCVS(cvsTag, "Name"); // CVS tag + psStringAppend(&version, " (cvs tag %s) %s, %s", tag, __DATE__, __TIME__); + psFree(tag); + return version; +} + Index: /tags/rel-0-1-0/pswarp/test/simple.image.op =================================================================== --- /tags/rel-0-1-0/pswarp/test/simple.image.op (revision 11633) +++ /tags/rel-0-1-0/pswarp/test/simple.image.op (revision 11633) @@ -0,0 +1,56 @@ + +$NX = 2048 +$NY = 2048 + +macro go + mkref + mksrc + exec pswarp -file src.fits out ref.fits + exec rm -f src.fits out.fits ref.fits +end + +macro mkref + mcreate ref $NX $NY + keyword ref CTYPE1 -w "RA---TAN" + keyword ref CTYPE2 -w "DEC--TAN" + keyword ref CRVAL1 -wf 0.0 + keyword ref CRVAL2 -wf 0.0 + keyword ref CRPIX1 -wf {$NX/2} + keyword ref CRPIX2 -wf {$NY/2} + + keyword ref CDELT1 -wf {1/3600} + keyword ref CDELT2 -wf {1/3600} + + keyword ref PC001001 -wf 1.0 + keyword ref PC001002 -wf 0.0 + keyword ref PC002001 -wf 0.0 + keyword ref PC002002 -wf 1.0 + + wd ref ref.fits -bitpix -32 -bzero 0 -bscale 1 +end + +macro mksrc + mcreate src $NX $NY + for i 0 $NX 64 + for j 0 $NY 64 + zap src $i $j 1 1 -v 1 + end + end + + keyword src CTYPE1 -w "RA---TAN" + keyword src CTYPE2 -w "DEC--TAN" + keyword src CRVAL1 -wf 0.0 + keyword src CRVAL2 -wf 0.0 + keyword src CRPIX1 -wf {$NX/2} + keyword src CRPIX2 -wf {$NY/2} + + keyword src CDELT1 -wf {1/3600} + keyword src CDELT2 -wf {1/3600} + + keyword src PC001001 -wf {dcos(20)} + keyword src PC001002 -wf {-1*dsin(20)} + keyword src PC002001 -wf {dsin(20)} + keyword src PC002002 -wf {dcos(20)} + + wd src src.fits -bitpix -32 -bzero 0 -bscale 1 +end