DM-11957: Cannot round-trip >7th degree Chebyshev photometry models

include/lsst/afw/formatters/Utils.h

@@ -32,8 +32,12 @@
 #ifndef LSST_AFW_FORMATTERS_UTILS_H
 #define LSST_AFW_FORMATTERS_UTILS_H
 
+#include <cstdint>
 #include <set>
 #include <string>
+#include <vector>
+
+#include "ndarray.h"
 
 #include "lsst/base.h"
 #include "lsst/daf/base.h"
@@ -167,6 +171,16 @@ std::string formatFitsProperties(lsst::daf::base::PropertyList const& prop,
                                  std::set<std::string> const& excludeNames = {});
 int countFitsHeaderCards(lsst::daf::base::PropertySet const& prop);
 
+/**
+ * Encode a std::string as a vector of uint8
+ */
+ndarray::Array<std::uint8_t, 1, 1> stringToBytes(std::string const &str);
+
+/**
+ * Decode a std::string from a vector of uint8 returned by stringToBytes
+ */
+std::string bytesToString(ndarray::Array<std::uint8_t const, 1, 1> const& bytes);
+
 }
 }
 }  // namespace lsst::afw::formatters

include/lsst/afw/image/ImageBase.h

@@ -385,15 +385,15 @@ class ImageBase : public lsst::daf::base::Persistable, public lsst::daf::base::C
      * @throws lsst::pex::exceptions::RuntimeError Argument `contiguous` is false, or the pixels are not in
      * fact contiguous
      */
-    fast_iterator begin(bool) const;
+    fast_iterator begin(bool contiguous) const;
     /** Return a fast STL compliant iterator to the end of the %image which must be contiguous
      *
      * @param contiguous Pixels are contiguous (must be true)
      *
      * @throws lsst::pex::exceptions::RuntimeError Argument `contiguous` is false, or the pixels are not in
      * fact contiguous
      */
-    fast_iterator end(bool) const;
+    fast_iterator end(bool contiguous) const;
 
     /** Return an `x_iterator` to the start of the `y`'th row
      *

include/lsst/afw/image/MaskedImage.h

@@ -1183,7 +1183,7 @@ class MaskedImage : public lsst::daf::base::Persistable, public lsst::daf::base:
      * @throws lsst::pex::exceptions::RuntimeError Argument `contiguous` is false, or the pixels are not in
      * fact contiguous
      */
-    fast_iterator begin(bool) const;
+    fast_iterator begin(bool contiguous) const;
     /** Return a fast `iterator` to the end of the %image, which must be contiguous
      * Note that the order in which pixels are visited is undefined.
      *
@@ -1192,7 +1192,7 @@ class MaskedImage : public lsst::daf::base::Persistable, public lsst::daf::base:
      * @throws lsst::pex::exceptions::RuntimeError Argument `contiguous` is false, or the pixels are not in
      * fact contiguous
      */
-    fast_iterator end(bool) const;
+    fast_iterator end(bool contiguous) const;
 
     /// Return an `x_iterator` to the start of the %image
     x_iterator row_begin(int y) const;

include/lsst/afw/table/types.h

@@ -25,20 +25,20 @@
  */
 
 // Scalar types: those that can serve as elements for other types, and use the default FieldBase template.
-#define AFW_TABLE_SCALAR_FIELD_TYPE_N 6
-#define AFW_TABLE_SCALAR_FIELD_TYPES RecordId, std::uint16_t, std::int32_t, float, double, Angle
+#define AFW_TABLE_SCALAR_FIELD_TYPE_N 7
+#define AFW_TABLE_SCALAR_FIELD_TYPES RecordId, std::uint16_t, std::int32_t, float, double, Angle, std::uint8_t
 #define AFW_TABLE_SCALAR_FIELD_TYPE_TUPLE BOOST_PP_LPAREN() AFW_TABLE_SCALAR_FIELD_TYPES BOOST_PP_RPAREN()
 
 // Arrays types: the types we allow for Array fields.
-#define AFW_TABLE_ARRAY_FIELD_TYPE_N 4
-#define AFW_TABLE_ARRAY_FIELD_TYPES std::uint16_t, int, float, double
+#define AFW_TABLE_ARRAY_FIELD_TYPE_N 5
+#define AFW_TABLE_ARRAY_FIELD_TYPES std::uint16_t, int, float, double, std::uint8_t
 #define AFW_TABLE_ARRAY_FIELD_TYPE_TUPLE BOOST_PP_LPAREN() AFW_TABLE_ARRAY_FIELD_TYPES BOOST_PP_RPAREN()
 
 // Field types: all the types we allow for fields.
-#define AFW_TABLE_FIELD_TYPE_N 12
+#define AFW_TABLE_FIELD_TYPE_N 14
 #define AFW_TABLE_FIELD_TYPES                                                                        \
     AFW_TABLE_SCALAR_FIELD_TYPES, Flag, std::string, Array<std::uint16_t>, Array<int>, Array<float>, \
-            Array<double>
+            Array<double>, Array<std::uint8_t>
 
 #define AFW_TABLE_FIELD_TYPE_TUPLE BOOST_PP_LPAREN() AFW_TABLE_FIELD_TYPES BOOST_PP_RPAREN()
 

python/lsst/afw/formatters/SConscript

@@ -1,3 +1,8 @@
 ## -*- python -*-
-#from lsst.sconsUtils import scripts
-#scripts.BasicSConscript.pybind11([])
+from lsst.sconsUtils import scripts
+scripts.BasicSConscript.pybind11(
+    [
+        "utils",
+    ],
+    addUnderscore=False,
+)

python/lsst/afw/formatters/__init__.py

@@ -22,4 +22,4 @@
 """lsst.afw.formatters
 """
 from __future__ import absolute_import
-from .detectionLib import *
+from .utils import *

python/lsst/afw/formatters/utils.cc

@@ -1,6 +1,6 @@
 /*
  * LSST Data Management System
- * Copyright 2008-2016  AURA/LSST.
+ * Copyright 2017  AURA/LSST.
  *
  * This product includes software developed by the
  * LSST Project (http://www.lsst.org/).
@@ -21,25 +21,38 @@
  */
 
 #include <pybind11/pybind11.h>
-//#include <pybind11/operators.h>
-//#include <pybind11/stl.h>
+#include <pybind11/stl.h>
 
-namespace py = pybind11;
-
-using namespace lsst::afw::formatters;
+#include "numpy/arrayobject.h"
+#include "ndarray/pybind11.h"
 
-PYBIND11_PLUGIN(_utils) {
-    py::module mod("_utils", "Python wrapper for afw _utils library");
+#include "lsst/afw/formatters/Utils.h"
 
-    /* Module level */
+namespace py = pybind11;
+using namespace py::literals;
 
-    /* Member types and enums */
+namespace py = pybind11;
+namespace lsst {
+namespace afw {
+namespace formatters {
+namespace {
 
-    /* Constructors */
+PYBIND11_PLUGIN(utils) {
+    py::module mod("utils");
 
-    /* Operators */
+    // Need to import numpy for ndarray and eigen conversions
+    if (_import_array() < 0) {
+        PyErr_SetString(PyExc_ImportError, "numpy.core.multiarray failed to import");
+        return nullptr;
+    }
 
-    /* Members */
+    mod.def("stringToBytes", stringToBytes);
+    mod.def("bytesToString", bytesToString);
 
     return mod.ptr();
-}
+}
+
+}  // namespace
+}  // namespace formatters
+}  // namespace afw
+}  // namespace lsst

python/lsst/afw/table/base/base.cc

@@ -103,11 +103,13 @@ PyBaseRecord declareBaseRecord(py::module &mod) {
     declareBaseRecordOverloads<double>(cls, "D");
     declareBaseRecordOverloads<float>(cls, "F");
     declareBaseRecordOverloads<lsst::afw::table::Flag>(cls, "Flag");
+    declareBaseRecordOverloads<std::uint8_t>(cls, "B");
     declareBaseRecordOverloads<std::uint16_t>(cls, "U");
     declareBaseRecordOverloads<std::int32_t>(cls, "I");
     declareBaseRecordOverloads<std::int64_t>(cls, "L");
     declareBaseRecordOverloads<std::string>(cls, "String");
     declareBaseRecordOverloads<lsst::afw::geom::Angle>(cls, "Angle");
+    declareBaseRecordArrayOverloads<std::uint8_t>(cls, "ArrayB");
     declareBaseRecordArrayOverloads<std::uint16_t>(cls, "ArrayU");
     declareBaseRecordArrayOverloads<int>(cls, "ArrayI");
     declareBaseRecordArrayOverloads<float>(cls, "ArrayF");

python/lsst/afw/table/baseColumnView/baseColumnView.cc

@@ -85,6 +85,7 @@ static void declareBaseColumnView(py::module &mod) {
     // _getBits supports a Python version of getBits that accepts None and field names as keys
     cls.def("_getBits", &BaseColumnView::getBits);
     cls.def("getAllBits", &BaseColumnView::getAllBits);
+    declareBaseColumnViewOverloads<std::uint8_t>(cls);
     declareBaseColumnViewOverloads<std::uint16_t>(cls);
     declareBaseColumnViewOverloads<std::int32_t>(cls);
     declareBaseColumnViewOverloads<std::int64_t>(cls);
@@ -94,6 +95,7 @@ static void declareBaseColumnView(py::module &mod) {
     // std::string columns are not supported, because numpy string arrays
     // do not have the same memory model as ours.
     declareBaseColumnViewOverloads<lsst::afw::geom::Angle>(cls);
+    declareBaseColumnViewArrayOverloads<std::uint8_t>(cls);
     declareBaseColumnViewArrayOverloads<std::uint16_t>(cls);
     declareBaseColumnViewArrayOverloads<int>(cls);
     declareBaseColumnViewArrayOverloads<float>(cls);

python/lsst/afw/table/schema/schema.cc

@@ -407,13 +407,15 @@ PYBIND11_PLUGIN(schema) {
     mod.attr("_Key") = py::dict();
     mod.attr("_SchemaItem") = py::dict();
 
+    declareSchemaType<std::uint8_t>(mod);
     declareSchemaType<std::uint16_t>(mod);
     declareSchemaType<std::int32_t>(mod);
     declareSchemaType<std::int64_t>(mod);
     declareSchemaType<float>(mod);
     declareSchemaType<double>(mod);
     declareSchemaType<std::string>(mod);
     declareSchemaType<geom::Angle>(mod);
+    declareSchemaType<Array<std::uint8_t>>(mod);
     declareSchemaType<Array<std::uint16_t>>(mod);
     declareSchemaType<Array<int>>(mod);
     declareSchemaType<Array<float>>(mod);

python/lsst/afw/table/schema/schemaContinued.py

@@ -42,6 +42,7 @@
 # Key/Field/SchemaItem types.
 _dtypes = {
     "String": str,
+    "B": np.uint8,
     "U": np.uint16,
     "I": np.int32,
     "L": np.int64,

python/lsst/afw/table/schemaMapper/schemaMapper.cc

@@ -60,6 +60,7 @@ PYBIND11_PLUGIN(schemaMapper) {
     cls.def("addMinimalSchema", &SchemaMapper::addMinimalSchema, "minimal"_a, "doMap"_a = true);
     cls.def_static("removeMinimalSchema", &SchemaMapper::removeMinimalSchema);
     cls.def_static("join", &SchemaMapper::join, "inputs"_a, "prefixes"_a = std::vector<std::string>());
+    declareSchemaMapperOverloads<std::uint8_t>(cls, "B");
     declareSchemaMapperOverloads<std::uint16_t>(cls, "U");
     declareSchemaMapperOverloads<std::int32_t>(cls, "I");
     declareSchemaMapperOverloads<std::int64_t>(cls, "L");
@@ -68,6 +69,7 @@ PYBIND11_PLUGIN(schemaMapper) {
     declareSchemaMapperOverloads<std::string>(cls, "String");
     declareSchemaMapperOverloads<lsst::afw::geom::Angle>(cls, "Angle");
     declareSchemaMapperOverloads<lsst::afw::table::Flag>(cls, "Flag");
+    declareSchemaMapperOverloads<lsst::afw::table::Array<std::uint8_t>>(cls, "ArrayB");
     declareSchemaMapperOverloads<lsst::afw::table::Array<std::uint16_t>>(cls, "ArrayU");
     declareSchemaMapperOverloads<lsst::afw::table::Array<int>>(cls, "ArrayI");
     declareSchemaMapperOverloads<lsst::afw::table::Array<float>>(cls, "ArrayF");

src/fits.cc

@@ -87,6 +87,10 @@ struct FitsType<char> {
     static int const CONSTANT = TSTRING;
 };
 template <>
+struct FitsType<signed char> {
+    static int const CONSTANT = TSBYTE;
+};
+template <>
 struct FitsType<unsigned char> {
     static int const CONSTANT = TBYTE;
 };
@@ -690,6 +694,15 @@ void writeKeyFromProperty(Fits &fits, daf::base::PropertySet const &metadata, st
         } else {
             writeKeyImpl(fits, key.c_str(), metadata.get<bool>(key), comment);
         }
+    } else if (valueType == typeid(std::uint8_t)) {
+        if (metadata.isArray(key)) {
+            std::vector<std::uint8_t> tmp = metadata.getArray<std::uint8_t>(key);
+            for (std::size_t i = 0; i != tmp.size(); ++i) {
+                writeKeyImpl(fits, key.c_str(), tmp[i], comment);
+            }
+        } else {
+            writeKeyImpl(fits, key.c_str(), metadata.get<std::uint8_t>(key), comment);
+        }
     } else if (valueType == typeid(int)) {
         if (metadata.isArray(key)) {
             std::vector<int> tmp = metadata.getArray<int>(key);

src/fitsCompression.cc

@@ -298,6 +298,7 @@ struct Bzero<T, typename std::enable_if<std::numeric_limits<T>::is_integer &&
 };
 
 
+#ifndef DOXYGEN  // suppress a bogus Doxygen complaint about an documented symbol
 template <typename T, int N>
 ImageScale ImageScalingOptions::determine(
     ndarray::Array<T const, N, N> const& image,
@@ -322,6 +323,7 @@ ImageScale ImageScalingOptions::determine(
         std::abort();  // should never get here
     }
 }
+#endif
 
 namespace {
 

src/formatters/Utils.cc

@@ -286,6 +286,24 @@ std::string formatFitsProperties(daf::base::PropertyList const& prop,
 int countFitsHeaderCards(lsst::daf::base::PropertySet const& prop) {
     return prop.paramNames(false).size();
 }
+
+ndarray::Array<std::uint8_t, 1, 1> stringToBytes(std::string const& str) {
+    auto nbytes = str.size() * sizeof(char) / sizeof(std::uint8_t);
+    std::uint8_t const* byteCArr = reinterpret_cast<std::uint8_t const*>(str.data());
+    auto shape = ndarray::makeVector(nbytes);
+    auto strides = ndarray::makeVector(1);
+    // Make an Array that shares memory with `str` (and does not free that memory when destroyed),
+    // then return a copy; this is simpler than manually copying the data into a newly allocated array
+    ndarray::Array<std::uint8_t const, 1, 1> localArray = ndarray::external(byteCArr, shape, strides);
+    return ndarray::copy(localArray);
+}
+
+std::string bytesToString(ndarray::Array<std::uint8_t const, 1, 1> const& bytes) {
+    auto nchars = bytes.size() * sizeof(std::uint8_t) / sizeof(char);
+    char const* charCArr = reinterpret_cast<char const*>(bytes.getData());
+    return std::string(charCArr, nchars);
+}
+
 }
 }
 }  // namespace lsst::afw::formatters

src/math/TransformBoundedField.cc

@@ -21,11 +21,13 @@
  * see <http://www.lsstcorp.org/LegalNotices/>.
  */
 
+#include <cstdint>
 #include <memory>
 #include <string>
 
 #include "ndarray/eigen.h"
 #include "astshim.h"
+#include "lsst/afw/formatters/Utils.h"
 #include "lsst/afw/math/TransformBoundedField.h"
 #include "lsst/afw/table/io/InputArchive.h"
 #include "lsst/afw/table/io/OutputArchive.h"
@@ -81,16 +83,17 @@ struct PersistenceHelper {
     table::Schema schema;
     table::PointKey<int> bboxMin;
     table::PointKey<int> bboxMax;
-    table::Key<std::string> frameSet;
+    // store the FrameSet as string encoded as a variable-length vector of bytes
+    table::Key<table::Array<std::uint8_t>> frameSet;
 
     PersistenceHelper()
             : schema(),
               bboxMin(table::PointKey<int>::addFields(schema, "bbox_min", "lower-left corner of bounding box",
                                                       "pixel")),
               bboxMax(table::PointKey<int>::addFields(schema, "bbox_max",
                                                       "upper-right corner of bounding box", "pixel")),
-              frameSet(schema.addField<std::string>("frameSet", "FrameSet contained in the Transform", 0,
-                                                    false)) {}
+              frameSet(schema.addField<table::Array<std::uint8_t>>(
+                      "frameSet", "FrameSet contained in the Transform", "", 0)) {}
 
     PersistenceHelper(table::Schema const& s)
             : schema(s), bboxMin(s["bbox_min"]), bboxMax(s["bbox_max"]), frameSet(s["frameSet"]) {}
@@ -109,7 +112,7 @@ class TransformBoundedFieldFactory : public table::io::PersistableFactory {
         PersistenceHelper const keys(record.getSchema());
         // NOTE: needed invert=false in case min=-1, max=0 (empty bbox). See RFC-324 and DM-10200
         geom::Box2I bbox(record.get(keys.bboxMin), record.get(keys.bboxMax), false);
-        auto frameSetStr = record.get(keys.frameSet);
+        auto frameSetStr = formatters::bytesToString(record.get(keys.frameSet));
         auto transform =
                 geom::Transform<geom::Point2Endpoint, geom::GenericEndpoint>::readString(frameSetStr);
         return std::make_shared<TransformBoundedField>(bbox, transform);
@@ -134,7 +137,7 @@ void TransformBoundedField::write(OutputArchiveHandle& handle) const {
     std::shared_ptr<table::BaseRecord> record = catalog.addNew();
     record->set(keys.bboxMin, getBBox().getMin());
     record->set(keys.bboxMax, getBBox().getMax());
-    record->set(keys.frameSet, getTransform().writeString());
+    record->set(keys.frameSet, formatters::stringToBytes(getTransform().writeString()));
     handle.saveCatalog(catalog);
 }
 

src/table/FieldBase.cc

@@ -19,6 +19,10 @@ namespace {
 template <typename T>
 struct TypeTraits;
 
+template <>
+struct TypeTraits<std::uint8_t> {
+    static char const *getName() { return "B"; }
+};
 template <>
 struct TypeTraits<std::uint16_t> {
     static char const *getName() { return "U"; }

src/table/io/FitsSchemaInputMapper.cc

@@ -591,6 +591,18 @@ std::unique_ptr<FitsColumnReader> makeColumnReader(Schema &schema, FitsSchemaIte
     }
     // switch code over FITS codes that correspond to different element types
     switch (code) {
+        case 'B':  // 8-bit unsigned integers -- can only be scalars or Arrays
+            if (size == 1) {
+                if (item.tccls == "Array") {
+                    return StandardReader<Array<std::uint8_t>>::make(schema, item, size);
+                }
+                return StandardReader<std::uint8_t>::make(schema, item);
+            }
+            if (size == 0) {
+                return VariableLengthArrayReader<std::uint8_t>::make(schema, item);
+            }
+            return StandardReader<Array<std::uint8_t>>::make(schema, item, size);
+
         case 'I':  // 16-bit integers - can only be scalars or Arrays (we assume they're unsigned, since
                    // that's all we ever write, and CFITSIO will complain later if they aren't)
             if (size == 1) {