Index: /branches/ccl_branches/ipponly-20191108/ppSub/Doxyfile.in =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/Doxyfile.in (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/Doxyfile.in (revision 41094) @@ -0,0 +1,1310 @@ +# Doxyfile 1.5.4 + +# This file describes the settings to be used by the documentation system +# doxygen (www.doxygen.org) for a project +# +# All text after a hash (#) is considered a comment and will be ignored +# The format is: +# TAG = value [value, ...] +# For lists items can also be appended using: +# TAG += value [value, ...] +# Values that contain spaces should be placed between quotes (" ") + +#--------------------------------------------------------------------------- +# Project related configuration options +#--------------------------------------------------------------------------- + +# This tag specifies the encoding used for all characters in the config file that +# follow. The default is UTF-8 which is also the encoding used for all text before +# the first occurrence of this tag. Doxygen uses libiconv (or the iconv built into +# libc) for the transcoding. See http://www.gnu.org/software/libiconv for the list of +# possible encodings. + +DOXYFILE_ENCODING = UTF-8 + +# The PROJECT_NAME tag is a single word (or a sequence of words surrounded +# by quotes) that should identify the project. + +PROJECT_NAME = ppSub + +# The PROJECT_NUMBER tag can be used to enter a project or revision number. +# This could be handy for archiving the generated documentation or +# if some version control system is used. + +PROJECT_NUMBER = ipp-2.7.dev + +# The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) +# base path where the generated documentation will be put. +# If a relative path is entered, it will be relative to the location +# where doxygen was started. If left blank the current directory will be used. + +OUTPUT_DIRECTORY = ./docs + +# If the CREATE_SUBDIRS tag is set to YES, then doxygen will create +# 4096 sub-directories (in 2 levels) under the output directory of each output +# format and will distribute the generated files over these directories. +# Enabling this option can be useful when feeding doxygen a huge amount of +# source files, where putting all generated files in the same directory would +# otherwise cause performance problems for the file system. + +CREATE_SUBDIRS = NO + +# The OUTPUT_LANGUAGE tag is used to specify the language in which all +# documentation generated by doxygen is written. Doxygen will use this +# information to generate all constant output in the proper language. +# The default language is English, other supported languages are: +# Afrikaans, Arabic, Brazilian, Catalan, Chinese, Chinese-Traditional, +# Croatian, Czech, Danish, Dutch, Finnish, French, German, Greek, Hungarian, +# Italian, Japanese, Japanese-en (Japanese with English messages), Korean, +# Korean-en, Lithuanian, Norwegian, Polish, Portuguese, Romanian, Russian, +# Serbian, Slovak, Slovene, Spanish, Swedish, and Ukrainian. + +OUTPUT_LANGUAGE = English + +# If the BRIEF_MEMBER_DESC tag is set to YES (the default) Doxygen will +# include brief member descriptions after the members that are listed in +# the file and class documentation (similar to JavaDoc). +# Set to NO to disable this. + +BRIEF_MEMBER_DESC = YES + +# If the REPEAT_BRIEF tag is set to YES (the default) Doxygen will prepend +# the brief description of a member or function before the detailed description. +# Note: if both HIDE_UNDOC_MEMBERS and BRIEF_MEMBER_DESC are set to NO, the +# brief descriptions will be completely suppressed. + +REPEAT_BRIEF = YES + +# This tag implements a quasi-intelligent brief description abbreviator +# that is used to form the text in various listings. Each string +# in this list, if found as the leading text of the brief description, will be +# stripped from the text and the result after processing the whole list, is +# used as the annotated text. Otherwise, the brief description is used as-is. +# If left blank, the following values are used ("$name" is automatically +# replaced with the name of the entity): "The $name class" "The $name widget" +# "The $name file" "is" "provides" "specifies" "contains" +# "represents" "a" "an" "the" + +ABBREVIATE_BRIEF = + +# If the ALWAYS_DETAILED_SEC and REPEAT_BRIEF tags are both set to YES then +# Doxygen will generate a detailed section even if there is only a brief +# description. + +ALWAYS_DETAILED_SEC = NO + +# If the INLINE_INHERITED_MEMB tag is set to YES, doxygen will show all +# inherited members of a class in the documentation of that class as if those +# members were ordinary class members. Constructors, destructors and assignment +# operators of the base classes will not be shown. + +INLINE_INHERITED_MEMB = NO + +# If the FULL_PATH_NAMES tag is set to YES then Doxygen will prepend the full +# path before files name in the file list and in the header files. If set +# to NO the shortest path that makes the file name unique will be used. + +FULL_PATH_NAMES = NO + +# If the FULL_PATH_NAMES tag is set to YES then the STRIP_FROM_PATH tag +# can be used to strip a user-defined part of the path. Stripping is +# only done if one of the specified strings matches the left-hand part of +# the path. The tag can be used to show relative paths in the file list. +# If left blank the directory from which doxygen is run is used as the +# path to strip. + +STRIP_FROM_PATH = + +# The STRIP_FROM_INC_PATH tag can be used to strip a user-defined part of +# the path mentioned in the documentation of a class, which tells +# the reader which header file to include in order to use a class. +# If left blank only the name of the header file containing the class +# definition is used. Otherwise one should specify the include paths that +# are normally passed to the compiler using the -I flag. + +STRIP_FROM_INC_PATH = + +# If the SHORT_NAMES tag is set to YES, doxygen will generate much shorter +# (but less readable) file names. This can be useful is your file systems +# doesn't support long names like on DOS, Mac, or CD-ROM. + +SHORT_NAMES = NO + +# If the JAVADOC_AUTOBRIEF tag is set to YES then Doxygen +# will interpret the first line (until the first dot) of a JavaDoc-style +# comment as the brief description. If set to NO, the JavaDoc +# comments will behave just like regular Qt-style comments +# (thus requiring an explicit @brief command for a brief description.) + +JAVADOC_AUTOBRIEF = NO + +# If the QT_AUTOBRIEF tag is set to YES then Doxygen will +# interpret the first line (until the first dot) of a Qt-style +# comment as the brief description. If set to NO, the comments +# will behave just like regular Qt-style comments (thus requiring +# an explicit \brief command for a brief description.) + +QT_AUTOBRIEF = NO + +# The MULTILINE_CPP_IS_BRIEF tag can be set to YES to make Doxygen +# treat a multi-line C++ special comment block (i.e. a block of //! or /// +# comments) as a brief description. This used to be the default behaviour. +# The new default is to treat a multi-line C++ comment block as a detailed +# description. Set this tag to YES if you prefer the old behaviour instead. + +MULTILINE_CPP_IS_BRIEF = NO + +# If the DETAILS_AT_TOP tag is set to YES then Doxygen +# will output the detailed description near the top, like JavaDoc. +# If set to NO, the detailed description appears after the member +# documentation. + +DETAILS_AT_TOP = NO + +# If the INHERIT_DOCS tag is set to YES (the default) then an undocumented +# member inherits the documentation from any documented member that it +# re-implements. + +INHERIT_DOCS = YES + +# If the SEPARATE_MEMBER_PAGES tag is set to YES, then doxygen will produce +# a new page for each member. If set to NO, the documentation of a member will +# be part of the file/class/namespace that contains it. + +SEPARATE_MEMBER_PAGES = NO + +# The TAB_SIZE tag can be used to set the number of spaces in a tab. +# Doxygen uses this value to replace tabs by spaces in code fragments. + +TAB_SIZE = 4 + +# This tag can be used to specify a number of aliases that acts +# as commands in the documentation. An alias has the form "name=value". +# For example adding "sideeffect=\par Side Effects:\n" will allow you to +# put the command \sideeffect (or @sideeffect) in the documentation, which +# will result in a user-defined paragraph with heading "Side Effects:". +# You can put \n's in the value part of an alias to insert newlines. + +ALIASES = + +# Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C +# sources only. Doxygen will then generate output that is more tailored for C. +# For instance, some of the names that are used will be different. The list +# of all members will be omitted, etc. + +OPTIMIZE_OUTPUT_FOR_C = NO + +# Set the OPTIMIZE_OUTPUT_JAVA tag to YES if your project consists of Java +# sources only. Doxygen will then generate output that is more tailored for Java. +# For instance, namespaces will be presented as packages, qualified scopes +# will look different, etc. + +OPTIMIZE_OUTPUT_JAVA = NO + +# If you use STL classes (i.e. std::string, std::vector, etc.) but do not want to +# include (a tag file for) the STL sources as input, then you should +# set this tag to YES in order to let doxygen match functions declarations and +# definitions whose arguments contain STL classes (e.g. func(std::string); v.s. +# func(std::string) {}). This also make the inheritance and collaboration +# diagrams that involve STL classes more complete and accurate. + +BUILTIN_STL_SUPPORT = NO + +# If you use Microsoft's C++/CLI language, you should set this option to YES to +# enable parsing support. + +CPP_CLI_SUPPORT = NO + +# Set the SIP_SUPPORT tag to YES if your project consists of sip sources only. +# Doxygen will parse them like normal C++ but will assume all classes use public +# instead of private inheritance when no explicit protection keyword is present. + +SIP_SUPPORT = NO + +# If member grouping is used in the documentation and the DISTRIBUTE_GROUP_DOC +# tag is set to YES, then doxygen will reuse the documentation of the first +# member in the group (if any) for the other members of the group. By default +# all members of a group must be documented explicitly. + +DISTRIBUTE_GROUP_DOC = NO + +# Set the SUBGROUPING tag to YES (the default) to allow class member groups of +# the same type (for instance a group of public functions) to be put as a +# subgroup of that type (e.g. under the Public Functions section). Set it to +# NO to prevent subgrouping. Alternatively, this can be done per class using +# the \nosubgrouping command. + +SUBGROUPING = YES + +# When TYPEDEF_HIDES_STRUCT is enabled, a typedef of a struct (or union) is +# documented as struct with the name of the typedef. So +# typedef struct TypeS {} TypeT, will appear in the documentation as a struct +# with name TypeT. When disabled the typedef will appear as a member of a file, +# namespace, or class. And the struct will be named TypeS. This can typically +# be useful for C code where the coding convention is that all structs are +# typedef'ed and only the typedef is referenced never the struct's name. + +TYPEDEF_HIDES_STRUCT = NO + +#--------------------------------------------------------------------------- +# Build related configuration options +#--------------------------------------------------------------------------- + +# If the EXTRACT_ALL tag is set to YES doxygen will assume all entities in +# documentation are documented, even if no documentation was available. +# Private class members and static file members will be hidden unless +# the EXTRACT_PRIVATE and EXTRACT_STATIC tags are set to YES + +EXTRACT_ALL = NO + +# If the EXTRACT_PRIVATE tag is set to YES all private members of a class +# will be included in the documentation. + +EXTRACT_PRIVATE = NO + +# If the EXTRACT_STATIC tag is set to YES all static members of a file +# will be included in the documentation. + +EXTRACT_STATIC = NO + +# If the EXTRACT_LOCAL_CLASSES tag is set to YES classes (and structs) +# defined locally in source files will be included in the documentation. +# If set to NO only classes defined in header files are included. + +EXTRACT_LOCAL_CLASSES = YES + +# This flag is only useful for Objective-C code. When set to YES local +# methods, which are defined in the implementation section but not in +# the interface are included in the documentation. +# If set to NO (the default) only methods in the interface are included. + +EXTRACT_LOCAL_METHODS = NO + +# If this flag is set to YES, the members of anonymous namespaces will be extracted +# and appear in the documentation as a namespace called 'anonymous_namespace{file}', +# where file will be replaced with the base name of the file that contains the anonymous +# namespace. By default anonymous namespace are hidden. + +EXTRACT_ANON_NSPACES = NO + +# If the HIDE_UNDOC_MEMBERS tag is set to YES, Doxygen will hide all +# undocumented members of documented classes, files or namespaces. +# If set to NO (the default) these members will be included in the +# various overviews, but no documentation section is generated. +# This option has no effect if EXTRACT_ALL is enabled. + +HIDE_UNDOC_MEMBERS = NO + +# If the HIDE_UNDOC_CLASSES tag is set to YES, Doxygen will hide all +# undocumented classes that are normally visible in the class hierarchy. +# If set to NO (the default) these classes will be included in the various +# overviews. This option has no effect if EXTRACT_ALL is enabled. + +HIDE_UNDOC_CLASSES = NO + +# If the HIDE_FRIEND_COMPOUNDS tag is set to YES, Doxygen will hide all +# friend (class|struct|union) declarations. +# If set to NO (the default) these declarations will be included in the +# documentation. + +HIDE_FRIEND_COMPOUNDS = NO + +# If the HIDE_IN_BODY_DOCS tag is set to YES, Doxygen will hide any +# documentation blocks found inside the body of a function. +# If set to NO (the default) these blocks will be appended to the +# function's detailed documentation block. + +HIDE_IN_BODY_DOCS = NO + +# The INTERNAL_DOCS tag determines if documentation +# that is typed after a \internal command is included. If the tag is set +# to NO (the default) then the documentation will be excluded. +# Set it to YES to include the internal documentation. + +INTERNAL_DOCS = NO + +# If the CASE_SENSE_NAMES tag is set to NO then Doxygen will only generate +# file names in lower-case letters. If set to YES upper-case letters are also +# allowed. This is useful if you have classes or files whose names only differ +# in case and if your file system supports case sensitive file names. Windows +# and Mac users are advised to set this option to NO. + +CASE_SENSE_NAMES = YES + +# If the HIDE_SCOPE_NAMES tag is set to NO (the default) then Doxygen +# will show members with their full class and namespace scopes in the +# documentation. If set to YES the scope will be hidden. + +HIDE_SCOPE_NAMES = NO + +# If the SHOW_INCLUDE_FILES tag is set to YES (the default) then Doxygen +# will put a list of the files that are included by a file in the documentation +# of that file. + +SHOW_INCLUDE_FILES = YES + +# If the INLINE_INFO tag is set to YES (the default) then a tag [inline] +# is inserted in the documentation for inline members. + +INLINE_INFO = YES + +# If the SORT_MEMBER_DOCS tag is set to YES (the default) then doxygen +# will sort the (detailed) documentation of file and class members +# alphabetically by member name. If set to NO the members will appear in +# declaration order. + +SORT_MEMBER_DOCS = YES + +# If the SORT_BRIEF_DOCS tag is set to YES then doxygen will sort the +# brief documentation of file, namespace and class members alphabetically +# by member name. If set to NO (the default) the members will appear in +# declaration order. + +SORT_BRIEF_DOCS = NO + +# If the SORT_BY_SCOPE_NAME tag is set to YES, the class list will be +# sorted by fully-qualified names, including namespaces. If set to +# NO (the default), the class list will be sorted only by class name, +# not including the namespace part. +# Note: This option is not very useful if HIDE_SCOPE_NAMES is set to YES. +# Note: This option applies only to the class list, not to the +# alphabetical list. + +SORT_BY_SCOPE_NAME = NO + +# The GENERATE_TODOLIST tag can be used to enable (YES) or +# disable (NO) the todo list. This list is created by putting \todo +# commands in the documentation. + +GENERATE_TODOLIST = YES + +# The GENERATE_TESTLIST tag can be used to enable (YES) or +# disable (NO) the test list. This list is created by putting \test +# commands in the documentation. + +GENERATE_TESTLIST = YES + +# The GENERATE_BUGLIST tag can be used to enable (YES) or +# disable (NO) the bug list. This list is created by putting \bug +# commands in the documentation. + +GENERATE_BUGLIST = YES + +# The GENERATE_DEPRECATEDLIST tag can be used to enable (YES) or +# disable (NO) the deprecated list. This list is created by putting +# \deprecated commands in the documentation. + +GENERATE_DEPRECATEDLIST= YES + +# The ENABLED_SECTIONS tag can be used to enable conditional +# documentation sections, marked by \if sectionname ... \endif. + +ENABLED_SECTIONS = + +# The MAX_INITIALIZER_LINES tag determines the maximum number of lines +# the initial value of a variable or define consists of for it to appear in +# the documentation. If the initializer consists of more lines than specified +# here it will be hidden. Use a value of 0 to hide initializers completely. +# The appearance of the initializer of individual variables and defines in the +# documentation can be controlled using \showinitializer or \hideinitializer +# command in the documentation regardless of this setting. + +MAX_INITIALIZER_LINES = 30 + +# Set the SHOW_USED_FILES tag to NO to disable the list of files generated +# at the bottom of the documentation of classes and structs. If set to YES the +# list will mention the files that were used to generate the documentation. + +SHOW_USED_FILES = YES + +# If the sources in your project are distributed over multiple directories +# then setting the SHOW_DIRECTORIES tag to YES will show the directory hierarchy +# in the documentation. The default is NO. + +SHOW_DIRECTORIES = YES + +# The FILE_VERSION_FILTER tag can be used to specify a program or script that +# doxygen should invoke to get the current version for each file (typically from the +# version control system). Doxygen will invoke the program by executing (via +# popen()) the command , where is the value of +# the FILE_VERSION_FILTER tag, and is the name of an input file +# provided by doxygen. Whatever the program writes to standard output +# is used as the file version. See the manual for examples. + +FILE_VERSION_FILTER = + +#--------------------------------------------------------------------------- +# configuration options related to warning and progress messages +#--------------------------------------------------------------------------- + +# The QUIET tag can be used to turn on/off the messages that are generated +# by doxygen. Possible values are YES and NO. If left blank NO is used. + +QUIET = NO + +# The WARNINGS tag can be used to turn on/off the warning messages that are +# generated by doxygen. Possible values are YES and NO. If left blank +# NO is used. + +WARNINGS = YES + +# If WARN_IF_UNDOCUMENTED is set to YES, then doxygen will generate warnings +# for undocumented members. If EXTRACT_ALL is set to YES then this flag will +# automatically be disabled. + +WARN_IF_UNDOCUMENTED = YES + +# If WARN_IF_DOC_ERROR is set to YES, doxygen will generate warnings for +# potential errors in the documentation, such as not documenting some +# parameters in a documented function, or documenting parameters that +# don't exist or using markup commands wrongly. + +WARN_IF_DOC_ERROR = YES + +# This WARN_NO_PARAMDOC option can be abled to get warnings for +# functions that are documented, but have no documentation for their parameters +# or return value. If set to NO (the default) doxygen will only warn about +# wrong or incomplete parameter documentation, but not about the absence of +# documentation. + +WARN_NO_PARAMDOC = NO + +# The WARN_FORMAT tag determines the format of the warning messages that +# doxygen can produce. The string should contain the $file, $line, and $text +# tags, which will be replaced by the file and line number from which the +# warning originated and the warning text. Optionally the format may contain +# $version, which will be replaced by the version of the file (if it could +# be obtained via FILE_VERSION_FILTER) + +WARN_FORMAT = "$file:$line: $text " + +# The WARN_LOGFILE tag can be used to specify a file to which warning +# and error messages should be written. If left blank the output is written +# to stderr. + +WARN_LOGFILE = + +#--------------------------------------------------------------------------- +# configuration options related to the input files +#--------------------------------------------------------------------------- + +# The INPUT tag can be used to specify the files and/or directories that contain +# documented source files. You may enter file names like "myfile.cpp" or +# directories like "/usr/src/myproject". Separate the files or directories +# with spaces. + +INPUT = ./src + +# This tag can be used to specify the character encoding of the source files that +# doxygen parses. Internally doxygen uses the UTF-8 encoding, which is also the default +# input encoding. Doxygen uses libiconv (or the iconv built into libc) for the transcoding. +# See http://www.gnu.org/software/libiconv for the list of possible encodings. + +INPUT_ENCODING = UTF-8 + +# If the value of the INPUT tag contains directories, you can use the +# FILE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp +# and *.h) to filter out the source-files in the directories. If left +# blank the following patterns are tested: +# *.c *.cc *.cxx *.cpp *.c++ *.java *.ii *.ixx *.ipp *.i++ *.inl *.h *.hh *.hxx +# *.hpp *.h++ *.idl *.odl *.cs *.php *.php3 *.inc *.m *.mm *.py *.f90 + +FILE_PATTERNS = *.h + +# The RECURSIVE tag can be used to turn specify whether or not subdirectories +# should be searched for input files as well. Possible values are YES and NO. +# If left blank NO is used. + +RECURSIVE = YES + +# The EXCLUDE tag can be used to specify files and/or directories that should +# excluded from the INPUT source files. This way you can easily exclude a +# subdirectory from a directory tree whose root is specified with the INPUT tag. + +EXCLUDE = + +# The EXCLUDE_SYMLINKS tag can be used select whether or not files or +# directories that are symbolic links (a Unix filesystem feature) are excluded +# from the input. + +EXCLUDE_SYMLINKS = NO + +# If the value of the INPUT tag contains directories, you can use the +# EXCLUDE_PATTERNS tag to specify one or more wildcard patterns to exclude +# certain files from those directories. Note that the wildcards are matched +# against the file with absolute path, so to exclude all test directories +# for example use the pattern */test/* + +EXCLUDE_PATTERNS = + +# The EXCLUDE_SYMBOLS tag can be used to specify one or more symbol names +# (namespaces, classes, functions, etc.) that should be excluded from the output. +# The symbol name can be a fully qualified name, a word, or if the wildcard * is used, +# a substring. Examples: ANamespace, AClass, AClass::ANamespace, ANamespace::*Test + +EXCLUDE_SYMBOLS = + +# The EXAMPLE_PATH tag can be used to specify one or more files or +# directories that contain example code fragments that are included (see +# the \include command). + +EXAMPLE_PATH = + +# If the value of the EXAMPLE_PATH tag contains directories, you can use the +# EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp +# and *.h) to filter out the source-files in the directories. If left +# blank all files are included. + +EXAMPLE_PATTERNS = + +# If the EXAMPLE_RECURSIVE tag is set to YES then subdirectories will be +# searched for input files to be used with the \include or \dontinclude +# commands irrespective of the value of the RECURSIVE tag. +# Possible values are YES and NO. If left blank NO is used. + +EXAMPLE_RECURSIVE = NO + +# The IMAGE_PATH tag can be used to specify one or more files or +# directories that contain image that are included in the documentation (see +# the \image command). + +IMAGE_PATH = + +# The INPUT_FILTER tag can be used to specify a program that doxygen should +# invoke to filter for each input file. Doxygen will invoke the filter program +# by executing (via popen()) the command , where +# is the value of the INPUT_FILTER tag, and is the name of an +# input file. Doxygen will then use the output that the filter program writes +# to standard output. If FILTER_PATTERNS is specified, this tag will be +# ignored. + +INPUT_FILTER = + +# The FILTER_PATTERNS tag can be used to specify filters on a per file pattern +# basis. Doxygen will compare the file name with each pattern and apply the +# filter if there is a match. The filters are a list of the form: +# pattern=filter (like *.cpp=my_cpp_filter). See INPUT_FILTER for further +# info on how filters are used. If FILTER_PATTERNS is empty, INPUT_FILTER +# is applied to all files. + +FILTER_PATTERNS = + +# If the FILTER_SOURCE_FILES tag is set to YES, the input filter (if set using +# INPUT_FILTER) will be used to filter the input files when producing source +# files to browse (i.e. when SOURCE_BROWSER is set to YES). + +FILTER_SOURCE_FILES = NO + +#--------------------------------------------------------------------------- +# configuration options related to source browsing +#--------------------------------------------------------------------------- + +# If the SOURCE_BROWSER tag is set to YES then a list of source files will +# be generated. Documented entities will be cross-referenced with these sources. +# Note: To get rid of all source code in the generated output, make sure also +# VERBATIM_HEADERS is set to NO. If you have enabled CALL_GRAPH or CALLER_GRAPH +# then you must also enable this option. If you don't then doxygen will produce +# a warning and turn it on anyway + +SOURCE_BROWSER = NO + +# Setting the INLINE_SOURCES tag to YES will include the body +# of functions and classes directly in the documentation. + +INLINE_SOURCES = NO + +# Setting the STRIP_CODE_COMMENTS tag to YES (the default) will instruct +# doxygen to hide any special comment blocks from generated source code +# fragments. Normal C and C++ comments will always remain visible. + +STRIP_CODE_COMMENTS = YES + +# If the REFERENCED_BY_RELATION tag is set to YES (the default) +# then for each documented function all documented +# functions referencing it will be listed. + +REFERENCED_BY_RELATION = YES + +# If the REFERENCES_RELATION tag is set to YES (the default) +# then for each documented function all documented entities +# called/used by that function will be listed. + +REFERENCES_RELATION = YES + +# If the REFERENCES_LINK_SOURCE tag is set to YES (the default) +# and SOURCE_BROWSER tag is set to YES, then the hyperlinks from +# functions in REFERENCES_RELATION and REFERENCED_BY_RELATION lists will +# link to the source code. Otherwise they will link to the documentstion. + +REFERENCES_LINK_SOURCE = YES + +# If the USE_HTAGS tag is set to YES then the references to source code +# will point to the HTML generated by the htags(1) tool instead of doxygen +# built-in source browser. The htags tool is part of GNU's global source +# tagging system (see http://www.gnu.org/software/global/global.html). You +# will need version 4.8.6 or higher. + +USE_HTAGS = NO + +# If the VERBATIM_HEADERS tag is set to YES (the default) then Doxygen +# will generate a verbatim copy of the header file for each class for +# which an include is specified. Set to NO to disable this. + +VERBATIM_HEADERS = YES + +#--------------------------------------------------------------------------- +# configuration options related to the alphabetical class index +#--------------------------------------------------------------------------- + +# If the ALPHABETICAL_INDEX tag is set to YES, an alphabetical index +# of all compounds will be generated. Enable this if the project +# contains a lot of classes, structs, unions or interfaces. + +ALPHABETICAL_INDEX = NO + +# If the alphabetical index is enabled (see ALPHABETICAL_INDEX) then +# the COLS_IN_ALPHA_INDEX tag can be used to specify the number of columns +# in which this list will be split (can be a number in the range [1..20]) + +COLS_IN_ALPHA_INDEX = 5 + +# In case all classes in a project start with a common prefix, all +# classes will be put under the same header in the alphabetical index. +# The IGNORE_PREFIX tag can be used to specify one or more prefixes that +# should be ignored while generating the index headers. + +IGNORE_PREFIX = + +#--------------------------------------------------------------------------- +# configuration options related to the HTML output +#--------------------------------------------------------------------------- + +# If the GENERATE_HTML tag is set to YES (the default) Doxygen will +# generate HTML output. + +GENERATE_HTML = YES + +# The HTML_OUTPUT tag is used to specify where the HTML docs will be put. +# If a relative path is entered the value of OUTPUT_DIRECTORY will be +# put in front of it. If left blank `html' will be used as the default path. + +HTML_OUTPUT = html + +# The HTML_FILE_EXTENSION tag can be used to specify the file extension for +# each generated HTML page (for example: .htm,.php,.asp). If it is left blank +# doxygen will generate files with .html extension. + +HTML_FILE_EXTENSION = .html + +# The HTML_HEADER tag can be used to specify a personal HTML header for +# each generated HTML page. If it is left blank doxygen will generate a +# standard header. + +HTML_HEADER = + +# The HTML_FOOTER tag can be used to specify a personal HTML footer for +# each generated HTML page. If it is left blank doxygen will generate a +# standard footer. + +HTML_FOOTER = + +# The HTML_STYLESHEET tag can be used to specify a user-defined cascading +# style sheet that is used by each HTML page. It can be used to +# fine-tune the look of the HTML output. If the tag is left blank doxygen +# will generate a default style sheet. Note that doxygen will try to copy +# the style sheet file to the HTML output directory, so don't put your own +# stylesheet in the HTML output directory as well, or it will be erased! + +HTML_STYLESHEET = + +# If the HTML_ALIGN_MEMBERS tag is set to YES, the members of classes, +# files or namespaces will be aligned in HTML using tables. If set to +# NO a bullet list will be used. + +HTML_ALIGN_MEMBERS = YES + +# If the GENERATE_HTMLHELP tag is set to YES, additional index files +# will be generated that can be used as input for tools like the +# Microsoft HTML help workshop to generate a compressed HTML help file (.chm) +# of the generated HTML documentation. + +GENERATE_HTMLHELP = NO + +# If the HTML_DYNAMIC_SECTIONS tag is set to YES then the generated HTML +# documentation will contain sections that can be hidden and shown after the +# page has loaded. For this to work a browser that supports +# JavaScript and DHTML is required (for instance Mozilla 1.0+, Firefox +# Netscape 6.0+, Internet explorer 5.0+, Konqueror, or Safari). + +HTML_DYNAMIC_SECTIONS = NO + +# If the GENERATE_HTMLHELP tag is set to YES, the CHM_FILE tag can +# be used to specify the file name of the resulting .chm file. You +# can add a path in front of the file if the result should not be +# written to the html output directory. + +CHM_FILE = + +# If the GENERATE_HTMLHELP tag is set to YES, the HHC_LOCATION tag can +# be used to specify the location (absolute path including file name) of +# the HTML help compiler (hhc.exe). If non-empty doxygen will try to run +# the HTML help compiler on the generated index.hhp. + +HHC_LOCATION = + +# If the GENERATE_HTMLHELP tag is set to YES, the GENERATE_CHI flag +# controls if a separate .chi index file is generated (YES) or that +# it should be included in the master .chm file (NO). + +GENERATE_CHI = NO + +# If the GENERATE_HTMLHELP tag is set to YES, the BINARY_TOC flag +# controls whether a binary table of contents is generated (YES) or a +# normal table of contents (NO) in the .chm file. + +BINARY_TOC = NO + +# The TOC_EXPAND flag can be set to YES to add extra items for group members +# to the contents of the HTML help documentation and to the tree view. + +TOC_EXPAND = NO + +# The DISABLE_INDEX tag can be used to turn on/off the condensed index at +# top of each HTML page. The value NO (the default) enables the index and +# the value YES disables it. + +DISABLE_INDEX = NO + +# This tag can be used to set the number of enum values (range [1..20]) +# that doxygen will group on one line in the generated HTML documentation. + +ENUM_VALUES_PER_LINE = 4 + +# If the GENERATE_TREEVIEW tag is set to YES, a side panel will be +# generated containing a tree-like index structure (just like the one that +# is generated for HTML Help). For this to work a browser that supports +# JavaScript, DHTML, CSS and frames is required (for instance Mozilla 1.0+, +# Netscape 6.0+, Internet explorer 5.0+, or Konqueror). Windows users are +# probably better off using the HTML help feature. + +GENERATE_TREEVIEW = NO + +# If the treeview is enabled (see GENERATE_TREEVIEW) then this tag can be +# used to set the initial width (in pixels) of the frame in which the tree +# is shown. + +TREEVIEW_WIDTH = 250 + +#--------------------------------------------------------------------------- +# configuration options related to the LaTeX output +#--------------------------------------------------------------------------- + +# If the GENERATE_LATEX tag is set to YES (the default) Doxygen will +# generate Latex output. + +GENERATE_LATEX = NO + +# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put. +# If a relative path is entered the value of OUTPUT_DIRECTORY will be +# put in front of it. If left blank `latex' will be used as the default path. + +LATEX_OUTPUT = latex + +# The LATEX_CMD_NAME tag can be used to specify the LaTeX command name to be +# invoked. If left blank `latex' will be used as the default command name. + +LATEX_CMD_NAME = latex + +# The MAKEINDEX_CMD_NAME tag can be used to specify the command name to +# generate index for LaTeX. If left blank `makeindex' will be used as the +# default command name. + +MAKEINDEX_CMD_NAME = makeindex + +# If the COMPACT_LATEX tag is set to YES Doxygen generates more compact +# LaTeX documents. This may be useful for small projects and may help to +# save some trees in general. + +COMPACT_LATEX = NO + +# The PAPER_TYPE tag can be used to set the paper type that is used +# by the printer. Possible values are: a4, a4wide, letter, legal and +# executive. If left blank a4wide will be used. + +PAPER_TYPE = a4wide + +# The EXTRA_PACKAGES tag can be to specify one or more names of LaTeX +# packages that should be included in the LaTeX output. + +EXTRA_PACKAGES = + +# The LATEX_HEADER tag can be used to specify a personal LaTeX header for +# the generated latex document. The header should contain everything until +# the first chapter. If it is left blank doxygen will generate a +# standard header. Notice: only use this tag if you know what you are doing! + +LATEX_HEADER = + +# If the PDF_HYPERLINKS tag is set to YES, the LaTeX that is generated +# is prepared for conversion to pdf (using ps2pdf). The pdf file will +# contain links (just like the HTML output) instead of page references +# This makes the output suitable for online browsing using a pdf viewer. + +PDF_HYPERLINKS = NO + +# If the USE_PDFLATEX tag is set to YES, pdflatex will be used instead of +# plain latex in the generated Makefile. Set this option to YES to get a +# higher quality PDF documentation. + +USE_PDFLATEX = NO + +# If the LATEX_BATCHMODE tag is set to YES, doxygen will add the \\batchmode. +# command to the generated LaTeX files. This will instruct LaTeX to keep +# running if errors occur, instead of asking the user for help. +# This option is also used when generating formulas in HTML. + +LATEX_BATCHMODE = NO + +# If LATEX_HIDE_INDICES is set to YES then doxygen will not +# include the index chapters (such as File Index, Compound Index, etc.) +# in the output. + +LATEX_HIDE_INDICES = NO + +#--------------------------------------------------------------------------- +# configuration options related to the RTF output +#--------------------------------------------------------------------------- + +# If the GENERATE_RTF tag is set to YES Doxygen will generate RTF output +# The RTF output is optimized for Word 97 and may not look very pretty with +# other RTF readers or editors. + +GENERATE_RTF = NO + +# The RTF_OUTPUT tag is used to specify where the RTF docs will be put. +# If a relative path is entered the value of OUTPUT_DIRECTORY will be +# put in front of it. If left blank `rtf' will be used as the default path. + +RTF_OUTPUT = rtf + +# If the COMPACT_RTF tag is set to YES Doxygen generates more compact +# RTF documents. This may be useful for small projects and may help to +# save some trees in general. + +COMPACT_RTF = NO + +# If the RTF_HYPERLINKS tag is set to YES, the RTF that is generated +# will contain hyperlink fields. The RTF file will +# contain links (just like the HTML output) instead of page references. +# This makes the output suitable for online browsing using WORD or other +# programs which support those fields. +# Note: wordpad (write) and others do not support links. + +RTF_HYPERLINKS = NO + +# Load stylesheet definitions from file. Syntax is similar to doxygen's +# config file, i.e. a series of assignments. You only have to provide +# replacements, missing definitions are set to their default value. + +RTF_STYLESHEET_FILE = + +# Set optional variables used in the generation of an rtf document. +# Syntax is similar to doxygen's config file. + +RTF_EXTENSIONS_FILE = + +#--------------------------------------------------------------------------- +# configuration options related to the man page output +#--------------------------------------------------------------------------- + +# If the GENERATE_MAN tag is set to YES (the default) Doxygen will +# generate man pages + +GENERATE_MAN = YES + +# The MAN_OUTPUT tag is used to specify where the man pages will be put. +# If a relative path is entered the value of OUTPUT_DIRECTORY will be +# put in front of it. If left blank `man' will be used as the default path. + +MAN_OUTPUT = man + +# The MAN_EXTENSION tag determines the extension that is added to +# the generated man pages (default is the subroutine's section .3) + +MAN_EXTENSION = .3 + +# If the MAN_LINKS tag is set to YES and Doxygen generates man output, +# then it will generate one additional man file for each entity +# documented in the real man page(s). These additional files +# only source the real man page, but without them the man command +# would be unable to find the correct page. The default is NO. + +MAN_LINKS = NO + +#--------------------------------------------------------------------------- +# configuration options related to the XML output +#--------------------------------------------------------------------------- + +# If the GENERATE_XML tag is set to YES Doxygen will +# generate an XML file that captures the structure of +# the code including all documentation. + +GENERATE_XML = NO + +# The XML_OUTPUT tag is used to specify where the XML pages will be put. +# If a relative path is entered the value of OUTPUT_DIRECTORY will be +# put in front of it. If left blank `xml' will be used as the default path. + +XML_OUTPUT = xml + +# The XML_SCHEMA tag can be used to specify an XML schema, +# which can be used by a validating XML parser to check the +# syntax of the XML files. + +XML_SCHEMA = + +# The XML_DTD tag can be used to specify an XML DTD, +# which can be used by a validating XML parser to check the +# syntax of the XML files. + +XML_DTD = + +# If the XML_PROGRAMLISTING tag is set to YES Doxygen will +# dump the program listings (including syntax highlighting +# and cross-referencing information) to the XML output. Note that +# enabling this will significantly increase the size of the XML output. + +XML_PROGRAMLISTING = YES + +#--------------------------------------------------------------------------- +# configuration options for the AutoGen Definitions output +#--------------------------------------------------------------------------- + +# If the GENERATE_AUTOGEN_DEF tag is set to YES Doxygen will +# generate an AutoGen Definitions (see autogen.sf.net) file +# that captures the structure of the code including all +# documentation. Note that this feature is still experimental +# and incomplete at the moment. + +GENERATE_AUTOGEN_DEF = NO + +#--------------------------------------------------------------------------- +# configuration options related to the Perl module output +#--------------------------------------------------------------------------- + +# If the GENERATE_PERLMOD tag is set to YES Doxygen will +# generate a Perl module file that captures the structure of +# the code including all documentation. Note that this +# feature is still experimental and incomplete at the +# moment. + +GENERATE_PERLMOD = NO + +# If the PERLMOD_LATEX tag is set to YES Doxygen will generate +# the necessary Makefile rules, Perl scripts and LaTeX code to be able +# to generate PDF and DVI output from the Perl module output. + +PERLMOD_LATEX = NO + +# If the PERLMOD_PRETTY tag is set to YES the Perl module output will be +# nicely formatted so it can be parsed by a human reader. This is useful +# if you want to understand what is going on. On the other hand, if this +# tag is set to NO the size of the Perl module output will be much smaller +# and Perl will parse it just the same. + +PERLMOD_PRETTY = YES + +# The names of the make variables in the generated doxyrules.make file +# are prefixed with the string contained in PERLMOD_MAKEVAR_PREFIX. +# This is useful so different doxyrules.make files included by the same +# Makefile don't overwrite each other's variables. + +PERLMOD_MAKEVAR_PREFIX = + +#--------------------------------------------------------------------------- +# Configuration options related to the preprocessor +#--------------------------------------------------------------------------- + +# If the ENABLE_PREPROCESSING tag is set to YES (the default) Doxygen will +# evaluate all C-preprocessor directives found in the sources and include +# files. + +ENABLE_PREPROCESSING = YES + +# If the MACRO_EXPANSION tag is set to YES Doxygen will expand all macro +# names in the source code. If set to NO (the default) only conditional +# compilation will be performed. Macro expansion can be done in a controlled +# way by setting EXPAND_ONLY_PREDEF to YES. + +MACRO_EXPANSION = NO + +# If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES +# then the macro expansion is limited to the macros specified with the +# PREDEFINED and EXPAND_AS_DEFINED tags. + +EXPAND_ONLY_PREDEF = NO + +# If the SEARCH_INCLUDES tag is set to YES (the default) the includes files +# in the INCLUDE_PATH (see below) will be search if a #include is found. + +SEARCH_INCLUDES = YES + +# The INCLUDE_PATH tag can be used to specify one or more directories that +# contain include files that are not input files but should be processed by +# the preprocessor. + +INCLUDE_PATH = + +# You can use the INCLUDE_FILE_PATTERNS tag to specify one or more wildcard +# patterns (like *.h and *.hpp) to filter out the header-files in the +# directories. If left blank, the patterns specified with FILE_PATTERNS will +# be used. + +INCLUDE_FILE_PATTERNS = + +# The PREDEFINED tag can be used to specify one or more macro names that +# are defined before the preprocessor is started (similar to the -D option of +# gcc). The argument of the tag is a list of macros of the form: name +# or name=definition (no spaces). If the definition and the = are +# omitted =1 is assumed. To prevent a macro definition from being +# undefined via #undef or recursively expanded use the := operator +# instead of the = operator. + +PREDEFINED = + +# If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then +# this tag can be used to specify a list of macro names that should be expanded. +# The macro definition that is found in the sources will be used. +# Use the PREDEFINED tag if you want to use a different macro definition. + +EXPAND_AS_DEFINED = + +# If the SKIP_FUNCTION_MACROS tag is set to YES (the default) then +# doxygen's preprocessor will remove all function-like macros that are alone +# on a line, have an all uppercase name, and do not end with a semicolon. Such +# function macros are typically used for boiler-plate code, and will confuse +# the parser if not removed. + +SKIP_FUNCTION_MACROS = YES + +#--------------------------------------------------------------------------- +# Configuration::additions related to external references +#--------------------------------------------------------------------------- + +# The TAGFILES option can be used to specify one or more tagfiles. +# Optionally an initial location of the external documentation +# can be added for each tagfile. The format of a tag file without +# this location is as follows: +# TAGFILES = file1 file2 ... +# Adding location for the tag files is done as follows: +# TAGFILES = file1=loc1 "file2 = loc2" ... +# where "loc1" and "loc2" can be relative or absolute paths or +# URLs. If a location is present for each tag, the installdox tool +# does not have to be run to correct the links. +# Note that each tag file must have a unique name +# (where the name does NOT include the path) +# If a tag file is not located in the directory in which doxygen +# is run, you must also specify the path to the tagfile here. + +TAGFILES = + +# When a file name is specified after GENERATE_TAGFILE, doxygen will create +# a tag file that is based on the input files it reads. + +GENERATE_TAGFILE = + +# If the ALLEXTERNALS tag is set to YES all external classes will be listed +# in the class index. If set to NO only the inherited external classes +# will be listed. + +ALLEXTERNALS = NO + +# If the EXTERNAL_GROUPS tag is set to YES all external groups will be listed +# in the modules index. If set to NO, only the current project's groups will +# be listed. + +EXTERNAL_GROUPS = YES + +# The PERL_PATH should be the absolute path and name of the perl script +# interpreter (i.e. the result of `which perl'). + +PERL_PATH = /usr/bin/perl + +#--------------------------------------------------------------------------- +# Configuration options related to the dot tool +#--------------------------------------------------------------------------- + +# If the CLASS_DIAGRAMS tag is set to YES (the default) Doxygen will +# generate a inheritance diagram (in HTML, RTF and LaTeX) for classes with base +# or super classes. Setting the tag to NO turns the diagrams off. Note that +# this option is superseded by the HAVE_DOT option below. This is only a +# fallback. It is recommended to install and use dot, since it yields more +# powerful graphs. + +CLASS_DIAGRAMS = YES + +# You can define message sequence charts within doxygen comments using the \msc +# command. Doxygen will then run the mscgen tool (see http://www.mcternan.me.uk/mscgen/) to +# produce the chart and insert it in the documentation. The MSCGEN_PATH tag allows you to +# specify the directory where the mscgen tool resides. If left empty the tool is assumed to +# be found in the default search path. + +MSCGEN_PATH = + +# If set to YES, the inheritance and collaboration graphs will hide +# inheritance and usage relations if the target is undocumented +# or is not a class. + +HIDE_UNDOC_RELATIONS = YES + +# If you set the HAVE_DOT tag to YES then doxygen will assume the dot tool is +# available from the path. This tool is part of Graphviz, a graph visualization +# toolkit from AT&T and Lucent Bell Labs. The other options in this section +# have no effect if this option is set to NO (the default) + +HAVE_DOT = NO + +# If the CLASS_GRAPH and HAVE_DOT tags are set to YES then doxygen +# will generate a graph for each documented class showing the direct and +# indirect inheritance relations. Setting this tag to YES will force the +# the CLASS_DIAGRAMS tag to NO. + +CLASS_GRAPH = YES + +# If the COLLABORATION_GRAPH and HAVE_DOT tags are set to YES then doxygen +# will generate a graph for each documented class showing the direct and +# indirect implementation dependencies (inheritance, containment, and +# class references variables) of the class with other documented classes. + +COLLABORATION_GRAPH = YES + +# If the GROUP_GRAPHS and HAVE_DOT tags are set to YES then doxygen +# will generate a graph for groups, showing the direct groups dependencies + +GROUP_GRAPHS = YES + +# If the UML_LOOK tag is set to YES doxygen will generate inheritance and +# collaboration diagrams in a style similar to the OMG's Unified Modeling +# Language. + +UML_LOOK = NO + +# If set to YES, the inheritance and collaboration graphs will show the +# relations between templates and their instances. + +TEMPLATE_RELATIONS = NO + +# If the ENABLE_PREPROCESSING, SEARCH_INCLUDES, INCLUDE_GRAPH, and HAVE_DOT +# tags are set to YES then doxygen will generate a graph for each documented +# file showing the direct and indirect include dependencies of the file with +# other documented files. + +INCLUDE_GRAPH = YES + +# If the ENABLE_PREPROCESSING, SEARCH_INCLUDES, INCLUDED_BY_GRAPH, and +# HAVE_DOT tags are set to YES then doxygen will generate a graph for each +# documented header file showing the documented files that directly or +# indirectly include this file. + +INCLUDED_BY_GRAPH = YES + +# If the CALL_GRAPH, SOURCE_BROWSER and HAVE_DOT tags are set to YES then doxygen will +# generate a call dependency graph for every global function or class method. +# Note that enabling this option will significantly increase the time of a run. +# So in most cases it will be better to enable call graphs for selected +# functions only using the \callgraph command. + +CALL_GRAPH = NO + +# If the CALLER_GRAPH, SOURCE_BROWSER and HAVE_DOT tags are set to YES then doxygen will +# generate a caller dependency graph for every global function or class method. +# Note that enabling this option will significantly increase the time of a run. +# So in most cases it will be better to enable caller graphs for selected +# functions only using the \callergraph command. + +CALLER_GRAPH = NO + +# If the GRAPHICAL_HIERARCHY and HAVE_DOT tags are set to YES then doxygen +# will graphical hierarchy of all classes instead of a textual one. + +GRAPHICAL_HIERARCHY = YES + +# If the DIRECTORY_GRAPH, SHOW_DIRECTORIES and HAVE_DOT tags are set to YES +# then doxygen will show the dependencies a directory has on other directories +# in a graphical way. The dependency relations are determined by the #include +# relations between the files in the directories. + +DIRECTORY_GRAPH = YES + +# The DOT_IMAGE_FORMAT tag can be used to set the image format of the images +# generated by dot. Possible values are png, jpg, or gif +# If left blank png will be used. + +DOT_IMAGE_FORMAT = png + +# The tag DOT_PATH can be used to specify the path where the dot tool can be +# found. If left blank, it is assumed the dot tool can be found in the path. + +DOT_PATH = + +# The DOTFILE_DIRS tag can be used to specify one or more directories that +# contain dot files that are included in the documentation (see the +# \dotfile command). + +DOTFILE_DIRS = + +# The MAX_DOT_GRAPH_MAX_NODES tag can be used to set the maximum number of +# nodes that will be shown in the graph. If the number of nodes in a graph +# becomes larger than this value, doxygen will truncate the graph, which is +# visualized by representing a node as a red box. Note that doxygen if the number +# of direct children of the root node in a graph is already larger than +# MAX_DOT_GRAPH_NOTES then the graph will not be shown at all. Also note +# that the size of a graph can be further restricted by MAX_DOT_GRAPH_DEPTH. + +DOT_GRAPH_MAX_NODES = 50 + +# The MAX_DOT_GRAPH_DEPTH tag can be used to set the maximum depth of the +# graphs generated by dot. A depth value of 3 means that only nodes reachable +# from the root by following a path via at most 3 edges will be shown. Nodes +# that lay further from the root node will be omitted. Note that setting this +# option to 1 or 2 may greatly reduce the computation time needed for large +# code bases. Also note that the size of a graph can be further restricted by +# DOT_GRAPH_MAX_NODES. Using a depth of 0 means no depth restriction. + +MAX_DOT_GRAPH_DEPTH = 0 + +# Set the DOT_TRANSPARENT tag to YES to generate images with a transparent +# background. This is disabled by default, which results in a white background. +# Warning: Depending on the platform used, enabling this option may lead to +# badly anti-aliased labels on the edges of a graph (i.e. they become hard to +# read). + +DOT_TRANSPARENT = NO + +# Set the DOT_MULTI_TARGETS tag to YES allow dot to generate multiple output +# files in one run (i.e. multiple -o and -T options on the command line). This +# makes dot run faster, but since only newer versions of dot (>1.8.10) +# support this, this feature is disabled by default. + +DOT_MULTI_TARGETS = NO + +# If the GENERATE_LEGEND tag is set to YES (the default) Doxygen will +# generate a legend page explaining the meaning of the various boxes and +# arrows in the dot generated graphs. + +GENERATE_LEGEND = YES + +# If the DOT_CLEANUP tag is set to YES (the default) Doxygen will +# remove the intermediate dot files that are used to generate +# the various graphs. + +DOT_CLEANUP = YES + +#--------------------------------------------------------------------------- +# Configuration::additions related to the search engine +#--------------------------------------------------------------------------- + +# The SEARCHENGINE tag specifies whether or not a search engine should be +# used. If set to NO the values of all tags below this one will be ignored. + +SEARCHENGINE = NO Index: /branches/ccl_branches/ipponly-20191108/ppSub/Makefile.am =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/Makefile.am (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/Makefile.am (revision 41094) @@ -0,0 +1,3 @@ +SUBDIRS = src + +CLEANFILES = *.pyc *~ core core.* Index: /branches/ccl_branches/ipponly-20191108/ppSub/TO_DO =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/TO_DO (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/TO_DO (revision 41094) @@ -0,0 +1,177 @@ +Some ideas for improving the image subtraction +============================================== + +Paul Price +15 August 2007 + + +1. Dual convolution + +The Alard & Lupton (1998) and Alard (2000) algorithms are unable to +produce useful subtractions when the PSFs of the input and reference +images have misaligned position angles (e.g., / vs \), since this +would involve a deconvolution along one of the axes, which does not +work well. One way to solve this problem is through a dual +convolution. Instead of solving + + I(x,y) = k(u,v) * R(x,y) + +for k(u,v), we can solve + + k1(u,v) * I(x,y) = k2(u,v) * R(x,y) + +for k1(u,v) and k2(u,v). If we write + + k1(u,v) = sum_i a_i f_i(u,v) + +and + + k2(u,v) = sum_i b_i g_i(u,v) + +where f_i(u,v) and g_i(u,v) are sets of basis functions, then the +solution boils down to two equations: + + sum_j a_j A_j A_i = sum_j b_j B_j A_i + +and + + sum_j a_j A_j B_i = sum_j b_j B_j B_i + +where + + A_i = sum_x,y f_i(u,v) * I(x,y) / sigma(x,y) + +and + + B_i = sum_x,y g_i(u,v) * R(x,y) / sigma(x,y) + +Noting that the matrix terms of the RHS of the first equation and the +LHS of the second equation (A_i B_j) are the same leads us to attempt +to solve this system by iteration: + + 1. Set a_i = 1 if i = 0, otherwise a_i = 0; f_0 = delta(u,v); and + solve for b in the first equation. This corresponds to doing the + usual Alard & Lupton solution. + + 2. Use b in the second equation, and solve for a. + + 3. Use a in the first equation, and solve for b. + + 4. Proceed in this manner until the change in b between iterations + is small. + +This is basically doing A&L to get the reference image to match the +input image in the usual manner, then doing A&L on the input image to +match the convolved reference image, then doing A&L on the reference +image to match the convolved input image, etc., until the system +settles down to the solution. A more refined method of solving the +equations may exist, but this should work. + +There is a need to normalise one of the kernels: the solution is +currently non-unique because a and b can be scaled by some arbitrary +value. If we write + + k1(u,v) = delta(u,v) + sum_i a_i [k_i - delta(u,v)] + +and the k_i are normalised to a sum of unity, then + + sum_u,v k1(u,v) = 1. + +So k1 has a sum of unity, always. This means that k2 supplies the +scaling to match the photometry (as happens in Alard & Lupton), and k1 +is merely used to broaden the input image as required for the best +subtraction. + +Note that the expense for this method over the usual Alard & Lupton +roughly amounts to the extra image convolution, since that is +typically the dominant factor. Accumulation of the sums is not much +more than Alard & Lupton, and the iteration should be fairly fast. + + +2. Data-based kernel selection + +The quality of the subtraction is highly sensitive to the choice of +basis functions. In the case of ISIS kernels (which seem to be the +most useful of those experimented with so far, because it takes a +small number of parameters to generate a large kernel), the choice of +the Gaussian widths is very important. Some recipes exist for +choosing these widths, but these are generally motivated by experience +(through much trial and error) rather than directly from the data. It +would be nice to be able to throw down a large number of widths and +allow the least-squares solution to choose the best for the data at +hand (i.e., the most dominant contributors), but this is prohibitively +expensive --- at least for a full solution. + +What might be possible is to do a quick and dirty solution by reducing +the dimensionality. Instead of working with the full two-dimensional +kernel, k(u,v), let's work in one dimension, k(u). For each of our +stamps, let's take a cut through them in a consistent direction (e.g., +the x direction), and solve + + I(x) = k(u) * R(x) + +This is not nearly as computationally expensive as the full solution, +so we can pack k(u) with multiple Gaussian widths, solve the +least-squares problem, and identify the most important kernel +contributions which we will use (in suitable two-dimensional versions) +in solving the full problem. This method could even be applied with +multiple direction cuts to ensure the full range of required Gaussian +widths is obtained. + +Fitting the normalisation of a single kernel component, k(u) and the +background, we obtain the least-squares matrix and vector: + +M = ( sum_x C(x)^2/sigma(x)^2 sum_x C(x)/sigma(x)^2 ) + ( sum_x C(x)/sigma(x)^2 sum_x 1/sigma(x)^2 ) + +v = ( sum_x I(x)C(x)/sigma(x)^2 sum_x I(x)/sigma(x)^2 ) + +where C(x) = R(x) * k(u) + +Inverting the matrix, multiplying by the vector and taking the +component corresponding to the kernel normalisation, we get: + +f = (ab - cd) / (ae - c^2) + +where + +a = sum_x 1/sigma(x)^2 +b = sum_x I(x)C(x)/sigma(x)^2 +c = sum_x C(x)/sigma(x)^2 +d = sum_x I(x)/simga(x)^2 +e = sum_x C(x)^2/sigma(x)^2 + +Note that a and d are independent of the kernel, and so may be +measured once only. Note that we are ignoring any spatial variation +of the kernel here --- we are only interested in which kernel +components should be in the final solution, not the details of that +solution. + +Assuming that the kernel components are normalised, the kernel +component normalisation, f, is a measure of how important that kernel +component is in the final solution of the full problem. The following +algorithm is suggested: + +(a) For each original kernel component, extract a vector at 0, 45 and + 90 degrees (x, y and x-y axes). Ignore one of these sub-kernels + in what follows if the standard deviation of the subkernel is + zero. We use these three extraction angles to cover the u, v and + uv terms. + +(b) For each stamp, extract a vector at 0, 45, and 90 degrees (x, y, + and x-y axes). These sub-stamps go with the corresponding + sub-kernel. + +(c) For each kernel component, measure b, c, e and therefore determine + the kernel normalisation, f. Use the sub-kernel appropriate for + the sub-stamp (matched by extraction angle), and accumulate the + values together (sub-kernels all contribute to the same original + kernel component). + +(d) Take the kernel component with the largest |f| as the 'winner' of + this iteration. Apply this kernel to the sub-stamps for R(x), and + mark this kernel component in the list. It is part of the end + solution, but should not be used any more in the iteration. + +(e) Repeat from step (c) until the largest |f| is small, say 10^-3 of + the first obtained value of |f|. Index: /branches/ccl_branches/ipponly-20191108/ppSub/autogen.sh =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/autogen.sh (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/autogen.sh (revision 41094) @@ -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=ppSub +TEST_TYPE=-f +# change this to be a unique filename in the top level dir +FILE=autogen.sh + +DIE=0 + +LIBTOOLIZE=libtoolize +ACLOCAL="aclocal $ACLOCAL_FLAGS" +AUTOHEADER=autoheader +AUTOMAKE=automake +AUTOCONF=autoconf + +($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $LIBTOOLIZE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/" + DIE=1 +} + +($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $ACLOCAL installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOHEADER installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOMAKE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOCONF installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +if test "$DIE" -eq 1; then + exit 1 +fi + +test $TEST_TYPE $FILE || { + echo "You must run this script in the top-level $PROJECT directory" + exit 1 +} + +if test -z "$*"; then + echo "I am going to run ./configure with no arguments - if you wish " + echo "to pass any to it, please specify them on the $0 command line." +fi + +$LIBTOOLIZE --copy --force || echo "$LIBTOOLIZE failed" +$ACLOCAL || echo "$ACLOCAL failed" +$AUTOHEADER || echo "$AUTOHEADER failed" +$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE failed" +$AUTOCONF || echo "$AUTOCONF failed" + +cd $ORIGDIR + +run_configure=true +for arg in $*; do + case $arg in + --no-configure) + run_configure=false + ;; + *) + ;; + esac +done + +if $run_configure; then + $srcdir/configure --enable-maintainer-mode "$@" + echo + echo "Now type 'make' to compile $PROJECT." +else + echo + echo "Now run 'configure' and 'make' to compile $PROJECT." +fi Index: /branches/ccl_branches/ipponly-20191108/ppSub/configure.ac =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/configure.ac (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/configure.ac (revision 41094) @@ -0,0 +1,45 @@ +dnl Process this file with autoconf to produce a configure script. +AC_PREREQ(2.61) + +AC_INIT([ppSub], [0.1.1], [ipp-support@ifa.hawaii.edu]) +AC_CONFIG_SRCDIR([src]) + +AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2]) +AM_CONFIG_HEADER([src/config.h]) +AM_MAINTAINER_MODE + +IPP_STDLDFLAGS + +AC_LANG(C) +AC_GNU_SOURCE +AC_PROG_CC_C99 +AC_PROG_INSTALL +AC_PROG_LIBTOOL +AC_SYS_LARGEFILE + +PKG_CHECK_MODULES([PSLIB], [pslib >= 1.0.0]) +PKG_CHECK_MODULES([PSMODULE], [psmodules >= 1.0.0]) +PKG_CHECK_MODULES([PPSTATS], [ppStats >= 1.0.0]) +PKG_CHECK_MODULES([PSPHOT], [psphot >= 0.9.0]) + +AC_PATH_PROG([ERRORCODES], [psParseErrorCodes], [missing]) +if test "$ERRORCODES" = "missing" ; then + AC_MSG_ERROR([psParseErrorCodes is required]) +fi + +dnl Set CFLAGS for build +IPP_STDOPTS +IPP_STDCFLAGS + +IPP_VERSION + +AC_SUBST([PPSUB_CFLAGS]) +AC_SUBST([PPSUB_LIBS]) + +AC_CONFIG_FILES([ + Makefile + src/Makefile + Doxyfile +]) + +AC_OUTPUT Index: /branches/ccl_branches/ipponly-20191108/ppSub/doc/outline.txt =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/doc/outline.txt (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/doc/outline.txt (revision 41094) @@ -0,0 +1,26 @@ + +bool ppSubMatchPSFs(ppSubData *data) + - setup huge number of recipe values + - load or find sources + - pmSubtractionMatch + - data normalization + +bool pmSubtractionMatch() + - setup + - loop over iso-kernel regions + - identify stamp sources + - extract stamp pixels + - generate kernel basis functions + - iterate to reject stamps + - pmSubtractionCalculateEquation + - pmSubtractionSolveEquation + - pmSubtractionCalculateDeviations + - pmSubtractionRejectStamps + +bool pmSubtractionCalculateEquationStamp() + - convolve stamp with each kernel basis function + - generate the polynomial terms for the basis functions + - sum the elements of the chisq equation + +bool pmSubtractionSolveEquation() + Index: /branches/ccl_branches/ipponly-20191108/ppSub/doc/psphot.txt =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/doc/psphot.txt (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/doc/psphot.txt (revision 41094) @@ -0,0 +1,36 @@ + +Notes on psphot function calls within ppsub: + +we have these psphotReadout varients: + +* psphotReadout : ppImagePhotom.c, ppSimPhotom*, +* psphotReadoutFindPSF : used by ppSubMatchPSFs.c, ppSubMakePSF.c, pswarpLoop.c +* psphotReadoutKnownSources : used by ppStackPhotometry.c +* psphotReadoutMinimal : used by ppSubReadoutPhotometry.c + +* psphotStackReadout : used only by psphotStack.c +* psphotForcedReadout : used only by psphotForced.c +* psphotMakePSFReadout : used only by psphotMakePSF.c + +psphotReadout: + * operates on pmFPAfile "PSPHOT.INPUT" -> convert to an input argument (filerule) + * expects a readout at the current view, populated with a image (can create mask and variance) + * XXX what about the headers? + + * adds a photcode to readout->analysis + * creates a mask and optionally modifies the variance (NOTE!) + * if needed, generates PSPHOT.BACKMDL & PSPHOT.BACKMDL.STDEV pmFPAfiles (internal) + * optionally copies PSF model from PSPHOT.PSF.LOAD to PSPHOT.INPUT + * detetections / sources saved on readout->analysis as PSPHOT.DETECTIONS + * result psf goes on readout->analyss as PSPHOT.PSF + * detection efficiency data goes on readout->analyss as DETEFF + * copy input sources from PSPHOT.SOURCES.CMF & PSPHOT.SOURCES.TXT -> PSPHOT.DETECTIONS + * header data goes into PSPHOT.HEADER + * psf clump stuff goes on readout->analysis + + + +--- + +in ppSubMatchPSFs.c : we call psphotReadoutFindPSF to measure the PSF using sources already loaded + (PPSUB.INPUT & PPSUB.REF) Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/Makefile.am =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/Makefile.am (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/Makefile.am (revision 41094) @@ -0,0 +1,74 @@ +bin_PROGRAMS = ppSub ppSubKernel ppSubConvolve + +# Force recompilation of ppSubVersion.c, since it gets the version information +ppSubVersion.c: ppSubVersionDefinitions.h +ppSubVersionDefinitions.h: ppSubVersionDefinitions.h.in FORCE + pslib-setsvnversion.pl PPSUB ppSubVersionDefinitions.h.in ppSubVersionDefinitions.h +FORCE: ; + +ppSub_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSTATS_CFLAGS) $(PSPHOT_CFLAGS) $(PPSUB_CFLAGS) +ppSub_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PPSTATS_LIBS) $(PSPHOT_LIBS) $(PPSUB_LIBS) + +ppSub_SOURCES = \ + ppSub.c \ + ppSubArguments.c \ + ppSubVersion.c \ + ppSubBackground.c \ + ppSubVarianceRescale.c \ + ppSubCamera.c \ + ppSubData.c \ + ppSubErrorCodes.c \ + ppSubExit.c \ + ppSubFiles.c \ + ppSubLoop.c \ + ppSubDefineOutput.c \ + ppSubInputDetections.c \ + ppSubExtras.c \ + ppSubFlagNeighbors.c \ + ppSubMakePSF.c \ + ppSubMatchPSFs.c \ + ppSubMaskStats.c \ + ppSubReadoutInverse.c \ + ppSubReadoutJpeg.c \ + ppSubReadoutPhotometry.c \ + ppSubReadoutForcedPhot.c \ + ppSubReadoutStats.c \ + ppSubReadoutSubtract.c \ + ppSubSetMasks.c \ + ppSubThreshold.c + +ppSubKernel_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSUB_CFLAGS) +ppSubKernel_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PPSUB_LIBS) + +ppSubKernel_SOURCES = \ + ppSubKernel.c + +ppSubConvolve_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSTATS_CFLAGS) $(PSPHOT_CFLAGS) $(PPSUB_CFLAGS) +ppSubConvolve_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PPSTATS_LIBS) $(PSPHOT_LIBS) $(PPSUB_LIBS) + +ppSubConvolve_SOURCES = \ + ppSubConvolve.c \ + ppSubExit.c \ + ppSubVersion.c + +noinst_HEADERS = \ + ppSub.h + + +### Error codes. +BUILT_SOURCES = ppSubErrorCodes.h ppSubErrorCodes.c +CLEANFILES = ppSubErrorCodes.h ppSubErrorCodes.c + +ppSubErrorCodes.h : ppSubErrorCodes.dat ppSubErrorCodes.h.in + $(ERRORCODES) --data=ppSubErrorCodes.dat --outdir=. ppSubErrorCodes.h + +ppSubErrorCodes.c : ppSubErrorCodes.dat ppSubErrorCodes.c.in ppSubErrorCodes.h + $(ERRORCODES) --data=ppSubErrorCodes.dat --outdir=. ppSubErrorCodes.c + + +clean-local: + -rm -f TAGS + +# Tags for emacs +tags: + etags `find . -name \*.[ch] -print` Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSub.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSub.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSub.c (revision 41094) @@ -0,0 +1,120 @@ +/** @file ppSub.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +int main(int argc, char *argv[]) +{ + psLibInit(NULL); + psTimerStart("ppSub"); + + pmErrorRegister(); + ppSubErrorRegister(); + psphotErrorRegister(); + + ppSubData *data = NULL; // Processing data + pmConfig *config = pmConfigRead(&argc, argv, PPSUB_RECIPE); // Configuration + if (!config) { + goto die; + } + + ppSubVersionPrint(); + + if (!pmModelClassInit()) { + psError(PPSUB_ERR_PROG, false, "Unable to initialise model classes."); + psFree(config); + goto die; + } + + if (!psphotInit()) { + psError(PPSUB_ERR_PROG, false, "Error initialising psphot."); + psFree(config); + goto die; + } + + data = ppSubDataAlloc(config); + + if (!ppSubArguments(argc, argv, data)) { + psError(psErrorCodeLast(), false, "Error reading arguments."); + goto die; + } + + if (!ppSubCamera(data)) { + goto die; + } + + if (!ppSubLoop(data)) { + goto die; + } + + die: + { + psExit exitValue = ppSubExitCode(PS_EXIT_SUCCESS); // Exit code + + if (data && data->stats && data->statsFile) { + psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_DIFF", 0, "Total time (sec)", psTimerMark("ppSub")); + psString stats = psMetadataConfigFormat(data->stats); // Statistics to output + if (!stats || strlen(stats) == 0) { + psError(PPSUB_ERR_IO, false, "Unable to format statistics file"); + } else if (fprintf(data->statsFile, "%s", stats) != strlen(stats)) { + psError(PPSUB_ERR_IO, true, "Unable to write statistics file"); + } + psFree(stats); + if (fclose(data->statsFile) == EOF) { + psError(PPSUB_ERR_IO, true, "Unable to close statistics file"); + } + data->statsFile = NULL; + pmConfigRunFilenameAddWrite(data->config, "STATS", data->statsName); + exitValue = ppSubExitCode(exitValue); + } + + if (config && !ppSubFilesIterateUp(config, PPSUB_FILES_ALL)) { + psError(psErrorCodeLast(), false, "Unable to close files."); + exitValue = ppSubExitCode(exitValue); + pmFPAfileFreeSetStrict(false); + } + + if (data) { + psString dump_file = psMetadataLookupStr(NULL, data->config->arguments, "DUMP_CONFIG"); + if (dump_file) { + if (!pmConfigDump(data->config, dump_file)) { + psError(PPSUB_ERR_IO, false, "Unable to dump configuration."); + exitValue = ppSubExitCode(exitValue); + } + } + psFree(data); + } + + psTrace("ppSub", 1, "Finished at %f sec\n", psTimerMark("ppSub")); + psLogMsg("ppSub", PS_LOG_INFO, "Complete ppSub run: %f sec\n", psTimerMark("ppSub")); + psTimerStop(); + + pmVisualClose(); //close plot windows, if -visual is set + pmModelClassCleanup(); + pmConfigDone(); + pmVisualCleanup (); + psLibFinalize(); + + fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "ppSub"); + exitValue = ppSubExitCode(exitValue); + exit(exitValue); + } +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSub.h =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSub.h (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSub.h (revision 41094) @@ -0,0 +1,193 @@ +/** @file ppSub.h + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef PP_SUB_H +#define PP_SUB_H + +#include +#include +#include + +#include "ppSubErrorCodes.h" + +/// @addtogroup ppSub +/// @{ + +#define PPSUB_RECIPE "PPSUB" /// Name of the recipe to use +#define WCS_TOLERANCE 0.001 // Tolerance for WCS + +// Output files, for activation/deactivation +typedef enum { + PPSUB_FILES_INPUT = 0x01, // Input files + PPSUB_FILES_CONV = 0x02, // Convolved files (output) + PPSUB_FILES_SUB = 0x04, // Subtracted files (output) + PPSUB_FILES_INV = 0x08, // Inverse subtracted files (output) + PPSUB_FILES_PSF = 0x10, // PSF files (output) + PPSUB_FILES_PHOT_SUB = 0x20, // Subtraction photometry files (output) + PPSUB_FILES_PHOT_INV = 0x40, // Inverse subtraction photometry files (output) + PPSUB_FILES_PHOT = 0x80, // General photometry files (internal) + PPSUB_FILES_ALL = 0xFF, // All files +} ppSubFiles; + +/// Data for processing +typedef struct { + pmConfig *config; // Configuration + psErrorCode quality; // Quality code; 0 for no problem + bool photometry; // Perform photometry? + bool inverse; // Output inverse subtraction as well? + bool forcedPhot1; // perform forced photometry? + bool forcedPhot2; // perform forced photometry? + bool saveInConv; // Save convolved input? + bool saveRefConv; // Save convolved reference? + psString stamps; // Stamps file + pmPSF *psf; // Point Spread Function + psString statsName; // Name of statistics file + FILE *statsFile; // Statistics file + psMetadata *stats; // Statistics +} ppSubData; + +/// Constructor +ppSubData *ppSubDataAlloc(pmConfig *config ///< Configuration + ); + +/// Setup the arguments parsing +bool ppSubArguments(int argc, char *argv[], ///< Command-line arguments + ppSubData *data ///< Processing data + ); + +/// Parse the camera input +bool ppSubCamera(ppSubData *data ///< Processing data + ); + +/// Loop over the FPA hierarchy +bool ppSubLoop(ppSubData *data ///< Processing data + ); + +/// Generate (if needed) and set or update the masks for input and reference images +bool ppSubSetMasks(pmConfig *config ///< Configuration + ); + +// Calculate mask statistics +bool ppSubMaskStats(pmConfig *config, pmFPAview *view, psMetadata *stats); + +/// Generate the PSF-matching kernel and convolve the images as needed. Most of this function involves +/// looking up the parameters in the recipe and supplying them to the function pmSubtractionMatch() +bool ppSubMatchPSFs(ppSubData *data ///< Processing data + ); + +/// Threshold low pixels in image +bool ppSubLowThreshold(ppSubData *data ///< Processing data + ); + +/// Generate the output readout and pass the kernel info to the header +bool ppSubDefineOutput(const char *name,///< Name of output to define + pmConfig *config ///< Configuration + ); + +/// Photometry stage 1: measure the PSF from the minuend image +bool ppSubMakePSF(ppSubData *data ///< Processing data + ); + +/// Perform the actual image subtraction, update output concepts +bool ppSubReadoutSubtract(pmConfig *config ///< Configuration + ); + + +/// Photometry stage 2: find and measure sources on the subtracted image +bool ppSubReadoutPhotometry(const char *name, ///< Name of file to photometer + ppSubData *data ///< Processing data + ); + +bool ppSubInputDetections (bool *foundDetections, const char *sourcesName, const char *imageName, ppSubData *data); +bool ppSubReadoutForcedPhot(const char *outputName, const char *targetName, const char *sourceName, ppSubData *data); +bool psphotCopyResults (bool *foundDetections, pmFPAfile *target, pmFPAfile *source, pmFPAview *view); + +/// Higher-order background subtraction +bool ppSubBackground(pmConfig *config ///< Configuration + ); + +/// Perform Variance correction (rescale within a modest range) +bool ppSubVarianceRescale(pmConfig *config, ///< Configuration + ppSubData *data ///< Processing data + ); + +/// Put the program version information into a header +bool ppSubVersionHeader(psMetadata *header ///< Header to populate + ); + +/// Print version information +void ppSubVersionPrint(void); + +/// write the version info to a string +psString ppSubVersionLong(void); + +/// Mark the data quality as bad and prepare to suspend processing +void ppSubDataQuality(ppSubData *data, ///< Processing data + psErrorCode error,///< Error code + ppSubFiles files ///< Files to deactivate + ); + + +/// Activate or deactivate files +void ppSubFilesActivate(pmConfig *config, // Configuration + ppSubFiles files, // File to activate/deactivate + bool state // Activation state + ); + +/// Generate a view suitable for a readout +/// +/// Assumes we're working with skycells +pmFPAview *ppSubViewReadout(void); + +/// Iterate down the FPA hierarchy, opening files +bool ppSubFilesIterateDown(pmConfig *config, // Configuration + ppSubFiles files // Files to open + ); + +/// Iterate up the FPA hierarchy, closing files +bool ppSubFilesIterateUp(pmConfig *config, // Configuration + ppSubFiles files // Files to open + ); + +/// Collect statistics +bool ppSubReadoutStats(ppSubData *data // Processing data + ); + +/// Generate JPEG images +bool ppSubReadoutJpeg(pmConfig *config // Configuration + ); + +/// Generate JPEG images +bool ppSubResidualSampleJpeg(pmConfig *config); + +/// Generate inverse subtraction +bool ppSubReadoutInverse(pmConfig *config // Configuration + ); + + +// Copy every instance of a single keyword from one metadata to another +bool psMetadataCopySingle(psMetadata *target, psMetadata *source, const char *name); + +bool ppSubCopyPSF (pmFPAfile *output, pmFPAfile *input, pmFPAview *view); + +/// Return appropriate exit code +psExit ppSubExitCode(psExit exitValue // Current exit value + ); + +bool ppSubFlagNeighbors(pmConfig *config, pmFPAview *view, psArray *sources, bool matchRef); +bool ppSubMatchSources (psArray *objects, psArray *sources, float RADIUS, float MIN_SN); +bool ppSubSetSourceImageIDs (psArray *sources, int imageID); + +void ppSubSetThreads (void); + +///@} +#endif Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubArguments.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubArguments.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubArguments.c (revision 41094) @@ -0,0 +1,181 @@ +/** @file ppSubArguments.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.59 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +// Print usage information and die +static void usage(const char *program, // Name of the program + psMetadata *arguments, // Command-line arguments + pmConfig *config // Configuration + ) +{ + fprintf(stderr, "\nPan-STARRS PSF-matched image subtraction\n\n"); + fprintf(stderr, "Usage: %s OUTPUT_ROOT \n" + "\t[-psf REFERENCE.psf.fits]\n\n" + "This subtracts the convolved REFERENCE from the INPUT, by default.\n", + program); + fprintf(stderr, "\n"); + psArgumentHelp(arguments); + psFree(config); + pmConfigDone(); + psLibFinalize(); + exit(PS_EXIT_CONFIG_ERROR); +} + +/** + * Add a single filename to the arguments as an array, so that it can be used with pmFPAfileBindFromArgs, etc + */ +static void fileList(const char *file, // The symbolic name for the file + const char *name, // The name of the file + const char *comment, // Description of the file + pmConfig *config // Configuration + ) +{ + psArray *files = psArrayAlloc(1); // Array with file names + files->data[0] = psStringCopy(name); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, file, 0, comment, files); + psFree(files); + return; +} + +void ppSubSetThreads (void) { + // ppSub does not have any of its own thread handlers + return; +} + +bool ppSubArguments(int argc, char *argv[], ppSubData *data) +{ + assert(data); + pmConfig *config = data->config; + assert(config); + + // generic arguments (version, dumpconfig) + PS_ARGUMENTS_GENERIC( ppSub, config, argc, argv ); + + // thread arguments + PS_ARGUMENTS_THREADS( ppSub, config, argc, argv ) + + int argNum = psArgumentGet(argc, argv, "-debug"); // Debugging argument number + if (argNum) { + psArgumentRemove(argNum, &argc, argv); + pmSubtractionRegions(true); + } + + + psMetadata *arguments = config->arguments; // Command-line arguments + psMetadataAddStr(arguments, PS_LIST_TAIL, "-inimage", 0, "Input image", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-inmask", 0, "Input mask image", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-invariance", 0, "Input variance image", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-insources", 0, "Input source list", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-refimage", 0, "Reference image", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-refmask", 0, "Reference mask image", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-refvariance", 0, "Reference variance image", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-refsources", 0, "Reference source list", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-kernel", 0, "Precalculated kernel to apply", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps filename; x,y on each line", NULL); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-convolve", 0, "Image to convolve [1 or 2]", 0); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Perform photometry?", NULL); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-forced-phot", 0, "Perform forced photometry?", NULL); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-forced-input1", 0, "Perform forced photometry?", NULL); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-forced-input2", 0, "Perform forced photometry?", NULL); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Zero point for photometry", NAN); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Generate inverse subtractions?", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-save-inconv", 0, "Save input convolved images?", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-save-refconv", 0, "Save reference convolved images?", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-visual", 0, "Show diagnostic plots", NULL); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-updatemode", 0, "update mode?", false); + + if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) { + usage(argv[0], arguments, config); + } + + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[1]); + + + const char *inImage = psMetadataLookupStr(NULL, arguments, "-inimage"); // Name of input image + if (inImage && strlen(inImage) > 0) { + fileList("INPUT", inImage, "Name of the input image", config); + } + const char *inMask = psMetadataLookupStr(NULL, arguments, "-inmask"); // Name of input mask + if (inMask && strlen(inMask) > 0) { + fileList("INPUT.MASK", inMask, "Name of the input mask image", config); + } + const char *inVariance = psMetadataLookupStr(NULL, arguments, "-invariance"); // Name of input variance + if (inVariance && strlen(inVariance) > 0) { + fileList("INPUT.VARIANCE", inVariance, "Name of the input variance image", config); + } + const char *inSources = psMetadataLookupStr(NULL, arguments, "-insources"); // Name of input source list + if (inSources && strlen(inSources) > 0) { + fileList("INPUT.SOURCES", inSources, "Name of the input source list", config); + } + + const char *refImage = psMetadataLookupStr(NULL, arguments, "-refimage"); // Name of reference image + if (refImage && strlen(refImage) > 0) { + fileList("REF", refImage, "Name of the reference image", config); + } + const char *refMask = psMetadataLookupStr(NULL, arguments, "-refmask"); // Name of reference mask + if (refMask && strlen(refMask) > 0) { + fileList("REF.MASK", refMask, "Name of the reference mask image", config); + } + const char *refVariance = psMetadataLookupStr(NULL, arguments, "-refvariance"); // Name of ref variance + if (refVariance && strlen(refVariance) > 0) { + fileList("REF.VARIANCE", refVariance, "Name of the reference variance image", config); + } + const char *refSources = psMetadataLookupStr(NULL, arguments, "-refsources"); // Name of ref source list + if (refSources && strlen(refSources) > 0) { + fileList("REF.SOURCES", refSources, "Name of the reference source list", config); + } + + const char *kernel = psMetadataLookupStr(NULL, arguments, "-kernel"); // Name of kernel + if (kernel && strlen(kernel) > 0) { + fileList("KERNEL", kernel, "Name of the kernel to apply", config); + } + + data->stamps = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-stamps")); + + data->statsName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-stats")); + if (data->statsName && strlen(data->statsName) > 0) { + psString resolved = pmConfigConvertFilename(data->statsName, config, true, true); // Resolved filename + if (!resolved) { + psError(psErrorCodeLast(), false, "Unable to resolve statistics file: %s", data->statsName); + return false; + } + data->statsFile = fopen(resolved, "w"); + if (!data->statsFile) { + psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved); + psFree(resolved); + return false; + } + psFree(resolved); + } + + data->saveInConv = psMetadataLookupBool(NULL, arguments, "-save-inconv"); + data->saveRefConv = psMetadataLookupBool(NULL, arguments, "-save-refconv"); + + if (psMetadataLookupBool(NULL, arguments, "-visual")) { + pmVisualSetVisual(true); + } + + psTrace("ppSub", 1, "Done reading command-line arguments\n"); + + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubBackground.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubBackground.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubBackground.c (revision 41094) @@ -0,0 +1,82 @@ +/** @file ppSubBackground.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +bool ppSubBackground(pmConfig *config) +{ + psAssert(config, "Require configuration"); + + bool mdok; // Status of metadata lookups + + psMetadata *ppSubRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSub + psAssert(ppSubRecipe, "Need PPSUB recipe"); + psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE); // Recipe for psphot + psAssert(psphotRecipe, "Need PSPHOT recipe for binning"); + + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask + + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output image + + // Generate the background model + if (!psphotModelBackground(config, view, "PPSUB.OUTPUT")) { + psError(psErrorCodeLast(), false, "Unable to model background"); + psFree(view); + return false; + } + + // select the model readout (should now exist) + pmReadout *modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL"); + if (!modelRO) { + psError(psErrorCodeLast(), false, "Unable to find background model"); + psFree(view); + return false; + } + psFree(view); + + psImageBinning *binning = psMetadataLookupPtr(&mdok, modelRO->analysis, "PSPHOT.BACKGROUND.BINNING"); // Binning for model + psImage *modelImage = modelRO->image; // Background model + psImage *image = outRO->image; // Image of interest + psImage *mask = outRO->mask; // Mask of interest + + // Do the background subtraction + int numCols = image->numCols, numRows = image->numRows; // Size of image + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + // special case 1: NAN the masked pixels + if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskBad) { + image->data.F32[y][x] = NAN; + continue; + } + // special case 1: NAN & mask pixels without a valid background model + float value = psImageUnbinPixel(x + 0.5, y + 0.5, modelImage, binning); // Background value + if (!isfinite(value)) { + image->data.F32[y][x] = NAN; + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskBad; + continue; + } + image->data.F32[y][x] -= value; + } + } + + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubCamera.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubCamera.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubCamera.c (revision 41094) @@ -0,0 +1,500 @@ +/** @file ppSubCamera.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +bool checkFileruleFileSave (pmFPAfile *file, pmConfig *config) { + + bool status; + + psMetadata *filerule = pmConfigFileRule(config, file->camera, file->name); // File rule + if (!filerule) return false; + + char *myString = psMetadataLookupStr(&status, filerule, "FILE.SAVE"); + if (!myString) return false; + + // do not change the value from the default unless TRUE or FALSE are found + if (!strcasecmp(myString, "TRUE")) { + file->save = true; + return true; + } + if (!strcasecmp(myString, "FALSE")) { + file->save = false; + return true; + } + + return false; +} + +// Define an input file +static pmFPAfile *defineInputFile(bool *success, + pmConfig *config,// Configuration + pmFPAfile *bind, // File to which to bind, or NULL + char *filerule, // Name of file rule + char *argname, // Argument name + pmFPAfileType fileType // Type of file + ) +{ + bool status; + + *success = false; + + pmFPAfile *file = NULL; + + // look for the file on the argument list + if (bind) { + file = pmFPAfileBindFromArgs(&status, bind, config, filerule, argname); + } else { + file = pmFPAfileDefineFromArgs(&status, config, filerule, argname); + } + + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return false; + } + if (!file) { + // look for the file on the RUN metadata + file = pmFPAfileDefineFromRun(&status, bind, config, filerule); // File to return + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return NULL; + } + } + + if (!file) { + // no file defined + *success = true; + return NULL; + } + + if (file->type != fileType) { + psError(PPSUB_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); + return NULL; + } + + *success = true; + return file; +} + +// Define an output file +static pmFPAfile *defineOutputFile(pmConfig *config, // Configuration + pmFPAfile *template, // File to use as basis for definition + bool source, // Is template a source (T), or for binding (F)? + char *filerule, // Name of file rule + pmFPAfileType fileType // Type of file + ) +{ + + pmFPAfile *file = source ? pmFPAfileDefineFromFile(config, template, 1, 1, filerule) : + pmFPAfileDefineOutput(config, template->fpa, filerule); + if (!file) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from %s"), filerule); + return NULL; + } + if (file->type != fileType) { + psError(PPSUB_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); + return NULL; + } + + return file; +} + + +// Define an output file that will be used in a calculation +// This means it might already be available in the RUN metadata +static pmFPAfile *defineCalcFile(pmConfig *config, // Configuration + pmFPAfile *bind, // File to which to bind, or NULL + char *filerule, // Name of file rule + const char *argname, // Argument name for file + pmFPAfileType fileType // Type of file + ) +{ + bool status; + + // look for the file on the RUN metadata + pmFPAfile *file = pmFPAfileDefineFromRun(&status, NULL, config, filerule); // File to return + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return NULL; + } + file = pmFPAfileBindFromArgs(&status, bind, config, filerule, argname); + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return NULL; + } + if (file) { + // It's an input + file->save = false; + } + + // define new version of file + if (!file) { + file = pmFPAfileDefineOutput(config, bind ? bind->fpa : NULL, filerule); + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return false; + } + if (file) { + // It's an output (respect filerules) + checkFileruleFileSave(file, config); + } + } + + if (!file) { + return NULL; + } + + if (file->type != fileType) { + psError(PPSUB_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); + return NULL; + } + + return file; +} + + +bool ppSubCamera(ppSubData *data) +{ + bool success = true; + + psAssert(data, "Require processing data"); + pmConfig *config = data->config; + psAssert(config, "Require configuration"); + + // Input image + pmFPAfile *input = defineInputFile(&success, config, NULL, "PPSUB.INPUT", "INPUT", PM_FPA_FILE_IMAGE); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.INPUT"); + return false; + } + + defineInputFile(&success, config, input, "PPSUB.INPUT.MASK", "INPUT.MASK", PM_FPA_FILE_MASK); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.INPUT.MASK"); + return false; + } + + pmFPAfile *inVar = defineInputFile(&success, config, input, "PPSUB.INPUT.VARIANCE", "INPUT.VARIANCE", PM_FPA_FILE_VARIANCE); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.INPUT.VARIANCE"); + return false; + } + + defineInputFile(&success, config, NULL, "PPSUB.INPUT.SOURCES", "INPUT.SOURCES", PM_FPA_FILE_CMF); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.INPUT.SOURCES"); + return false; + } + + // Reference image + pmFPAfile *ref = defineInputFile(&success, config, NULL, "PPSUB.REF", "REF", PM_FPA_FILE_IMAGE); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.REF"); + return false; + } + + defineInputFile(&success, config, ref, "PPSUB.REF.MASK", "REF.MASK", PM_FPA_FILE_MASK); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.REF.MASK"); + return false; + } + + pmFPAfile *refVar = defineInputFile(&success, config, ref, "PPSUB.REF.VARIANCE", "REF.VARIANCE", PM_FPA_FILE_VARIANCE); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.REF.VARIANCE"); + return false; + } + + defineInputFile(&success, config, NULL, "PPSUB.REF.SOURCES", "REF.SOURCES", PM_FPA_FILE_CMF); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.REF.SOURCES"); + return false; + } + + // Now that the camera has been determined, we can read the recipe + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + if (!recipe) { + psError(PPSUB_ERR_CONFIG, false, "Unable to find recipe %s", PPSUB_RECIPE); + return false; + } + if (psMetadataLookupBool(NULL, config->arguments, "-photometry")) { + psMetadataAddBool(recipe, PS_LIST_TAIL, "PHOTOMETRY", PS_META_REPLACE, + "Perform photometry?", true); + } + if (psMetadataLookupBool(NULL, config->arguments, "-inverse")) { + psMetadataAddBool(recipe, PS_LIST_TAIL, "INVERSE", PS_META_REPLACE, + "Generate inverse subtractions?", true); + } + if (psMetadataLookupBool(NULL, config->arguments, "-forced-phot")) { + psMetadataAddBool(recipe, PS_LIST_TAIL, "FORCED.PHOTOMETRY.BOTH", PS_META_REPLACE, "Perform forced photometry?", true); + } + if (psMetadataLookupBool(NULL, config->arguments, "-forced-input1")) { + psMetadataAddBool(recipe, PS_LIST_TAIL, "FORCED.PHOTOMETRY.INPUT1", PS_META_REPLACE, "Perform forced photometry?", true); + } + if (psMetadataLookupBool(NULL, config->arguments, "-forced-input2")) { + psMetadataAddBool(recipe, PS_LIST_TAIL, "FORCED.PHOTOMETRY.INPUT2", PS_META_REPLACE, "Perform forced photometry?", true); + } + + data->inverse = psMetadataLookupBool(NULL, recipe, "INVERSE"); + data->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY"); + data->forcedPhot1 = psMetadataLookupBool(NULL, recipe, "FORCED.PHOTOMETRY.BOTH") || psMetadataLookupBool(NULL, recipe, "FORCED.PHOTOMETRY.INPUT1"); + data->forcedPhot2 = psMetadataLookupBool(NULL, recipe, "FORCED.PHOTOMETRY.BOTH") || psMetadataLookupBool(NULL, recipe, "FORCED.PHOTOMETRY.INPUT2"); + + // Convolved input image + pmFPAfile *inConvImage = defineOutputFile(config, input, true, "PPSUB.INPUT.CONV", PM_FPA_FILE_IMAGE); + pmFPAfile *inConvMask = defineOutputFile(config, inConvImage, false, "PPSUB.INPUT.CONV.MASK", + PM_FPA_FILE_MASK); + if (!inConvImage || !inConvMask) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + // these ->save values below are set by command-line arguments (in ppSubArguments.c) : -save-inconv, -save-refconv + checkFileruleFileSave (inConvImage, config); + inConvImage->save = inConvImage->save || data->saveInConv; + + checkFileruleFileSave (inConvMask, config); + inConvMask->save = inConvMask->save || data->saveInConv; + if (inVar) { + pmFPAfile *inConvVar = defineOutputFile(config, inConvImage, false, "PPSUB.INPUT.CONV.VARIANCE", + PM_FPA_FILE_VARIANCE); + if (!inConvVar) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave (inConvVar, config); + inConvVar->save = inConvVar->save || data->saveInConv; + } + + // Convolved ref image + pmFPAfile *refConvImage = defineOutputFile(config, input, true, "PPSUB.REF.CONV", PM_FPA_FILE_IMAGE); + pmFPAfile *refConvMask = defineOutputFile(config, refConvImage, false, "PPSUB.REF.CONV.MASK", + PM_FPA_FILE_MASK); + if (!refConvImage || !refConvMask) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave (refConvImage, config); + refConvImage->save = refConvImage->save || data->saveRefConv; + + checkFileruleFileSave (refConvMask, config); + refConvMask->save = refConvMask->save || data->saveRefConv; + if (refVar) { + pmFPAfile *refConvVar = defineOutputFile(config, refConvImage, false, "PPSUB.REF.CONV.VARIANCE", + PM_FPA_FILE_VARIANCE); + if (!refConvVar) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave (refConvVar, config); + refConvVar->save = refConvVar->save || data->saveRefConv; + } + + // Output image + pmFPAfile *output = defineOutputFile(config, inConvImage, true, "PPSUB.OUTPUT", PM_FPA_FILE_IMAGE); + pmFPAfile *outMask = defineOutputFile(config, output, false, "PPSUB.OUTPUT.MASK", PM_FPA_FILE_MASK); + if (!output || !outMask) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave (output, config); + checkFileruleFileSave (outMask, config); + pmFPAfile *outVar = NULL; + if (inVar && refVar) { + outVar = defineOutputFile(config, output, false, "PPSUB.OUTPUT.VARIANCE", + PM_FPA_FILE_VARIANCE); + if (!outVar) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave (outVar, config); + } + // If we are in update mode unconditionally save the output files + bool updateMode = psMetadataLookupBool(NULL, config->arguments, "-updatemode"); + if (updateMode) { + output->save = true; + outMask->save = true; + if (outVar) { + outVar->save = true; + } + } + + pmFPAfile *inverse = NULL; // Inverse output image + if (data->inverse) { + // Inverse output image + inverse = defineOutputFile(config, output, true, "PPSUB.INVERSE", PM_FPA_FILE_IMAGE); + pmFPAfile *invMask = defineOutputFile(config, inverse, false, "PPSUB.INVERSE.MASK", + PM_FPA_FILE_MASK); + if (!inverse || !invMask) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave (inverse, config); + checkFileruleFileSave (invMask, config); + if (inVar && refVar) { + pmFPAfile *invVar = defineOutputFile(config, inverse, false, "PPSUB.INVERSE.VARIANCE", + PM_FPA_FILE_VARIANCE); + if (!invVar) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave(invVar, config); + } + } + + + // Output JPEGs + pmFPAfile *jpeg1 = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT.JPEG1"); + if (!jpeg1) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSUB.OUTPUT.JPEG1")); + return false; + } + if (jpeg1->type != PM_FPA_FILE_JPEG) { + psError(psErrorCodeLast(), true, "PPSUB.OUTPUT.JPEG1 is not of type JPEG"); + return false; + } + checkFileruleFileSave(jpeg1, config); + pmFPAfile *jpeg2 = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT.JPEG2"); + if (!jpeg2) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSUB.OUTPUT.JPEG2")); + return false; + } + if (jpeg2->type != PM_FPA_FILE_JPEG) { + psError(psErrorCodeLast(), true, "PPSUB.OUTPUT.JPEG2 is not of type JPEG"); + return false; + } + checkFileruleFileSave(jpeg2, config); + + // Output residual JPEG + pmFPAfile *jpeg3 = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT.RESID.JPEG"); + if (!jpeg3) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSUB.OUTPUT.RESID.JPEG")); + return false; + } + if (jpeg3->type != PM_FPA_FILE_JPEG) { + psError(psErrorCodeLast(), true, "PPSUB.OUTPUT.RESID.JPEG is not of type JPEG"); + return false; + } + checkFileruleFileSave(jpeg3, config); + + // Output subtraction kernel + pmFPAfile *kernel = defineCalcFile(config, output, "PPSUB.OUTPUT.KERNELS", "KERNEL", + PM_FPA_FILE_SUBKERNEL); + if (!kernel) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to define file PPSUB.OUTPUT.KERNELS"); + return false; + } + + // psPhot input + if (data->photometry || 1) { + psphotModelClassInit(); // load implementation-specific models + + pmFPAfile *psphot = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT"); + if (!psphot) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.INPUT"); + return false; + } + if (psphot->type != PM_FPA_FILE_IMAGE) { + psError(psErrorCodeLast(), true, "PSPHOT.INPUT is not of type IMAGE"); + return false; + } + pmFPAfileActivate(config->files, false, "PSPHOT.INPUT"); + + // Internal file for getting the PSF from the minuend + pmFPAfile *psf = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.PSF.LOAD"); + if (!psf) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.PSF.LOAD"); + return false; + } + if (psf->type != PM_FPA_FILE_PSF) { + psError(psErrorCodeLast(), true, "PSPHOT.PSF.LOAD is not of type PSF"); + return false; + } + pmFPAfileActivate(config->files, false, "PSPHOT.PSF.LOAD"); + + if (!psphotDefineFiles(config, psphot)) { + psError(psErrorCodeLast(), false, "Unable to set up psphot files."); + return false; + } + + // Deactivate psphot output sources --- we want to define output source files of our own + pmFPAfile *psphotOutput = pmFPAfileSelectSingle(config->files, "PSPHOT.OUTPUT", 0); + psphotOutput->save = false; // this one should NOT be set + + pmFPAfile *outSources = defineOutputFile(config, output, false, "PPSUB.OUTPUT.SOURCES", + PM_FPA_FILE_CMF); + if (!outSources) { + psError(psErrorCodeLast(), false, "Unable to set up output source file."); + return false; + } + checkFileruleFileSave(outSources, config); + + if (data->inverse) { + pmFPAfile *invSources = defineOutputFile(config, inverse, false, "PPSUB.INVERSE.SOURCES", + PM_FPA_FILE_CMF); + if (!invSources) { + psError(psErrorCodeLast(), false, "Unable to set up inverse source file."); + return false; + } + checkFileruleFileSave(invSources, config); + } + + // files need to do the forced photometry on the positions of sources in the positive images + if (data->forcedPhot1) { + // this pmFPAfile is used to carry sources detected in the positive image #1 + pmFPAfile *posSources1 = defineOutputFile(config, input, true, "PPSUB.POS1.SOURCES", PM_FPA_FILE_CMF); + if (!posSources1) { + psError(psErrorCodeLast(), false, "Unable to set up forced source file."); + return false; + } + checkFileruleFileSave(posSources1, config); + + // this pmFPAfile is used to carry sources detected in the diff image @ the positions from positive image #1 + pmFPAfile *frcSources1 = defineOutputFile(config, input, true, "PPSUB.FORCED1.SOURCES", PM_FPA_FILE_CMF); + if (!frcSources1) { + psError(psErrorCodeLast(), false, "Unable to set up forced source file."); + return false; + } + checkFileruleFileSave(frcSources1, config); + } + + if (data->forcedPhot2) { + // this pmFPAfile is used to carry sources detected in the positive image #2 + pmFPAfile *posSources2 = defineOutputFile(config, ref, true, "PPSUB.POS2.SOURCES", PM_FPA_FILE_CMF); + if (!posSources2) { + psError(psErrorCodeLast(), false, "Unable to set up forced source file."); + return false; + } + checkFileruleFileSave(posSources2, config); + + // this pmFPAfile is used to carry sources detected in the diff image @ the positions from positive image #2 + pmFPAfile *frcSources2 = defineOutputFile(config, ref, true, "PPSUB.FORCED2.SOURCES", PM_FPA_FILE_CMF); + if (!frcSources2) { + psError(psErrorCodeLast(), false, "Unable to set up forced source file."); + return false; + } + checkFileruleFileSave(frcSources2, config); + } + } + + return true; +} + + Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubConvolve.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubConvolve.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubConvolve.c (revision 41094) @@ -0,0 +1,338 @@ +#include +#include +#include + +#include "ppSub.h" + +// Print usage information and die +static void usage(const char *program, // Name of the program + psMetadata *arguments, // Command-line arguments + pmConfig *config // Configuration + ) +{ + fprintf(stderr, "\nPan-STARRS image convolution\n\n"); + fprintf(stderr, "Usage: %s OUTPUT_ROOT\n\n", + program); + fprintf(stderr, "\n"); + psArgumentHelp(arguments); + psFree(config); + pmConfigDone(); + psLibFinalize(); + exit(PS_EXIT_CONFIG_ERROR); +} + +// Put filename in it's own list +static void fileList(const char *file, // The symbolic name for the file + const char *name, // The name of the file + const char *comment, // Description of the file + pmConfig *config // Configuration + ) +{ + if (!name) { + return; + } + psArray *files = psArrayAlloc(1); // Array with file names + files->data[0] = psStringCopy(name); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, file, 0, comment, files); + psFree(files); + return; +} + +void ppSubSetThreads (void) { + // ppSub does not have any of its own thread handlers + return; +} + +int main(int argc, char *argv[]) +{ + pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy + + pmConfig *config = pmConfigRead(&argc, argv, PPSUB_RECIPE); // Configuration + if (!config) { + psError(psErrorCodeLast(), false, "Unable to read configuration"); + goto die; + } + + // generic arguments (version, dumpconfig) + PS_ARGUMENTS_GENERIC( ppSub, config, argc, argv ); + + // thread arguments + PS_ARGUMENTS_THREADS( ppSub, config, argc, argv ) + + bool reference = false; // Input is actually the reference image? + int threads = 0; // Number of threads + { + psMetadata *arguments = config->arguments; // Command-line arguments + psMetadataAddStr(arguments, PS_LIST_TAIL, "-image", 0, "Input image", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask", 0, "Input mask", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-variance", 0, "Input variance", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-kernel", 0, "Convolution kernel", NULL); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-reference", 0, "Input is actually reference?", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-save-all", 0, "Save all outputs?", false); + + if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) { + usage(argv[0], arguments, config); + } + + const char *inImage = psMetadataLookupStr(NULL, arguments, "-image"); // Input image + const char *inMask = psMetadataLookupStr(NULL, arguments, "-mask"); // Input mask + const char *inVariance = psMetadataLookupStr(NULL, arguments, "-variance"); // Input variance + const char *inKernel = psMetadataLookupStr(NULL, arguments, "-kernel"); // Input kernel + if (!inImage || !inKernel) { + usage(argv[0], arguments, config); + } + + reference = psMetadataLookupBool(NULL, arguments, "-reference"); + threads = psMetadataLookupS32(NULL, arguments, "NTHREADS"); + fileList("PPSUB.INPUT", inImage, "Input image", config); + fileList("PPSUB.INPUT.MASK", inMask, "Input mask", config); + fileList("PPSUB.INPUT.VARIANCE", inVariance, "Input variance", config); + fileList("PPSUB.INPUT.KERNEL", inKernel, "Input kernel", config); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[1]); + + bool status; // Status of file definition + pmFPAfile *image = pmFPAfileDefineFromArgs(&status, config, "PPSUB.INPUT", "PPSUB.INPUT"); + if (!image || !status || image->type != PM_FPA_FILE_IMAGE) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define input image file."); + goto die; + } + if (inMask) { + pmFPAfile *mask = pmFPAfileBindFromArgs(&status, image, config, "PPSUB.INPUT.MASK", + "PPSUB.INPUT.MASK"); + if (!mask || !status || mask->type != PM_FPA_FILE_MASK) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define input mask file."); + goto die; + } + } + if (inVariance) { + pmFPAfile *variance = pmFPAfileBindFromArgs(&status, image, config, "PPSUB.INPUT.VARIANCE", + "PPSUB.INPUT.VARIANCE"); + if (!variance || !status || variance->type != PM_FPA_FILE_VARIANCE) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define input variance file."); + goto die; + } + } + + pmFPAfile *kernel = pmFPAfileBindFromArgs(&status, image, config, "PPSUB.INPUT.KERNEL", + "PPSUB.INPUT.KERNEL"); + if (!kernel || !status || kernel->type != PM_FPA_FILE_SUBKERNEL) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define kernel file."); + goto die; + } + + pmFPAfile *output = pmFPAfileDefineFromFile(config, image, 1, 1, "PPSUB.INPUT.CONV"); + if (!output || output->type != PM_FPA_FILE_IMAGE) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define output file."); + goto die; + } + output->save = true; + + if (psMetadataLookupBool(NULL, arguments, "-save-all")) { + pmFPAfile *outMask = pmFPAfileDefineOutput(config, output->fpa, "PPSUB.INPUT.CONV.MASK"); + if (!outMask || outMask->type != PM_FPA_FILE_MASK) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define output mask file."); + goto die; + } + outMask->save = true; + + pmFPAfile *outVar = pmFPAfileDefineOutput(config, output->fpa, "PPSUB.INPUT.CONV.VARIANCE"); + if (!outVar || outVar->type != PM_FPA_FILE_VARIANCE) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define output variance file."); + goto die; + } + outVar->save = true; + } + } + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + if (!recipe) { + psError(PPSUB_ERR_CONFIG, false, "Unable to find recipe %s", PPSUB_RECIPE); + goto die; + } + + // Read everything + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "File checks failed."); + goto die; + } + view->chip = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "File checks failed."); + goto die; + } + view->cell = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "File checks failed."); + goto die; + } + view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "File checks failed."); + goto die; + } + + + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout + pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, inRO->analysis, + PM_SUBTRACTION_ANALYSIS_KERNEL); // Convolution kernel + pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Output cell + pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout + pmConceptsCopyFPA(outRO->parent->parent->parent, inRO->parent->parent->parent, true, true); + pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Input header + pmHDU *outHDU = pmHDUFromCell(outCell); // Output header + outHDU->header = psMetadataCopy(outHDU->header, inHDU->header); + ppSubVersionHeader(outHDU->header); + + pmReadout *input, *ref, *inConv, *refConv; // Input and output readouts + if (reference) { + ref = inRO; + refConv = outRO; + input = NULL; + inConv = NULL; + } else { + input = inRO; + inConv = outRO; + ref = NULL; + refConv = NULL; + } + + bool convolve = false; // Do we need to convolve? + float norm = 1.0; // Normalisation to apply + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_1: + if (!reference) { + convolve = true; + } + norm = 1.0 / kernels->solution1->data.F64[PM_SUBTRACTION_INDEX_NORM(kernels)]; + break; + case PM_SUBTRACTION_MODE_2: + if (reference) { + convolve = true; + } + // No normalisation; should already be normalised appropriately + break; + case PM_SUBTRACTION_MODE_DUAL: { + // We only have a single input, so hack together a single-mode kernel + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index of normalisation term + norm = 1.0 / kernels->solution1->data.F64[normIndex]; + convolve = true; + if (reference) { + for (int i = 0; i < kernels->solution2->n; i++) { + kernels->solution1->data.F64[i] = -kernels->solution2->data.F64[i]; + } + kernels->solution1->data.F64[normIndex] = 1.0; + kernels->solution1->data.F64[PM_SUBTRACTION_INDEX_BG(kernels)] = 0.0; + kernels->mode = PM_SUBTRACTION_MODE_2; + } else { + kernels->mode = PM_SUBTRACTION_MODE_1; + } + psFree(kernels->solution2); + kernels->solution2 = NULL; + break; + } + default: + psAbort("Unrecognised kernel mode: %x", kernels->mode); + } + + if (convolve) { + int stride = psMetadataLookupS32(NULL, recipe, "STRIDE"); // Size of convolution patches + float kernelErr = psMetadataLookupF32(NULL, recipe, "KERNEL.ERR"); // Relative sys error in kernel + float covarFrac = psMetadataLookupF32(NULL, recipe, "COVAR.FRAC"); // Fraction for covar calculation + float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction + float poorFrac = psMetadataLookupF32(NULL, recipe, "POOR.FRACTION"); // Fraction for "poor" + + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs + psImageMaskType maskPoor = pmConfigMaskGet("CONV.POOR", config); // Bits to mask for poor pixels + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + + if (threads > 0) { + pmSubtractionThreadsInit(); + } + + // Mask the NAN values (USE BLANK instead of SAT?) + if (!pmReadoutMaskInvalid(input, maskVal, maskBad)) { + psError(PPSUB_ERR_DATA, false, "Unable to mask non-finite pixels in input."); + return false; + } + if (!pmReadoutMaskInvalid(ref, maskVal, maskBad)) { + psError(PPSUB_ERR_DATA, false, "Unable to mask non-finite pixels in reference."); + return false; + } + + if (!pmSubtractionMatchPrecalc(inConv, refConv, input, ref, inRO->analysis, + stride, kernelErr, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac)) { + psError(psErrorCodeLast(), false, "Unable to convolve images."); + if (threads > 0) { + pmSubtractionThreadsFinalize(); + } + goto die; + } + if (threads > 0) { + pmSubtractionThreadsFinalize(); + } + } else { + // Subtract background, since that's what ppSub will do with it + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs + psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + psStatsInit(bg); + if (!psImageBackground(bg, NULL, inRO->image, inRO->mask, maskVal, rng)) { + psError(PPSUB_ERR_DATA, false, "Unable to measure background statistics."); + psFree(bg); + psFree(rng); + goto die; + } + psFree(rng); + outRO->image = (psImage*)psBinaryOp(outRO->image, inRO->image, "-", + psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32)); + psFree(bg); + outRO->data_exists = outRO->parent->data_exists = outRO->parent->parent->data_exists = true; + } + + if (norm != 1.0) { + psBinaryOp(outRO->image, outRO->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); + } + + psFree(outRO); + + die: + { + psExit exitValue = ppSubExitCode(PS_EXIT_SUCCESS); // Exit value + + // Write everything + view->chip = view->cell = view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "File checks failed."); + exitValue = ppSubExitCode(exitValue); + pmFPAfileFreeSetStrict(false); + } + view->readout = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "File checks failed."); + exitValue = ppSubExitCode(exitValue); + pmFPAfileFreeSetStrict(false); + } + view->cell = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "File checks failed."); + exitValue = ppSubExitCode(exitValue); + pmFPAfileFreeSetStrict(false); + } + view->chip = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "File checks failed."); + exitValue = ppSubExitCode(exitValue); + pmFPAfileFreeSetStrict(false); + } + psFree(view); + + psFree(config); + pmConceptsDone(); + pmConfigDone(); + psLibFinalize(); + exit(exitValue); + } + + psAbort("Should never reach here."); +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubData.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubData.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubData.c (revision 41094) @@ -0,0 +1,62 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + + +static void subDataFree(ppSubData *data) +{ + psAssert(!data->statsFile, "Statistics file still open."); + psFree(data->statsName); + psFree(data->stamps); + psFree(data->psf); + psFree(data->stats); + + psFree(data->config); + + return; +} + +ppSubData *ppSubDataAlloc(pmConfig *config) +{ + ppSubData *data = psAlloc(sizeof(ppSubData)); // Processing data, to return + psMemSetDeallocator(data, (psFreeFunc)subDataFree); + + data->config = config; + data->quality = 0; + data->photometry = false; + data->inverse = false; + data->saveInConv = false; + data->saveRefConv = false; + data->stamps = NULL; + data->psf = NULL; + data->statsName = NULL; + data->statsFile = NULL; + data->stats = psMetadataAlloc(); + psMetadataAddS32(data->stats, PS_LIST_TAIL, "QUALITY", 0, "Data quality", 0); + + return data; +} + + +void ppSubDataQuality(ppSubData *data, psErrorCode error, ppSubFiles files) +{ + psAssert(data, "Require processing data"); + + if (psMetadataLookupS32(NULL, data->stats, "QUALITY") == 0) { + data->quality = error; + psMetadataAddS32(data->stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, "Data quality", error); + } + + ppSubFilesActivate(data->config, files, false); + + psErrorClear(); + + return; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubDefineOutput.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubDefineOutput.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubDefineOutput.c (revision 41094) @@ -0,0 +1,283 @@ +/** @file ppSubDefineOutput.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +bool ppSubCopyWarpToChip (psMetadata *tgtHeader, psMetadata *srcHeader); + +bool ppSubDefineOutput(const char *name, pmConfig *config) +{ + psAssert(name, "Require name"); + psAssert(config, "Require configuration"); + + bool status = false; + + pmFPAview *view = ppSubViewReadout(); // View to readout + pmCell *outCell = pmFPAfileThisCell(config->files, view, name); // Output cell + pmFPA *outFPA = outCell->parent->parent; // Output FPA + pmHDU *outHDU = outFPA->hdu; // Output HDU + if (!outHDU->header) { + outHDU->header = psMetadataAlloc(); + } + + // The output readout may already be present if we read in the convolution kernel + pmReadout *outRO = NULL; // Output readout + if (outCell->readouts && outCell->readouts->n > 0 && outCell->readouts->data[0]) { + outRO = psMemIncrRefCounter(outCell->readouts->data[0]); + } else { + outRO = pmReadoutAlloc(outCell); + } + + // Convolved input images + psMetadata *recipe = psMetadataLookupPtr(&status, config->recipes, PPSUB_RECIPE); + bool noConvolve = psMetadataLookupBool(&status, recipe, "NOCONVOLVE"); // Do not use convolved images. + bool reverse = psMetadataLookupBool(&status, config->arguments, "REVERSE"); // Reverse sense of subtraction? + bool mapFromPositive = psMetadataLookupBool(&status, config->arguments, "CHIP.MAP.FROM.POSITIVE"); + // mapFromPositive = true means "always grab the warp->chip map from the positive image (eg, A for A - B, B for B - A) + // mapFromPositive = false means "always grab the warp->chip map from the first image (eg, PPSUB.INPUT if not 'reverse') + + pmReadout *inConv; + if (noConvolve) { + inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout + } else { + inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input readout + } + pmReadout *refConv; + if (noConvolve) { + refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout + } else { + refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference readout + } + + // Add kernel descrption to header. + // We don't know which readout has the kernels because it depends on which PSF is larger + bool mdok; // Status of MD lookup + psMetadata *analysis = inConv->analysis; // Analysis metadata with kernel information + pmHDU *hdu = pmHDUFromCell(inConv->parent); + pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel + + if (!kernels) { + hdu = pmHDUFromCell(refConv->parent); + analysis = refConv->analysis; + kernels = psMetadataLookupPtr(&mdok, analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); + } + if (!kernels && !noConvolve) { + psError(PPSUB_ERR_UNKNOWN, true, "Unable to find SUBTRACTION.KERNEL"); + psFree(outRO); + psFree(view); + return false; + } + psAssert (hdu, "unable to find HDU"); + if (!noConvolve) { + psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, "Subtraction kernel", kernels->description); + } + outHDU->header = psMetadataCopy(outHDU->header, hdu->header); + + // in warp-stack mode, we should always use the warp + // in warp-warp mode, we should use the positive warp for the positive subtraction + + bool normalDiff = true; + if (!strcmp (name, "PPSUB.INVERSE")) { + normalDiff = false; + } + + if (reverse) { + // normal = PPSUB.REF - PPSUB.INPUT + // inverse = PPSUB.INPUT - PPSUB.REF + if (mapFromPositive) { + pmCell *cell_Pos = normalDiff ? pmFPAfileThisCell(config->files, view, "PPSUB.REF") : pmFPAfileThisCell(config->files, view, "PPSUB.INPUT"); + pmHDU *hdu_Pos = pmHDUFromCell(cell_Pos); + ppSubCopyWarpToChip (outHDU->header, hdu_Pos->header); + } else { + pmCell *cell_Pos = pmFPAfileThisCell(config->files, view, "PPSUB.REF"); + pmHDU *hdu_Pos = pmHDUFromCell(cell_Pos); + ppSubCopyWarpToChip (outHDU->header, hdu_Pos->header); + } + } else { + // normal = PPSUB.INPUT - PPSUB.REF + // inverse = PPSUB.REF - PPSUB.INPUT + if (mapFromPositive) { + pmCell *cell_Pos = normalDiff ? pmFPAfileThisCell(config->files, view, "PPSUB.INPUT") : pmFPAfileThisCell(config->files, view, "PPSUB.REF"); + pmHDU *hdu_Pos = pmHDUFromCell(cell_Pos); + ppSubCopyWarpToChip (outHDU->header, hdu_Pos->header); + } else { + pmCell *cell_Pos = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT"); + pmHDU *hdu_Pos = pmHDUFromCell(cell_Pos); + ppSubCopyWarpToChip (outHDU->header, hdu_Pos->header); + } + } + + // Add additional data to the header + pmFPAfile *refFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF"); // Reference file + pmFPAfile *inFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.INPUT"); // Input file + + // save the names of the input and reference image in the header + psString refBase = psStringFileBasename(refFile->origname); // Basename of reference + psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.REFERENCE", 0, "Subtraction reference", refBase); + psFree(refBase); + + psString inBase = psStringFileBasename(inFile->origname); // Basename of input + psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.INPUT", 0, "Subtraction input", inBase); + psFree(inBase); + + ppSubVersionHeader(outHDU->header); + + outRO->analysis = psMetadataCopy(outRO->analysis, analysis); + + psFree(outRO); + psFree(view); + + return true; +} + +// we have 4 sets of header keywords to copy: +// SRC_nnnn, SEC_nnnn, MPX_nnnn, MPY_nnnn + +bool ppSubCopyWarpToChip (psMetadata *tgtHeader, psMetadata *srcHeader) { + + char keyword[80]; + + bool status = false; + + int Nchip = 0; + while (true) { + snprintf (keyword, 80, "SRC_%04d", Nchip); + char *string = psMetadataLookupStr (&status, srcHeader, keyword); + if (!status) { + break; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, keyword, PS_META_REPLACE, "input image", string); + Nchip ++; + } + + for (int i = 0; i < Nchip; i++) { + snprintf (keyword, 80, "SEC_%04d", i); + char *string = psMetadataLookupStr (&status, srcHeader, keyword); + if (!status) { + psWarning ("cannot find keyword %s\n", keyword); + continue; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, keyword, PS_META_REPLACE, "input image", string); + } + + for (int i = 0; i < Nchip; i++) { + snprintf (keyword, 80, "MPX_%04d", i); + char *string = psMetadataLookupStr (&status, srcHeader, keyword); + if (!status) { + psWarning ("cannot find keyword %s\n", keyword); + continue; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, keyword, PS_META_REPLACE, "input image", string); + } + + for (int i = 0; i < Nchip; i++) { + snprintf (keyword, 80, "MPY_%04d", i); + char *string = psMetadataLookupStr (&status, srcHeader, keyword); + if (!status) { + psWarning ("cannot find keyword %s\n", keyword); + continue; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, keyword, PS_META_REPLACE, "input image", string); + } + + // Copy PSREFCAT from the source to the target as well + snprintf(keyword, 80, "PSREFCAT"); + char *string = psMetadataLookupStr(&status, srcHeader, keyword); + if (status) { + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, keyword, PS_META_REPLACE, "input image reference catalog", string); + } + + return true; +} + +bool ppSubCopyWarpToChip_Alt (psMetadata *tgtHeader, psMetadata *srcHeader, bool isPositive) { + + char srcKeyword[80], tgtKeyword[80]; + + bool status = false; + + int Nchip = 0; + while (true) { + snprintf (srcKeyword, 80, "SRC_%04d", Nchip); + if (isPositive) { + snprintf (tgtKeyword, 80, "SRCP_%03d", Nchip); + } else { + snprintf (tgtKeyword, 80, "SRCM_%03d", Nchip); + } + + char *string = psMetadataLookupStr (&status, srcHeader, srcKeyword); + if (!status) { + break; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, srcKeyword, PS_META_REPLACE, "input image", string); + Nchip ++; + } + + for (int i = 0; i < Nchip; i++) { + snprintf (srcKeyword, 80, "SEC_%04d", i); + if (isPositive) { + snprintf (tgtKeyword, 80, "SECP_%03d", Nchip); + } else { + snprintf (tgtKeyword, 80, "SECM_%03d", Nchip); + } + + char *string = psMetadataLookupStr (&status, srcHeader, srcKeyword); + if (!status) { + psWarning ("cannot find keyword %s\n", srcKeyword); + continue; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, srcKeyword, PS_META_REPLACE, "input image", string); + } + + for (int i = 0; i < Nchip; i++) { + snprintf (srcKeyword, 80, "MPX_%04d", i); + if (isPositive) { + snprintf (tgtKeyword, 80, "MPXP_%03d", Nchip); + } else { + snprintf (tgtKeyword, 80, "MPXM_%03d", Nchip); + } + + char *string = psMetadataLookupStr (&status, srcHeader, srcKeyword); + if (!status) { + psWarning ("cannot find keyword %s\n", srcKeyword); + continue; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, srcKeyword, PS_META_REPLACE, "input image", string); + } + + for (int i = 0; i < Nchip; i++) { + snprintf (srcKeyword, 80, "MPY_%04d", i); + if (isPositive) { + snprintf (tgtKeyword, 80, "MPYP_%03d", Nchip); + } else { + snprintf (tgtKeyword, 80, "MPYM_%03d", Nchip); + } + + char *string = psMetadataLookupStr (&status, srcHeader, srcKeyword); + if (!status) { + psWarning ("cannot find keyword %s\n", srcKeyword); + continue; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, srcKeyword, PS_META_REPLACE, "input image", string); + } + return true; +} + Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubErrorCodes.c.in =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubErrorCodes.c.in (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubErrorCodes.c.in (revision 41094) @@ -0,0 +1,37 @@ +/** @file ppSubErrorCodes.c.in + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-05 20:44:04 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#include "pslib.h" +#include "ppSubErrorCodes.h" + +/* + * The line + { PPSUB_ERR_$X{ErrorCode}, "$X{ErrorDescription}"}, + * (without the Xs) + * will be replaced by values from errorCodes.dat + */ +void ppSubErrorRegister(void) +{ + static psErrorDescription errors[] = { + { PPSUB_ERR_BASE, "First value we use; lower values belong to psLib" }, + { PPSUB_ERR_${ErrorCode}, "${ErrorDescription}"}, + }; + static int nerror = PPSUB_ERR_NERROR - PPSUB_ERR_BASE; ///< number of values in enum + + for (int i = 0; i < nerror; i++) { + psErrorDescription *tmp = psAlloc(sizeof(psErrorDescription)); + *tmp = errors[i]; + psErrorRegister(tmp, 1); + psFree(tmp); /* it's on the internal list */ + } + nerror = 0; // don't register more than once +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubErrorCodes.dat =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubErrorCodes.dat (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubErrorCodes.dat (revision 41094) @@ -0,0 +1,13 @@ +# +# This file is used to generate pswarpErrorClasses.h +# +BASE = 14000 First value we use; lower values belong to psLib +UNKNOWN Unknown ppSub error code +NOT_IMPLEMENTED Desired feature is not yet implemented +ARGUMENTS Incorrect arguments +CONFIG Problem in configure files +IO Problem in FITS I/O +DATA Problem in data values +NO_OVERLAP No overlap between input and skycell +PROG Programming error +VARIANCE Variance renormalisation out of bounds Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubErrorCodes.h.in =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubErrorCodes.h.in (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubErrorCodes.h.in (revision 41094) @@ -0,0 +1,30 @@ +/** @file ppSubErrorCodes.h.in + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-05 20:44:04 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#if !defined(PPSUB_ERROR_CODES_H) +#define PPSUB_ERROR_CODES_H +/* + * The line + * PPSUB_ERR_$X{ErrorCode}, + * (without the X) + * + * will be replaced by values from errorCodes.dat + */ +typedef enum { + PPSUB_ERR_BASE = 14000, + PPSUB_ERR_${ErrorCode}, + PPSUB_ERR_NERROR +} ppSubErrorCode; + +void ppSubErrorRegister(void); + +#endif Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubExit.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubExit.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubExit.c (revision 41094) @@ -0,0 +1,72 @@ +#include +#include +#include + +#include "ppSub.h" + +psExit ppSubExitCode(psExit exitValue) +{ + if (exitValue != PS_EXIT_SUCCESS) { + return exitValue; + } + + // gcc -Wswitch complains here if err is declared as type psErrorCode + // the collection of ps*ErrorCode values are enums defined separately for + // each module (psphot, pswarp, etc). the lowest type, psErrorCode is only the base set and does + // not include the possible psphot values + + // for now, to get around this, we just use an int for the switch + + // psErrorCode errorCode = psErrorCodeLast(); // Error code + int errorCode = psErrorCodeLast(); // Error code + if (errorCode != PS_ERR_NONE) { + pmFPAfileFreeSetStrict(false); + psErrorStackPrint(stderr, "Error in subtraction:"); + switch (errorCode) { + case PPSUB_ERR_UNKNOWN: + case PS_ERR_UNKNOWN: + psLogMsg("ppSub", PS_LOG_WARN, "Unknown error code: %x", errorCode); + exitValue = PS_EXIT_UNKNOWN_ERROR; + break; + case PS_ERR_IO: + case PS_ERR_DB_CLIENT: + case PS_ERR_DB_SERVER: + case PS_ERR_BAD_FITS: + case PS_ERR_OS_CALL_FAILED: + case PM_ERR_SYS: + case PPSUB_ERR_IO: + psLogMsg("ppSub", PS_LOG_WARN, "I/O error code: %x", errorCode); + exitValue = PS_EXIT_SYS_ERROR; + break; + case PS_ERR_BAD_PARAMETER_VALUE: + case PS_ERR_BAD_PARAMETER_TYPE: + case PS_ERR_BAD_PARAMETER_NULL: + case PS_ERR_BAD_PARAMETER_SIZE: + case PPSUB_ERR_ARGUMENTS: + case PPSUB_ERR_CONFIG: + psLogMsg("ppSub", PS_LOG_WARN, "Configuration error code: %x", errorCode); + exitValue = PS_EXIT_CONFIG_ERROR; + break; + case PSPHOT_ERR_PSF: + case PPSUB_ERR_DATA: + case PPSUB_ERR_NO_OVERLAP: + psLogMsg("ppSub", PS_LOG_WARN, "Data error code: %x", errorCode); + exitValue = PS_EXIT_DATA_ERROR; + break; + case PS_ERR_UNEXPECTED_NULL: + case PS_ERR_PROGRAMMING: + case PPSUB_ERR_NOT_IMPLEMENTED: + case PPSUB_ERR_PROG: + psLogMsg("ppSub", PS_LOG_WARN, "Programming error code: %x", errorCode); + exitValue = PS_EXIT_PROG_ERROR; + break; + default: + // It's a programming error if we're not dealing with the error correctly + psLogMsg("ppSub", PS_LOG_WARN, "Unrecognised error code: %x", errorCode); + exitValue = PS_EXIT_PROG_ERROR; + break; + } + } + + return exitValue; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubExtras.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubExtras.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubExtras.c (revision 41094) @@ -0,0 +1,38 @@ +/** @file ppSubExtras.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 01:37:17 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#include "ppSub.h" + +// This file contains minor functions used in ppSub. some of these should be pushed into +// psModules or psLib. + +/** + * Copy every instance of a single keyword from one metadata to another + */ +bool psMetadataCopySingle(psMetadata *target, psMetadata *source, const char *name) +{ + PS_ASSERT_METADATA_NON_NULL(target, false); + PS_ASSERT_METADATA_NON_NULL(source, false); + PS_ASSERT_STRING_NON_EMPTY(name, false); + + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", name); + psMetadataIterator *iter = psMetadataIteratorAlloc(source, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + psMetadataItem *item; // Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + psMetadataAddItem(target, item, PS_LIST_TAIL, PS_META_DUPLICATE_OK); + } + psFree(iter); + + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubFiles.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubFiles.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubFiles.c (revision 41094) @@ -0,0 +1,235 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + + +// Input files +static const char *inputFiles[] = { "PPSUB.INPUT", "PPSUB.INPUT.MASK", "PPSUB.INPUT.VARIANCE", + "PPSUB.INPUT.SOURCES", "PPSUB.REF", "PPSUB.REF.MASK", + "PPSUB.REF.VARIANCE", "PPSUB.REF.SOURCES", NULL }; + +// Convolved files +static const char *convFiles[] = { "PPSUB.INPUT.CONV", "PPSUB.INPUT.CONV.MASK", "PPSUB.INPUT.CONV.VARIANCE", + "PPSUB.REF.CONV", "PPSUB.REF.CONV.MASK", "PPSUB.REF.CONV.VARIANCE", + NULL }; + +// Subtraction files +static const char *subFiles[] = { "PPSUB.OUTPUT", "PPSUB.OUTPUT.MASK", "PPSUB.OUTPUT.VARIANCE", + "PPSUB.OUTPUT.JPEG1", "PPSUB.OUTPUT.JPEG2", "PPSUB.OUTPUT.RESID.JPEG", + NULL }; + +// Subtraction photometry +static const char *subPhotFiles[] = { "PPSUB.OUTPUT.SOURCES", NULL }; + +// Inverse subtraction files +static const char *invFiles[] = { "PPSUB.INVERSE", "PPSUB.INVERSE.MASK", "PPSUB.INVERSE.VARIANCE", NULL }; + +// Inverse subtraction photometry +static const char *invPhotFiles[] = { "PPSUB.INVERSE.SOURCES", NULL }; + +// PSF files +static const char *psfFiles[] = { "PSPHOT.PSF.SAVE", NULL }; + +// Calculation (may be either input or output) files +static const char *calcFiles[] = { "PPSUB.OUTPUT.KERNELS", NULL }; + + +// Activate/deactivate a list of files +static void filesActivate(pmConfig *config, // Configuration + const char **files, // List of files + bool state // Activation status to set + ) +{ + for (int i = 0; files[i]; i++) { + pmFPAfileActivate(config->files, state, files[i]); + } + return; +} + +// Activate/deactivate a list of files depending on their 'save' boolean. +// This is so we can activate/deactivate the 'calculation' files, which may be either input or output, which +// is indicated by their 'save' boolean. +static void filesActivateSave(pmConfig *config, // Configuration + const char **files, // List of files + bool save, // Activate when this save state is set + bool state // Activation status to set + ) +{ + for (int i = 0; files[i]; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, files[i], 0); + if (file && file->save == save) { + pmFPAfileActivate(config->files, state, files[i]); + } + } + return; +} + +void ppSubFilesActivate(pmConfig *config, ppSubFiles files, bool state) +{ + psAssert(config, "Require configuration"); + + if (files & PPSUB_FILES_INPUT) { + filesActivate(config, inputFiles, state); + filesActivateSave(config, calcFiles, false, state); + } + if (files & PPSUB_FILES_CONV) { + filesActivate(config, convFiles, state); + } + if (files & PPSUB_FILES_SUB) { + filesActivate(config, subFiles, state); + filesActivateSave(config, calcFiles, true, state); + } + if (files & PPSUB_FILES_INV) { + filesActivate(config, invFiles, state); + } + if (files & PPSUB_FILES_PHOT_SUB) { + filesActivate(config, subPhotFiles, state); + } + if (files & PPSUB_FILES_PHOT_INV) { + filesActivate(config, invPhotFiles, state); + } + if (files & PPSUB_FILES_PSF) { + filesActivate(config, psfFiles, state); + } + if (files & PPSUB_FILES_PHOT) { + psphotFilesActivate(config, state); + } + + return; +} + + +pmFPAview *ppSubViewReadout(void) +{ + pmFPAview *view = pmFPAviewAlloc(0); + view->chip = view->cell = view->readout = 0; + return view; +} + +bool ppSubFilesIterateDown(pmConfig *config, ppSubFiles files) +{ + psAssert(config, "Require configuration"); + + ppSubFilesActivate(config, PPSUB_FILES_ALL, false); + ppSubFilesActivate(config, files, true); + + pmFPAview *view = pmFPAviewAlloc(0);// View to FPA top + + // FPA + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psFree (view); + return false; + } + + // Chip + view->chip = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psFree (view); + return false; + } + + // Cell + view->cell = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psFree (view); + return false; + } + + // Readout + view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psFree (view); + return false; + } + + ppSubFilesActivate(config, PPSUB_FILES_ALL, false); + + psFree (view); + return true; +} + +bool ppSubFilesIterateUp(pmConfig *config, ppSubFiles files) +{ + psAssert(config, "Require configuration"); + + ppSubFilesActivate(config, PPSUB_FILES_ALL, false); + ppSubFilesActivate(config, files, true); + + pmFPAview *view = ppSubViewReadout(); // View to readout + + // Readout + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psFree (view); + return false; + } + + // Cell + view->readout = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psFree (view); + return false; + } + + // Chip + view->cell = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psFree (view); + return false; + } + + // FPA + view->chip = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psFree (view); + return false; + } + + psFree(view); + + // Disable writing of data now that we've written things out + if (files & PPSUB_FILES_CONV) { + pmFPAview *view = ppSubViewReadout(); // View to readout + { + pmReadout *ro = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); + if (ro) { + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + } + { + pmReadout *ro = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); + if (ro) { + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + } + psFree(view); + } + if (files & PPSUB_FILES_SUB) { + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *ro = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); + if (ro) { + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + psFree(view); + } + // Disable writing of data now that we've written things out + if (files & PPSUB_FILES_INV) { + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *ro = pmFPAfileThisReadout(config->files, view, "PPSUB.INVERSE"); + if (ro) { + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + psFree(view); + } + + ppSubFilesActivate(config, PPSUB_FILES_ALL, false); + + psFree (view); + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubFlagNeighbors.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubFlagNeighbors.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubFlagNeighbors.c (revision 41094) @@ -0,0 +1,292 @@ +/** @file ppSubFlagNeighbors.c + * + * @ingroup ppSub + * @author EAM, IfA + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +# define MIN_ERROR 0.02 /* padding for self-match flux test */ +# define MAX_OFFSET 5.0 /* delta mag greater than this is not a self match */ + +// match the diff sources to the detections from the positive or negative input images +// matchRef : false if this is the standard diff image, true if this is the inverse (in-ref vs ref-in) +// this is needed to choose which value is positive and which is negative relative to the difference +bool ppSubFlagNeighbors(pmConfig *config, pmFPAview *view, psArray *sources, bool subInverse) { + + bool status; + + psAssert(config, "Require config"); + psAssert(view, "Require config"); + psAssert (sources, "missing sources?"); + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSUB_RECIPE); // should this be psphot? + psAssert(recipe, "We checked this earlier, so it should be here."); + + float radius = psMetadataLookupF32(&status, recipe, "INPUT.MATCH.RADIUS"); + psAssert (status, "cannot find INPUT.MATCH.RADIUS"); + + float minSN = psMetadataLookupF32(&status, recipe, "INPUT.MATCH.MIN.SN"); + psAssert (status, "cannot find INPUT.MATCH.MIN.SN"); + + // Input sources + pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES"); + pmDetections *inDetections = inSourceRO ? psMetadataLookupPtr(NULL, inSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Input sources + if (!inDetections) { + psWarning("Unable to filter detections based on image A."); + } + + pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES"); + pmDetections *refDetections = refSourceRO ? psMetadataLookupPtr(NULL, refSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Ref sources + if (!refDetections) { + psWarning("Unable to filter detections based on image B."); + } + if (!inDetections && !refDetections) { + return true; + } + + psArray *inSources = inDetections->allSources; + psAssert (inSources, "missing inSources?"); + + psArray *refSources = refDetections->allSources; + psAssert (refSources, "missing refSources?"); + + ppSubSetSourceImageIDs(sources, 0); + ppSubSetSourceImageIDs(inSources, 1); + ppSubSetSourceImageIDs(refSources, 2); + + // define the objects from the diff sources (1-to-1) + psArray *objects = psArrayAllocEmpty(sources->n); + for (int i = 0; i < sources->n; i++) { + pmSource *src = sources->data[i]; + pmPhotObj *obj = pmPhotObjAlloc(); + pmPhotObjAddSource(obj, src); + psArrayAdd (objects, 100, obj); + psFree (obj); + } + + // Match the sources to both inSources and refSources. We want to find the closest source + // in the positive source lists to each source in the diff source list. + // XXX we need to either (a) apply S/N and other filters before we pass the in/ref sources + // to this function, or (b) add those filters in the function + ppSubMatchSources (objects, inSources, radius, minSN); + ppSubMatchSources (objects, refSources, radius, minSN); + + // now mark or flag the sources with matches that meet some criterion + for (int i = 0; i < objects->n; i++) { + pmPhotObj *obj = objects->data[i]; + + // the first entry should be from image 0 (the diff image). + // we should have only 1 of 4 possible results: + // a) only one entry -- the diff source (no positive matches) + // b) 2 matches -- diff and a source from image 1 + // c) 2 matches -- diff and a source from image 2 + // d) 3 matches -- diff and a source from image 1 & 2 + + switch (obj->sources->n) { + case 1: { + pmSource *source = obj->sources->data[0]; + psAssert (source->imageID == 0, "error in pmPhotObj construction: first entry is not image ID 0"); + break; + } + case 2: { + pmSource *source = obj->sources->data[0]; + psAssert (source->imageID == 0, "error in pmPhotObj construction: first entry not 0 (case 2)"); + + pmSource *sourceM1 = obj->sources->data[1]; + psAssert ((sourceM1->imageID == 1) || (sourceM1->imageID == 2), "error in pmPhotObj construction (second entry is not from 1 or 2"); + + // which image gives the match, the positive or negative one? + bool positive = false; + if (subInverse) { + positive = (sourceM1->imageID == 2); + } else { + positive = (sourceM1->imageID == 1); + } + + // check if source and sourceM1 are likely to be the same flux (ie, isolated detection in one image) + if (positive && isfinite(sourceM1->psfMag) && isfinite(source->psfMag)) { + float dMag = sourceM1->psfMag - source->psfMag; + float nSig = fabs(dMag) / hypot(MIN_ERROR, source->psfMagErr); + if (nSig < MAX_OFFSET) { + source->mode2 |= PM_SOURCE_MODE2_DIFF_SELF_MATCH; + } + } + + // only one match. set a flag? + source->mode2 |= PM_SOURCE_MODE2_DIFF_WITH_SINGLE; + if (source->diffStats) { + float SN1 = isfinite(sourceM1->psfMagErr) ? 1.0 / sourceM1->psfMagErr : NAN; + if (positive) { + source->diffStats->SNp = SN1; + source->diffStats->Rp = hypot(source->peak->xf - sourceM1->peak->xf, source->peak->yf - sourceM1->peak->yf); + } else { + source->diffStats->SNm = SN1; + source->diffStats->Rm = hypot(source->peak->xf - sourceM1->peak->xf, source->peak->yf - sourceM1->peak->yf); + } + } + break; + } + case 3: { + pmSource *source = obj->sources->data[0]; + psAssert (source->imageID == 0, "error in pmPhotObj construction: first entry not 0 (case 3)"); + + pmSource *sourceM1 = obj->sources->data[1]; + pmSource *sourceM2 = obj->sources->data[2]; + psAssert ((sourceM1->imageID == 1) || (sourceM1->imageID == 2), "error in pmPhotObj construction (case 3.1)"); + psAssert ((sourceM2->imageID == 1) || (sourceM2->imageID == 2), "error in pmPhotObj construction (case 3.2)"); + psAssert (sourceM1->imageID != sourceM2->imageID, "error in pmPhotObj construction (case 3.3)"); + + // we don't know which matched source (sourceM1 or sourceM2) is the positive detection, + // and the answer depends on both the ID and the subtraction order (in-ref vs ref-in) + // positive is true if sourceM1 is the positive detection + bool positive = false; + if (subInverse) { + positive = (sourceM1->imageID == 2); + } else { + positive = (sourceM1->imageID == 1); + } + + // check for self-detection + float magPSF1 = positive ? sourceM1->psfMag : sourceM2->psfMag; + if (isfinite(magPSF1) && isfinite(source->psfMag)) { + float dMag = magPSF1 - source->psfMag; + float nSig = fabs(dMag) / hypot(MIN_ERROR, source->psfMagErr); + if (nSig < MAX_OFFSET) { + source->mode2 |= PM_SOURCE_MODE2_DIFF_SELF_MATCH; + } + } + + source->mode2 |= PM_SOURCE_MODE2_DIFF_WITH_DOUBLE; + if (source->diffStats) { + float SN1 = isfinite(sourceM1->psfMagErr) ? 1.0 / sourceM1->psfMagErr : NAN; + float SN2 = isfinite(sourceM2->psfMagErr) ? 1.0 / sourceM2->psfMagErr : NAN; + if (positive) { + source->diffStats->SNp = SN1; + source->diffStats->Rp = hypot(source->peak->xf - sourceM1->peak->xf, source->peak->yf - sourceM1->peak->yf); + source->diffStats->SNm = SN2; + source->diffStats->Rm = hypot(source->peak->xf - sourceM2->peak->xf, source->peak->yf - sourceM2->peak->yf); + } else { + source->diffStats->SNp = SN2; + source->diffStats->Rp = hypot(source->peak->xf - sourceM2->peak->xf, source->peak->yf - sourceM2->peak->yf); + source->diffStats->SNm = SN1; + source->diffStats->Rm = hypot(source->peak->xf - sourceM1->peak->xf, source->peak->yf - sourceM1->peak->yf); + } + } + break; + } + default: + psAbort ("error in pmPhotObj construction (unexpected matches)"); + } + } + psFree (objects); + return true; +} + +// XXX based on psphotMatchSourcesToObjects +# define NEXT1 { i++; continue; } +# define NEXT2 { j++; continue; } +bool ppSubMatchSources (psArray *objects, psArray *sources, float RADIUS, float MIN_SN) { + + float dx, dy; + + float RADIUS2 = RADIUS*RADIUS; + + // sort the source list by X + sources = psArraySort (sources, pmSourceSortByX); + objects = psArraySort (objects, pmPhotObjSortByX); + + // match sources to existing objects + + psLogMsg ("psphot", PS_LOG_DETAIL, "attempt to match sources (%ld vs %ld) (SN > %f)", sources->n, objects->n, MIN_SN); + + int nMatch = 0; + int i = 0; + int j = 0; + for (i = j = 0; (i < sources->n) && (j < objects->n); ) { + + pmSource *src = sources->data[i]; + pmPhotObj *obj = objects->data[j]; + + // fprintf (stderr, "%f, %f vs %f, %f\n", src->peak->xf, src->peak->yf, obj->x, obj->y); + + if (!src) NEXT1; + if (!src->peak) NEXT1; + if (!isfinite(src->peak->xf)) NEXT1; + if (!isfinite(src->peak->yf)) NEXT1; + if (!isfinite(src->peak->detValue)) NEXT1; + if (sqrt(src->peak->detValue) < MIN_SN) NEXT1; + + if (!obj) NEXT2; + if (!isfinite(obj->x)) NEXT2; + if (!isfinite(obj->y)) NEXT2; + + dx = src->peak->xf - obj->x; + if (dx < -1.02*RADIUS) NEXT1; + if (dx > +1.02*RADIUS) NEXT2; + + // we are within match range, look for closest match to object j + int Imin = -1; + float Rmin = RADIUS2; + for (int I = i; (dx < +1.02*RADIUS) && (I < sources->n); I++) { + + src = sources->data[I]; + + if (!src) continue; + if (!src->peak) continue; + if (!isfinite(src->peak->xf)) continue; + if (!isfinite(src->peak->yf)) continue; + if (!isfinite(src->peak->detValue)) continue; + if (sqrt(src->peak->detValue) < MIN_SN) continue; + + dx = src->peak->xf - obj->x; + dy = src->peak->yf - obj->y; + + float dr = dx*dx + dy*dy; + if (dr > RADIUS2) continue; + if (dr > Rmin) continue; + Rmin = dr; + Imin = I; + // fprintf (stderr, "j: %d, I: %d, Imin: %d, dr: %f, dx: %f, dy: %f\n", j, I, Imin, dr, dx, dy); + } + + // no match, try next object + if (Imin == -1) { + // fprintf (stderr, "*** missed j: %d, Imin: %d, obj x,y: %f, %f\n", j, Imin, obj->x, obj->y); + j++; + continue; + } + src = sources->data[Imin]; + + // fprintf (stderr, "j: %d, Imin: %d, obj x,y: %f, %f src x,y: %f, %f, SN: %f, ID: %d\n", j, Imin, obj->x, obj->y, src->peak->xf, src->peak->yf, sqrt(src->peak->detValue), src->id); + + // add source to object + pmPhotObjAddSource (obj, src); + nMatch ++; + j++; + } + + psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%d matches)", nMatch); + return true; +} + +bool ppSubSetSourceImageIDs (psArray *sources, int imageID) { + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + source->imageID = imageID; + } + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubInputDetections.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubInputDetections.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubInputDetections.c (revision 41094) @@ -0,0 +1,193 @@ +/** @file ppSubInputDetections.c + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +pmHDU *getHDUandLevel (pmFPALevel *level, pmFPA *fpa, pmFPAview *view) { + + pmHDU *hdu = NULL; + *level = PM_FPA_LEVEL_NONE; + + pmCell *cell = pmFPAviewThisCell(view, fpa); + if (cell) { + hdu = cell->hdu; // cell + if (hdu) { + *level = PM_FPA_LEVEL_CELL; + return hdu; + } + } + + pmChip *chip = pmFPAviewThisChip(view, fpa); + if (chip) { + hdu = chip->hdu; // chip + if (hdu) { + *level = PM_FPA_LEVEL_CHIP; + return hdu; + } + } + + hdu = fpa->hdu; // fpa + if (hdu) { + *level = PM_FPA_LEVEL_FPA; + return hdu; + } + return NULL; +} + +pmHDU *setHDUatLevel (pmFPALevel level, pmFPA *fpa, pmFPAview *view, char *extname) { + + switch (level) { + case (PM_FPA_LEVEL_FPA): { + if (!fpa->hdu) { + fpa->hdu = pmHDUAlloc(extname); + } + return fpa->hdu; + } + case (PM_FPA_LEVEL_CHIP): { + pmChip *chip = pmFPAviewThisChip(view, fpa); + if (!chip->hdu) { + chip->hdu = pmHDUAlloc(extname); + } + return chip->hdu; + } + case (PM_FPA_LEVEL_CELL): { + pmCell *cell = pmFPAviewThisCell(view, fpa); + if (!cell->hdu) { + cell->hdu = pmHDUAlloc(extname); + } + return cell->hdu; + } + default: + return NULL; + } + return NULL; +} + +bool ppSubInputDetections (bool *foundDetections, const char *sourcesName, const char *imageName, ppSubData *data) { + + bool mdok = false; + + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // select the view of interest + pmFPAview *view = ppSubViewReadout(); // View to readout + + // point PSPHOT.INPUT at the positive image + pmFPAfile *imageFile = psMetadataLookupPtr(&mdok, config->files, imageName); // Image to photometer + pmFPAfile *sourcesFile = psMetadataLookupPtr(&mdok, config->files, sourcesName); // Place results her + pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // holder for the working image + + // create a psphot working file, so we don't damage the image or affect existing sources + if (!pmFPACopy(photFile->fpa, imageFile->fpa)) { + psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry"); + psFree(view); + return false; + } + + pmFPALevel hduLevel = PM_FPA_LEVEL_NONE; + pmHDU *imageHDU = getHDUandLevel (&hduLevel, imageFile->fpa, view); + + if (imageHDU) { + pmFPALevel outLevel = PM_FPA_LEVEL_NONE; + pmHDU *sourcesHDU = getHDUandLevel (&outLevel, sourcesFile->fpa, view); + if (!sourcesHDU) { + sourcesHDU = setHDUatLevel (hduLevel, sourcesFile->fpa, view, imageHDU->extname); + } + + if (!sourcesHDU->header) { + sourcesHDU->header = psMetadataAlloc(); + } + if (!psMetadataCopy(sourcesHDU->header, imageHDU->header)) { + psError(PPSUB_ERR_PROG, false, "Unable to copy header"); + psFree(view); + return false; + } + } + + if (!psphotReadout(config, view, "PSPHOT.INPUT")) { + psErrorStackPrint(stderr, "Unable to perform photometry on image"); + psWarning("Unable to perform photometry on image --- suspect bad data quality."); + ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); + } + // save the outputs on PPSUB.POS1.SOURCES + if (!psphotCopyResults (foundDetections, sourcesFile, photFile, view)) { + psError(PPSUB_ERR_PROG, false, "Unable to copy psphot outputs"); + psFree(view); + return false; + } + // if no sources were found here, we report that back and let them handle it + + psFree(view); + return true; +} + +bool psphotCopyResults (bool *foundDetections, pmFPAfile *target, pmFPAfile *source, pmFPAview *view) { + + pmReadout *sourceRO = pmFPAviewThisReadout(view, source->fpa); + psAssert(sourceRO, "programming error: source readout not defined"); + + pmReadout *targetRO = pmFPAviewThisReadout(view, target->fpa); + if (!targetRO) { + pmCell *targetCell = pmFPAviewThisCell(view, target->fpa); + targetRO = pmReadoutAlloc(targetCell); + psAssert(targetRO, "programming error: could not make target readout"); + } + + // If no sources were found, there's no error, but we need to inform the calling function + pmDetections *detections = psMetadataLookupPtr(NULL, sourceRO->analysis, "PSPHOT.DETECTIONS"); // Sources + if (!detections) { + *foundDetections = false; + return true; + } else { + *foundDetections = true; + } + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!psMetadataCopySingle(targetRO->analysis, sourceRO->analysis, "PSPHOT.DETECTIONS")) { + psError(PPSUB_ERR_PROG, false, "Unable to copy PSPHOT.DETECTIONS"); + return false; + } + if (!psMetadataCopySingle(targetRO->analysis, sourceRO->analysis, "PSPHOT.HEADER")) { + psError(PPSUB_ERR_PROG, false, "Unable to copy PSPHOT.HEADER"); + return false; + } + if (!psMetadataCopySingle(targetRO->analysis, sourceRO->analysis, PM_DETEFF_ANALYSIS)) { + psError(PPSUB_ERR_PROG, false, "Unable to copy Detection Efficiency"); + return false; + } + + // Ensure photometry information is put in the header + // XXX create one if it does not exist? + pmHDU *hdu = pmHDUFromReadout(targetRO); // HDU for readout + if (hdu) { + psMetadata *header = psMetadataLookupMetadata(NULL, targetRO->analysis, "PSPHOT.HEADER"); // Header + hdu->header = psMetadataCopy(hdu->header, header); + } + + targetRO->data_exists = true; + targetRO->parent->data_exists = true; + targetRO->parent->parent->data_exists = true; + + return true; +} + Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubKernel.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubKernel.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubKernel.c (revision 41094) @@ -0,0 +1,130 @@ +/** @file ppSubKernel.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 01:37:17 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#define KERNEL_MOSAIC 2 ///< Half-number of kernel instances in the mosaic image + + +int main(int argc, char *argv[]) +{ + if (argc != 3) { + fprintf(stderr, "Usage: %s INPUT_KERNEL OUTPUT_IMAGE\n\n", argv[0]); + exit(PS_EXIT_SYS_ERROR); + } + + (void) psTraceSetLevel("psModules.imcombine", 7); + + const char *inName = argv[1]; // Input file name + const char *outName = argv[2]; // Output file name + + psFits *inFile = psFitsOpen(inName, "r"); // Input file + if (!inFile) { + psErrorStackPrint(stderr, "Unable to open input kernel %s", inName); + exit(PS_EXIT_DATA_ERROR); + } + + pmReadout *ro = pmReadoutAlloc(NULL); // Readout to hold kernel + if (!pmReadoutReadSubtractionKernels(ro, inFile)) { + psErrorStackPrint(stderr, "Unable to read input kernel"); + psFree(ro); + psFitsClose(inFile); + exit(PS_EXIT_DATA_ERROR); + } + psFitsClose(inFile); + + pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, ro->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); + if (!psMemIncrRefCounter(kernels)) { + psErrorStackPrint(stderr, "Unable to find input kernel"); + psFree(ro); + exit(PS_EXIT_PROG_ERROR); + } + psFree(ro); + + pmSubtractionMode mode = kernels->mode; // Mode for subtraction + int fullSize = 2 * kernels->size + 1 + 1; // Full size of kernel + int imageSize = (2 * KERNEL_MOSAIC + 1) * fullSize; // Size of image + psImage *image = psImageAlloc((mode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) * imageSize - 1 + + (mode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0), imageSize - 1, PS_TYPE_F32); + psImageInit(image, NAN); + for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) { + for (int i = -KERNEL_MOSAIC; i <= KERNEL_MOSAIC; i++) { + psImage *kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC, + (float)j / (float)KERNEL_MOSAIC, + false); // Image of the kernel + if (!kernel) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); + psFree(image); + psFree(kernel); + psFree(kernels); + exit(PS_EXIT_SYS_ERROR); + } + + if (psImageOverlaySection(image, kernel, (i + KERNEL_MOSAIC) * fullSize, + (j + KERNEL_MOSAIC) * fullSize, "=") == 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image."); + psFree(kernel); + psFree(image); + psFree(kernels); + exit(PS_EXIT_PROG_ERROR); + } + psFree(kernel); + + if (mode == PM_SUBTRACTION_MODE_DUAL) { + kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC, + (float)j / (float)KERNEL_MOSAIC, + true); // Image of the kernel + if (!kernel) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); + psFree(image); + psFree(kernels); + exit(PS_EXIT_SYS_ERROR); + } + + if (psImageOverlaySection(image, kernel, + (2 * KERNEL_MOSAIC + 1 + i + KERNEL_MOSAIC) * fullSize + 4, + (j + KERNEL_MOSAIC) * fullSize, "=") == 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image."); + psFree(kernel); + psFree(image); + psFree(kernels); + exit(PS_EXIT_PROG_ERROR); + } + psFree(kernel); + } + } + } + psFree(kernels); + + psFits *outFile = psFitsOpen(outName, "w"); // Output file + if (!outFile) { + psErrorStackPrint(stderr, "Unable to open output file %s", outName); + psFree(image); + exit(PS_EXIT_SYS_ERROR); + } + if (!psFitsWriteImage(outFile, NULL, image, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write output file"); + psFree(image); + psFitsClose(outFile); + exit(PS_EXIT_SYS_ERROR); + } + psFitsClose(outFile); + psFree(image); + + exit(PS_EXIT_SUCCESS); +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubLoop.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubLoop.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubLoop.c (revision 41094) @@ -0,0 +1,356 @@ +/** @file ppSubLoop.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppSub.h" + +bool dumpout(pmConfig *config, char *name) +{ + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *out = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); + psphotSaveImage (NULL, out->image, name); + psFree(view); + return true; +} + + +bool ppSubLoop(ppSubData *data) +{ + bool mdok = false; + bool success = true; + + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration."); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images. + + pmConfigCamerasCull(config, NULL); + pmConfigRecipesCull(config, "PPSUB,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG"); + + pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSUB.INPUT"); + pmFPAfile *reference = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF"); + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSUB.OUTPUT"); + psAssert(input && reference && output, "Require files"); + + if (!ppSubFilesIterateDown(config, PPSUB_FILES_INPUT | PPSUB_FILES_CONV)) { + psError(PPSUB_ERR_IO, false, "Unable to load files."); + return false; + } + + psTimerStart("PPSUB_MATCH"); + + if (!ppSubSetMasks(config)) { + psError(psErrorCodeLast(), false, "Unable to set masks."); + return false; + } + + if (data->forcedPhot1) { + bool foundDetections = false; + if (!ppSubInputDetections(&foundDetections, "PPSUB.POS1.SOURCES", "PPSUB.INPUT", data)) { + psError(psErrorCodeLast(), false, "Unable to measure positive detections (1)"); + success = false; + goto ESCAPE; + } + // if nothing was found, don't bother doing the forced photometry below + if (!foundDetections) { + psWarning ("no sources found in positive image 1, skipping forced photometry"); + data->forcedPhot1 = false; + } + } + if (data->forcedPhot2) { + // Change the recipe to use a higher nsigma limit and quit after pass1 + psMetadata *psphotRecipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + + psF32 nsigma_peak_save = psMetadataLookupF32 (NULL, psphotRecipe, "PEAKS_NSIGMA_LIMIT"); + char *breakPt_save = psMetadataLookupStr (NULL, psphotRecipe, "BREAK_POINT"); + + psF32 pos2_nsigma_peak = psMetadataLookupF32 (&mdok, psphotRecipe, "PEAKS_POS2_NSIGMA_LIMIT"); + if (!mdok) { + psWarning("PEAKS_POS2_NSIGMA_LIMIT not found in psphotRecipe. Will use 25.\n"); + pos2_nsigma_peak = 25.; + } + psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", pos2_nsigma_peak); + psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", "PASS1"); + + bool foundDetections = false; + if (!ppSubInputDetections(&foundDetections, "PPSUB.POS2.SOURCES", "PPSUB.REF", data)) { + psError(psErrorCodeLast(), false, "Unable to measure positive detections (2)"); + psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save); + psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save); + success = false; + goto ESCAPE; + } + psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save); + psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save); + // if nothing was found, don't bother doing the forced photometry below + if (!foundDetections) { + psWarning ("no sources found in positive image 2, skipping forced photometry"); + data->forcedPhot2 = false; + } + } + + // XXX if it exists, use the POS1, POS2 successs for the FWHMs + if (!ppSubMatchPSFs(data)) { + psError(psErrorCodeLast(), false, "Unable to match PSFs."); + success = false; + goto ESCAPE; + } + + if (data->quality) { + // Can't do anything at all + success = false; + goto ESCAPE; + } + // generate the residual stamp grid for visualization + if (!ppSubResidualSampleJpeg(config)) { + psError(psErrorCodeLast(), false, "Unable to update."); + success = false; + goto ESCAPE; + } + + // XXX add in a positive image detection step here (if needed) + + + psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs", + psTimerClear("PPSUB_MATCH")); + + // Close input files (freeing up space) : for the noConvolve case, we cannot close the + // inputs since we will use them for subtraction below. wait until later to do the work + if (!noConvolve) { + if (!ppSubFilesIterateUp(config, PPSUB_FILES_INPUT)) { + psError(PPSUB_ERR_IO, false, "Unable to close input files."); + success = false; + goto ESCAPE; + } + } + if (!ppSubLowThreshold(data)) { + psError(psErrorCodeLast(), false, "Unable to threshold images."); + success = false; + goto ESCAPE; + } + + // Set up subtraction files + if (!ppSubFilesIterateDown(config, PPSUB_FILES_SUB)) { + psError(PPSUB_ERR_IO, false, "Unable to set up subtraction files."); + success = false; + goto ESCAPE; + } + + if (!ppSubDefineOutput("PPSUB.OUTPUT", config)) { + psError(psErrorCodeLast(), false, "Unable to define output."); + success = false; + goto ESCAPE; + } + + if (!data->quality && !ppSubMakePSF(data)) { + psError(psErrorCodeLast(), false, "Unable to generate PSF."); + success = false; + goto ESCAPE; + } + + // Now we've got a PSF, blow away detections in case they're confused with real output detections + { + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *out = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); + if (psMetadataLookup(out->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(out->analysis, "PSPHOT.DETECTIONS"); + } + psFree(view); + } + + if (!ppSubReadoutSubtract(config)) { + psError(psErrorCodeLast(), false, "Unable to subtract images."); + success = false; + goto ESCAPE; + } + // dumpout(config, "diff.1.fits"); + + // Close convolved files + if (!ppSubFilesIterateUp(config, PPSUB_FILES_PSF | PPSUB_FILES_CONV)) { + psError(PPSUB_ERR_IO, false, "Unable to close input files."); + success = false; + goto ESCAPE; + } + // dumpout(config, "diff.2a.fits"); + if (noConvolve) { + if (!ppSubFilesIterateUp(config, PPSUB_FILES_INPUT)) { + psError(PPSUB_ERR_IO, false, "Unable to close input files."); + success = false; + goto ESCAPE; + } + } + + // Higher order background subtraction using psphot + if (!ppSubBackground(config)) { + psError(psErrorCodeLast(), false, "Unable to subtract background."); + success = false; + goto ESCAPE; + } + // dumpout(config, "diff.2b.fits"); + + // Perform Variance correction (rescale within a modest range) + if (!ppSubVarianceRescale(config, data)) { + psError(psErrorCodeLast(), false, "Unable to rescale variance."); + success = false; + goto ESCAPE; + } + // dumpout(config, "diff.2c.fits"); + + if (data->quality) { + // Done all we can do up to this point + success = false; + goto ESCAPE; + } + + if (!ppSubFilesIterateDown(config, PPSUB_FILES_PHOT_SUB)) { + psError(PPSUB_ERR_IO, false, "Unable to set up photometry files."); + success = false; + goto ESCAPE; + } + + if (!data->quality && !ppSubReadoutPhotometry("PPSUB.OUTPUT", data)) { + psError(psErrorCodeLast(), false, "Unable to perform photometry."); + success = false; + goto ESCAPE; + } + // dumpout(config, "diff.3.fits"); + + // forced photometry for positive image 1 + if (data->forcedPhot1 && !data->quality) { + if (!ppSubReadoutForcedPhot("PPSUB.FORCED1.SOURCES", "PPSUB.OUTPUT", "PPSUB.POS1.SOURCES", data)) { + psError(psErrorCodeLast(), false, "Unable to perform photometry."); + success = false; + goto ESCAPE; + } + } + + // forced photometry for positive image 2 + if (data->forcedPhot2 && !data->quality) { + if (!ppSubReadoutForcedPhot("PPSUB.FORCED2.SOURCES", "PPSUB.OUTPUT", "PPSUB.POS2.SOURCES", data)) { + psError(psErrorCodeLast(), false, "Unable to perform photometry."); + success = false; + goto ESCAPE; + } + } + + if (!ppSubFilesIterateUp(config, PPSUB_FILES_PHOT_SUB)) { + psError(PPSUB_ERR_IO, false, "Unable to set up photometry files."); + success = false; + goto ESCAPE; + } + + // Perform statistics on the cell + if (!ppSubReadoutStats(data)) { + psError(psErrorCodeLast(), false, "Unable to collect statistics"); + success = false; + goto ESCAPE; + } + // Do Mask Stats + { + pmFPAview *view = ppSubViewReadout(); // View to readout + if (!ppSubMaskStats(config, view,data->stats)) { + psError(psErrorCodeLast(), false, "Unable to generate mask statistics"); + success = false; + psFree(view); + goto ESCAPE; + } + psFree(view); + } + // dumpout(config, "diff.4.fits"); + + // generate the binned image used to write the jpeg + if (!ppSubReadoutJpeg(config)) { + psError(psErrorCodeLast(), false, "Unable to update."); + success = false; + goto ESCAPE; + } + + if (data->inverse) { + // Set up inverse subtraction files + if (!ppSubFilesIterateDown(config, PPSUB_FILES_INV)) { + psError(PPSUB_ERR_IO, false, "Unable to set up inverse files."); + success = false; + goto ESCAPE; + } + + if (data->inverse && !ppSubDefineOutput("PPSUB.INVERSE", config)) { + psError(psErrorCodeLast(), false, "Unable to define inverse."); + success = false; + goto ESCAPE; + } + + if (!ppSubReadoutInverse(config)) { + psError(psErrorCodeLast(), false, "Unable to invert images."); + success = false; + goto ESCAPE; + } + + //MEH -- need to also clear out detections or inv.cmf corrupted on update + pmFPAview *viewinv = ppSubViewReadout(); // View to readout + pmReadout *inv = pmFPAfileThisReadout(config->files, viewinv, "PPSUB.INVERSE"); + if (psMetadataLookup(inv->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(inv->analysis, "PSPHOT.DETECTIONS"); + } + psFree(viewinv); + + // Close subtraction files and open inverse photometry files + if (!ppSubFilesIterateUp(config, PPSUB_FILES_SUB)) { + psError(PPSUB_ERR_IO, false, "Unable to close subtraction files."); + success = false; + goto ESCAPE; + } + + if (!ppSubFilesIterateDown(config, PPSUB_FILES_PHOT_INV)) { + psError(PPSUB_ERR_IO, false, "Unable to set up inverse files."); + success = false; + goto ESCAPE; + } + + if (!data->quality && !ppSubReadoutPhotometry("PPSUB.INVERSE", data)) { + psError(psErrorCodeLast(), false, "Unable to perform photometry."); + success = false; + goto ESCAPE; + } + + // Close inverse subtraction files + if (!ppSubFilesIterateUp(config, PPSUB_FILES_INV | PPSUB_FILES_PHOT_INV)) { + psError(PPSUB_ERR_IO, false, "Unable to close subtraction files."); + success = false; + goto ESCAPE; + } + } else { + // Close subtraction files + // dumpout(config, "diff.5.fits"); + if (!ppSubFilesIterateUp(config, PPSUB_FILES_SUB)) { + psError(PPSUB_ERR_IO, false, "Unable to close subtraction files."); + success = false; + goto ESCAPE; + } + } + + ESCAPE: + pmFPAfileDropInternal(config->files, "PSPHOT.BACKGND"); + pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL"); + pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV"); + + return success; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubMakePSF.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubMakePSF.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubMakePSF.c (revision 41094) @@ -0,0 +1,203 @@ +/** @file ppSubMakePSF.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +psArray *ppSubSelectPSFSources(psArray *sources); + +bool ppSubMakePSF(ppSubData *data) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + if (!data->photometry) { + return true; + } + + psTimerStart("PPSUB_PHOT"); + + bool reverse = psMetadataLookupBool(NULL, config->arguments, "REVERSE"); // Reverse sense of subtraction? + + bool mdok = false; // Status of MD lookup + pmReadout *minuend = NULL; // Image that will be positive following subtraction + pmFPAfile *minuendFile = NULL; // File for minuend image + pmFPAview *view = ppSubViewReadout(); // View to readout + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSub + bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images. + + if (noConvolve) { + // if we do not convolve, we need to copy the detections to the image for analysis + pmReadout *psfSourcesRO = NULL; // readout containing loaded sources for psf model + psWarning("Not using Convolved images because NOCONVOLVE is TRUE\n"); + if (reverse) { + minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); + minuendFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.REF"); + psfSourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES"); + } else { + minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); + minuendFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.INPUT"); + psfSourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES"); + } + psAssert (psfSourcesRO, "missing readout with sources"); + pmDetections *psfDetections = psMetadataLookupPtr(&mdok, psfSourcesRO->analysis, "PSPHOT.DETECTIONS"); + psMetadataAddPtr(minuend->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "Merged source list", psfDetections); + } else { + printf("Using Convolved images because NOCONVOLVE is FALSE\n"); + if (reverse) { + minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); + minuendFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.REF.CONV"); + } else { + minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); + minuendFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.INPUT.CONV"); + } + } + + pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file + if (!pmFPACopy(photFile->fpa, minuendFile->fpa)) { + psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry"); + psFree(view); + return false; + } + + pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer + + // we need to remove any existing PSPHOT.DETECTIONS (why not do this in psphot?) + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) { + psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF"); + } + +# ifdef TESTING + // XXX for testing, dump these images: + psphotSaveImage (NULL, photRO->image, "findpsf.im.fits"); + psphotSaveImage (NULL, photRO->variance, "findpsf.wt.fits"); + psphotSaveImage (NULL, photRO->mask, "findpsf.mk.fits"); +# endif + + // Extract the loaded sources from the associated readout, and generate PSF + // Here, we assume the image is background-subtracted + pmDetections *detections = psMetadataLookupPtr(&mdok, minuend->analysis, "PSPHOT.DETECTIONS"); + if (!detections || !detections->allSources) { + psError(PPSUB_ERR_CONFIG, true, "No sources from which to determine PSF."); + psFree(view); + return false; + } + psArray *sources = detections->allSources; + + // XXX filter sources? limit the total number and return only brighter objects? + // use flags to toss totally bogus entries? + psArray *goodSources = ppSubSelectPSFSources (sources); + + if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", goodSources)) { + // This is likely a data quality issue + // XXX Split into multiple cases using error codes? + psErrorStackPrint(stderr, "Unable to determine PSF"); + psWarning("Unable to determine PSF --- suspect bad data quality."); + ppSubDataQuality(data, PSPHOT_ERR_PSF, PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); + psFree(view); + psFree(goodSources); + return true; + } + + // save the resulting PSF information on the pmFPAfile PSPHOT.PSF.LOAD + pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file + if (!ppSubCopyPSF(psfFile, photFile, view)) { + psErrorStackPrint(stderr, "PSF was not generated"); + psWarning("PSF was not generated --- suspect bad data quality."); + ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); + } + + // Get rid of the generated header; it will be regenerated by the real photometry run + psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER"); + + psFree(goodSources); + psFree(view); + + return true; +} + +bool ppSubCopyPSF(pmFPAfile *output, pmFPAfile *input, pmFPAview *view) +{ + pmChip *inputChip = pmFPAviewThisChip(view, input->fpa); // Chip with PSF info + pmChip *outputChip = pmFPAviewThisChip(view, output->fpa); // Chip to store PSF info + + pmReadout *inputRO = pmFPAviewThisReadout(view, input->fpa); // Readout with PSF info + pmReadout *outputRO = pmFPAviewThisReadout(view, output->fpa); // Readout to store PSF info + + if (!outputRO) { + pmCell *outputCell = pmFPAviewThisCell(view, output->fpa); + outputRO = pmReadoutAlloc(outputCell); + outputRO->image = psMemIncrRefCounter(inputRO->image); + psFree(outputRO); // I have a copy on the outputCell + } + + // Copy the PSF-related data + psMetadataIterator *iter = psMetadataIteratorAlloc(inputRO->analysis, PS_LIST_HEAD, "^PSF.*"); + psMetadataItem *item; // Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + psMetadataAddItem(outputRO->analysis, item, PS_LIST_TAIL, PS_META_REPLACE); + } + psFree(iter); + + // copy the PSF model data + pmPSF *psf = psMetadataLookupPtr(NULL, inputChip->analysis, "PSPHOT.PSF"); // PSF for photometry + if (!psf) { + psErrorStackPrint(stderr, "No PSF available"); + return false; + } + + psMetadataAddPtr(outputChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE, + "PSF from ppSubMakePSF", psf); + + return true; +} + + +// XXX hardwired MAX for now +# define MAX_NPSF 500 + +psArray *ppSubSelectPSFSources(psArray *sources){ + + sources = psArraySort (sources, pmSourceSortByFlux); + + psArray *subset = psArrayAllocEmpty(MAX_NPSF); + + int nPSF = 0; + for (int i = 0; (nPSF < MAX_NPSF) && (i < sources->n); i++) { + + pmSource *source = sources->data[i]; + if (!source) continue; + + // skip non-astronomical objects (very likely defects) + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + + psArrayAdd (subset, 100, source); + nPSF++; + } + + return subset; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubMaskStats.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubMaskStats.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubMaskStats.c (revision 41094) @@ -0,0 +1,55 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppSub.h" + + +bool ppSubMaskStats(pmConfig *config, pmFPAview *view, psMetadata *stats) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(config, false); + + pmReadout *readout = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); + if (!stats || !readout || !readout->data_exists) { + // Nothing to process + return(true); + } + + bool status; + + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSUB_RECIPE); + psImageMaskType staticMaskVal = psMetadataLookupImageMask(&status, recipe, "MASKSTAT.STATIC"); + psImageMaskType magicMaskVal = psMetadataLookupImageMask(&status, recipe, "MASKSTAT.MAGIC"); + psImageMaskType dynamicMaskVal = psMetadataLookupImageMask(&status, recipe, "MASKSTAT.DYNAMIC"); + psImageMaskType advisoryMaskVal = psMetadataLookupImageMask(&status, recipe, "MASKSTAT.ADVISORY"); + + psS32 Npix_valid = 0; + psS32 Npix_static = 0; + psS32 Npix_magic = 0; + psS32 Npix_dynamic = 0; + psS32 Npix_advisory = 0; + + psImage *mask = readout->mask; // Mask of interest; + if (!pmSingleImageMaskStats(mask,&Npix_valid,&Npix_static,&Npix_magic, + &Npix_dynamic,&Npix_advisory, + staticMaskVal,magicMaskVal, + dynamicMaskVal,advisoryMaskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate masks for readout."); + return(false); + } + psMetadataAddS32(stats, PS_LIST_TAIL,"MASKFRAC_NPIX", 0, + "Number of valid pixels", Npix_valid); + psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_STATIC", 0, + "Fraction of pixels statically masked", (float) Npix_static / Npix_valid); + psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_DYNAMIC", 0, + "Fraction of pixels dynamically masked", (float) Npix_dynamic / Npix_valid); + psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_MAGIC", 0, + "Fraction of pixels magically masked", (float) Npix_magic / Npix_valid); + psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_ADVISORY", 0, + "Fraction of pixels masked as an advisory", (float) Npix_advisory / Npix_valid); + return(true); +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubMatchPSFs.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubMatchPSFs.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubMatchPSFs.c (revision 41094) @@ -0,0 +1,538 @@ +/** @file ppSubMatchPSFs.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +#define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix + +// Normalise a region on an image +static void normaliseRegion(psImage *image, // Image to normalise + const psRegion *region, // Region of image to normalise + float norm // Normalisation + ) +{ + if (!image) { + return; + } + psAssert(region, "Expect region"); + psImage *subImage = psImageSubset(image, *region); // Sub-image + psBinaryOp(subImage, subImage, "*", psScalarAlloc(norm, PS_TYPE_F32)); + psFree(subImage); + return; +} + +// Measure the PSF for an image +static float subImagePSF(ppSubData *data, // Processing data + const pmReadout *ro, // Readout for which to measure PSF + psArray *sources // Sources with positions at which to measure PSF + ) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + psAssert(ro, "Need readout"); + psAssert(sources, "Need sources."); + + pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // Photometry file + psAssert(photFile, "Need photometry file."); + if (!pmFPACopy(photFile->fpa, ro->parent->parent->parent)) { + psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry"); + return NAN; + } + + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer + + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) { + psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF"); + } + + // Extract the loaded sources from the associated readout, and generate PSF + // Here, we assume the image is background-subtracted + if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", sources)) { + psErrorStackPrint(stderr, "Unable to determine PSF"); + psWarning("Unable to determine PSF."); + psFree(view); + return NAN; + } + psFree(view); + + psMetadata *header = psMetadataLookupMetadata(NULL, photRO->analysis, "PSPHOT.HEADER"); + psAssert(header, "Require header."); + float fwhm = psMetadataLookupF32(NULL, header, "FWHM_MAJ"); + + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) { + psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF"); + } + + return fwhm; +} + +// Scale the kernel parameters according to the PSFs +static bool subScaleKernel(ppSubData *data, // Processing data + psVector *kernelWidths, // Widths for kernel + int *kernelSize, // Size of kernel + int *stampSize // Size of stamps (footprint) + ) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // Nothing to do if pre-calculated kernel exists + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT.KERNELS"); // RO with kernel + if (kernelRO) { + psFree(view); + return true; + } + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + // Input images + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout + pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout + + // Input sources + pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES"); + pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES"); + if (!inSourceRO || !refSourceRO) { + psError(PPSUB_ERR_DATA, false, "Unable to scale kernel, since no sources were provided."); + psFree(view); + return false; + } + + pmDetections *inDetections = psMetadataLookupPtr(NULL, inSourceRO->analysis, "PSPHOT.DETECTIONS"); // Input sources + pmDetections *refDetections = psMetadataLookupPtr(NULL, refSourceRO->analysis, "PSPHOT.DETECTIONS"); // Ref sources + + psFree(view); + + if (!inDetections || !refDetections) { + psError(PPSUB_ERR_DATA, false, "Unable to set FWHM or scale kernel, since no sources were provided."); + return false; + } + + psArray *inSources = inDetections->allSources; + psAssert (inSources, "missing inSources?"); + + psArray *refSources = refDetections->allSources; + psAssert (refSources, "missing refSources?"); + + float inFWHM = psMetadataLookupF32(NULL, inRO->parent->parent->concepts, "CHIP.SEEING"); // FWHM for input + if (!isfinite(inFWHM) || inFWHM == 0.0) { + inFWHM = subImagePSF(data, inRO, inSources); + } + float refFWHM = psMetadataLookupF32(NULL, refRO->parent->parent->concepts, "CHIP.SEEING"); // FWHM for ref + if (!isfinite(refFWHM) || refFWHM == 0.0) { + refFWHM = subImagePSF(data, refRO, refSources); + } + psLogMsg("ppSub", PS_LOG_INFO, "Input FWHM: %f\nReference FWHM: %f\n", inFWHM, refFWHM); + if (!isfinite(inFWHM) || !isfinite(refFWHM)) { +#ifdef SET_QUALITY_INSTEAD_OF_ERROR + psErrorStackPrint(stderr, "Cannot determine FHWM for images, giving up."); + int error = psErrorCodeLast(); // Error code + ppSubDataQuality(data, error, PPSUB_FILES_ALL); + return true; +#else + psError(PPSUB_ERR_DATA, false, "Cannot determine FHWM for images, giving up."); + return false; +#endif + } + + // we need to register the FWHM values for use downstream + pmSubtractionSetFWHMs(inFWHM, refFWHM); + + // is auto-scaling needed? + bool scale = psMetadataLookupBool(NULL, recipe, "SCALE"); + float scaleRef = psMetadataLookupF32(NULL, recipe, "SCALE.REF"); // Reference for scaling + float scaleMin = psMetadataLookupF32(NULL, recipe, "SCALE.MIN"); // Minimum for scaling + float scaleMax = psMetadataLookupF32(NULL, recipe, "SCALE.MAX"); // Maximum for scaling + + if (scale && (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax))) { + psError(PPSUB_ERR_ARGUMENTS, false, + "auto-scale selected but scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in recipe.", + scaleRef, scaleMin, scaleMax); + return false; + } + + pmSubtractionParamScaleOptions(scale, scaleRef, scaleMin, scaleMax); + + // if (!pmSubtractionParamsScale(kernelSize, stampSize, kernelWidths)) { + // psError(PPSUB_ERR_DATA, false, "Unable to scale parameters."); + // return false; + // } + + return true; +} + +pmSubtractionMode subModeFromString (char *string) { + + if (!strcasecmp(string, "AUTO")) return PM_SUBTRACTION_MODE_UNSURE; + if (!strcasecmp(string, "DUAL")) return PM_SUBTRACTION_MODE_DUAL; + if (!strcasecmp(string, "SINGLE_AUTO")) return PM_SUBTRACTION_MODE_SINGLE_AUTO; + if (!strcasecmp(string, "SINGLE1")) return PM_SUBTRACTION_MODE_1; + if (!strcasecmp(string, "SINGLE2")) return PM_SUBTRACTION_MODE_2; + if (!strcasecmp(string, "1")) return PM_SUBTRACTION_MODE_1; + if (!strcasecmp(string, "2")) return PM_SUBTRACTION_MODE_2; + + return PM_SUBTRACTION_MODE_UNSURE; +} + +bool ppSubMatchPSFs(ppSubData *data) +{ + bool mdok = false; + + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images. + if (noConvolve) { + psWarning("not matching PSFs because NOCONVOLVE is TRUE\n"); + return true; + } + + pmFPAview *view = ppSubViewReadout(); // View to readout + + // Input images + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout + pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout + + // Output image holders + pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved + if (!inConv) { + pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Cell for convolved input + inConv = pmReadoutAlloc(cell); + psFree(inConv); + } + pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference convolved + if (!refConv) { + pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.REF.CONV"); // Cell for convolved ref. + refConv = pmReadoutAlloc(cell); + psFree(refConv); + } + + // Load pre-calculated kernel, if available + pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT.KERNELS"); // RO with kernel + + // Sources in image, used for stamps: these must be loaded from previous analysis stages + pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES"); + pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES"); + + pmDetections *inDetections = inSourceRO ? psMetadataLookupPtr(&mdok, inSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Input sources + pmDetections *refDetections = refSourceRO ? psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Ref sources + + psFree(view); + + pmDetections *detections = NULL; // Merged detection set + if (inDetections && refDetections) { + psArray *inSources = inDetections->allSources; + psArray *refSources = refDetections->allSources; + + psAssert (inSources, "missing in sources?"); + psAssert (refSources, "missing ref sources?"); + + detections = pmDetectionsAlloc(); + float radius = psMetadataLookupF32(NULL, recipe, "SOURCE.RADIUS"); // Matching radius + psArray *lists = psArrayAlloc(2); // Source lists + lists->data[0] = psMemIncrRefCounter(inSources); + lists->data[1] = psMemIncrRefCounter(refSources); + // XXX MEH changed to get only match (true to cull single), no apparent need of unmatched sources and can cause trouble + detections->allSources = pmSourceMatchMerge(lists, radius, true); + psFree(lists); + if (!detections->allSources) { + psFree(detections); + int error = psErrorCodeLast(); // Error code + if (error == PM_ERR_OBJECTS) { + psErrorStackPrint(stderr, "Unable to match source lists"); + psWarning("Unable to match source lists --- suspect bad data quality."); + ppSubDataQuality(data, error, PPSUB_FILES_ALL); + return true; + } else { + psError(error, false, "Unable to merge source lists"); + return false; + } + } + } + if (!detections && inDetections) { + detections = psMemIncrRefCounter(inDetections); + } + if (!detections && refDetections) { + detections = psMemIncrRefCounter(refDetections); + } + + psMetadataAddPtr(inConv->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "Merged source list", detections); + psMetadataAddPtr(refConv->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "Merged source list", detections); + psFree(detections); // Drop reference + + int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size + int stride = psMetadataLookupS32(NULL, recipe, "STRIDE"); // Size of convolution patches + float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs + float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing + float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps + + const char *typeStr = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE"); // Kernel type + psAssert(typeStr, "We put it here in ppSubArguments.c"); + pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Type of kernel + if (type == PM_SUBTRACTION_KERNEL_NONE) { + psError(PPSUB_ERR_ARGUMENTS, true, "Unrecognised kernel type: %s", typeStr); + return false; + } + + int size = psMetadataLookupS32(NULL, recipe, "KERNEL.SIZE"); // Kernel half-size + int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order + psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths + psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders + + int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius + int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order + int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel + float penalty = psMetadataLookupF32(NULL, recipe, "PENALTY"); // Penalty for wideness + + int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations + float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold + float kernelErr = psMetadataLookupF32(NULL, recipe, "KERNEL.ERR"); // Relative systematic error in kernel + float normFrac = psMetadataLookupF32(NULL, recipe, "NORM.FRAC"); // Fraction of window for normalisn windw + float sysErr = psMetadataLookupF32(NULL, recipe, "SYS.ERR"); // Relative systematic error in images + float skyErr = psMetadataLookupF32(NULL, recipe, "SKY.ERR"); // Additional error in sky + float covarFrac = psMetadataLookupF32(NULL, recipe, "COVAR.FRAC"); // Fraction for covariance calculation + + float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction + float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor" + + // Options which control the generation of optimal parameters + bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters? + float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search + float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search + float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search + psVector *optWidths = NULL; // Vector with FWHMs for optimum search + if (optimum) { + optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32); + } + int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search + float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search + + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs + psImageMaskType maskPoor = pmConfigMaskGet("CONV.POOR", config); // Bits to mask for poor pixels + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + + // we have three things to control the convolution choice: + + // 1) in the recipe, the keyword DUAL can be TRUE or FALSE. if TRUE, DUAL is used. If not + // DUAL, then we are doing SINGLE convolution (but see CONVOLVE.TARGET below). + + // 2) if the -convolve option is given on the command line, the specified input (1 or 2) is + // the one convolved. Allowed values are 1 or 2. + + // 3) if the -convolve option is NOT given, then the value of the recipe keyword + // CONVOLVE.TARGET is used. this may have the value DUAL, AUTO, SINGLE1, SINGLE2, 1, 2. + // thus, DUAL convolution can be turned on with either DUAL = T or CONVOLVE.TARGET = DUAL + + bool dual = psMetadataLookupBool(&mdok, recipe, "DUAL"); // Dual convolution? + pmSubtractionMode subMode; // Subtraction mode + if (dual) { + subMode = PM_SUBTRACTION_MODE_DUAL; + } else { + char *convolveName = psMetadataLookupStr(&mdok, recipe, "CONVOLVE.TARGET"); // recipe value for target + int convolve = psMetadataLookupS32(&mdok, config->arguments, "-convolve"); // override with command-line option + switch (convolve) { + case 0: + // convolve is 0 if it is not supplied on the command line + subMode = subModeFromString(convolveName); + break; + case 1: + subMode = PM_SUBTRACTION_MODE_1; + break; + case 2: + subMode = PM_SUBTRACTION_MODE_2; + break; + default: + psError(PPSUB_ERR_ARGUMENTS, false, "Invalid value for -convolve"); + return false; + } + } + + if (!subScaleKernel(data, widths, &size, &footprint)) { + psError(PPSUB_ERR_DATA, false, "Unable to scale kernel parameters"); + return false; + } + + int threads = psMetadataLookupS32(NULL, config->arguments, "NTHREADS"); // Number of threads + if (threads > 0) { + pmSubtractionThreadsInit(); + } + + if (inRO->covariance) { + psKernel *truncated = psImageCovarianceTruncate(inRO->covariance, COVAR_FRAC); + psFree(inRO->covariance); + inRO->covariance = truncated; + } + if (refRO->covariance) { + psKernel *truncated = psImageCovarianceTruncate(refRO->covariance, COVAR_FRAC); + psFree(refRO->covariance); + refRO->covariance = truncated; + } + + // Data doesn't exist until it's created in pmSubtraction... + inConv->data_exists = inConv->parent->data_exists = inConv->parent->parent->data_exists = false; + refConv->data_exists = refConv->parent->data_exists = refConv->parent->parent->data_exists = false; + + // Match the PSFs + bool success = false; // Operation was successful? + if (kernelRO) { + success = pmSubtractionMatchPrecalc(inConv, refConv, inRO, refRO, kernelRO->analysis, + stride, kernelErr, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac); + } else { + success = pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize, + spacing, threshold, detections ? detections->allSources : NULL, + data->stamps, type, size, order, widths, orders, inner, ringsOrder, + binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, + normFrac, sysErr, skyErr, kernelErr, covarFrac, + maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode); + } + +# ifdef TESTING + // XXX for testing + psphotSaveImage (NULL, refRO->image, "refRO.im.fits"); + psphotSaveImage (NULL, refRO->variance, "refRO.wt.fits"); + psphotSaveImage (NULL, refRO->mask, "refRO.mk.fits"); + + psphotSaveImage (NULL, inRO->image, "inRO.im.fits"); + psphotSaveImage (NULL, inRO->variance, "inRO.wt.fits"); + psphotSaveImage (NULL, inRO->mask, "inRO.mk.fits"); + + psphotSaveImage (NULL, inConv->image, "inConv.im.fits"); + psphotSaveImage (NULL, inConv->variance, "inConv.wt.fits"); + psphotSaveImage (NULL, inConv->mask, "inConv.mk.fits"); + + psphotSaveImage (NULL, refConv->image, "refConv.im.fits"); + psphotSaveImage (NULL, refConv->variance, "refConv.wt.fits"); + psphotSaveImage (NULL, refConv->mask, "refConv.mk.fits"); +# endif + + psFree(optWidths); + pmSubtractionThreadsFinalize(); + + if (!success) { + int error = psErrorCodeLast(); // Error code + if (error == PM_ERR_STAMPS) { + psErrorStackPrint(stderr, "Unable to find stamps"); + psWarning("Unable to find stamps --- suspect bad data quality."); + ppSubDataQuality(data, error, PPSUB_FILES_ALL); + return true; + } else if (error == PM_ERR_SMALL_AREA) { + psErrorStackPrint(stderr, "Insufficient area for PSF matching"); + psWarning("Insufficient area for PSF matching --- suspect bad data quality."); + ppSubDataQuality(data, error, PPSUB_FILES_ALL); + return true; + // XX } else if (error == PM_ERR_DATA) { + // XX psErrorStackPrint(stderr, "Unable to solve for the kernel to match images"); + // XX psWarning("Failed to find PSF match kernel --- suspect bad data quality."); + // XX ppSubDataQuality(data, error, PPSUB_FILES_ALL); + // XX return true; + } else { + psError(PPSUB_ERR_DATA, false, "Unable to match images."); + return false; + } + } + + // Need to be careful with the normalisation + // We will normalise everything to the normalisation of the *input* image + { + // Since the entries are MULTI, we have to retrieve them differently + psMetadataIterator *regIter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, + "^" PM_SUBTRACTION_ANALYSIS_REGION "$"); + psMetadataIterator *modeIter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, + "^" PM_SUBTRACTION_ANALYSIS_MODE "$"); + psMetadataIterator *normIter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, + "^" PM_SUBTRACTION_ANALYSIS_NORM "$"); + psMetadataItem *regItem; // Item with region + while ((regItem = psMetadataGetAndIncrement(regIter))) { + psAssert(regItem->type == PS_DATA_REGION && regItem->data.V, "Expect region type"); + psRegion *region = regItem->data.V; // Region of interest + psMetadataItem *modeItem = psMetadataGetAndIncrement(modeIter); // Item with mode + psAssert(modeItem && modeItem->type == PS_DATA_S32, "Expect subtraction mode"); + pmSubtractionMode mode = modeItem->data.S32; // Subtraction mode + psMetadataItem *normItem = psMetadataGetAndIncrement(normIter); // Item with normalisation + psAssert(normItem && normItem->type == PS_DATA_F32 && isfinite(normItem->data.F32), + "Expect normalisation"); + float norm = normItem->data.F32; // Normalisation + + switch (mode) { + case PM_SUBTRACTION_MODE_1: // Convolved the input to match template + case PM_SUBTRACTION_MODE_DUAL: // Convolved both; template should have flux conserved + psLogMsg("ppSub", PS_LOG_INFO, "Correcting image for normalisation of %f\n", norm); + normaliseRegion(inConv->image, region, 1.0 / norm); + normaliseRegion(refConv->image, region, 1.0 / norm); + normaliseRegion(inConv->variance, region, 1.0 / PS_SQR(norm)); + normaliseRegion(refConv->variance, region, 1.0 / PS_SQR(norm)); + break; + case PM_SUBTRACTION_MODE_2: // Convolved the template to match input + // We're already happy! + psLogMsg("ppSub", PS_LOG_INFO, "Image normalisation is correct\n"); + break; + default: + psAbort("Invalid subtraction mode: %x", mode); + } + } + psFree(regIter); + psFree(modeIter); + psFree(normIter); + } + + + pmConceptsCopyFPA(inConv->parent->parent->parent, inRO->parent->parent->parent, true, true); + pmConceptsCopyFPA(refConv->parent->parent->parent, refRO->parent->parent->parent, true, true); + + if (inConv->covariance) { + psKernel *truncated = psImageCovarianceTruncate(inConv->covariance, COVAR_FRAC); + psFree(inConv->covariance); + inConv->covariance = truncated; + } + if (refConv->covariance) { + psKernel *truncated = psImageCovarianceTruncate(refConv->covariance, COVAR_FRAC); + psFree(refConv->covariance); + refConv->covariance = truncated; + } + + if (inConv->variance) { + psImageCovarianceTransfer(inConv->variance, inConv->covariance); + } + if (refConv->variance) { + psImageCovarianceTransfer(refConv->variance, refConv->covariance); + } + + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutForcedPhot.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutForcedPhot.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutForcedPhot.c (revision 41094) @@ -0,0 +1,103 @@ +/** @file ppSubReadoutForcedphot.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +// run forced photometry on the image in 'targetName' at the positions of the sources in 'sourceName' +// this function is only called if (data->forcedPhot1 || data->forcedPhot2) (ppSubLoop.c:186,194) +bool ppSubReadoutForcedPhot(const char *outputName, const char *targetName, const char *sourceName, ppSubData *data) +{ + bool foundDetections = false; + bool mdok = false; + + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // select the view of interest + pmFPAview *view = ppSubViewReadout(); // View to readout + + pmFPAfile *sourceFile = psMetadataLookupPtr(&mdok, config->files, sourceName); // positive image sources + pmFPAfile *targetFile = psMetadataLookupPtr(&mdok, config->files, targetName); // diff image pixels + pmFPAfile *outputFile = psMetadataLookupPtr(&mdok, config->files, outputName); // result file + + psAssert (sourceFile, "failed to construct sourceName %s", sourceName); + psAssert (targetFile, "failed to construct targetName %s", targetName); + psAssert (outputFile, "failed to construct outputName %s", outputName); + + // copy the image data to PSPHOT.INPUT so that psphotReadoutForcedKnownSources does + // not affect the prior results on targetName (ie, diff detections) + pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file + if (!pmFPACopy(photFile->fpa, targetFile->fpa)) { + psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry"); + psFree(view); + return false; + } + + // select the reference position sources from sourceName + pmReadout *sourceRO = pmFPAviewThisReadout(view, sourceFile->fpa); + psAssert(sourceRO, "programming error: source readout not defined"); + + pmDetections *sourceDet = psMetadataLookupPtr(NULL, sourceRO->analysis, "PSPHOT.DETECTIONS"); // Sources + if (!sourceDet) { + // XXX remove the pixels from photFile? + // XXX other cleanup operations? + psFree(view); + return true; + } + psArray *sources = sourceDet->allSources; + psAssert (sources, "missing sources?"); + + // grab the PSF information from the pmFPAfile PSPHOT.PSF.LOAD + pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file + ppSubCopyPSF (photFile, psfFile, view); + + // psphotSaveImage (photFile->fpa->hdu->header, photRO->image, "findsrc.im.fits"); + // psphotSaveImage (photFile->fpa->hdu->header, photRO->variance, "findsrc.wt.fits"); + // psphotSaveImage (photFile->fpa->hdu->header, photRO->mask, "findsrc.mk.fits"); + + // erase the overlays from a previous psphot-related step + if (pmVisualIsVisual()) { + // psphotVisualEraseOverlays (1, "all"); + } + + if (!psphotReadoutForcedKnownSources(config, view, "PSPHOT.INPUT", sources)) { + // This is likely a data quality issue + // XXX Split into multiple cases using error codes? + psErrorStackPrint(stderr, "Unable to perform photometry on image"); + psWarning("Unable to perform photometry on image --- suspect bad data quality."); + ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); + } + // save the outputs on PPSUB.POS1.SOURCES + if (!psphotCopyResults (&foundDetections, outputFile, photFile, view)) { + psError(PPSUB_ERR_PROG, false, "Unable to copy psphot outputs"); + psFree(view); + return false; + } + + float newTime = psTimerClear("PPSUB_PHOT"); // Time for photometry + float oldTime = psMetadataLookupF32(&mdok, data->stats, "TIME_PHOT"); // Previous time for photometry + float elapsed = isfinite(oldTime) ? oldTime + newTime : newTime; + psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_PHOT", PS_META_REPLACE, "Time to do photometry", elapsed); + + psFree(view); + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutInverse.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutInverse.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutInverse.c (revision 41094) @@ -0,0 +1,59 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppSub.h" + +bool ppSubReadoutInverse(pmConfig *config) +{ + psAssert(config, "Require configuration"); + // Check configuration for convolve option + psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE); + bool noConvolve = psMetadataLookupBool(NULL, recipe, "NOCONVOLVE"); // Do not use convolved images. + + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); + pmReadout *invRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INVERSE"); + + invRO->image = (psImage*)psBinaryOp(invRO->image, outRO->image, "*", psScalarAlloc(-1.0, PS_TYPE_F32)); + invRO->mask = psMemIncrRefCounter(outRO->mask); + invRO->variance = psMemIncrRefCounter(outRO->variance); + invRO->covariance = psMemIncrRefCounter(outRO->covariance); + invRO->analysis = psMetadataCopy(invRO->analysis, outRO->analysis); + + invRO->data_exists = invRO->parent->data_exists = invRO->parent->parent->data_exists = true; + + // Get concepts from reference + pmFPAfile *refFile; + if (noConvolve) { + refFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF"); // File with concepts + } + else { + refFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF.CONV"); // File with concepts + } + pmFPA *invFPA = invRO->parent->parent->parent; // Inverse FPA + pmConceptsCopyFPA(invFPA, refFile->fpa, true, true); + + // Get astrometry from (forward) subtraction + pmChip *outChip = outRO->parent->parent; // Output chip + pmFPA *outFPA = outChip->parent; // Output FPA + pmChip *invChip = invRO->parent->parent; // Inverse chip + pmHDU *invHDU = invFPA->hdu; // Inverse HDU + if (!pmAstromWriteWCS(invHDU->header, outFPA, outChip, WCS_TOLERANCE)) { + psError(psErrorCodeLast(), false, "Unable to write WCS astrometry to PPSUB.INVERSE."); + psFree(view); + return false; + } + // Read from newly written astrometry so that it exists in the "inverse" FPA (for sources) + if (!pmAstromReadWCS(invFPA, invChip, invHDU->header, 1.0)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry."); + psFree(view); + return false; + } + + psFree(view); + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutJpeg.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutJpeg.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutJpeg.c (revision 41094) @@ -0,0 +1,143 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppSub.h" + +bool ppSubReadoutJpeg(pmConfig *config) +{ + psAssert(config, "Require configuration"); + + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Readout to jpeg + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + // Generate binned JPEGs + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + + int bin1 = psMetadataLookupS32(NULL, recipe, "BIN1"); // First binning level + int bin2 = psMetadataLookupS32(NULL, recipe, "BIN2"); // Second binning level + + // Target cells + pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG1"); // Rebinned cell once + pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG2"); // Rebinned cell twice + psFree(view); + + pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts + if (!pmReadoutRebin(ro1, outRO, maskBad, bin1, bin1)) { + psError(PPSUB_ERR_DATA, false, "Unable to bin output (1st binning)"); + psFree(ro1); + psFree(ro2); + return false; + } + if (!pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) { + psError(PPSUB_ERR_DATA, false, "Unable to bin output (2nd binning)"); + psFree(ro1); + psFree(ro2); + return false; + } + psFree(ro1); + psFree(ro2); + + return true; +} + +bool ppSubResidualSampleJpeg(pmConfig *config) +{ + return true; + + pmFPAview *view = ppSubViewReadout(); // View to readout + + // we save sample difference stamps on the convolved image readout->analysis metadata + pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved + + psArray *samples = psArrayAllocEmpty(16); // Array of sample stamp images + + // we may have multiple entries with the same name; pull them off into a separate array: + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", "SUBTRACTION.SAMPLE.STAMP.SET"); + psMetadataIterator *iter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_ARRAY); + psArray *sampleStamps = item->data.V; + samples = psArrayAdd(samples, 16, sampleStamps); + } + psFree(iter); + psAssert (samples, "no sample stamps?"); + psAssert (samples->n, "no sample stamps?"); + + // get the kernel sizes + psArray *kernels = samples->data[0]; + psAssert (kernels, "no valid kernel?"); + psAssert (kernels->n, "no valid kernel?"); + + psImage *kernel = kernels->data[0]; + psAssert (kernel, "missing valid kernel?"); + + int DX = kernel->numCols; + int DY = kernel->numRows; + + // each array contains up to 9 sample stamps. generate an image mosaicking these together in 3x3 blocks. + int innerBorder = 1; + int outerBorder = 2; + int NXblock = sqrt(samples->n); + int NYblock = samples->n / NXblock; + if (samples->n % NXblock) NYblock ++; + + int NXpix = NXblock * (3 * (DX + innerBorder) + outerBorder); + int NYpix = NYblock * (3 * (DY + innerBorder) + outerBorder); + + // output cell + pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.RESID.JPEG"); + pmReadout *readout = pmReadoutAlloc(cell); + readout->image = psImageAlloc(NXpix, NYpix, PS_TYPE_F32); + + cell->data_exists = true; + cell->parent->data_exists = true; + + psImageInit (readout->image, 0.0); + + for (int i = 0; i < samples->n; i++) { + + int xBlock = i % NXblock; + int yBlock = i / NYblock; + + psArray *kernels = samples->data[i]; + + for (int j = 0; j < kernels->n; j++) { + + psImage *kernel = kernels->data[j]; + + int xSub = j % 3; + int ySub = j / 3; + + int xPix = xBlock * (3 * (DX + innerBorder) + outerBorder) + xSub * (DX + innerBorder); + int yPix = yBlock * (3 * (DX + innerBorder) + outerBorder) + ySub * (DY + innerBorder); + + for (int y = 0; y < kernel->numRows; y++) { + for (int x = 0; x < kernel->numCols; x++) { + readout->image->data.F32[y + yPix][x + xPix] = kernel->data.F32[y][x]; + } + } + } + } + + { + psFits *fits = psFitsOpen ("resid.stamps.fits", "w"); + psFitsWriteImage (fits, NULL, readout->image, 0, NULL); + psFitsClose (fits); + } + + psFree(view); + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutPhotometry.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutPhotometry.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutPhotometry.c (revision 41094) @@ -0,0 +1,188 @@ +/** @file ppSubReadoutPhotometry.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +bool ppSubReadoutPhotometry(const char *name, ppSubData *data) +{ + bool mdok = false; + + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + if (!data->photometry) { + return true; + } + + // select the view of interest + pmFPAview *view = ppSubViewReadout(); // View to readout + + // psphotReadoutMinimal performs the photometry analysis on PSPHOT.INPUT; we need to move + // around the pointers so PSPHOT.INPUT corresponds to the output image of interest (on one + // pass this is the subtraction image, in another it is negative of the subtraction + pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file + pmFPAfile *inFile = psMetadataLookupPtr(&mdok, config->files, name); // Input file + if (!pmFPACopy(photFile->fpa, inFile->fpa)) { + psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry"); + psFree(view); + return false; + } + + // psphot may need the PHU header + if (!photFile->fpa->hdu) { + photFile->fpa->hdu = psMemIncrRefCounter (inFile->fpa->hdu); + } + + // drop references to PSPHOT.DETECTIONS on both of these (why is this needed for both??) + pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, name); // Readout with image and sources + if (psMetadataLookup(inRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(inRO->analysis, "PSPHOT.DETECTIONS"); + } + + // grab the PSF information from the pmFPAfile PSPHOT.PSF.LOAD + pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file + + ppSubCopyPSF (photFile, psfFile, view); + + // psphotSaveImage (photFile->fpa->hdu->header, photRO->image, "findsrc.im.fits"); + // psphotSaveImage (photFile->fpa->hdu->header, photRO->variance, "findsrc.wt.fits"); + // psphotSaveImage (photFile->fpa->hdu->header, photRO->mask, "findsrc.mk.fits"); + + // erase the overlays from a previous psphot-related step + if (pmVisualIsVisual()) { + // psphotVisualEraseOverlays (1, "all"); + } + + if (!psphotReadoutMinimal(config, view, "PSPHOT.INPUT")) { + // This is likely a data quality issue + // XXX Split into multiple cases using error codes? + psErrorStackPrint(stderr, "Unable to perform photometry on image"); + psWarning("Unable to perform photometry on image --- suspect bad data quality."); + ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); + } + + // If no sources were found, there's no error, but we want to trigger 'bad quality' + psWarning("no sources found: why is this being set to bad quality??"); + pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // Sources + if (!detections) { + ppSubDataQuality(data, PSPHOT_ERR_DATA, PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); + } + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + // a likely source of false positives is poor subtractions. this results in + // detections in the wings (or cores) of bright(er) stars found in both images. + // flag detections based on their distance from the bright(er) input sources. + bool subInverse = !strcasecmp(name, "PPSUB.INVERSE"); + ppSubFlagNeighbors (config, view, sources, subInverse); + + if (data->stats) { + bool mdok; + int numSources = psMetadataLookupS32(&mdok, data->stats, "NUM_SOURCES"); // Number of sources + numSources += sources ? sources->n : 0; + psMetadataAddS32(data->stats, PS_LIST_TAIL, "NUM_SOURCES", PS_META_REPLACE, + "Total number of sources detected", numSources); + float newTime = psTimerClear("PPSUB_PHOT"); // Time for photometry + float oldTime = psMetadataLookupF32(&mdok, data->stats, "TIME_PHOT"); // Previous time for photometry + psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_PHOT", PS_META_REPLACE, "Time to do photometry", + isfinite(oldTime) ? oldTime + newTime : newTime); + } + + if (!data->quality) { + if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.DETECTIONS")) { + psError(PPSUB_ERR_PROG, false, "Unable to copy PSPHOT.DETECTIONS"); + psFree(view); + return false; + } + if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.HEADER")) { + psError(PPSUB_ERR_PROG, false, "Unable to copy PSPHOT.HEADER"); + psFree(view); + return false; + } + if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, PM_DETEFF_ANALYSIS)) { + psError(PPSUB_ERR_PROG, false, "Unable to copy Detection Efficiency"); + psFree(view); + return false; + } + + // Ensure photometry information is put in the header + pmHDU *hdu = pmHDUFromReadout(inRO); // HDU for readout + if (hdu) { + psMetadata *photHeader = psMetadataLookupMetadata(NULL, inRO->analysis, "PSPHOT.HEADER"); // Header + hdu->header = psMetadataCopy(hdu->header, photHeader); + } + } + + psFree(view); + return true; +} + +#ifdef TESTING +// Record data about sources: not everything gets into the output CMF files + { + pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources + pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // Sources + psArray *sources = detections->allSources; + FILE *sourceFile = fopen("sources.dat", "w"); // File for sources + fprintf(sourceFile, + "# x y mag mag_err psf_chisq cr_nsigma ext_nsigma psf_qf flags m_x m_y m_xx m_xy m_yy\n"); + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (!source) { + continue; + } + + float x, y; // Position of source + float chi2; // chi^2 for source + if (source->modelPSF) { + x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + chi2 = source->modelPSF->chisq; + } else if (source->peak) { + x = source->peak->xf; + y = source->peak->yf; + chi2 = NAN; + } else { + psWarning("No position available for source."); + continue; + } + + float xMoment = NAN, yMoment = NAN, xxMoment = NAN, xyMoment = NAN, yyMoment = NAN; + if (source->moments) { + xMoment = source->moments->Mx; + yMoment = source->moments->My; + xxMoment = source->moments->Mxx; + xyMoment = source->moments->Mxy; + yyMoment = source->moments->Myy; + } + + fprintf(sourceFile, "%f %f %f %f %f %f %f %f %d %f %f %f %f %f\n", + x, y, source->psfMag, source->psfMagErr, chi2, source->crNsigma, source->extNsigma, + source->pixWeight, source->mode, xMoment, yMoment, xxMoment, xyMoment, yyMoment); + } + fclose(sourceFile); + } +#endif Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutStats.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutStats.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutStats.c (revision 41094) @@ -0,0 +1,51 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppSub.h" + + +bool ppSubReadoutStats(ppSubData *data) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + if (!data->statsFile) { + // Nothing to do + return true; + } + + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSUB.OUTPUT"); // Output file + if (!output) { + psError(PPSUB_ERR_PROG, true, "Unable to find file PPSUB.OUTPUT.\n"); + return false; + } + psImageMaskType maskValue = pmConfigMaskGet("MASK.VALUE", config); + pmFPAview *view = ppSubViewReadout(); // View to readout + ppStatsFPA(data->stats, output->fpa, view, maskValue, config); + + pmReadout *outRO = pmFPAviewThisReadout(view, output->fpa); // Readout of interest + psFree(view); + + // Statistics on the matching + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MODE); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_STAMPS); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_MEAN); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_RMS); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_NORM); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_BGDIFF); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MX); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MY); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXX); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXY); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MYY); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_CONVOL_MAX); + + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutSubtract.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutSubtract.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubReadoutSubtract.c (revision 41094) @@ -0,0 +1,155 @@ +/** @file ppSubReadoutSubtract.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppSub.h" + +bool ppSubReadoutSubtract(pmConfig *config) +{ + psAssert(config, "Require configuration"); + + // Look up recipe values + bool mdok = false; // Status of MD lookup + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSub + psAssert(recipe, "We checked this earlier, so it should be here."); + + bool reverse = psMetadataLookupBool(&mdok, config->arguments, "REVERSE"); // Reverse sense of subtraction? + bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images. + bool addPair = psMetadataLookupBool(&mdok, recipe, "ADD.NOT.SUBTRACT"); // add instead of subtracting + + pmFPAview *view = ppSubViewReadout(); // View to readout + + // Subtraction is: minuend - subtrahend + pmReadout *minuend = NULL; // Positive image + pmReadout *subtrahend = NULL; // Negative image + if (!noConvolve) { + printf("Using Convolved images because NOCONVOLVE is FALSE\n"); + if (reverse) { + minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); + subtrahend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); + } else { + minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); + subtrahend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); + } + } + else { + psWarning("Not using Convolved images because NOCONVOLVE is TRUE\n"); + if (reverse) { + minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); + subtrahend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); + } else { + minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); + subtrahend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); + } + } + + // Do the actual subtraction + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); + + if (addPair) { + outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "+", subtrahend->image); + } else { + outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image); + } + outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask); + outRO->variance = (psImage*)psBinaryOp(outRO->variance, minuend->variance, "+", subtrahend->variance); + + // NAN the masked pixels in the diff image (pixels masked in A are not yet NAN'ed in B) + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config) | pmConfigMaskGet("BLANK", config); // Bits to mask in inputs + for (int iy = 0; iy < outRO->image->numRows; iy++) { + for (int ix = 0; ix < outRO->image->numCols; ix++) { + if ((outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) == 0) continue; + outRO->image->data.F32[iy][ix] = NAN; + } + } + + // Measure the variance scales + psStats *varStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance images + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + psImageBackground(varStats, NULL, minuend->variance, minuend->mask, maskVal, rng); + float minuendVar = varStats->robustMedian; // Mean variance for minuend + psImageBackground(varStats, NULL, subtrahend->variance, subtrahend->mask, maskVal, rng); + float subtrahendVar = varStats->robustMedian; // Mean variance for subtrahend + psFree(varStats); + psFree(rng); + + // XXX EAM : I suspect that the weighted averaging is giving the wrong answer for + // single-direction PSF matching (and maybe DUAL as well). I am testing the difference by + // generating A+B images instead of A-B images and then checking the significance of the + // sources in the image + + // Combine the covariances + // These are weighted by the appropriate mean variance. This is probably not perfectly correct, but it + // does seem to reproduce the correct magnitude limit in psphot. + psArray *covars = psArrayAlloc(2); // Covariance pseudo-matrices + psVector *covarWeights = psVectorAlloc(2, PS_TYPE_F32); // Weights for covariances + covars->data[0] = psMemIncrRefCounter(minuend->covariance); + covars->data[1] = psMemIncrRefCounter(subtrahend->covariance); + covarWeights->data.F32[0] = minuendVar; + covarWeights->data.F32[1] = subtrahendVar; +# if (0) + outRO->covariance = psImageCovarianceAverageWeighted(covars, covarWeights); +# else + outRO->covariance = psImageCovarianceAverage(covars); +# endif + psFree(covars); + psFree(covarWeights); + + psImageCovarianceTransfer(outRO->variance, outRO->covariance); + + outRO->data_exists = true; + outRO->parent->data_exists = true; + outRO->parent->parent->data_exists = true; + + pmSubtractionVisualShowSubtraction(minuend->image, subtrahend->image, outRO->image); + + // Copy concepts from the input to the output + pmFPAfile *inFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.INPUT"); // Input file + pmFPA *inFPA = inFile->fpa; // Input FPA + pmFPAfile *outFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.OUTPUT"); // Output file + pmFPA *outFPA = outFile->fpa; // Output FPA + if (!pmConceptsCopyFPA(outFPA, inFPA, true, true)) { + psError(PPSUB_ERR_CONFIG, false, "Unable to copy concepts from input to output."); + psFree(outRO); + psFree(view); + return false; + } + + // Copy astrometry over + // It should find its way into the output images and photometry + pmHDU *inHDU = inFPA->hdu; // Input HDU + pmHDU *outHDU = outFPA->hdu; // Output HDU + pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSUB.OUTPUT"); // Output chip + psFree(view); + + if (!outHDU || !inHDU) { + psError(PPSUB_ERR_PROG, true, "Unable to find HDU at FPA level to copy astrometry."); + return false; + } + if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry from input FPA."); + return false; + } + if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { + psError(psErrorCodeLast(), false, "Unable to write WCS astrometry to output FPA."); + return false; + } + + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubSetMasks.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubSetMasks.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubSetMasks.c (revision 41094) @@ -0,0 +1,122 @@ +/** @file ppSubSetMasks.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +bool ppSubSetMasks(pmConfig *config) +{ + psAssert(config, "Require configuration"); + + psImageMaskType maskValue, markValue; // Mask values + if (!pmConfigMaskSetBits(&maskValue, &markValue, config)) { + psError(PPSUB_ERR_CONFIG, false, "Unable to determine mask value."); + return false; + } + + // Set the mask bits needed by psphot (in psphot recipe) + psphotSetMaskRecipe(config, maskValue, markValue); + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + psImageMaskType satValue = pmConfigMaskGet("SAT", config); + psAssert(satValue, "SAT must be non-zero"); + + psImageMaskType lowValue = pmConfigMaskGet("LOW", config); + if (!lowValue) { + // Look up old name for backward compatability + lowValue = pmConfigMaskGet("BAD", config); + } + psAssert(lowValue, "LOW or BAD must be non-zero"); + + // Input images + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout + pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout + + psImage *input = inRO->image; // Input image + psImage *reference = refRO->image; // Reference image + PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false); + int numCols = input->numCols, numRows = input->numRows; // Size of image + + // Generate masks if they don't exist + if (!inRO->mask) { + if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) { + pmReadoutSetMask(inRO, satValue, lowValue); + } else { + inRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + psImageInit(inRO->mask, 0); + } + } + if (!refRO->mask) { + if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) { + pmReadoutSetMask(refRO, satValue, lowValue); + } else { + refRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + psImageInit(refRO->mask, 0); + } + } + + // Mask the NAN values (USE BLANK instead of SAT?) + if (!pmReadoutMaskInvalid(inRO, maskValue, satValue)) { + psError(PPSUB_ERR_DATA, false, "Unable to mask non-finite pixels in input."); + psFree(view); + return false; + } + if (!pmReadoutMaskInvalid(refRO, maskValue, satValue)) { + psError(PPSUB_ERR_DATA, false, "Unable to mask non-finite pixels in reference."); + psFree(view); + return false; + } + +#if 0 + // Interpolation over bad pixels: this takes a while! + bool mdok = false; + psString interpModeStr = psMetadataLookupStr(&mdok, recipe, "INTERPOLATION"); // Interpolation mode + psImageInterpolateMode interpMode = psImageInterpolateModeFromString(interpModeStr); // Interp + if (interpMode == PS_INTERPOLATE_NONE) { + psError(PPSUB_ERR_CONFIG, false, "Unknown interpolation mode: %s", interpModeStr); + psFree(view); + return false; + } + float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor" + + psImageMaskType maskPoor = pmConfigMaskGet("CONV.POOR", config); // Bits to mask for poor pixels + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + + // Interpolate over bad pixels, so the bad pixels don't explode + if (!pmReadoutInterpolateBadPixels(inRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) { + psError(PPSUB_ERR_DATA, false, "Unable to interpolate bad pixels for input image."); + psFree(view); + return false; + } + if (!pmReadoutInterpolateBadPixels(refRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) { + psError(PPSUB_ERR_DATA, false, "Unable to interpolate bad pixels for reference image."); + psFree(view); + return false; + } + maskVal |= maskBad; +#endif + + psFree(view); + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubThreshold.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubThreshold.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubThreshold.c (revision 41094) @@ -0,0 +1,147 @@ +/** @file ppSubThreshold.c + * + * @brief Mask pixels below threshold + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.3 $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppSub.h" + +// Apply low threshold to image +static bool lowThreshold( + pmReadout *ro, // Readout to threshold + float thresh, // Threshold to apply + psImageMaskType maskIgnore, // Ignore pixels with this mask + psImageMaskType maskThresh, // Give pixels this mask if below threshold + psRegion *region, // Region of interest + const char *description // Description of image + ) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_IMAGE(ro, false); + PM_ASSERT_READOUT_MASK(ro, false); + PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false); + + psImage *image = ro->image; // Image + psImage *mask = ro->mask; // Mask + + psImage *subImage = psImageSubset(image, *region); // Image with region of interest + psImage *subMask = psImageSubset(mask, *region); // Maks with region of interest + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + if (!psImageBackground(stats, NULL, subImage, subMask, maskIgnore, rng)) { + psError(PPSUB_ERR_DATA, false, "Unable to determine threshold."); + psFree(rng); + psFree(stats); + return false; + } + psFree(rng); + + psFree(subImage); + psFree(subMask); + + float threshold = stats->robustMedian - thresh * stats->robustStdev; // Threshold below which to clip + psFree(stats); + psLogMsg("ppSub", PS_LOG_INFO, "Masking pixels below %f in %s", threshold, description); + + for (int y = region->y0; y < region->y1; y++) { + for (int x = region->x0; x < region->x1; x++) { + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskIgnore) { + continue; + } + if (image->data.F32[y][x] < threshold) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskThresh; + } + } + } + + return true; +} + + +bool ppSubLowThreshold(ppSubData *data) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs + psImageMaskType maskThresh = pmConfigMaskGet("LOW", config); // Bits to mask for low pixels + + float thresh = psMetadataLookupF32(NULL, recipe, "LOW.THRESHOLD"); // Threshold for masking + if (!isfinite(thresh)) { + return true; + } + + pmFPAview *view = ppSubViewReadout(); // View to readout + + // Input images + bool noConvolve = psMetadataLookupBool(NULL, recipe, "NOCONVOLVE"); + + pmReadout *in; + if (noConvolve) { + in = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input image + } + else { + in = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input image + } + if (!in) { + psError(PPSUB_ERR_UNKNOWN, false, "Unable to find readout."); + psFree(view); + return false; + } + + pmReadout *ref; + if (noConvolve) { + ref = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference image + } + else { + ref = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference image + } + if (!ref) { + psError(PPSUB_ERR_UNKNOWN, false, "Unable to find readout."); + psFree(view); + return false; + } + + psMetadataIterator *regIter = psMetadataIteratorAlloc(in->analysis, PS_LIST_HEAD, + "^" PM_SUBTRACTION_ANALYSIS_REGION "$"); + psMetadataItem *regItem; // Item with region + while ((regItem = psMetadataGetAndIncrement(regIter))) { + psAssert(regItem->type == PS_DATA_REGION && regItem->data.V, "Expect region type"); + psRegion *region = regItem->data.V; // Region of interest + if (!lowThreshold(in, thresh, maskVal, maskThresh, region, "input convolved image")) { + psError(psErrorCodeLast(), false, "Unable to threshold input image."); + psFree(view); + return false; + } + if (!lowThreshold(ref, thresh, maskVal, maskThresh, region, "reference convolved image")) { + psError(psErrorCodeLast(), false, "Unable to threshold input image."); + psFree(view); + return false; + } + } + psFree(regIter); + + psFree(view); + + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubVarianceRescale.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubVarianceRescale.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubVarianceRescale.c (revision 41094) @@ -0,0 +1,74 @@ +/** @file ppSubVarianceRescale.c + * + * @brief measure the background signal/noise distribution and rescale the variance if needed + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppSub.h" + +bool ppSubVarianceRescale(pmConfig *config, ppSubData *data) +{ + psAssert(config, "Require configuration"); + + bool mdok; // Status of metadata lookups + + psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSub + psAssert(recipe, "Need PPSUB recipe"); + + if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true; + + int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); + if (!mdok) { + psError(PPSUB_ERR_ARGUMENTS, true, "RENORM.NUM is not set in the recipe"); + return false; + } + float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN"); + if (!mdok) { + psError(PPSUB_ERR_ARGUMENTS, true, "RENORM.MIN is not set in the recipe"); + return false; + } + float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX"); + if (!mdok) { + psError(PPSUB_ERR_ARGUMENTS, true, "RENORM.MAX is not set in the recipe"); + return false; + } + + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask + + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *readout = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output image + + if (!readout->variance) { + // Nothing to renormalise + psWarning("Renormalisation of the variance requested, but no variance provided."); + psFree(view); + return true; + } + + if (!pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid)) { + psErrorStackPrint(stderr, "Unable to renormalise variances"); + psWarning("Unable to renormalise variances --- suspect bad data quality."); + // Allow the convolved and subtracted images to be written + ppSubDataQuality(data, PPSUB_ERR_VARIANCE, + PPSUB_FILES_INPUT | PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV | + PPSUB_FILES_PSF | PPSUB_FILES_PHOT); + } + + psFree(view); + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubVersion.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubVersion.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubVersion.c (revision 41094) @@ -0,0 +1,129 @@ +/** @file ppSubVersion.c + * + * @brief + * + * @ingroup ppSub + * + * @author IfA + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#include "ppSub.h" +#include "ppSubVersionDefinitions.h" + +#ifndef PPSUB_VERSION +#error "PPSUB_VERSION is not set" +#endif +#ifndef PPSUB_BRANCH +#error "PPSUB_BRANCH is not set" +#endif +#ifndef PPSUB_SOURCE +#error "PPSUB_SOURCE is not set" +#endif + +psString ppSubVersion(void) +{ + char *value = NULL; + psStringAppend(&value, "%s@%s", PPSUB_BRANCH, PPSUB_VERSION); + return value; +} + +psString ppSubSource(void) +{ + return psStringCopy(PPSUB_SOURCE); +} + +psString ppSubVersionLong(void) +{ + psString version = ppSubVersion(); // Version, to return + psString source = ppSubSource(); // Source + + psStringPrepend(&version, "ppSub "); + psStringAppend(&version, " from %s, built %s, %s", source, __DATE__, __TIME__); + psFree(source); + +#ifdef __OPTIMIZE__ + psStringAppend(&version, " optimised"); +#else + psStringAppend(&version, " unoptimised"); +#endif + + return version; +}; + + +bool ppSubVersionHeader(psMetadata *header) +{ + PS_ASSERT_METADATA_NON_NULL(header, false); + + psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now + psString timeString = psTimeToISO(time); // The time in an ISO string + psFree(time); + psString history = NULL; // History string + psStringAppend(&history, "ppSub at %s", timeString); + psFree(timeString); + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history); + psFree(history); + + psLibVersionHeader(header); + psModulesVersionHeader(header); + psphotVersionHeader(header); + ppStatsVersionHeader(header); + + psString version = ppSubVersion(); // Software version + psString source = ppSubSource(); // Software source + + psMetadataAddStr(header, PS_LIST_TAIL, "PPSUB_V", PS_META_REPLACE, NULL, PPSUB_VERSION); + + psStringPrepend(&version, "ppSub version: "); + psStringPrepend(&source, "ppSub source: "); + + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, version); + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, source); + + psFree(version); + psFree(source); + + return true; +} + +void ppSubVersionPrint(void) +{ + psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now + psString timeString = psTimeToISO(time); // The time in an ISO string + psFree(time); + psLogMsg("ppSub", PS_LOG_INFO, "ppSub at %s", timeString); + psFree(timeString); + + psString pslib = psLibVersionLong();// psLib version + psString psmodules = psModulesVersionLong(); // psModules version + psString psphot = psphotVersionLong(); // psphot version + psString ppStats = ppStatsVersionLong(); // ppStats version + psString ppSub = ppSubVersionLong(); // ppSub version + + psLogMsg("ppSub", PS_LOG_INFO, "%s", pslib); + psLogMsg("ppSub", PS_LOG_INFO, "%s", psmodules); + psLogMsg("ppSub", PS_LOG_INFO, "%s", psphot); + psLogMsg("ppSub", PS_LOG_INFO, "%s", ppStats); + psLogMsg("ppSub", PS_LOG_INFO, "%s", ppSub); + + psFree(pslib); + psFree(psmodules); + psFree(psphot); + psFree(ppStats); + psFree(ppSub); + + return; +} Index: /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubVersionDefinitions.h.in =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubVersionDefinitions.h.in (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/src/ppSubVersionDefinitions.h.in (revision 41094) @@ -0,0 +1,8 @@ +#ifndef PPSUB_VERSION_DEFINITIONS_H +#define PPSUB_VERSION_DEFINITIONS_H + +#define PPSUB_VERSION @PPSUB_VERSION@ // SVN version +#define PPSUB_BRANCH @PPSUB_BRANCH@ // SVN branch +#define PPSUB_SOURCE @PPSUB_SOURCE@ // SVN source + +#endif Index: /branches/ccl_branches/ipponly-20191108/ppSub/test/fake.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSub/test/fake.c (revision 41094) +++ /branches/ccl_branches/ipponly-20191108/ppSub/test/fake.c (revision 41094) @@ -0,0 +1,182 @@ +#include +#include + +// Quick and dirty program to generate cross-directed PSFs for testing ppSub. +// To build: gcc fake.c -o fake -g -std=gnu99 `pslib-config --cflags --libs` + + +// PSF for images +#define X_AXIS_1 5.0 // Length of x axis (FWHM) for image 1 +#define Y_AXIS_1 3.0 // Length of y axis (FWHM) for image 1 +#define X_AXIS_2 5.0 // Length of x axis (FWHM) for image 1 +#define Y_AXIS_2 5.0 // Length of y axis (FWHM) for image 1 + +// Image parameters +#define NUMCOLS 1024 // Number of columns +#define NUMROWS 1024 // Number of rows +#define IMAGE_SCALE 0.2 // Arcsec per pixel + +// Source parameters +#define SOURCES_LUM 0.4 // Luminosity function +#define SOURCES_MAG 14.5 // Magnitude for density +#define SOURCES_DENSITY 30000.0 // Source density at nominated magnitude (per deg^2) +#define SOURCES_ZP 25.0 // Magnitude ZP +#define SOURCES_GAUSSIAN 0 // Gaussian sources (boolean)? + +// Noise properties +#define NOISE_1 10.0 // Noise in image 1 +#define NOISE_2 10.0 // Noise in image 2 +#define NOISE_FRAC 0.1 // Go out to this fraction of the noise + +// Conversion factor by which to multiply sigma to get FWHM. Roughly 2.35 +#define SIGMA_FWHM (2.0*sqrt((double)2.0*(log((double)2.0)))) + +// Add a star to an image +static void addstar(float x, float y, // Coordinates of source + float flux, // Total (integrated) flux of source + float xSigma, float ySigma, // Gaussian width of source + float noise, // Lowest noise level to consider + psImage *image // Image to which to add source + ) +{ + float peak = flux / 2.0 / M_PI / sqrtf(PS_SQR(xSigma) + PS_SQR(ySigma)); // Peak flux + if (peak < noise) { + return; + } + +#if SOURCES_GAUSSIAN + // Gaussian + float radius = sqrtf(2.0 * (PS_SQR(xSigma) + PS_SQR(ySigma)) * logf(peak / noise)); // Radius for source +#else + // Basic Waussian + float radius; + if (peak / noise > 10.0) { + // Need to treat the bright limit differently because the wings are much more noticeable + radius = sqrtf(2.0 * (PS_SQR(xSigma) + PS_SQR(ySigma)) * peak / noise); + } else { + radius = sqrtf(2.0) * hypot(1.0 / xSigma, 1.0 / ySigma) * sqrtf(2.0 * logf(peak / noise)); + } +#endif + + // Bounds of source + int xMin = PS_MAX(floorf(x - radius), 0); + int xMax = PS_MIN(ceilf(x + radius), NUMCOLS - 1); + int yMin = PS_MAX(floorf(y - radius), 0); + int yMax = PS_MIN(ceilf(y + radius), NUMROWS - 1); + + for (int j = yMin; j <= yMax; j++) { + float dy = j - y; + for (int i = xMin; i <= xMax; i++) { + float dx = i - x; +#if SOURCES_GAUSSIAN + // Gaussian + float value = peak * expf(- PS_SQR(dx) / 2.0 / PS_SQR(xSigma) - + PS_SQR(dy) / 2.0 / PS_SQR(ySigma)); +#else + // Basic Waussian + float z = 0.5 * (PS_SQR(dx / xSigma) + PS_SQR(dy / ySigma)); + float value = peak / (1.0 + z + PS_SQR(z)); +#endif + + image->data.F32[j][i] += value; + } + } + + return; +} + + +int main(int argc, char *argv[]) +{ + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + + // Images to generate + psImage *image1 = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32); + psImage *image2 = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32); + + psImageInit(image1, 0.0); + psImageInit(image2, 0.0); + + float magNorm = SOURCES_DENSITY * NUMCOLS * NUMROWS * PS_SQR(IMAGE_SCALE / 3600.0) / + powf(10.0, (SOURCES_LUM * SOURCES_MAG)); // Normalisation for magnitudes + + float faintMag = -2.5 * log10(PS_MIN(NOISE_1 * sqrtf(X_AXIS_1 * Y_AXIS_1), + NOISE_2 * sqrtf(X_AXIS_2 * Y_AXIS_2)) * + NOISE_FRAC) + SOURCES_ZP; // Faintest magnitude + printf("Faintest mag: %f\n", faintMag); + + long numStars = magNorm * powf(10.0, SOURCES_LUM * faintMag); // Number of stars + printf("%ld stars\n", numStars); + + // Gaussian widths, in sigma + float xSigma1 = X_AXIS_1 / SIGMA_FWHM; + float ySigma1 = Y_AXIS_1 / SIGMA_FWHM; + float xSigma2 = X_AXIS_2 / SIGMA_FWHM; + float ySigma2 = Y_AXIS_2 / SIGMA_FWHM; + + for (int i = 0; i < numStars; i++) { + float x = psRandomUniform(rng) * NUMCOLS; // x coordinate + float y = psRandomUniform(rng) * NUMROWS; // y coordinate + float mag = log10((i + 1.0) / magNorm) / SOURCES_LUM; // Magnitude of source + float flux = powf(10.0, -0.4 * (mag - SOURCES_ZP)); // Total flux of source + + addstar(x, y, flux, xSigma1, ySigma1, NOISE_1 * NOISE_FRAC, image1); + addstar(x, y, flux, xSigma2, ySigma2, NOISE_2 * NOISE_FRAC, image2); + } + +#if 1 + psImage *var1 = (psImage*)psBinaryOp(NULL, image1, "+", psScalarAlloc(PS_SQR(NOISE_1), PS_TYPE_F32)); + psImage *var2 = (psImage*)psBinaryOp(NULL, image2, "+", psScalarAlloc(PS_SQR(NOISE_2), PS_TYPE_F32)); +#else + psImage *var1 = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32); + psImage *var2 = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32); + psImageInit(var1, PS_SQR(NOISE_1)); + psImageInit(var2, PS_SQR(NOISE_2)); +#endif + + for (int y = 0; y < NUMROWS; y++) { + for (int x = 0; x < NUMCOLS; x++) { + image1->data.F32[y][x] += psRandomGaussian(rng) * sqrtf(var1->data.F32[y][x]); + image2->data.F32[y][x] += psRandomGaussian(rng) * sqrtf(var2->data.F32[y][x]); + } + } +#if 1 + var1 = (psImage*)psBinaryOp(var1, image1, "+", psScalarAlloc(PS_SQR(NOISE_1), PS_TYPE_F32)); + var2 = (psImage*)psBinaryOp(var2, image2, "+", psScalarAlloc(PS_SQR(NOISE_2), PS_TYPE_F32)); +#else + psImageInit(var1, PS_SQR(NOISE_1)); + psImageInit(var2, PS_SQR(NOISE_2)); +#endif + + { + psFits *fits = psFitsOpen("image1.fits", "w"); + psFitsWriteImage(fits, NULL, image1, 0, NULL); + psFitsClose(fits); + } + { + psFits *fits = psFitsOpen("image2.fits", "w"); + psFitsWriteImage(fits, NULL, image2, 0, NULL); + psFitsClose(fits); + } + + { + psFits *fits = psFitsOpen("variance1.fits", "w"); + psFitsWriteImage(fits, NULL, var1, 0, NULL); + psFitsClose(fits); + } + { + psFits *fits = psFitsOpen("variance2.fits", "w"); + psFitsWriteImage(fits, NULL, var2, 0, NULL); + psFitsClose(fits); + } + + psFree(image1); + psFree(image2); + + psFree(var1); + psFree(var2); + + psFree(rng); + + return PS_EXIT_SUCCESS; +}