Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Makefile =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Makefile (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Makefile (revision 21655) @@ -0,0 +1,30 @@ +CC = cc -std=c99 +CFLAGS = -Wall -g $(INCLUDES) + +# +# N.b. Utils must come first: +# +DIRS = Utils Metadata + +.PHONY : init +init : + @for d in $(DIRS); do \ + echo "cd $$d; $(MAKE)"; \ + (cd $$d; $(MAKE)); \ + done +# +_tags :; +tags : _tags + etags */*.[ch] +# +clean : + $(RM) $(PROGS) *.o *.a *.log *.dvi *.aux *.toc *.log *.out *~ core TAGS + @for d in $(DIRS); do \ + echo "cd $$d; $(MAKE) clean"; \ + (cd $$d; $(MAKE) clean); \ + done +empty : clean + @for d in $(DIRS); do \ + echo "cd $$d; $(MAKE) empty"; \ + (cd $$d; $(MAKE) empty); \ + done Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/Makefile =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/Makefile (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/Makefile (revision 21655) @@ -0,0 +1,23 @@ +CC = cc -std=c99 +CFLAGS = -Wall -g $(INCLUDES) + +INCLUDES = -I$(UTILS_DIR) +LIBS = -L$(UTILS_DIR) -lUtils + +UTILS_DIR = ../Utils + +PROGS = metadata + +.PHONY: init +init : $(PROGS) + +metadata : metadata.o + $(CC) -o metadata metadata.o $(LIBS) +# +_tags :; +tags : _tags + etags *.[ch] +# +clean : + $(RM) $(PROGS) *.o *.a *.log *.dvi *.aux *.toc *.log *.out *~ core TAGS +empty : clean Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/metadata.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/metadata.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/metadata.c (revision 21655) @@ -0,0 +1,209 @@ +/* + * A strawman implementation of the Pan-STARRS metadata system + */ +#include +#include +#include +#include +#include "psUtils.h" +#include "psMetaData.h" +/* + * Return a properly filled-out psMetaData + * + * N.b. For types PS_META_FLOAT, PS_META_INT, and PS_META_STR, a COPY + * of the data will be kept. In all other cases, only a POINTER + * to the data will be kept. + * + * An unused argument, copy, is provided to handle this better + * once we have copy constructors available. + */ +psMetaData *psMetaDataNew( + const char *name, // name of new item of metadata + psMetaDataType type, // type of this piece of metadata + const void *data, // value of new item + const char *comment, // comment associated with item; + // may be NULL + int copy) // always copy data; NOTUSED +{ + static int id = 0; // unique ID for an item of metadata + psMetaData *ms = psAlloc(sizeof(psMetaData)); + + *(int *)&ms->id = id++; // need to cast away const + ms->type = type; + + assert (name != NULL); + ms->name = psAlloc(strlen(name) + 1); + strcpy(ms->name, name); + + if(comment == NULL) { + comment = ""; + } + ms->comment = psAlloc(strlen(comment) + 1); + strcpy(ms->comment, comment); + + assert (type >= 0 && type < PS_META_NTYPE); + switch (type) { + case PS_META_FLOAT: + *(float *)&ms->val.f = *(float *)data; + break; + + case PS_META_INT: + *(int *)&ms->val.i = *(int *)data; + break; + + case PS_META_STR: + *(void **)&ms->val.v = psAlloc(strlen(data) + 1); + strcpy(ms->val.v, data); + break; + + case PS_META_IMG: + case PS_META_JPEG: + case PS_META_PNG: + case PS_META_ASTROM: + case PS_META_UNKNOWN: + *(void **)ms->val.v = (void *)data; + break; + case PS_META_NTYPE: + psAbort(__FILE__, "You cannot get here (%d)\n", __LINE__); + break; // NOTREACHED + } + + return ms; +} + +void psMetaDataDel(psMetaData *ms) // piece of metadata to destroy +{ + if(ms == NULL) { + return; + } + + psFree(ms->name); + psFree(ms->comment); + + assert (ms->type >= 0 && ms->type < PS_META_NTYPE); + switch (ms->type) { + case PS_META_FLOAT: + case PS_META_INT: + break; // nothing to do + + case PS_META_STR: + psFree(ms->val.v); + break; + + case PS_META_IMG: + case PS_META_JPEG: + case PS_META_PNG: + case PS_META_ASTROM: + case PS_META_UNKNOWN: + break; // we don't own the memory + case PS_META_NTYPE: + psAbort(__FILE__, "You cannot get here (%d)\n", __LINE__); + break; // NOTREACHED + } + + psFree(ms); + +} + +/*****************************************************************************/ + +void psMetaDataPrint(FILE *fd, // file descriptor to write to + const psMetaData *ms) // item of metadata to print +{ + if(ms == NULL) { + return; + } + + fprintf(fd, "Name: %-40s (ID: %d)\n", ms->name, ms->id); + if (*ms->comment != '\0') { + fprintf(fd, "Comment: %s\n", ms->comment); + } + fprintf(fd, "Value: "); + + assert (ms->type >= 0 && ms->type < PS_META_NTYPE); + switch (ms->type) { + case PS_META_FLOAT: + fprintf(fd, "%g\n", ms->val.f); + break; + + case PS_META_INT: + fprintf(fd, "%d\n", ms->val.i); + break; + + case PS_META_STR: + fprintf(fd, "%s\n", (char *)ms->val.v); + break; + + case PS_META_IMG: + fprintf(fd, "(IMG)\n"); + break; + + case PS_META_JPEG: + fprintf(fd, "(JPG)\n"); + break; + + case PS_META_PNG: + fprintf(fd, "(PNG)\n"); + break; + + case PS_META_ASTROM: + fprintf(fd, "(ASTROM)\n"); + break; + + case PS_META_UNKNOWN: + fprintf(fd, "(UNKNOWN)\n"); + break; + + case PS_META_NTYPE: + psAbort(__FILE__, "You cannot get here (%d)\n", __LINE__); + break; // NOTREACHED + } +} + +/*****************************************************************************/ +/* + * Test code + */ +#if 1 +#include + +int main(void) +{ + psDlist *list = psDlistNew(NULL); // list of metadata + + const int ten = 10; + psDlistAppend(list, psMetaDataNew("const.ten", PS_META_INT, &ten, NULL,0)); + const float sqrt2 = sqrt(2); + psDlistAppend(list, psMetaDataNew("const.sqrt2", PS_META_FLOAT, &sqrt2, + "square root of two", 0)); + psDlistAppend(list, psMetaDataNew("lang.hello", PS_META_STR, "Bonjour", + "French", 0)); + /* + * Print that data + */ + for (psDlistElem *ptr = list->head; ptr != NULL; ptr = ptr->next) { + psMetaDataPrint(stdout, ptr->data); + } + + fprintf(stdout, "------\n"); + /* + * Print it again + */ + (void)psDlistGet(list, PS_DLIST_TAIL); // initialise iterator + psMetaData *ms = NULL; + while ((ms = psDlistGet(list, PS_DLIST_PREV)) != NULL) { + psMetaDataPrint(stdout, ms); + } + /* + * Clean up + */ + psDlistDel(list, (void (*)(void *))psMetaDataDel); + /* + * Check for memory leaks + */ + fprintf(stderr,"Checking for memory leaks:\n"); + psMemCheckLeaks(0, NULL, stderr); + + return 0; +} +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/psAstrom.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/psAstrom.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/psAstrom.h (revision 21655) @@ -0,0 +1,10 @@ +#if !defined(PS_ASTROM_H) +#define PS_ASTROM_H +/* + * This is not a full implementation... + */ +typedef struct { + const int id; +} psAstrom; + +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/psImage.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/psImage.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/psImage.h (revision 21655) @@ -0,0 +1,10 @@ +#if !defined(PS_IMAGE_H) +#define PS_IMAGE_H +/* + * This is not a full implementation... + */ +typedef struct { + const int id; +} psImage; + +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/psMetaData.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/psMetaData.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Metadata/psMetaData.h (revision 21655) @@ -0,0 +1,60 @@ +#if !defined(PS_META_DATA_H) +#define PS_META_DATA_H 1 + +#include "psImage.h" // for psImage +#include "psAstrom.h" // for psAstrom +/* + * Possible types of metadata. + */ +typedef enum { // type of val is: + PS_META_FLOAT = 0, // float (.f) + PS_META_INT, // int (.i) + PS_META_STR, // string (.v) + PS_META_IMG, // image (.v) + PS_META_JPEG, // JPEG (.v) + PS_META_PNG, // PNG (.v) + PS_META_ASTROM, // astrometric coefficients (.v) + PS_META_UNKNOWN, // other (.v) + PS_META_NTYPE // Number of types; must be last +} psMetaDataType; + +/* + * A struct to define a single item of metadata + */ +typedef struct { + const int id; // unique ID for this item + + char *name; // Name of item + psMetaDataType type; // type of this item + const union { + float f; // floating value + int i; // integer value + void *v; // other type + } val; // value of metadata + char *comment; // optional comment ("", not NULL) +} psMetaData; + +/*****************************************************************************/ +/* + * Constructors and destructors + */ +psMetaData *psMetaDataNew( + const char *name, // name of new item of metadata + psMetaDataType type, // type of this piece of metadata + const void *val, // value of new item + // N.b. a pointer even if the item + // is of type e.g. int + const char *comment, // comment associated with item + int copy); // always copy data; NOTUSED + +void psMetaDataDel(psMetaData *ms); // piece of metadata to destroy + +/*****************************************************************************/ +/* + * Utilities + */ +void psMetaDataPrint(FILE *fd, // file descriptor to write to + const psMetaData *ms); // item of metadata to print + + +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/.gdb_history =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/.gdb_history (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/.gdb_history (revision 21655) @@ -0,0 +1,477 @@ +rr +up +rr +b 180 +rr +p froot +p *froot +p *froot->subfacil[0] +p *froot->subfacil[0]->subfacil[0] +p *froot->subfacil[0]->subfacil[0]->subfacil[0] +rr +rr +dis +c +rr +d +b 101 +rr +c +p2 facil->name, cpt0 +p2 facil->name cpt0 +rr +dis +c +rr +rr +rr +rr +rr +rr +rr +b 258 +rr +call psGetTraceLevel("aaa") +call psGetTraceLevel("aaa.bbb") +call psGetTraceLevel("aaa.bbb.ccc") +call psGetTraceLevel("aaa.BBB") +call psGetTraceLevel("aaa.BBB.aa") +call psGetTraceLevel("") +rr +q +b 254 +rr +call psGetTraceLevel("aaa.bbb.ddd") +p highest_level +q +rr +rr +rr +ww +rr +ww +up +p fmt +q +rr +ww +up +p ap +whatis ap +list +p va_list +whatis va_list +rr +b 216 +rr +p name +p cachedName +b 212 +rr +c +p name +p cachedName +n +p name +p level +c +rr +cond 2 name[0] == 'u' +c +p2 name cachedName +n +b 108 +rr +p cpt0 +p rest +p *facil +rr +dis +c +c +rr +ena +rr +p *facil +n +rr +p *facil +rr +dis +c +c +rr +rr +up +l +p facil +p level +l +where +call psGetTraceLevel("aaa.bbb.ddd") +b 311 +rr +call psGetTraceLevel("aaa.bbb.ddd") +call psGetTraceLevel("aaa.bbb.ccc") +call psGetTraceLevel("aaa.bbb") +call psGetTraceLevel("") +call psGetTraceLevel("aaa.bbb.ccc") +b 207 +call psGetTraceLevel("aaa.bbb.ccc") +rr +dis 2 +c +ena +c +up +p *facil +down +p *facil +p level +rr +dis +c +ena 1 +rr +n +call psGetTraceLevel("aaa.bbb.ccc") +call psGetTraceLevel("") +p froot->level +rr +rr +rr +rr +rr +q +b 533 +b memory.c:533 +rr +rr +b memAllocateCB0 +b main +call psMemSetAllocateID(1) +r +call psMemSetAllocateID(1) +c +up +b psDlistDel +f0 +p *ptr +c +p *list +p *(psMemBlock *)list +p *((psMemBlock *)list - 1) +n +rr +c +q +p/x 512 +p/x 511 +4b +rr +b psHashDel +rr +n +p *table +p *table->buckets +p table->buckets[0]@16 +b 172 +rr +s +n +p hash +n +p ptr +n +rr +c +p table->buckets[0]@16 +c +rr +dis +c +rr +q +rr +b psHashDel +rr +n +l +b 93 +c +p *(psMemBlock *)(bucket - 1) +p *(psMemBlock *)(ptr - 1) +p *((psMemBlock *)ptr - 1) +p *ptr +c +b psFree +rr +c +ww +up +l +c +q +rr +rr +rr +ww +p *((psMemBlock *)ptr - 1) +p *((psMemBlock *)id - 1) +up +p table->buckets[0]@16 +up +p table->buckets[0]@16 +p i +b psHashDel +rr +p table->buckets[0]@16 +p *table->buckets[5] +p *table->buckets[4] +p *table->buckets[8] +p *table->buckets[4]->data +p table->buckets[4]->data +p *((psMemBlock*)table->buckets[4]->data - 1) +s +l +p i +n +p i +l +b 92 +c +n +p tmp +s +n +p *((psMemBlock*)table->buckets[4]->data - 1) +up +p *((psMemBlock*)table->buckets[4]->data - 1) +down +n +dis +rr +q +b psHashLookup +r +s +l +b 130 +c +p hash +p ptr +n +p ptr->key +p key +n +p ptr->data +q +q +rr +q +rr +b 217 +rr +n +p *id +p name +n +rr +n +p *name +p id +p id->name +c +dis +c +ww +q +rr +rr +q +rr +up +p ptr +p table->buckets[0]@16 +b 149 +rr +p *((psMemBlock *)ptr - 1) +n +p *((psMemBlock *)data - 1) +rr +c +q +rr +up +p fd +q +q +rr +q +rr +rrrr +r +ww +p2 i j +rr +q +rr +rr +b psLogMsg +rr +n +p hostname +p name +n +p name +n +p head +rr +dis +c +b p_psRealloc +rr +n +p *ptr +n +p *ptr +rr +n +p2 ptr *ptr +n +p2 ptr *ptr +c +rr +n +p *ptr +whatis ptr +p2 ptr *ptr +n +p2 ptr *ptr +n +n +p t +p *((psMemBlock *)t - 1) +p *((psMemBlock *)t->arr - 1) +dis +rr +q +rr +p *m +up +rr +ww +up +down +rr +q +rr +p ptrt +p ptr +p *ptr +p i +up +ww +b 287] +b 287 +rr +p list +p *list +p *list->head +p *list->head->next +p *list->head->next->next +p *list->head->next->next->data +p *(X *)(list->head->next->next->data) +n +p *ptr +rr +dis +c +q +rr +ww +b psDlistToArray +rr +p *dlist +n +p *arr +p *arr->arr[0] +p arr->arr[0] +n +p i +p n +p *ptr +n +p *ptr +n +p arr->arr[0]@3 +rr +c +ww +f 5 +p elem +p *elem +rr +dis +c +rr +rr +rr +dis +d +b 307 +rr +p arr +whatis arr +ptype *arr +ptype *arr->arr +ptype *arr->arr[0] +ptype psVoidPtr +whatis arr->arr +whatis arr->arr[0] +p (X *)(arr->arr[0])->i +p (X *)(arr->arr[0]) +rr +c +d +rr +rr +rr +rr +rr +up +ww +b psVoidPtrArrayDel +rr +p *arr +n +rr +dis +c +p arr +rr +up +down +p *vptr +p *ptr +ww +where +rr +q +b psLogMsg +rr +up +l +p file +p *ptr +up +l +rr +up +down +p *ptr +p sizeof(ptr) +p sizeof(*ptr) +q +rr +q Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/.log =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/.log (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/.log (revision 21655) @@ -0,0 +1,34 @@ +This is pdfTeX, Version 3.14159-1.10a-devel (Web2C 7.3.8) (format=pdflatex 2003.1.31) 2 MAR 2004 11:42 +**utils.tex +(/usr/local/tetex/share/texmf/tex/latex/tools/.tex{/usr/local/tetex/share/texmf +.local/pdftex/config/pdftex.cfg +Warning: pdftex (file /usr/local/tetex/share/texmf.local/pdftex/config/pdftex.c +fg): invalid line in config file: `pdf_minorversion 4' +} +LaTeX2e <2001/06/01> +Babel and hyphenation patterns for american, british, nohyphenation, lo +aded. +File ignored) +* +! Interruption. +<*> + +? +! Emergency stop. +<*> + +End of file on the terminal! + + +Here is how much of TeX's memory you used: + 6 strings out of 60983 + 162 string characters out of 1189764 + 44491 words of memory out of 2500001 + 3092 multiletter control sequences out of 10000+50000 + 3640 words of font info for 14 fonts, out of 500000 for 1000 + 22 hyphenation exceptions out of 1000 + 5i,0n,1p,89b,7s stack positions out of 1500i,500n,5000p,200000b,5000s + 0 PDF objects out of 300000 + 0 named destinations out of 131072 + 0 words of extra memory for PDF output out of 65536 +No pages of output. Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/Makefile =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/Makefile (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/Makefile (revision 21655) @@ -0,0 +1,26 @@ +CC = cc -std=c99 +CFLAGS = -Wall -g + +libUtils.a : array.o dlist.o hash.o logmsg.o memory.o misc.o trace.o + ar r libUtils.a $? + ranlib libUtils.a +# +# Test code (change #if 0 to #if 1 in ONE source file to build) +# +array : libUtils.a array.o + $(CC) -o array array.o libUtils.a +dlist : libUtils.a dlist.o + $(CC) -o dlist dlist.o libUtils.a +logmsg : libUtils.a logmsg.o + $(CC) -o logmsg logmsg.o libUtils.a +memory : libUtils.a + $(CC) -o memory memory.o libUtils.a +# +# +_tags :; +tags : _tags + etags *.[ch] +# +clean : + $(RM) $(PROGS) *.o *.a *.log *.dvi *.aux *.toc *.log *.out *~ core TAGS +empty : clean Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/Private/p_psMemory.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/Private/p_psMemory.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/Private/p_psMemory.h (revision 21655) @@ -0,0 +1,17 @@ +#if !defined(P_PS_MEMORY_H) +#define P_PS_MEMORY_H + +#define P_PS_MEMMAGIC (void *)0xdeadbeef // Magic number in psMemBlock header + +/*****************************************************************************/ +/* + * When to call the allocate/free callbacks. Set using + * psMemSetAllocateID/psMemSetFreeID, but it can be + * convenient to see these directly from the debugger + */ +extern long p_psMemAllocateID; // notify user this block is allocated +extern long p_psMemFreeID; // notify user this block is freed + +/*****************************************************************************/ + +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/array.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/array.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/array.c (revision 21655) @@ -0,0 +1,92 @@ +/* + * Arrays. Almost all the array support is in fact in psArray.h, + * but we directly support arrays of (void *) here, so as to be + * able to free each element properly + */ +#include "psUtils.h" + +typedef void *psVoidPtr; + +P_PS_CREATE_ARRAY_TYPE(static, my_, psVoidPtr); + +psVoidPtrArray *psVoidPtrArrayNew(int n, int s) +{ + return my_psVoidPtrArrayNew(n, s); +} + +psVoidPtrArray *psVoidPtrArrayRealloc(psVoidPtrArray *arr, int n) +{ + return my_psVoidPtrArrayRealloc(arr, n); +} + +void psVoidPtrArrayDel(psVoidPtrArray *arr, + void (*elemDel)(void *)) // destructor for data on array +{ + if (arr == NULL) { + return; + } + + for (int i = 0; i < arr->n; i++) { + if (elemDel == NULL) { + psMemDecrRefCounter(arr->arr[i]); + } else { + elemDel(psMemDecrRefCounter(arr->arr[i])); + } + } + my_psVoidPtrArrayDel(arr); +} + +/*****************************************************************************/ +#if 0 // testing + +typedef struct { + int x, y; +} psXY; + +psXY *psXYNew(void) +{ + return psAlloc(sizeof(psXY)); +} + +void psXYDel(psXY *xy) +{ + psFree(xy); +} + +PS_DECLARE_ARRAY_TYPE(psXY); +PS_CREATE_ARRAY_TYPE(psXY); + +PS_DECLARE_ARRAY_PTR_TYPE(psXY); +PS_CREATE_ARRAY_PTR_TYPE(psXY); + +int main(void) +{ + psXYArray *t = psXYArrayNew(10, 15); + psXYPtrArray *pt = psXYPtrArrayNew(10, 10); + + for (int i = 0; i < t->n; i++) { + t->arr[i].x = i; + pt->arr[i]->y = 10*i; + } + + t = psXYArrayRealloc(t, 5); + t = psXYArrayRealloc(t, 8); + + for (int i = 0; i < t->n; i++) { + printf("%d %d ", t->arr[i].x, pt->arr[i]->y); + } + printf("\n"); + + psXYArrayDel(t); + + psXY *xy = psMemDecrRefCounter(pt->arr[0]); + pt->arr[0] = NULL; + psXYDel(xy); + + psXYPtrArrayDel(pt); + + psMemCheckLeaks(0, NULL, stderr); + + return 0; +} +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/dlist.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/dlist.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/dlist.c (revision 21655) @@ -0,0 +1,360 @@ +/* + * Minimal support for doubly linked lists + */ +#include +#include +#include +#include "psUtils.h" + +static psDlistElem *dlistElemNew(void) +{ + return(psAlloc(sizeof(psDlistElem))); +} + +static void dlistElemDel(psDlistElem *elem) +{ + psFree(elem); +} + +/*****************************************************************************/ + +psDlist *psDlistNew(void *data) // initial data item; may be NULL +{ + psDlist *list = psAlloc(sizeof(psDlist)); + + list->n = 0; + list->head = list->tail = NULL; + list->iter = NULL; + if (data != NULL) { + psDlistAdd(list, data, PS_DLIST_TAIL); + } + + return list; +} + +void psDlistDel(psDlist *list, // list to destroy + void (*elemDel)(void *)) // destructor for data on list +{ + if (list == NULL) { + return; + } + + psDlistElem *ptr = list->head; + while (ptr != NULL) { + psDlistElem *next = ptr->next; + + if (elemDel == NULL) { + psMemDecrRefCounter(ptr->data); + } else { + elemDel(psMemDecrRefCounter(ptr->data)); + } + dlistElemDel(ptr); + + ptr = next; + } + + psFree(list); +} + +/*****************************************************************************/ +/* + * Add an element to a list + */ +psDlist *psDlistAdd(psDlist *list, // list to add to; may be NULL + void *data, // data to add + int where) // where to add data. PS_DLIST_HEAD, + // PS_DLIST_TAIL, or an integer +{ + psDlistElem *elem = dlistElemNew(); + + if (list == NULL) { + list = psDlistNew(NULL); + } + + if (where > list->n) { // XXX need better diagnostic here + fprintf(stderr, "Invalid index %d (only %d elements)\n", + where, list->n); + dlistElemDel(elem); // cleanup + + return list; + } + + if (where == PS_DLIST_HEAD) { // easy + elem->prev = NULL; + elem->next = list->head; + elem->data = psMemIncrRefCounter(data); + + if (list->head != NULL) { + list->head->prev = elem; + } + + list->head = elem; + if (list->tail == NULL) { + list->tail = elem; + } + list->n++; + } else if (where == PS_DLIST_TAIL) { // also easy + elem->prev = list->tail; + elem->next = NULL; + elem->data = psMemIncrRefCounter(data); + + if (list->tail != NULL) { + list->tail->next = elem; + } + + list->tail = elem; + if (list->head == NULL) { + list->head = elem; + } + list->n++; + } else { // XXX + fprintf(stderr, "Insertion into psDlists is not yet supported\n"); + return psDlistAdd(list, data, PS_DLIST_TAIL); + } + + return list; +} + +/*****************************************************************************/ + +psDlist *psDlistAppend(psDlist *list, // list to add to; may be NULL + void *data) // data to add at end +{ + return psDlistAdd(list, data, PS_DLIST_TAIL); +} + +/*****************************************************************************/ +/* + * Remove an element from a list + */ +void *psDlistRemove(psDlist *list, // list to remove element from + void *data, // data item to remove (or NULL) + int which) // index of item, if known. + // PS_DLIST_UNKNOWN, or PS_DLIST_HEAD, + // or PS_DLIST_TAIL +{ + assert (list != NULL); + + if (data != NULL) { + assert (which == PS_DLIST_UNKNOWN); + + fprintf(stderr,"psDlistRemove does not yet support removal by data"); + return NULL; + } + + if (which == PS_DLIST_HEAD) { + if (list->head == NULL) { + return NULL; + } else { + data = list->head->data; + if (list->head->next != NULL) { + list->head->next->prev = NULL; + } + list->head = list->head->next; + list->n--; + + if (list->head == NULL) { + list->tail = NULL; + } + } + } else if (which == PS_DLIST_TAIL) { + if (list->tail == NULL) { + return NULL; + } else { + data = list->tail->data; + if (list->tail->prev != NULL) { + list->tail->prev->next = NULL; + } + list->tail = list->tail->prev; + list->n--; + + if (list->tail == NULL) { + list->head = NULL; + } + } + } else if (which != PS_DLIST_UNKNOWN) { + fprintf(stderr,"psDlistRemove does not yet support removal by index"); + return NULL; + } + + return psMemDecrRefCounter(data); +} + +/*****************************************************************************/ +/* + * Retrieve an element from the list, or iterate over list. + * + * If which is psDlist{Next,Prev} return the next/previous + * data and move the iteration cursor + * + * If which is psDlist{Head,Tail} initialise the iteration pointer + * and return the head/tail + * + * Otherwise, return the desired element + */ +#define ITER_INIT_HEAD ((void *)1) // next iteration should return head +#define ITER_INIT_TAIL ((void *)2) // next iteration should return tail + +void *psDlistGet( + const psDlist *list, // list to retrieve element from + int which) // index of item, or PS_DLIST_NEXT, + // or PS_DLIST_PREV +{ + assert (list != NULL); + + switch (which) { + case PS_DLIST_HEAD: // easy + ((psDlist *)list)->iter = ITER_INIT_HEAD; + return (list->head == NULL) ? NULL : list->head->data; + case PS_DLIST_TAIL: // also easy + ((psDlist *)list)->iter = ITER_INIT_TAIL; + return (list->tail == NULL) ? NULL : list->tail->data; + case PS_DLIST_PREV: // use iterator + case PS_DLIST_NEXT: + if (list->iter == NULL) { + return NULL; + } else if (list->iter == ITER_INIT_HEAD) { + ((psDlist *)list)->iter = list->head; + } else if (list->iter == ITER_INIT_TAIL) { + ((psDlist *)list)->iter = list->tail; + } else if (which == PS_DLIST_PREV) { + ((psDlist *)list)->iter = list->iter->prev; + } else if (which == PS_DLIST_NEXT) { + ((psDlist *)list)->iter = list->iter->next; + } else { + fprintf(stderr,"You cannot get here\n"); + abort(); + } + + if (list->iter == NULL) { + return NULL; + } + return list->iter->data; // XXX what if data is NULL? + + break; + default: // XXX + fprintf(stderr, + "Retrieval from psDlists by index is not yet supported\n"); + return NULL; + } +} + +/*****************************************************************************/ +/* + * Convert a psDlist to/from a psVoidPtrArray + */ +psVoidPtrArray *psDlistToArray(psDlist *restrict dlist) +{ + if (dlist == NULL) { + return NULL; + } + + psVoidPtrArray *restrict arr = psVoidPtrArrayNew(dlist->n, dlist->n); + + psDlistElem *ptr = dlist->head; + for (int i = 0, n = dlist->n; i < n; i++) { + arr->arr[i] = ptr->data; + + ptr->data = NULL; + ptr = ptr->next; + } + + psDlistDel(dlist, NULL); + + return arr; +} + +psDlist *psArrayToDlist(psVoidPtrArray *arr) +{ + psDlist *list = psDlistNew(NULL); // list of elements + + for (int i = 0, n = arr->n; i < n; i++) { + psDlistAppend(list, + psMemDecrRefCounter(arr->arr[i])); // it's already Incr + arr->arr[i] = NULL; + } + + psVoidPtrArrayDel(arr, NULL); + + return list; +} + +/*****************************************************************************/ +#if 0 // testing + +typedef struct { int i; } X; +static X *xNew(int i) { X *x = psAlloc(sizeof(X)); x->i = i; return x; } +static void xDel(X *x) { psFree(x); } + +int main(void) +{ + psDlist *list = psDlistNew(NULL); // list of metadata + + psDlistAppend(list, xNew(1)); + psDlistAppend(list, xNew(2)); + psDlistAppend(list, xNew(3)); + /* + * Print that data + */ + for (psDlistElem *ptr = list->head; ptr != NULL; ptr = ptr->next) { + printf("%d ", ((X *)(ptr->data))->i); + } + printf("\n"); + /* + * Print it again + */ + (void)psDlistGet(list, PS_DLIST_TAIL); // initialise iterator + X *x = NULL; + while ((x = psDlistGet(list, PS_DLIST_PREV)) != NULL) { + printf("%d ", x->i); + } + printf("\n"); +#if 1 + /* + * Convert to an array + */ + psVoidPtrArray *arr = psDlistToArray(list); + list = NULL; // it's been freed + + for (int i = 0; i < arr->n; i++) { + printf("%d ", ((X *)(arr->arr[i]))->i); + } + printf("\n"); +#if 1 + /* + * Convert back to a list + */ + list = psArrayToDlist(arr); + arr = NULL; // it's been freed + + (void)psDlistGet(list, PS_DLIST_TAIL); // initialise iterator + while ((x = psDlistGet(list, PS_DLIST_PREV)) != NULL) { + printf("%d ", x->i); + } + printf("\n"); +#endif + psVoidPtrArrayDel(arr, (void (*)(void *))xDel); +#endif + +#if 1 + list = psArrayToDlist(psDlistToArray(list)); + + (void)psDlistGet(list, PS_DLIST_TAIL); // initialise iterator + while ((x = psDlistGet(list, PS_DLIST_PREV)) != NULL) { + printf("%d ", x->i); + } + printf("\n"); +#endif + /* + * Clean up + */ + psDlistDel(list, (void (*)(void *))xDel); + /* + * Check for memory leaks + */ + fprintf(stderr,"Checking for memory leaks:\n"); + psMemCheckLeaks(0, NULL, stderr); + + return 0; +} +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/hash.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/hash.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/hash.c (revision 21655) @@ -0,0 +1,245 @@ +/* + * An interface to hash tables for Pan-STARRS + * + * Collisions are handled by simple linked lists + */ +#include +#include +#include +#include "psUtils.h" + +/*****************************************************************************/ +/* + * A bucket that holds an item of data + */ +typedef struct HashBucket { + char *key; // key for this item of data + void *data; // the data itself + struct HashBucket *next; // list of other possible keys +} HashBucket; + +/* + * An entire hash table + */ +struct HashTable { + int nbucket; // number of buckets + HashBucket **buckets; // the buckets themselves +}; + +/*****************************************************************************/ +/* + * Con/Destruct buckets + */ +static HashBucket *hashBucketNew(const char *key, + void *data, + HashBucket *next) +{ + HashBucket *bucket = psAlloc(sizeof(HashBucket *)); + + bucket->key = psStringCopy(key); + bucket->data = psMemIncrRefCounter(data); + bucket->next = next; + + return bucket; +} + +static void hashBucketDel( + HashBucket *bucket, // bucket to free + void (*itemDel)(void *item)) // how to free hashed data; + // or NULL +{ + if (bucket == NULL) { + return; + } + + psFree(bucket->key); + psMemDecrRefCounter(bucket->data); + + if (itemDel != NULL) { + itemDel(bucket->data); + } + + psFree(bucket); +} + +/*****************************************************************************/ +/* + * Con/Destruct psHash tables + */ +psHash *psHashNew(int nbucket) // initial number of buckets +{ + psHash *table = psAlloc(sizeof(psHash)); + table->buckets = psAlloc(nbucket*sizeof(HashBucket *)); + table->nbucket = nbucket; + + psTrace("utils.hash", 1, "Creating %d-element hash table\n", nbucket); + + for (int i = 0; i < nbucket; i++) { + table->buckets[i] = NULL; + } + + return table; +} + +void psHashDel(psHash *table, // hash table to be freed + void (*itemDel)(void *item)) // how to free hashed data; or NULL +{ + if (table == NULL) { + return; + } +/* + * Release data interned on the list, and delete the buckets + */ + for (int i = 0; i < table->nbucket; i++) { + if (table->buckets[i] != NULL) { + HashBucket *ptr = table->buckets[i]; + while (ptr != NULL) { + HashBucket *tmp = ptr->next; + hashBucketDel(ptr, itemDel); + ptr = tmp; + } + } + } + + psFree(table->buckets); + psFree(table); +} + +/*****************************************************************************/ +/* + * Here's the routine to do the work of insertion/retrieval. + * + * We need itemDel if the key is already in use and we're inserting + * + * N.b. this is NOT a good hash function! See Knuth for some better ones + */ +static void *doHashWork( + psHash *table, // table to insert in + const char *key, // key to use + void *data, // data to insert, + // or (if NULL) retrieve + void (*itemDel)(void *item)) // how to free hashed data + // or NULL +{ + long int hash = 1; + + for (int i = 0, len = strlen(key); i < len; i++) { + hash = (hash << 1) ^ key[i]; + } + + hash &= (table->nbucket - 1); + + HashBucket *ptr = table->buckets[hash]; + /* + * We've found the correct hash bucket, now we need to know what to do + */ + if (data == NULL) { // retrieve + while (ptr != NULL) { + if (strcmp(key, ptr->key) == 0) { // found it! + return ptr->data; + } + + ptr = ptr->next; + } + + return NULL; // not in hash + } else { // insert + while (ptr != NULL) { + if (strcmp(key, ptr->key) == 0) { // found it! + if (itemDel == NULL) { + return NULL; // can't insert + } + + psTrace("utils.hash.insert", 3, "Replacing data for %s\n", + key); + + itemDel(psMemDecrRefCounter(ptr->data)); + ptr->data = psMemIncrRefCounter(data); + + return data; + } + } + /* + * Not found, so insert at the front of the list + */ + table->buckets[hash] = hashBucketNew(key, data, table->buckets[hash]); + + return data; + } +} + +/*****************************************************************************/ +/* + * Insert a value into a hash table + */ +void *psHashInsert(psHash *table, // table to insert in + const char *key, // key to use + void *data, // data to insert + void (*itemDel)(void *item)) // how to free hashed data; + // or NULL +{ + return doHashWork(table, key, data, itemDel); +} + +/*****************************************************************************/ +/* + * Lookup a value in a hash table + */ +void *psHashLookup(psHash *table, // table to lookup key in + const char *key) // key to lookup +{ + return doHashWork(table, key, NULL, NULL); +} + +/*****************************************************************************/ +#if 0 // testing + +typedef struct { char *name; } ID; +static ID *IdNew(const char *name) +{ + ID *id = psAlloc(sizeof(ID)); + id->name = psStringCopy(name); + + return id; +} + +static void IdDel(ID *id) { + if (psMemGetRefCounter(id) > 1) { + psMemDecrRefCounter(id); + return; + } + + psFree(id->name); + psFree(id); +} + +int main(void) +{ + psSetTraceLevel("utils.hash", 3); + long memId0 = psMemGetId(); + + psHash *hash = psHashNew(16); + + psHashInsert(hash, "Lynda", IdNew("Lee"), (void (*)(void *))IdDel); + psHashInsert(hash, "Lynda", IdNew("Lee"), (void (*)(void *))IdDel); + psHashInsert(hash, "Robert", IdNew("Lupton"), (void (*)(void *))IdDel); + psHashInsert(hash, "Robert", IdNew("the Good"), (void (*)(void *))IdDel); + + char *names[] = {"Robert", "Lynda", "Rowan", NULL}; + for (char **name = names; *name != NULL; name++) { + ID *id = psHashLookup(hash, *name); + printf("%s's surname is:\t%s\n", *name, + ((id == NULL) ? "(unknown)" : id->name)); + } + + psHashDel(hash, (void (*)(void *))IdDel); + /* + * Check for memory leaks + */ + fprintf(stderr,"Checking for memory leaks:\n"); + psMemCheckLeaks(memId0, NULL, stderr); + + return 0; +} + +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/logmsg.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/logmsg.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/logmsg.c (revision 21655) @@ -0,0 +1,168 @@ +/* + * A simple implementation of a logging facility for Pan-STARRS + */ +#include +#include +#include +#include +#include +#include +#include +#include "psUtils.h" + +/*****************************************************************************/ + +static int logDest = PS_LOG_TO_STDERR; // where to log messages +static int logLevel = PS_LOG_INFO; // log all messages at this or above + +/*****************************************************************************/ +/* + * Set the current log level and return old level + */ +int psSetLogLevel(int level) // new level +{ + int oldLevel = logLevel; + + if (level < 0 || level > 9) { + psLogMsg("logmsg", PS_LOG_WARN, + "Attempt to set invalid logMsg level: %d", level); + level = (level < 0) ? 0 : 9; + } + + logLevel = level; + + return oldLevel; +} + +/*****************************************************************************/ + +int psSetLogDestination(int dest) // new destination +{ + int old = logDest; + + switch (dest) { + case PS_LOG_TO_STDOUT: + case PS_LOG_TO_STDERR: + logDest = dest; + break; + + default: + fprintf(stderr,"Unknown logDestination: %d (ignored)\n", dest); + break; + } + + return old; +} + +/*****************************************************************************/ +#if !defined(HOST_NAME_MAX) // should be in limits.h +# define HOST_NAME_MAX 256 +#endif + +void p_psVLogMsg(const char *name, int level, const char *fmt, va_list ap) +{ + static int first = 1; + static char hostname[HOST_NAME_MAX + 1]; + + if (level > logLevel) { + return; + } + + if (first) { + first = 0; + gethostname(hostname, HOST_NAME_MAX); + } + + char clevel; // letter-name for level + switch (level) { + case PS_LOG_ABORT: + clevel = 'A'; + break; + + case PS_LOG_ERROR: + clevel = 'E'; + break; + + case PS_LOG_WARN: + clevel = 'W'; + break; + + case PS_LOG_INFO: + clevel = 'I'; + break; + + case 4: + case 5: + case 6: + case 7: + case 8: + case 9: + clevel = level + '0'; + break; + + default: + psTrace("utils.logMsg", 2, "Invalid logMsg level: %d (%s)\n", + level, fmt); + level = (level < 0) ? 0 : 9; + break; + } + + char head[HOST_NAME_MAX + 40]; // yes, this is long enough + time_t clock = time(NULL); + struct tm *utc = gmtime(&clock); + + sprintf(head, "%4d:%02d:%02d %02d:%02d:%02dZ|%-20s|%c|%-15s|", + utc->tm_year + 1900, utc->tm_mon + 1, utc->tm_mday, + utc->tm_hour, utc->tm_min, utc->tm_sec, + hostname, clevel, name); + + switch (logDest) { + case PS_LOG_TO_STDOUT: + puts(head); + vprintf(fmt, ap); + if (fmt[strlen(fmt) - 1] != '\n') { + putc('\n', stdout); + } + break; + + case PS_LOG_TO_STDERR: + fputs(head, stderr); + vfprintf(stderr, fmt, ap); + if (fmt[strlen(fmt) - 1] != '\n') { + putc('\n', stderr); + } + break; + + default: + fprintf(stderr, "psLogMsg: You cannot get here (%s:%d)\n", + __FILE__, __LINE__); + abort(); + } + + va_end(ap); +} + +void psLogMsg(const char *name, int level, const char *fmt, ...) +{ + va_list ap; + va_start(ap, fmt); + + p_psVLogMsg(name, level, fmt, ap); +} + +/*****************************************************************************/ +#if 0 // testing + +int main(void) +{ + psSetTraceLevel("utils.logMsg", 10); + + psSetLogDestination(PS_LOG_TO_STDERR); + psSetLogLevel(PS_LOG_INFO); + + psLogMsg("aaa", PS_LOG_INFO, "Hello World\n"); + psAbort("bbb", "I'm in Hawai`i"); + + return 0; +} +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/memory.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/memory.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/memory.c (revision 21655) @@ -0,0 +1,572 @@ +#include +#include +#include +#include "psUtils.h" +#include "Private/p_psMemory.h" + +#undef psAlloc // don't get the __FILE__ versions +#undef psRealloc +#undef psFree + +static int bad_memblock(const psMemBlock *m, int quiet); +/* + * This strawman allocation package only handles up to NMEMBLOCK allocations + */ +#define NMEMBLOCK 10000 // max. number of calls to psAlloc +static psMemBlock *memBlocks[NMEMBLOCK]; + +/*****************************************************************************/ +/* + * Unique ID for allocated blocks + */ +static unsigned long memid = 0; + +/*****************************************************************************/ + +void *psAlloc(size_t size) +{ + return p_psAlloc(size, "(unknown)", 0); +} + +void *psRealloc(void *ptr, size_t size) +{ + return p_psRealloc(ptr, size, "(unknown)", 0); +} + +void psFree(void *ptr) +{ + p_psFree(ptr, "(unknown)", 0); +} + +/*****************************************************************************/ +/* + * Callbacks + */ +/* + * First the I-can't-get-the-memory callback + */ +static void *memExhaustedCB0(size_t size) +{ + return NULL; +} +static psMemExhaustedCallback memExhaustedCB = memExhaustedCB0; + +psMemExhaustedCallback psMemExhaustedSetCB(psMemExhaustedCallback func) +{ + psMemExhaustedCallback old = memExhaustedCB; + + memExhaustedCB = (func != NULL) ? func : memExhaustedCB0; + + return old; +} +/* + * then the I-have-detected-a-problem callback + */ +static void memProblemCB0(const psMemBlock *ptr, + const char *file, int lineno) +{ + if (bad_memblock(ptr, 0)) { + fprintf(stderr, "Memory corruption or an attempt to manipulate memory " + "not allocated with psAlloc at %s", file); + if (lineno > 0) { + fprintf(stderr, ":%d", lineno); + } + fprintf(stderr, "\n"); + } else if(ptr->refCounter <= 0) { + psError(__func__, "Block %ld allocated at %s:%d freed more than once at %s:%d\n", + ptr->id, ptr->file, ptr->lineno, file, lineno); + } else if(ptr->refCounter > 1) { + psError(__func__, "Block %ld allocated at %s:%d freed while still referenced at %s:%d\n", + ptr->id, ptr->file, ptr->lineno, file, lineno); + } + + if (lineno > 0) { + psAbort(__func__, "Detected memory corruption"); + } +} +static psMemProblemCallback memProblemCB = memProblemCB0; + +psMemProblemCallback psMemProblemSetCB(psMemProblemCallback func) +{ + psMemProblemCallback old = memProblemCB; + + memProblemCB = (func != NULL) ? func : memProblemCB0; + + return old; +} +/* + * And now the I-want-to-be-informed callbacks + * + * Call the callbacks when these IDs are allocated/freed + */ +long p_psMemAllocateID = 0; // notify user this block is allocated +long p_psMemFreeID = 0; // notify user this block is freed + +long psMemSetAllocateID(long id) // set p_psMemAllocateID to id +{ + long old = p_psMemAllocateID; + p_psMemAllocateID = id; + + return old; +} + +long psMemSetFreeID(long id) // set p_psMemFreeID to id +{ + long old = p_psMemFreeID; + p_psMemFreeID = id; + + return old; +} + +/* + * Default callback for both allocate and free. Note that the + * value of p_psMemAllocateID/p_psMemFreeID is incremented + * by the return value (so returning 0 means that the callback + * isn't resignalled) + */ +static int memAllocateCB0(const psMemBlock *ptr) +{ + static int incr = 0; // "p_psMemAllocateID += incr" + + assert (ptr != NULL); // prevent compiler warnings + + return incr; +} + +/* + * The default callbacks, and the routines to change them + */ +static psMemCallback memAllocateCB = memAllocateCB0; +static psMemCallback memFreeCB = memAllocateCB0; + +psMemCallback psMemAllocateSetCB(psMemCallback func) +{ + psMemCallback old = memAllocateCB; + + memAllocateCB = (func != NULL) ? func : memAllocateCB0; + + return old; +} + +psMemCallback psMemFreeSetCB(psMemCallback func) +{ + psMemCallback old = memFreeCB; + + memFreeCB = (func != NULL) ? func : memAllocateCB0; + + return old; +} + +/*****************************************************************************/ +/* + * Return memory ID counter for next block to be allocated + */ +int psMemGetId(void) +{ + return memid + 1; +} + +/*****************************************************************************/ +/* + * Routines to check the consistency of the allocated and/or free memory arena + * + * N.b. If the block wasn't allocated by psAlloc, it will appear corrupted + */ +#define ALIGNED(P) ((void *)((long)(P) & ~03) == (P)) + +static int +bad_memblock(const psMemBlock *m, // block to check + int quiet) // be quiet? +{ + if(m == NULL) { + if (!quiet) { + psError(__func__, "psMemCheckCorruption: NULL memory block\n"); + } + return(1); + } + + if(!ALIGNED(m)) { + if (!quiet) { + psError(__func__, "psMemCheckCorruption: non-aligned memory block\n"); + } + return(1); + } + + if(m->magic0 != P_PS_MEMMAGIC || m->magic != P_PS_MEMMAGIC) { + if (!quiet) { + psError(__func__, + "psMemCheckCorruption: memory block %ld is corrupted\n", m->id); + } + return(1); + } + + return(0); +} + +int psMemCheckCorruption(int abort_on_error) +{ + int nbad = 0; // number of bad blocks + + for (int i = 1; i <= memid; i++) { + if (memBlocks[i] != NULL) { + if (bad_memblock(memBlocks[i], 1)) { + nbad++; + + memProblemCB(memBlocks[i], "(psMemCheckCorruption)", 0); + + if (abort_on_error) { + psAbort(__func__, "Detected memory corruption"); + } + } + } + } + + return nbad; +} + +/*****************************************************************************/ + +void *p_psAlloc(size_t size, const char *file, int lineno) +{ + psMemBlock *ptr = malloc(sizeof(psMemBlock) + size); + + if (ptr == NULL) { + ptr = memExhaustedCB(size); + if (ptr == NULL) { + psAbort(__func__, "Failed to allocate %ld bytes at %s:%d", + size, file, lineno); + } + } + + *(unsigned long *)&ptr->id = ++memid; + ptr->file = file; + *(int *)&ptr->lineno = lineno; + *(void **)&ptr->magic0 = *(void **)&ptr->magic = P_PS_MEMMAGIC; + + ptr->refCounter = 1; // one user so far + /* + * Did the user ask to be informed about this allocation? + */ + if(ptr->id == p_psMemAllocateID) { + p_psMemAllocateID += memAllocateCB(ptr); + } + /* + * Register the allocation. This is not production quality! + */ + if (memid == NMEMBLOCK) { + psAbort(__func__, "I'm sorry, but I can only handle %d allocations", + NMEMBLOCK); + } + memBlocks[ptr->id] = ptr; + /* + * And return the user the memory that they allocated + */ + return ptr + 1; // usr memory +} + +void *p_psRealloc(void *vptr, size_t size, const char *file, int lineno) +{ + if (vptr == NULL) { + return p_psAlloc(size, file, lineno); + } else { + psMemBlock *ptr = ((psMemBlock *)vptr) - 1; + ptr = (psMemBlock *)realloc(ptr, sizeof(psMemBlock) + size); + memBlocks[ptr->id] = ptr; + + return ptr + 1; // usr memory + } +} + +void p_psFree(void *vptr, const char *file, int lineno) +{ + if (vptr == NULL) { + return; + } + + psMemBlock *ptr = ((psMemBlock *)vptr) - 1; + + if (bad_memblock(ptr, 0) || ptr->refCounter != 1) { + memProblemCB(ptr, file, lineno); + } + /* + * Did the user ask to be informed about this allocation? + */ + if(ptr->id == p_psMemFreeID) { + p_psMemFreeID += memFreeCB(ptr); + } + + ptr->refCounter--; + memBlocks[ptr->id] = NULL; // XXX sub-optimal! Can't check + // freed blocks for corruption + free(ptr); +} + +/*****************************************************************************/ +/* + * Check for memory leaks. Not production quality code + */ +int psMemCheckLeaks( + int id0, // don't list blocks with id < id0 + psMemBlock ***arr, // pointer to array of pointers to + // leaked blocks, or NULL + FILE *fd) // print list of leaks to fd (or NULL) +{ + int nleak = 0; + + for (int i = id0; i <= memid; i++) { + if (memBlocks[i] != NULL) { + if (memBlocks[i]->refCounter > 0) { + nleak++; + + if (fd != NULL) { + if (nleak == 1) { + fprintf(fd, " %20s:line ID\n", "file"); + } + + fprintf(fd, " %20s:%-4d %ld\n", + memBlocks[i]->file, memBlocks[i]->lineno, + memBlocks[i]->id); + } + } + } + } + + if (nleak == 0 || arr == NULL) { + return nleak; + } + + assert (*arr == NULL); // Don't generate a memory leak! + + *arr = p_psAlloc(nleak*sizeof(psMemBlock), __FILE__, __LINE__); + for (int i = id0, j = 0; i < memid; i++) { + if (memBlocks[i] != NULL) { + if (memBlocks[i]->refCounter > 0) { + (*arr)[j++] = memBlocks[i]; + if (j == nleak) { // found them all + break; + } + } + } + } + + return nleak; +} + +/*****************************************************************************/ +/* + * Reference counting APIs + */ +int psMemGetRefCounter(void *vptr) // return refCounter +{ + psMemBlock *ptr = ((psMemBlock *)vptr) - 1; + + if (bad_memblock(ptr, 0)) { + memProblemCB(ptr, "(psMemGetRefCounter)", -1); + } + + return ptr->refCounter; +} + +void *psMemIncrRefCounter(void *vptr) // increment and return refCounter +{ + if (vptr == NULL) { + return vptr; + } + + psMemBlock *ptr = ((psMemBlock *)vptr) - 1; + + if (bad_memblock(ptr, 0)) { + memProblemCB(ptr, "(psMemIncrRefCounter)", -1); + } + + ptr->refCounter++; + + return vptr; +} + +void *psMemDecrRefCounter(void *vptr) // decrement and return refCounter +{ + if (vptr == NULL) { + return vptr; + } + + psMemBlock *ptr = ((psMemBlock *)vptr) - 1; + + if (bad_memblock(ptr, 0)) { + memProblemCB(ptr, "(psMemDecrRefCounter)", -1); + } + + ptr->refCounter--; + + return vptr; +} + +/*****************************************************************************/ +#if 0 // Testing +#include +/* + * These are usually defined in psMemory.h, but in this file (only!) + * they are undefined explicitly, so let's get them back + */ +#define psAlloc(size) p_psAlloc(size, __FILE__, __LINE__) +#define psRealloc(ptr, size) p_psRealloc(ptr, size, __FILE__, __LINE__) +#define psFree(size) p_psFree(size, __FILE__, __LINE__) + +/*****************************************************************************/ +/* + * My callbacks + */ +static int my_MemAllocateCB(const psMemBlock *ptr) +{ + static int incr = 0; // "p_psMemAllocateID += incr" + + assert (ptr != NULL); // prevent compiler warnings + + fprintf(stderr, "Allocating block %ld (%s:%d)\n", + ptr->id, ptr->file, ptr->lineno); + + incr = 2; // See every other allocation + + return incr; +} + +static void *my_MemExhaustedCB(size_t size) +{ + psError(__func__, "Bye bye sweet world\n"); + + return NULL; +} + +static psMemProblemCallback default_MemProblemCB = NULL; + +static void my_MemProblemCB(const psMemBlock *ptr, + const char *file, int lineno) +{ + fprintf(stderr, "I foresee trouble:\n"); + + default_MemProblemCB(ptr, file, lineno); +} + +/*****************************************************************************/ + +int main(void) +{ + /* + * Survive any aborts that are generated + * This doesn't seem to work on mac os/x 10.2.8 + */ + if (signal(SIGABRT, SIG_IGN) == SIG_ERR) { + perror("Installing a SIGABRT handler"); + } + /* + * Install my callbacks + */ + (void)psMemAllocateSetCB(my_MemAllocateCB); + default_MemProblemCB = psMemProblemSetCB(my_MemProblemCB); + (void)psMemExhaustedSetCB(my_MemExhaustedCB); + /* + * Request callback on first allocation + */ + (void)psMemSetAllocateID(1); + /* + * Start allocating + */ + char *str = psAlloc(10); + double *x = psAlloc(1); + int *i = psAlloc(1); + *i = 1; + *x = 1.234; + + psFree(i); + + str = psRealloc(str, 20); + psFree(str); +#if 1 + /* + * An illegal double-free + */ + psFree(str); +#endif + +#if 1 + /* + * Free memory not allocated with psAlloc + */ + psFree(malloc(1)); +#endif + +#if 0 + /* + * Ask for more memory than we can get + */ + psAlloc((size_t)-1 - sizeof(psMemBlock)); +#endif + +#if 0 + /* + * Increment refCounter on non-malloced pointer + */ + psMemGetRefCounter(malloc(1)); +#endif + +#if 0 + /* + * Mess with the refCounter + */ + (void)psMemIncrRefCounter(x); + (void)psMemIncrRefCounter(x); + (void)psMemDecrRefCounter(x); // now == 2 +#endif + +#if 0 + /* + * A memory leak + */ + { + char *str2 = psAlloc(20); + + sprintf(str2, "XXX"); + } +#else + psFree(x); +#endif + +#if 1 + /* + * Corrupt some memory + */ + char *c = psAlloc(1); + c[-4] = 0; +#endif + + /*************************************************************************/ + /* + * Check for memory corruption + */ + if(psMemCheckCorruption(1)) { + psError(__func__, "Detected corrupted memory\n\n"); + } + /* + * Check for memory leaks + */ + psMemBlock **leaks = NULL; + int nleak = psMemCheckLeaks(0, &leaks, NULL); + + if (nleak > 0) { + fprintf(stderr, "Memory is leaking (%d blocks):\n", nleak); + fprintf(stderr, " %20s:line ID\n", "file"); + for (int i = 0; i < nleak; i++) { + fprintf(stderr, " %20s:%-4d %ld\n", + leaks[i]->file, leaks[i]->lineno, leaks[i]->id); + } + psFree(leaks); + } + + if ((nleak = psMemCheckLeaks(0, NULL, NULL)) > 0) { + fprintf(stderr, "Memory is still leaking (%d blocks):\n", nleak); + psMemCheckLeaks(0, NULL, stderr); + } + + return 0; +} +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/misc.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/misc.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/misc.c (revision 21655) @@ -0,0 +1,41 @@ +/* + * A set of miscellaneous utilities for Pan-STARRS + */ +#include +#include +#include +#include "psUtils.h" + +/*****************************************************************************/ +/* + * Return a copy of a string + */ +char *psStringCopy(const char *str) +{ + return strcpy(psAlloc(strlen(str) + 1), str); +} + +/*****************************************************************************/ +/* + * Send a psLogFatal message and abort + */ +void psAbort(const char *name, const char *fmt, ...) +{ + va_list ap; + va_start(ap, fmt); + + p_psVLogMsg(name, PS_LOG_ABORT, fmt, ap); + abort(); +} + +/*****************************************************************************/ +/* + * Print an error message + */ +void psError(const char *name, const char *fmt, ...) +{ + va_list ap; + va_start(ap, fmt); + + p_psVLogMsg(name, PS_LOG_ERROR, fmt, ap); +} Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psArray.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psArray.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psArray.h (revision 21655) @@ -0,0 +1,137 @@ +#if !defined(PS_ARRAY_H) +#define PS_ARRAY_H + +#include +/* + * Declare TYPEArray + */ +#define PS_DECLARE_ARRAY_TYPE(TYPE) \ + typedef struct { \ + int size; /* number of allocated elements */ \ + int n; /* number of elements in use */ \ + \ + TYPE *arr; /* the actual data */ \ + } PS_CONCAT(TYPE, Array); \ + \ +PS_CONCAT(TYPE, Array) *PS_CONCAT3(TYPE,Array,New)(int n, int s); \ +PS_CONCAT(TYPE, Array) *PS_CONCAT3(TYPE,Array,Realloc)(PS_CONCAT(TYPE, Array) *arr, int n); \ +void PS_CONCAT3(TYPE,Array,Del)(PS_CONCAT(TYPE, Array) *arr) + +/*****************************************************************************/ +/* + * Generate the code TYPEArray's constructors/destructors + */ +#define PS_CREATE_ARRAY_TYPE(TYPE) \ + P_PS_CREATE_ARRAY_TYPE(,,TYPE) + +#define P_PS_CREATE_ARRAY_TYPE(QUAL,PREFIX,TYPE) \ +\ +QUAL PS_CONCAT(TYPE, Array) *PS_CONCAT4(PREFIX,TYPE,Array,New)(int n, /* initial dimension */\ + int s) /* initial number of elements */ \ +{ \ + PS_CONCAT(TYPE, Array) *arr = psAlloc(sizeof(PS_CONCAT(TYPE, Array))); \ + \ + if (n > s) { \ + psLogMsg(PS_STRING(TYPE) "ArrayNew", PS_LOG_ERROR, "Too few elements requested in array (%d < %d)\n", s, n); \ + s = n; \ + } \ + \ + arr->size = s; \ + arr->n = n; \ + \ + if (s == 0) { \ + arr->arr = NULL; \ + } else { \ + arr->arr = psAlloc(s*sizeof(TYPE)); \ + } \ + \ + return arr; \ +} \ +\ +QUAL PS_CONCAT(TYPE, Array) *PS_CONCAT4(PREFIX,TYPE,Array,Realloc)(PS_CONCAT(TYPE, Array) *arr, int n) \ +{ \ + if (arr == NULL) { \ + return PS_CONCAT4(PREFIX,TYPE,Array,New)(n, n); \ + } \ + \ + if (n <= arr->size) { \ + if (arr->n < n) { \ + arr->n = n; \ + } \ + } else { \ + arr->arr = psRealloc(arr->arr, n*sizeof(TYPE)); \ + arr->size = n; \ + } \ + \ + return arr; \ +} \ +\ +QUAL void PS_CONCAT4(PREFIX,TYPE,Array,Del)(PS_CONCAT(TYPE,Array) *arr) \ +{ \ + if (arr == NULL) { \ + return; \ + } \ + \ + psFree(arr->arr); \ + psFree(arr); \ +} + +/*****************************************************************************/ +/* + * Support for pointer types + */ +#define PS_DECLARE_ARRAY_PTR_TYPE(TYPE) \ + typedef TYPE *P_PS_CONCAT(TYPE, Ptr); \ + PS_DECLARE_ARRAY_TYPE(PS_CONCAT(TYPE, Ptr)) + +#define PS_CREATE_ARRAY_PTR_TYPE(TYPE); \ + P_PS_CREATE_ARRAY_TYPE(static, my_, P_PS_CONCAT(TYPE, Ptr)) \ + \ + PS_CONCAT(TYPE, PtrArray) *PS_CONCAT3(TYPE,PtrArray,New)(int n, int s) \ + { \ + PS_CONCAT(TYPE, PtrArray) *arr = PS_CONCAT4(my_,TYPE,PtrArray,New)(n, s); \ + for (int i = 0; i < arr->n; i++) { \ + arr->arr[i] = psMemIncrRefCounter,(PS_CONCAT(TYPE, New)());\ + } \ + \ + return arr; \ + } \ + \ + PS_CONCAT(TYPE, PtrArray) *PS_CONCAT3(TYPE,PtrArray,Realloc)(PS_CONCAT(TYPE, PtrArray) *arr, int n) \ + { \ + for (int i = n; i < arr->n; i++) { \ + PS_CONCAT(TYPE, Del)(psMemDecrRefCounter(arr->arr[i]));\ + } \ + \ + return PS_CONCAT4(my_,TYPE,PtrArray,Realloc)(arr, n); \ + } \ + \ + void PS_CONCAT3(TYPE,PtrArray,Del)(PS_CONCAT2(TYPE,PtrArray) *arr) \ + { \ + if (arr == NULL) { return; } \ + for (int i = 0; i < arr->size; i++) { \ + PS_CONCAT(TYPE, Del)(psMemDecrRefCounter(arr->arr[i])); \ + } \ + PS_CONCAT4(my_, TYPE, PtrArray, Del)(arr); \ + } + +/*****************************************************************************/ +/* + * Declare some common types of arrays + * + * psVoidPtrArray is special, as it needs to have a destructor that + * accepts a destructor for the elements, so we cannot simply use + * PS_DECLARE_ARRAY_TYPE(psVoidPtr); + */ +typedef struct { + int n, size; + void **arr; +} psVoidPtrArray; + +psVoidPtrArray *psVoidPtrArrayNew(int n, int s); +psVoidPtrArray *psVoidPtrArrayRealloc(psVoidPtrArray *arr, int n); +void psVoidPtrArrayDel(psVoidPtrArray *arr, + void (*elemDel)(void *)); // destructor for array data + + +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psDlist.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psDlist.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psDlist.h (revision 21655) @@ -0,0 +1,56 @@ +#if !defined(PS_DLIST_H) +#define PS_DLIST_H +/* + * Support for doubly linked lists + */ +typedef struct psDlistElem { + struct psDlistElem *prev; // previous link in list + struct psDlistElem *next; // next link in list + void *data; // real data item +} psDlistElem; + +typedef struct { + int n; // number of elements on list + psDlistElem *head; // first element on list (may be NULL) + psDlistElem *tail; // last element on list (may be NULL) + psDlistElem *iter; // iteration cursor +} psDlist; + +enum { // Special values of index into list + PS_DLIST_HEAD = 0, // at head + PS_DLIST_TAIL = -1, // at tail + PS_DLIST_UNKNOWN = -2, // unknown position + PS_DLIST_PREV = -3, // previous element + PS_DLIST_NEXT = -4 // next element +}; + +/*****************************************************************************/ +/* + * Constructors and Destructors + */ +psDlist *psDlistNew(void *data); // initial data item; may be NULL + +void psDlistDel(psDlist *list, // list to destroy + void (*elemDel)(void *)); // destructor for data on list + +/* + * List maintainence functions + */ +psDlist *psDlistAdd(psDlist *list, // list to add to (may be NULL) + void *data, // data item to add + int where); // index, PS_DLIST_HEAD, or PS_DLIST_TAIL +psDlist *psDlistAppend(psDlist *list, // list to append to (may be NULL) + void *data); // data item to add +void *psDlistRemove(psDlist *list, // list to remove element from + void *data, // data item to remove + int which); // index of item, or PS_DLIST_UNKNOWN, + // or PS_DLIST_NEXT, or PS_DLIST_PREV +void *psDlistGet(const psDlist *list, // list to retrieve element from + int which); // index of item, or PS_DLIST_NEXT, + // or PS_DLIST_PREV +/* + * Conversions to/from arrays + */ +psVoidPtrArray *psDlistToArray(psDlist *dlist); +psDlist *psArrayToDlist(psVoidPtrArray *arr); +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psHash.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psHash.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psHash.h (revision 21655) @@ -0,0 +1,19 @@ +#if !defined(PS_HASH_H) +#define PS_HASH_H + +typedef struct HashTable psHash; + +psHash *psHashNew(int nbucket); // initial number of buckets +void psHashDel(psHash *table, // hash table to be freed + void (*itemDel)(void *item)); // how to free hashed data; + // or NULL + +void *psHashInsert(psHash *table, // table to insert in + const char *key, // key to use + void *data, // data to insert + void (*itemDel)(void *item)); // how to free hashed data; + // or NULL +void *psHashLookup(psHash *table, // table to lookup key in + const char *key); // key to lookup + +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psLogMsg.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psLogMsg.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psLogMsg.h (revision 21655) @@ -0,0 +1,14 @@ +#if !defined(PS_LOG_MSG_H) +#define PS_LOG_MSG_H +#include + +enum { PS_LOG_ABORT = 0, PS_LOG_ERROR, PS_LOG_WARN, PS_LOG_INFO }; + +enum { PS_LOG_TO_STDERR, PS_LOG_TO_STDOUT }; + +int psSetLogDestination(int dest); +int psSetLogLevel(int level); +void psLogMsg(const char *name, int level, const char *fmt, ...); +void p_psVLogMsg(const char *name, int level, const char *fmt, va_list ap); + +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psMemory.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psMemory.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psMemory.h (revision 21655) @@ -0,0 +1,71 @@ +#if !defined(PS_MEMORY_H) +#define PS_MEMORY_H +#include +/* + * Book-keeping data for storage allocator + * + * N.b. sizeof(psMemBlock) must be chosen such that if ptr is a pointer + * returned by malloc, then ((char *)ptr + sizeof(psMemBlock)) is properly + * aligned for all storage types. + */ +typedef struct { + const void *magic0; // initialised to p_psMEMMAGIC + const unsigned long id; // a unique ID for this allocation + const char *file; // set from __FILE__ in e.g. p_psAlloc + const int lineno; // set from __LINE__ in e.g. p_psAlloc + int refCounter; // how many times pointer is referenced + const void *magic; // initialised to p_psMEMMAGIC +} psMemBlock; + +typedef int (*psMemCallback)(const psMemBlock *ptr); +typedef void (*psMemProblemCallback)(const psMemBlock *ptr, + const char *file, int lineno); +typedef void *(*psMemExhaustedCallback)(size_t size); + +/*****************************************************************************/ + +void *psAlloc(size_t size); +void *psRealloc(void *ptr, size_t size); +void psFree(void *ptr); + +void *p_psAlloc(size_t size, const char *file, int lineno); +void *p_psRealloc(void *ptr, size_t size, const char *file, int lineno); +void p_psFree(void *ptr, const char *file, int lineno); + +#define psAlloc(size) p_psAlloc(size, __FILE__, __LINE__) +#define psRealloc(ptr, size) p_psRealloc(ptr, size, __FILE__, __LINE__) +#define psFree(size) p_psFree(size, __FILE__, __LINE__) + +/*****************************************************************************/ +/* + * Checks of memory system + */ +int psMemCheckLeaks( + int id0, // don't list blocks with id < id0 + psMemBlock ***arr, // pointer to array of pointers to + // leaked blocks, or NULL + FILE *fd); // print list of leaks to fd (or NULL) +int psMemCheckCorruption(int abort_on_error); + +/*****************************************************************************/ +/* + * Manipulate reference counter + */ +int psMemGetRefCounter(void *vptr); // return refCounter +void *psMemIncrRefCounter(void *vptr); // increment refCounter and return vptr +void *psMemDecrRefCounter(void *vptr); // decrement refCounter and return vptr + +/*****************************************************************************/ +/* + * Functions to set and control callbacks + */ +psMemProblemCallback psMemProblemSetCB(psMemProblemCallback func); +psMemExhaustedCallback psMemExhaustedSetCB(psMemExhaustedCallback func); + +psMemCallback psMemAllocateSetCB(psMemCallback func); +psMemCallback psMemFreeSetCB(psMemCallback func); + +int psMemGetId(void); // get next memory ID +long psMemSetAllocateID(long id); // set p_psMemAllocateID to id +long psMemSetFreeID(long id); // set p_psMemFreeID to id +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psMisc.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psMisc.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psMisc.h (revision 21655) @@ -0,0 +1,19 @@ +#if !defined(PS_MISC_H) +#define PS_MISC_H +/* + * Concatenate two macro arguments + */ +#define P_PS_CONCAT(A, B) A ## B // Expands to AB +#define PS_CONCAT(A, B) A ## B // Also Expands to AB +#define PS_CONCAT2(A, B) PS_CONCAT(A, B) // Also expands to AB +#define PS_CONCAT3(A, B, C) A ## B ## C // Expands to ABC +#define PS_CONCAT4(A, B, C, D) A ## B ## C ## D // Expands to ABCD + +#define PS_STRING(S) #S // converts argument to string + +void psAbort(const char *name, const char *fmt, ...); +void psError(const char *name, const char *fmt, ...); + +char *psStringCopy(const char *str); + +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psTrace.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psTrace.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psTrace.h (revision 21655) @@ -0,0 +1,20 @@ +#if !defined(PS_TRACE_H) +#define PS_TRACE_H 1 + +//#define PS_NTRACE 1 /* to turn off all tracing */ +#if defined(PS_NTRACE) +# define psTrace(facil, level, ...) /* do nothing */ +#else +# define psTrace(facil, level, ...) \ + p_psTrace(facil, level, __VA_ARGS__) +#endif + +void p_psTrace(const char *facil, int level, ...); + +int psSetTraceLevel(const char *facil, // facilty of interest + int level); // desired trace level +int psGetTraceLevel(const char *name); // facilty of interest + +void psPrintTraceLevels(void); // print trace levels + +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psUtils.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psUtils.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/psUtils.h (revision 21655) @@ -0,0 +1,13 @@ +#if !defined(PS_UTILS) +#define PS_UTILS + +#include "psArray.h" +#include "psDlist.h" +#include "psHash.h" +#include "psLogMsg.h" +#include "psMemory.h" +#include "psTrace.h" + +#include "psMisc.h" // must go last + +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/trace.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/trace.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/pslib/src/Utils/trace.c (revision 21655) @@ -0,0 +1,326 @@ +/* + * A simple implementation of a tracing facility for Pan-STARRS + * + * Tracing is controlled on a per "component" basis, where a "component" is a + * name of the form aaa.bbb.ccc where aaa is the Most significant part; for + * example, the utilities library might be called "utils", the doubly-linked + * list "utils.dlist", and the code to destroy a list "utils.dlist.del" + */ +#include +#include +#include +#include +#include "psUtils.h" + +/*****************************************************************************/ +/* + * A component is a string of the form aaa.bbb.ccc, and may itself + * contain further subcomponents. The Component structure doesn't + * in fact contain it's full name, but only the last part. + */ +typedef struct Component { + const char *name; // last part of name of this component + int level; // trace level for this component + int dimen; // dimension of subcomp + int n; // number of subcomponents + struct Component **subcomp; // next level of subcomponents +} Component; + +static int highest_level = 0; // highest level requested +/* + * tracelevel cache + */ +#define CACHE_NAME_LEN 50 +#define NO_CACHE "\a" // no name is cached + +#define UNKNOWN_LEVEL -9999 // we don't know this name's level + +static char cachedName[CACHE_NAME_LEN + 1] = NO_CACHE; // last name looked up +static int cachedLevel; // level of last looked up name + +/*****************************************************************************/ + +static Component *componentNew(const char *name, + int level) +{ + Component *comp = psAlloc(sizeof(Component)); + + comp->name = psStringCopy(name); + comp->level = level; + comp->dimen = comp->n = 0; + comp->subcomp = NULL; + + return comp; +} + +static void componentDel(Component *comp) +{ + if (comp == NULL) { return; } + + if (comp->subcomp != NULL) { + for (int i = 0; comp->subcomp[i] != NULL; i++) { + componentDel(comp->subcomp[i]); + } + psFree(comp->subcomp); + } + + psFree(comp); +} + +/*****************************************************************************/ +/* + * The root of the trace component tree + */ +static Component *croot = NULL; + +static void initTrace(void) { + if (croot == NULL) { + croot = componentNew("", UNKNOWN_LEVEL); + } +} + +/*****************************************************************************/ +/* + * Add a new component to the tree + */ +static void componentAdd(Component *comp, + const char *aname, + int level) +{ + char name[strlen(aname) + 1]; // need a writeable copy + strcpy(name, aname); + + char *cpt0; // first component of name + char *rest = name; // rest of name + cpt0 = strsep(&rest, "."); + /* + * Does cpt0 match this level? + */ + if (strcmp(comp->name, cpt0) == 0) { // a match + if (rest == NULL) { + comp->level = level; + } else { + componentAdd(comp, rest, level); + } + + return; + } + /* + * Look for a match for cpt0 in this level's subcomps + */ + for (int i = 0; i < comp->n; i++) { + if (strcmp(comp->subcomp[i]->name, cpt0) == 0) { // a match + if (rest == NULL) { + comp->subcomp[i]->level = level; + } else { + componentAdd(comp->subcomp[i], rest, level); + } + + return; + } + } + /* + * No match; add cpt0 to this level. Ensure that subcomp is sorted + */ + if (comp->n == comp->dimen) { + comp->dimen += 5; + comp->subcomp = psRealloc(comp->subcomp, comp->dimen); + } + + Component *fcpt0 = componentNew(cpt0, UNKNOWN_LEVEL); + int i; + for (i = 0; i < comp->n; i++) { + if (strcmp(cpt0, comp->subcomp[i]->name) > 0) { + for (int j = comp->n; j > i; j--) { + comp->subcomp[j] = comp->subcomp[j - 1]; + } + } + } + comp->subcomp[i] = fcpt0; + comp->n++; + + if (rest == NULL) { + fcpt0->level = level; + } else { + componentAdd(fcpt0, rest, level); + } +} + +/*****************************************************************************/ +/* + * Find the highest level present in the tree + */ +static void setHighestLevel(const Component *comp) +{ + if (comp->level > highest_level) { + highest_level = comp->level; + } + + for (int i = 0; i < comp->n; i++) { + setHighestLevel(comp->subcomp[i]); + } +} + +/*****************************************************************************/ + +int psSetTraceLevel(const char *comp, // component of interest + int level) // desired trace level +{ + if (croot == NULL) { + initTrace(); + } + + strncpy(cachedName, NO_CACHE, CACHE_NAME_LEN); // invalidate cache + + componentAdd(croot, comp, level); + + if (level > highest_level) { + highest_level = level; + } else { + highest_level = UNKNOWN_LEVEL; + setHighestLevel(croot); + } + + return 0; +} + +/*****************************************************************************/ +/* + * Return a trace level given a name + */ +static int doGetTraceLevel(const Component *comp, // end of component to search + const char *aname) // name to find +{ + char name[strlen(aname) + 1]; // need a writeable copy + strcpy(name, aname); + + char *cpt0; // first component of name + char *rest = name; // rest of name + cpt0 = strsep(&rest, "."); + /* + * Look for a match for cpt0 in this level's subcomps + */ + for (int i = 0; i < comp->n; i++) { + if (strcmp(comp->subcomp[i]->name, cpt0) == 0) { // a match + if (rest == NULL) { // we're there + int level = comp->subcomp[i]->level; + + return (level == UNKNOWN_LEVEL) ? comp->level : level; + } else { + return doGetTraceLevel(comp->subcomp[i], rest); + } + } + } + /* + * No match. This is as far as she goes + */ + return comp->level; +} + +int psGetTraceLevel(const char *name) // component of interest +{ + if (croot == NULL) { + initTrace(); + } + + if (strcmp(name, cachedName) == 0) { + return cachedLevel; + } + + int level = doGetTraceLevel(croot, name); + + strncpy(cachedName, name, CACHE_NAME_LEN); + cachedLevel = level; + + return level; +} + +/*****************************************************************************/ +/* + * Print a tree of trace levels + */ +static void doPrintTraceLevels(const Component *comp, int depth) +{ + if (comp->name[0] == '\0') { + printf("%*s%-*s %d\n", depth, "", 20 - depth, + "(root)", (comp->level == UNKNOWN_LEVEL) ? 0 : comp->level); + } else { + if (comp->level == UNKNOWN_LEVEL) { + printf("%*s%-*s %s\n", depth, "", 20 - depth, + comp->name, "."); + } else { + printf("%*s%-*s %d\n", depth, "", 20 - depth, + comp->name, comp->level); + } + } + + for (int i = 0; i < comp->n; i++) { + doPrintTraceLevels(comp->subcomp[i], depth + 1); + } +} + +void psPrintTraceLevels(void) +{ + if (croot == NULL) { + initTrace(); + } + + doPrintTraceLevels(croot, 0); +} + +/*****************************************************************************/ +/* + * Actually generate the trace message. Note that psTrace dealt with + * prepending appropriate indentation + */ +void p_psTrace(const char *comp, // component being traced + int level, // desired trace level + ...) // arguments +{ + if (level <= highest_level && psGetTraceLevel(comp) >= level) { + va_list ap; + va_start(ap, level); + for (int i = 0; i < level; i++) { + putchar(' '); + } + char *fmt = va_arg(ap, char *); + vprintf(fmt, ap); + va_end(ap); + } +} + +/*****************************************************************************/ +#if 0 // testing + +int main(void) +{ +#if 0 + psSetTraceLevel("", 10); +#else + psSetTraceLevel("aaa.bbb.ccc.ddd", 9); + psSetTraceLevel("aaa.bbb.ccc", 3); + psSetTraceLevel("aaa.bbb.ddd", 4); + psSetTraceLevel("aaa.BBB", 2); + psSetTraceLevel("BBB", 2); + psSetTraceLevel("BBB.XXX.YYY.ZZZ", 5); + psSetTraceLevel("aaa.bbb", 2); + psSetTraceLevel("aaa.bbb.ccc", 9); + psSetTraceLevel("aaa", 1); +#endif + + psPrintTraceLevels(); + printf("\n"); + + psTrace("aaa.bbb.ddd", 2, "(aaa.bbb.ddd) Hello %s\n", "world"); + psTrace("aaa", 2, "aaa\n"); + psTrace("", 2, "\"\"\n"); + psTrace("aaa.bbb", 2, "aaa.bbb\n"); + psTrace("aaa.bbb.ddd", 2, "aaa.bbb.ddd\n"); + psTrace("aaa.bbb.ddd", 4, "aaa.bbb.ddd\n"); + psTrace("aaa.bbb.ddd", 2, "aaa.bbb.ddd\n"); + psTrace("aaa.bbb.ccc", 1, "aaa.bbb.ccc\n"); + psTrace("aaa.bbb.eee", 2, "aaa.bbb.eee\n"); + + return 0; +} +#endif Index: /branches/ipp-1-X-branches/v1_0/archive/utils/check-namespace =================================================================== --- /branches/ipp-1-X-branches/v1_0/archive/utils/check-namespace (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/archive/utils/check-namespace (revision 21655) @@ -0,0 +1,301 @@ +#!/usr/bin/perl -w +# +# Check that a set of .h files and .o and .a files satisfy the Pan-STARRS +# rules about namespaces. +# +# Robert Lupton (rhl@astro.princeton.edu) February 2004 +# +# Parse arguments +# +$nm_opts = ""; +$prefix = "ps"; +$pprefix = "p_$prefix"; +$verbose = 0; + +if (!@ARGV) { + warn "Please specify some files to check! E.g. *.h *.o\n"; + push @ARGV, "-h"; +} + +while ($ARGV[0] =~ /^-/) { + if ($ARGV[0] eq "-h" || $ARGV[0] eq "-?") { + warn "\ +Usage: check-namespace [options] files.h files.o + +Check that a set of .h and .o (and/or .a) files obey the pan-STARRS +naming convention, namely that iff a symbol is externally visible, it +should have a name starting \"ps\" or \"p_ps\", and be declared in a +.h file. + +E.g. check-namespace *.o *.h + +The treatment of comments and typedefs is sleazy and non-robust. + +Options: + -h, -? Print this message + -p prefix Change prefix from ps to \"prefix\" + -v Be chatty (repeat for even more output) +"; + exit(0); + } elsif($ARGV[0] eq "-p") { + shift(@ARGV); + $arg = $ARGV[0]; + if(!$arg) { + warn "-p expects an argument\n"; + } else { + $prefix = "$arg"; + $pprefix = "p_$prefix"; + } + } elsif($ARGV[0] eq "-v") { + $verbose++; + } else { + warn "Unknown argument $ARGV[0]\n"; + } + + shift(@ARGV); +} +# +$uprefix = uc($prefix); +$upprefix = uc($pprefix); +# +# Get list of files to process +# +foreach (@ARGV) { + if(/\.h$/) { + push(@hfiles, $_); + } elsif(/\.[ao]$/) { + push(@ofiles, $_); + } else { + warn "I don't know what to do with $_; ignoring\n"; + } +} +# +# Read all the .h files looking for suitable symbols +# +foreach $file (@hfiles) { + $in_comment = 0; # in a multiline comment? + $in_typedef = 0; # parsing a typedef? + + if ($verbose > 1) { + warn " $file\n"; + } + + open(FD, $file) or die "Failed to open $file:$! \n"; + while () { + chomp; + # Deal with comments, at least approximately + if($in_comment) { + while (!m|\*/|) { + $_ = ; + } + $in_comment = 0; + } + + s|//.*||; + if(m|(.*)/\*(.*)|) { + $first = $1; $rest = $2; + if ($rest =~ m|\*/(.*)|) { + $_ = "$first /**/ $1"; + } else { + $_ = $first; + $in_comment = 1; + } + } + # + # OK, how about typedefs? + # + if (/typedef.*\W\(\*(($pprefix|$prefix)[a-zA-Z_0-9]+)\)/o) { + $name = $1; + + $typedefs{$name} = $file; + $in_typedef--; + } elsif (/typedef\s+struct/) { + if(/([a-zA-Z_0-9]+)\s*;/) { + $name = $1; + $typedefs{$name} = $file; + next; + } + $in_typedef = 1; + } elsif ($in_typedef) { + if (/\{/) { + $in_typedef++; + } elsif ($in_typedef > 1 && /\}/) { + $in_typedef--; + } + } + + if ($in_typedef) { + if (/\}\s*(($prefix|$pprefix)[a-zA-Z_0-9]+)\s*;/) { + $in_typedef--; + if ($in_typedef == 0) { + $name = $1; + $typedefs{$name} = $file; + } + } + } + # + # And #defines? + # + if (/^\s*\#\s*define\s+([a-zA-Z_0-9]+)(\()?/) { + $name = $1; + $func_macro = (defined($2)) ? 1 : 0; + + $defines{$name} = $file; + + if((!$func_macro && $name !~ /^($uprefix|$upprefix)/o) || + ($func_macro && $name !~ /^($prefix|$pprefix)/io)) { + printf STDERR + "\#define %s (%s) does not have a permitted prefix\n", + $name, $file; + } + } + # + # OK, how about enums? Their members look like external names, + # and should obey the same rules + # + if (s/((typedef\s+)?enum)//) { + $typedef = $1; + $is_enum_typedef = ($typedef =~ /typedef/) ? 1 : 0; + + while (!/\}/) { + $line = ; + chomp; + s|//.*||; + s|/\*[^\*]*\*/||; + + $_ .= " $line"; + } + + if ($is_enum_typedef) { + if ($line =~ /\}\s*([a-zA-Z_0-9]+)/) { + $name = $1; + $typedefs{$name} = $file; + } + } + + s/^\s*\{\s*//; + while (s/([a-zA-Z_0-9]+)(\s*=\s*-?[a-zA-Z_0-9]+)?(\s*,\s*)?\s*//) { + if(/^\}/) { + last; + } + + $name = $1; + + if ($name eq "enum") { + next; + } + + #$enums{$1} = $file; + + if($name !~ /^($uprefix|$upprefix)/o) { + printf STDERR + "enum element %s (%s) does not have a permitted prefix\n", + $name, $file; + } + } + } + # + # Include filenames can look like external names too + # + if (/^\s*\#\s*include\s/) { + next; + } + # Look for possible external symbols + # + while(s/(^|\W)(($pprefix|$prefix)[a-zA-Z_0-9]+)//o) { + $name = $2; + + $declarations{$name} = $file; + } + } + close(FD); +} + +if($verbose) { + warn "Typedefs found in header files:\n"; + foreach (sort keys %typedefs) { + printf STDERR "%-30s %s\n", $_, $typedefs{$_}; + } + + warn "External symbols found in header files:\n"; + foreach (sort keys %declarations) { + printf STDERR "%-30s %s\n", $_, $declarations{$_}; + } +} +# +# OK, now run nm on the .o and .a files and see if there are problems +# +foreach $file (@ofiles) { + if ($verbose) { + warn " $file\n"; + } + + open(FD, "nm $nm_opts $file|") or die "Failed to open $file:$! \n"; + while () { + if (/^\s*$/) { + next; + } + + if (/^(lib[A-Za-z0-9]+\.a\(\w+\.o\))/) { + $file = $1; + next; + } + + if (/^\s+U /) { + next; # an undefined symbol + } + + my($addr, $type, $name) = split; + + if (!defined($name)) { + next; + } + $name =~ s/^_//; + + if ($name eq "main") { + next; + } + if ($name =~ /\./) { + next; # an internal name to the compiler + } + + if ($type eq "b") { $type = "d"; } # don't distinguish BSS and Data + if ($type eq "B") { $type = "D"; } + + if ($type !~ /[DT]/i) { # Data or Text section + die "Unknown symbol type in $file: $_\n"; + } + + if ($verbose > 1) { + printf "%-20s %-40s %s\n", $file, $name, $type; + } + + if ($type =~ /[dt]/) { + if($name =~ /^($prefix|$pprefix)/o) { + warn "Non-global symbol $name in $file is in reserved namespace ($prefix|$pprefix)\n"; + } + } elsif ($type =~ /[DT]/) { + if ($declarations{$name}) { + $used_declarations{$name}++; + } elsif ($name !~ /^($pprefix|$prefix)/) { + if($name !~ /^($prefix|$pprefix)/io) { + printf STDERR + "global symbol %s (%s) does not have a permitted prefix\n", + $name, $file; + } + } else { + printf STDERR "External symbol %s (%s) appears in no .h file\n", + $name, $file; + } + } + } + close(FD); +} +# +# Now check for unused declarations +# +foreach $name (sort keys %declarations) { + if (!$used_declarations{$name} && !$typedefs{$name} && !$defines{$name}) { + warn "$name was declared in $declarations{$name} but not defined\n"; + } +} Index: /branches/ipp-1-X-branches/v1_0/pois/Makefile.am =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/Makefile.am (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/Makefile.am (revision 21655) @@ -0,0 +1,3 @@ +SUBDIRS = src swig + +CLEANFILES = *.pyc *~ core core.* Index: /branches/ipp-1-X-branches/v1_0/pois/NOTES =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/NOTES (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/NOTES (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/autogen.sh =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/autogen.sh (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/autogen.sh (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/configure.ac =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/configure.ac (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/configure.ac (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/python/ds9.py =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/python/ds9.py (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/python/ds9.py (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/python/pois.py =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/python/pois.py (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/python/pois.py (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/python/ps.py =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/python/ps.py (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/python/ps.py (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/Makefile.am =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/Makefile.am (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/Makefile.am (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/pois.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/pois.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/pois.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/pois.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/pois.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/pois.h (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/poisAddPenalty.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/poisAddPenalty.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/poisAddPenalty.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/poisCalculateDeviations.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/poisCalculateDeviations.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/poisCalculateDeviations.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/poisCalculateEquation.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/poisCalculateEquation.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/poisCalculateEquation.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/poisConvolveImage.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/poisConvolveImage.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/poisConvolveImage.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/poisExtractKernel.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/poisExtractKernel.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/poisExtractKernel.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/poisFindStamps.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/poisFindStamps.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/poisFindStamps.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/poisKernelBasisFunctions.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/poisKernelBasisFunctions.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/poisKernelBasisFunctions.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/poisMakeMask.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/poisMakeMask.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/poisMakeMask.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/poisParseConfig.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/poisParseConfig.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/poisParseConfig.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/poisRejectStamps.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/poisRejectStamps.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/poisRejectStamps.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/poisSolveEquation.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/poisSolveEquation.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/poisSolveEquation.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/src/poisStamp.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/src/poisStamp.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/src/poisStamp.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/swig/.cvsignore =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/swig/.cvsignore (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/swig/.cvsignore (revision 21655) @@ -0,0 +1,7 @@ +*.pyc +*_wrap.c +psLib.py +swigrun.py +xpa.py +.Makefile.depend +*.pm Index: /branches/ipp-1-X-branches/v1_0/pois/swig/Makefile.am =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/swig/Makefile.am (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/swig/Makefile.am (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/swig/p_psSwig.i =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/swig/p_psSwig.i (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/swig/p_psSwig.i (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/swig/poisSwig.i =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/swig/poisSwig.i (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/swig/poisSwig.i (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/swig/psLibForSwig.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/swig/psLibForSwig.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/swig/psLibForSwig.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/swig/psLibForSwig.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/swig/psLibForSwig.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/swig/psLibForSwig.h (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/swig/psLibSwig.i =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/swig/psLibSwig.i (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/swig/psLibSwig.i (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/swig/simpleFits.c =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/swig/simpleFits.c (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/swig/simpleFits.c (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/swig/simpleFits.h =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/swig/simpleFits.h (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/swig/simpleFits.h (revision 21655) @@ -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/ipp-1-X-branches/v1_0/pois/swig/xpaSwig.i =================================================================== --- /branches/ipp-1-X-branches/v1_0/pois/swig/xpaSwig.i (revision 21655) +++ /branches/ipp-1-X-branches/v1_0/pois/swig/xpaSwig.i (revision 21655) @@ -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"