Index: /branches/pap_branches/Price/pois/Makefile.am =================================================================== --- /branches/pap_branches/Price/pois/Makefile.am (revision 22256) +++ /branches/pap_branches/Price/pois/Makefile.am (revision 22256) @@ -0,0 +1,3 @@ +SUBDIRS = src swig + +CLEANFILES = *.pyc *~ core core.* Index: /branches/pap_branches/Price/pois/NOTES =================================================================== --- /branches/pap_branches/Price/pois/NOTES (revision 22256) +++ /branches/pap_branches/Price/pois/NOTES (revision 22256) @@ -0,0 +1,104 @@ +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 Index: /branches/pap_branches/Price/pois/autogen.sh =================================================================== --- /branches/pap_branches/Price/pois/autogen.sh (revision 22256) +++ /branches/pap_branches/Price/pois/autogen.sh (revision 22256) @@ -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-1.7 +AUTOHEADER=autoheader-2.59 +AUTOMAKE=automake-1.7 +AUTOCONF=autoconf-2.59 + +($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 --copy --foreign || 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: /branches/pap_branches/Price/pois/configure.ac =================================================================== --- /branches/pap_branches/Price/pois/configure.ac (revision 22256) +++ /branches/pap_branches/Price/pois/configure.ac (revision 22256) @@ -0,0 +1,27 @@ +dnl Process this file with autoconf to produce a configure script. + +AC_INIT([pois], [0.0.1], [pan-starrs.ifa.hawaii.edu]) +AC_CONFIG_SRCDIR([pois.pc.in]) + +AM_INIT_AUTOMAKE([dist-bzip2]) +AM_CONFIG_HEADER([config.h]) +AM_MAINTAINER_MODE + +AC_LANG(C) +AC_PROG_CC +AC_PROG_INSTALL +#AC_LIBTOOL_DLOPEN +AC_PROG_LIBTOOL + +PKG_CHECK_MODULES(PSLIB, pslib >= 0.5.0) + +dnl is this the best was to setup recursive CFLAGS? +pois_CFLAGS="-Wall -g -O3 -std=c99" +AC_SUBST([pois_CFLAGS]) + +AC_PROG_SWIG(1.3.24) +SWIG_MULTI_MODULE_SUPPORT +SWIG_PYTHON + +AC_CONFIG_FILES(pois.pc Makefile src/Makefile swig/Makefile python/Makefile) +AC_OUTPUT Index: /branches/pap_branches/Price/pois/python/ds9.py =================================================================== --- /branches/pap_branches/Price/pois/python/ds9.py (revision 22256) +++ /branches/pap_branches/Price/pois/python/ds9.py (revision 22256) @@ -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: /branches/pap_branches/Price/pois/python/pois.py =================================================================== --- /branches/pap_branches/Price/pois/python/pois.py (revision 22256) +++ /branches/pap_branches/Price/pois/python/pois.py (revision 22256) @@ -0,0 +1,315 @@ +# +# 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="testout.fits", + 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.""" + + 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'] + + 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 + config.outFile = fileDir + outFile + + # Read inputs + imageRegion = psLib.psRegionAlloc(0, 0, 0, 0) + reference = psLib.psFitsAlloc(config.refFile) + refImage = psLib.psFitsReadImage(None, reference, imageRegion, 0) + + psTrace("pois", 3, ("Read %s: %dx%d\n" % (config.refFile, refImage.numCols, refImage.numRows))) + del reference + + input = psLib.psFitsAlloc(config.inFile) + #inHeader = psLib.psFitsReadHeader(None, input) + inImage = psLib.psFitsReadImage(None, input, imageRegion, 0) + psTrace("pois", 3, ("Read %s: %dx%d\n" % (config.inFile, inImage.numCols, inImage.numRows))) + del input; 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.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? + if not pois.poisRejectStamps(stamps, mask, deviations, config): + break + + if display['stamps']: + for s in range(0, stamps.n): + stamp = pois.poisStamp_Cast(stamps.get_data(s)) # Stamp of interest + if stamp.status == pois.POIS_STAMP_RESET: + 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: + # 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 True: # start block scope + 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))) Index: /branches/pap_branches/Price/pois/python/ps.py =================================================================== --- /branches/pap_branches/Price/pois/python/ps.py (revision 22256) +++ /branches/pap_branches/Price/pois/python/ps.py (revision 22256) @@ -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: /branches/pap_branches/Price/pois/src/Makefile.am =================================================================== --- /branches/pap_branches/Price/pois/src/Makefile.am (revision 22256) +++ /branches/pap_branches/Price/pois/src/Makefile.am (revision 22256) @@ -0,0 +1,37 @@ +AM_CFLAGS = @pois_CFLAGS@ + +include_HEADERS = pois.h + +AM_CFLAGS += $(PSLIB_CFLAGS) +LIBS = $(PSLIB_LIBS) $(AM_LIBS) +AM_LDFLAGS = -release $(PACKAGE_VERSION) + +#noinst_PROGRAMS = pois +lib_LTLIBRARIES = libpois.la + +pois_SOURCES = \ + pois.c + +# +# Build a library +# +libpois_la_SOURCES = \ + pois.h \ + poisParseConfig.c \ + poisKernelBasisFunctions.c \ + poisCalculateEquation.c \ + poisConvolveImage.c \ + poisExtractKernel.c \ + poisMakeMask.c \ + poisFindStamps.c \ + poisCalculateDeviations.c \ + poisAddPenalty.c \ + poisStamp.c \ + poisRejectStamps.c \ + poisSolveEquation.c +# +# test program +# +test: pois test1.fits test2.fits + -$(RM) testout.fits stamp*.fits kernel*.fits + ./pois -v test1.fits test2.fits testout.fits Index: /branches/pap_branches/Price/pois/src/pois.c =================================================================== --- /branches/pap_branches/Price/pois/src/pois.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/pois.c (revision 22256) @@ -0,0 +1,266 @@ +/* + * 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. + * + */ + +#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 = getTime(); + + // Set trace levels + psTraceSetLevel(".", 0); + 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); + + // Set logging level + psLogSetLevel(9); + + // Parse the command line + poisConfig *config = poisParseConfig(argc, argv); // Configuration values + + // Read inputs + psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); + psFits *reference = psFitsAlloc(config->refFile); + //psMetadata *refHeader = psFitsReadHeader(NULL, reference); + psImage *refImage = psFitsReadImage(NULL, 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); + psFree(reference); + + psFits *input = psFitsAlloc(config->inFile); + //psMetadata *inHeader = psFitsReadHeader(NULL, input); + psImage *inImage = psFitsReadImage(NULL, 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); + psFree(input); + psFree(imageRegion); + + 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 = psFitsAlloc("mask.fits"); + 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"); + psFree(maskFile); +#endif + + psArray *stamps = NULL; // Array of stamps + psVector *stampMask = NULL; // Mask for stamps + psVector *solution = NULL; // Solution vector + + // 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 + 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 = psFitsAlloc(kernelName); + 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); + psFree(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); + } + + // 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_OK; + } + } + } + + // 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 = psFitsAlloc(convName); + if (!psFitsWriteImage(convFile, NULL, convImage, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write convolved image: %s\n", convName); + } + psFree(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); + (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; + } + } + } + + psFits *subFile = psFitsAlloc(config->outFile); + if (!psFitsWriteImage(subFile, NULL, subImage, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write subtracted image: %s\n", config->outFile); + } + psFree(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 + (void)psFree(kernelBasisFunctions); + (void)psFree(refImage); + (void)psFree(inImage); + (void)psFree(convImage); + (void)psFree(subImage); + (void)psFree(config); +} Index: /branches/pap_branches/Price/pois/src/pois.h =================================================================== --- /branches/pap_branches/Price/pois/src/pois.h (revision 22256) +++ /branches/pap_branches/Price/pois/src/pois.h (revision 22256) @@ -0,0 +1,202 @@ +#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 bad; // Value below which pixels are bad + 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 +} 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); + +/************************************************************************************************************/ + +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 + ); + +/************************************************************************************************************/ + +#endif Index: /branches/pap_branches/Price/pois/src/poisAddPenalty.c =================================================================== --- /branches/pap_branches/Price/pois/src/poisAddPenalty.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/poisAddPenalty.c (revision 22256) @@ -0,0 +1,57 @@ +#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: /branches/pap_branches/Price/pois/src/poisCalculateDeviations.c =================================================================== --- /branches/pap_branches/Price/pois/src/poisCalculateDeviations.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/poisCalculateDeviations.c (revision 22256) @@ -0,0 +1,100 @@ +#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 + + 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) { + psImage *refStamp = psImageSubset(refImage, x - xSize, y - ySize, x + xSize, y + ySize); + psImage *inStamp = psImageSubset(inImage, x - xSize, y - ySize, x + xSize, y + ySize); + psImage *maskStamp = psImageSubset(mask, x - xSize, y - ySize, x + xSize, y + ySize); + psImage *convRefStamp = poisConvolveImage(refStamp, maskStamp, solution, kernelParams, config); + // Calculate chi^2 + (void)psBinaryOp(subStamp, inStamp, "-", convRefStamp); + (void)psBinaryOp(subStamp, subStamp, "/", inStamp); + (void)psBinaryOp(subStamp, subStamp, "*", subStamp); + psImage *subStampTrim = psImageSubset(subStamp, config->xKernel, config->yKernel, + config->xKernel + 2 * footprint, + config->yKernel + 2 * footprint); + psImage *maskStampTrim = psImageSubset(maskStamp, config->xKernel, config->yKernel, + config->xKernel + 2 * footprint, + config->yKernel + 2 * footprint); + // Copy image to workaround bug 305 + psImage *tempImage = psImageCopy(NULL, subStampTrim, PS_TYPE_F32); + psImage *tempMask = psImageCopy(NULL, maskStampTrim, PS_TYPE_U8); + + (void)psImageStats(stats, tempImage, tempMask, POIS_MASK_BAD | POIS_MASK_NEAR_BAD); + + psFree(tempImage); + psFree(tempMask); + + deviations->data.F32[s] = stats->sampleMean * (float)footprint * (float)footprint * 4.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 = psFitsAlloc(stampName); + if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write stamp: %s\n", stampName); + } + psFree(stampFile); +#endif + +#if 0 + psFree(convRefStamp); + psFree(maskStampTrim); + psFree(subStampTrim); + psFree(maskStamp); + psFree(refStamp); + psFree(inStamp); +#endif + } + } + + psFree(stats); + psFree(subStamp); + + return deviations; +} Index: /branches/pap_branches/Price/pois/src/poisCalculateEquation.c =================================================================== --- /branches/pap_branches/Price/pois/src/poisCalculateEquation.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/poisCalculateEquation.c (revision 22256) @@ -0,0 +1,181 @@ +#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) + int xKernel = config->xKernel; // Half-size of kernel in x. + int yKernel = config->yKernel; // Half-size of kernel in y. + 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 + + psDPolynomial2D *poly = psDPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1, + PS_POLYNOMIAL_ORD); // 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); + psImage *matrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64); + for (int j = 0; j < numSolveParams; j++) { + for (int i = 0; i < numSolveParams; i++) { + matrix->data.F64[j][i] = 0.0; + } + } + stamp->matrix = matrix; + } + if (! stamp->vector) { + psTrace("pois.calculateEquation", 3, "Allocating vector: %d\n", numSolveParams); + psVector *vector = psVectorAlloc(numSolveParams, PS_TYPE_F64); + for (int i = 0; i < numSolveParams; i++) { + vector->data.F64[i] = 0.0; + } + stamp->vector = vector; + } + + // 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); + + // 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] = psDPolynomial2DEval(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: /branches/pap_branches/Price/pois/src/poisConvolveImage.c =================================================================== --- /branches/pap_branches/Price/pois/src/poisConvolveImage.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/poisConvolveImage.c (revision 22256) @@ -0,0 +1,103 @@ +#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 + + psDPolynomial2D *poly = psDPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1, + PS_POLYNOMIAL_ORD); // 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; + +#if 0 + // Evaluate the polynomial + poly->mask[xOrder][yOrder] = 0; + double polyVal = psDPolynomial2DEval(poly, imageX, imageY); + poly->mask[xOrder][yOrder] = 1; +#else + double polyVal = 1.0; +#endif + + 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: /branches/pap_branches/Price/pois/src/poisExtractKernel.c =================================================================== --- /branches/pap_branches/Price/pois/src/poisExtractKernel.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/poisExtractKernel.c (revision 22256) @@ -0,0 +1,70 @@ +#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; + } + } + + psDPolynomial2D *poly = psDPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1, + PS_POLYNOMIAL_ORD); // 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] * psDPolynomial2DEval(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: /branches/pap_branches/Price/pois/src/poisFindStamps.c =================================================================== --- /branches/pap_branches/Price/pois/src/poisFindStamps.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/poisFindStamps.c (revision 22256) @@ -0,0 +1,95 @@ +#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); + + float threshold = config->threshold;// Stamp threshold + 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 + psU8 **maskpix = mask->data.U8; // The mask 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 && pixels[y][x] >= threshold) { + // Check the footprint + int ok = 1; + for (int v = y - footprint; v <= y + footprint && ok; v++) { + for (int u = x - footprint; u <= x + footprint && ok; u++) { + if (maskpix[v][u] & + (POIS_MASK_BAD | POIS_MASK_NEAR_BAD | POIS_MASK_STAMP)) { + ok = 0; + } + } + } + if (ok) { + max = pixels[y][x]; + bestx = x; + besty = y; + } + } // Done checking the candidate pixel + } + } // 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: /branches/pap_branches/Price/pois/src/poisKernelBasisFunctions.c =================================================================== --- /branches/pap_branches/Price/pois/src/poisKernelBasisFunctions.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/poisKernelBasisFunctions.c (revision 22256) @@ -0,0 +1,77 @@ +#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: /branches/pap_branches/Price/pois/src/poisMakeMask.c =================================================================== --- /branches/pap_branches/Price/pois/src/poisMakeMask.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/poisMakeMask.c (revision 22256) @@ -0,0 +1,63 @@ +#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->bad || + 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->bad || + isnan(inImage->data.F32[y][x])) { + mask->data.U8[y][x] |= POIS_MASK_BAD; + } + } + } + + return mask; +} Index: /branches/pap_branches/Price/pois/src/poisParseConfig.c =================================================================== --- /branches/pap_branches/Price/pois/src/poisParseConfig.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/poisParseConfig.c (revision 22256) @@ -0,0 +1,165 @@ +#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 BAD] [-B BG] [-o N] [-n NSTAMPS] [-i ITER] [-r REJ] 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 BAD 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-i ITER Number of rejection iterations (10)\n" + "\t-r REJ Rejection level in standard deviations (2.5)\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->bad = 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; + + /* Parse command-line arguments using getopt */ + while ((opt = getopt(argc, argv, "hvk:f:q:t:s:o:b:B:n:i:r:p:")) != -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(optarg, "%f", &config->bad) != 1) { + 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; + 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: /branches/pap_branches/Price/pois/src/poisRejectStamps.c =================================================================== --- /branches/pap_branches/Price/pois/src/poisRejectStamps.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/poisRejectStamps.c (revision 22256) @@ -0,0 +1,66 @@ +#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? + psVector *devMask = psVectorAlloc(stamps->n, PS_TYPE_U8); // Mask for statistics + + // Don't want the standard deviation, but rather the deviation from zero + double meanDev = 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) { + meanDev += SQUARE(deviations->data.F32[i]); + numDev++; + } + } + float rmsDev = sqrt(meanDev / (double)(numDev - 1)); + float limit = rmsDev * config->sigmaRej; + psTrace("pois.rejectStamps", 2, "RMS deviation: %f\n", rmsDev); + 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) { + 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: /branches/pap_branches/Price/pois/src/poisSolveEquation.c =================================================================== --- /branches/pap_branches/Price/pois/src/poisSolveEquation.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/poisSolveEquation.c (revision 22256) @@ -0,0 +1,70 @@ +#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) { + +#ifdef TESTING + printf("Matrix %d:\n", s); + 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 + + (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]); + + // Clean up + psFree(luMatrix); + psFree(permutation); + psFree(matrix); + psFree(vector); + + return solution; +} Index: /branches/pap_branches/Price/pois/src/poisStamp.c =================================================================== --- /branches/pap_branches/Price/pois/src/poisStamp.c (revision 22256) +++ /branches/pap_branches/Price/pois/src/poisStamp.c (revision 22256) @@ -0,0 +1,27 @@ +#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; + p_psMemSetDeallocator(stamp, (psFreeFcn)poisStampFree); + + return stamp; +} + +void poisStampFree(poisStamp *stamp) +{ + if (stamp->matrix) { + psFree(stamp->matrix); + } + if (stamp->vector) { + psFree(stamp->vector); + } + psFree(stamp); +} Index: /branches/pap_branches/Price/pois/swig/.cvsignore =================================================================== --- /branches/pap_branches/Price/pois/swig/.cvsignore (revision 22256) +++ /branches/pap_branches/Price/pois/swig/.cvsignore (revision 22256) @@ -0,0 +1,7 @@ +*.pyc +*_wrap.c +psLib.py +swigrun.py +xpa.py +.Makefile.depend +*.pm Index: /branches/pap_branches/Price/pois/swig/Makefile.am =================================================================== --- /branches/pap_branches/Price/pois/swig/Makefile.am (revision 22256) +++ /branches/pap_branches/Price/pois/swig/Makefile.am (revision 22256) @@ -0,0 +1,83 @@ +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 += $(XPA_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_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 +# +#_xpa.so : $(XPA_OBJS) xpa.py +# cc -o _xpa.so -bundle -undefined suppress -flat_namespace $(XPA_OBJS) -L$(XPA_DIR)/lib -lxpa $(PS_LIBS) +XPA_DIR = /Users/rhl/ACT/products/Darwin/xpa/v2_1_5 +XPA_CFLAGS = -I$(XPA_DIR)/include +XPA_LIBS = -L$(XPA_DIR)/lib -lxpa + +libxpaSwig_la_LIBADD = $(PSLIB_LIBS) +libxpaSwig_la_LIBADD += $(XPA_LIBS) + +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) -I$(XPA_DIR)/include 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 lib$$f.so _$$f.so; done + +uninstall-local: + for f in $(SWIG_LIB_NAMES); do $(RM) $(prefix)/lib/_$$f.so; done + +CLEANFILES = *_wrap.c *.pyc poisSwig.py psLibSwig.py xpaSwig.py *~ core core.* Index: /branches/pap_branches/Price/pois/swig/p_psSwig.i =================================================================== --- /branches/pap_branches/Price/pois/swig/p_psSwig.i (revision 22256) +++ /branches/pap_branches/Price/pois/swig/p_psSwig.i (revision 22256) @@ -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: /branches/pap_branches/Price/pois/swig/poisSwig.i =================================================================== --- /branches/pap_branches/Price/pois/swig/poisSwig.i (revision 22256) +++ /branches/pap_branches/Price/pois/swig/poisSwig.i (revision 22256) @@ -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(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: /branches/pap_branches/Price/pois/swig/psLibForSwig.c =================================================================== --- /branches/pap_branches/Price/pois/swig/psLibForSwig.c (revision 22256) +++ /branches/pap_branches/Price/pois/swig/psLibForSwig.c (revision 22256) @@ -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: /branches/pap_branches/Price/pois/swig/psLibForSwig.h =================================================================== --- /branches/pap_branches/Price/pois/swig/psLibForSwig.h (revision 22256) +++ /branches/pap_branches/Price/pois/swig/psLibForSwig.h (revision 22256) @@ -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: /branches/pap_branches/Price/pois/swig/psLibSwig.i =================================================================== --- /branches/pap_branches/Price/pois/swig/psLibSwig.i (revision 22256) +++ /branches/pap_branches/Price/pois/swig/psLibSwig.i (revision 22256) @@ -0,0 +1,353 @@ +// -*- 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. + PS_TYPE_PTR = 0x0000 ///< Something else that's not supported for arithmetic. +} 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: /branches/pap_branches/Price/pois/swig/simpleFits.c =================================================================== --- /branches/pap_branches/Price/pois/swig/simpleFits.c (revision 22256) +++ /branches/pap_branches/Price/pois/swig/simpleFits.c (revision 22256) @@ -0,0 +1,684 @@ +/* + * 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" + +#if !defined(PS_LITTLE_ENDIAN) +# if defined(Darwin) +# define PS_LITTLE_ENDIAN 1 +# endif +#endif + +#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_TYPE_F64: + card.type = DOUBLE_CARD; + card.val.d = meta->data.F64; + break; + case PS_TYPE_F32: + card.type = FLOAT_CARD; + card.val.f = meta->data.F32; + break; + case PS_TYPE_S32: + card.type = INTEGER_CARD; + card.val.i = meta->data.S32; + break; + case PS_TYPE_PTR: + 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 defined(PS_LITTLE_ENDIAN) /* 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: /branches/pap_branches/Price/pois/swig/simpleFits.h =================================================================== --- /branches/pap_branches/Price/pois/swig/simpleFits.h (revision 22256) +++ /branches/pap_branches/Price/pois/swig/simpleFits.h (revision 22256) @@ -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: /branches/pap_branches/Price/pois/swig/xpaSwig.i =================================================================== --- /branches/pap_branches/Price/pois/swig/xpaSwig.i (revision 22256) +++ /branches/pap_branches/Price/pois/swig/xpaSwig.i (revision 22256) @@ -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"