Add support for subblocks to handle +2GiB compressed images

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.14)
 
-project(libXISF VERSION 0.2.7 LANGUAGES CXX C
+project(libXISF VERSION 0.2.11 LANGUAGES CXX C
     HOMEPAGE_URL https://gitea.nouspiro.space/nou/libXISF
     DESCRIPTION "LibXISF is C++ library that can read and write XISF files produced by PixInsight.")
 
@@ -92,10 +92,14 @@ else(USE_BUNDLED_LIBS)
     list(APPEND PC_LIBS_REQUIRE lz4 pugixml zlib)
 endif(USE_BUNDLED_LIBS)
 
-target_link_libraries(XISF PUBLIC zstd)
-
 set_target_properties(XISF PROPERTIES VERSION ${PROJECT_VERSION} SOVERSION ${PROJECT_VERSION_MAJOR})
 
+pkg_check_modules(ZSTD libzstd IMPORTED_TARGET)
+if(ZSTD_FOUND)
+    target_compile_definitions(XISF PRIVATE HAVE_ZSTD)
+    target_link_libraries(XISF PUBLIC PkgConfig::ZSTD)
+endif(ZSTD_FOUND)
+
 if(BUILD_SHARED_LIBS)
     target_compile_definitions(XISF PRIVATE LIBXISF_LIBRARY)
 else(BUILD_SHARED_LIBS)

README.md

@@ -9,7 +9,7 @@ To compile simply run these commands
 ```
 cmake -B build -S .
 cmake --build build --parallel
-cmake --install .
+cmake --install build
 ```
 
 By default it use bundled libraries. If you wish to use external libraries you will may add

libxisf.cpp

@@ -27,6 +27,9 @@
 #include <lz4hc.h>
 #include <pugixml.hpp>
 #include <zlib.h>
+#ifdef HAVE_ZSTD
+#include <zstd.h>
+#endif
 #include "streambuffer.h"
 
 namespace LibXISF
@@ -46,6 +49,7 @@ static std::unordered_map<Image::ColorSpace, String> colorSpaceToString;
 static DataBlock::CompressionCodec compressionCodecOverride = DataBlock::None;
 static bool byteShuffleOverride = false;
 static int compressionLevelOverride = -1;
+const size_t GiB = 1073741824;
 
 static const std::unordered_map<String, std::pair<String, Variant::Type>> fitsNameToPropertyIdTypeConvert = {
     {"OBSERVER", {"Observer:Name", Variant::Type::String}},
@@ -83,12 +87,12 @@ static void byteShuffle(ByteArray &data, int itemSize)
     {
         ByteArray &input = data;
         ByteArray output(input.size());
-        int num = input.size() / itemSize;
+        size_t num = input.size() / itemSize;
         char *s = output.data();
         for(int i=0; i<itemSize; i++)
         {
             const char *u = input.constData() + i;
-            for(int o=0; o<num; o++, s++, u += itemSize)
+            for(size_t o=0; o<num; o++, s++, u += itemSize)
                 *s = *u;
         }
         memcpy(s, input.constData() + num * itemSize, input.size() % itemSize);
@@ -102,12 +106,12 @@ static void byteUnshuffle(ByteArray &data, int itemSize)
     {
         ByteArray &input = data;
         ByteArray output(input.size());
-        int num = input.size() / itemSize;
+        size_t num = input.size() / itemSize;
         const char *s = input.constData();
         for(int i=0; i<itemSize; i++)
         {
             char *u = output.data() + i;
-            for(int o=0; o<num; o++, s++, u += itemSize)
+            for(size_t o=0; o<num; o++, s++, u += itemSize)
                 *u = *s;
         }
         memcpy(output.data() + num * itemSize, s, input.size() % itemSize);
@@ -124,6 +128,9 @@ void DataBlock::decompress(const ByteArray &input, const String &encoding)
     else if(encoding == "base16")
         tmp.decodeHex();
 
+    if(subblocks.size() == 0)
+        subblocks.push_back({tmp.size(), uncompressedSize});
+
     switch(codec)
     {
     case None:
@@ -132,17 +139,54 @@ void DataBlock::decompress(const ByteArray &input, const String &encoding)
     case Zlib:
     {
         data.resize(uncompressedSize);
-        uLongf size = uncompressedSize;
+        const char *srcPtr = tmp.constData();
-        ::uncompress((Bytef*)data.data(), &size, (Bytef*)tmp.data(), tmp.size());
+        char *dstPtr = data.data();
+        for(auto &block : subblocks)
+        {
+            uLongf size = block.second;
+            ::uncompress((Bytef*)dstPtr, &size, (const Bytef*)srcPtr, block.first);
+            srcPtr += block.first;
+            dstPtr += block.second;
+        }
         break;
     }
     case LZ4:
     case LZ4HC:
+    {
         data.resize(uncompressedSize);
-        if(LZ4_decompress_safe(tmp.constData(), data.data(), tmp.size(), data.size()) < 0)
+        const char *srcPtr = tmp.constData();
-            throw Error("LZ4 decompression failed");
+        char *dstPtr = data.data();
+        for(auto &block : subblocks)
+        {
+            if(LZ4_decompress_safe(srcPtr, dstPtr, block.first, block.second) < 0)
+                throw Error("LZ4 decompression failed");
+            srcPtr += block.first;
+            dstPtr += block.second;
+        }
         break;
     }
+    case ZSTD:
+#ifdef HAVE_ZSTD
+    {
+        data.resize(uncompressedSize);
+        const char *srcPtr = tmp.constData();
+        char *dstPtr = data.data();
+        for(auto &block : subblocks)
+        {
+            if(ZSTD_isError(ZSTD_decompress(dstPtr, block.second, srcPtr, block.first)))
+                throw Error("ZSTD decompression failed");
+            srcPtr += block.first;
+            dstPtr += block.second;
+        }
+
+#else
+        throw Error("ZSTD support not compiled");
+#endif
+        break;
+    }
+    }
+
+    subblocks.clear();
 
     byteUnshuffle(data, byteShuffling);
     attachmentPos = 0;
@@ -169,30 +213,77 @@ void DataBlock::compress(int sampleFormatSize)
         break;
     case Zlib:
     {
-        data.resize(compressBound(uncompressedSize));
+        int64_t size = tmp.size();
-        uLongf compressedSize = data.size();
+        int64_t compSize = 0;
-        if(::compress2((Bytef*)data.data(), &compressedSize, (Bytef*)tmp.data(), tmp.size(), compressLevel) != Z_OK)
+        int64_t inPtr = 0;
-            throw Error("Zlib compression failed");
+        while(inPtr < size)
-        data.resize(compressedSize);
+        {
+            int64_t inSize = UINT32_MAX < size - inPtr ? UINT32_MAX : size - inPtr;
+            data.resize(compSize + compressBound(inSize));
+            uLongf outSize = data.size() - compSize;
+            if(::compress2((Bytef*)data.data() + compSize, &outSize, (const Bytef*)tmp.constData() + inPtr, inSize, compressLevel) != Z_OK)
+                throw Error("Zlib compression failed");
+
+            compSize += outSize;
+            inPtr += inSize;
+            subblocks.push_back({outSize, inSize});
+        }
+        data.resize(compSize);
         break;
     }
     case LZ4:
     case LZ4HC:
     {
-        int compSize = 0;
+        int64_t size = tmp.size();
-        data.resize(LZ4_compressBound(tmp.size()));
+        int64_t compSize = 0;
-        if(codec == LZ4)
+        int64_t inPtr = 0;
-            compSize = LZ4_compress_default(tmp.constData(), data.data(), tmp.size(), data.size());
+        while(inPtr < size)
-        else
+        {
-            compSize = LZ4_compress_HC(tmp.constData(), data.data(), tmp.size(), data.size(), compressLevel < 0 ? LZ4HC_CLEVEL_DEFAULT : compressLevel);
+            int64_t inSize = LZ4_MAX_INPUT_SIZE < size - inPtr ? LZ4_MAX_INPUT_SIZE : size - inPtr;
+            data.resize(compSize + LZ4_compressBound(inSize));
+            int outSize = 0;
 
-        if(compSize <= 0)
+            if(codec == LZ4)
-            throw Error("LZ4 compression failed");
+                outSize = LZ4_compress_default(tmp.constData() + inPtr, data.data() + compSize, inSize, data.size() - compSize);
+            else
+                outSize = LZ4_compress_HC(tmp.constData() + inPtr, data.data() + compSize, inSize, data.size() - compSize, compressLevel < 0 ? LZ4HC_CLEVEL_DEFAULT : compressLevel);
 
+            if(outSize <= 0)
+                throw Error("LZ4 compression failed");
+
+            compSize += outSize;
+            inPtr += inSize;
+            subblocks.push_back({outSize, inSize});
+        }
         data.resize(compSize);
         break;
     }
+    case ZSTD:
+    {
+#ifdef HAVE_ZSTD
+        size_t compSize = 0;
+        data.resize(ZSTD_compressBound(uncompressedSize));
+        compSize = ZSTD_compress(data.data(), data.size(), tmp.data(), tmp.size(), compressLevel < 0 ? ZSTD_CLEVEL_DEFAULT : compressLevel);
+        if(ZSTD_isError(compSize))
+            throw Error("ZSTD compression failed");
+
+        data.resize(compSize);
+#else
+        throw Error("ZSTD support not compiled");
+#endif
+        break;
     }
+    }
+}
+
+bool DataBlock::CompressionCodecSupported(CompressionCodec codec)
+{
+    (void)codec;
+#ifndef HAVE_ZSTD
+    if(codec == ZSTD)
+        return false;
+#endif
+    return true;
 }
 
 Property::Property(const String &_id, const char *_value) :
@@ -352,6 +443,10 @@ bool Image::addFITSKeywordAsProperty(const String &name, const String &value)
     {
         auto &c = fitsNameToPropertyIdTypeConvert.at(name);
         Property prop(c.first, variantFromString(c.second, value));
+
+        if(name == "APTDIA" || name == "FOCALLEN")
+            prop.value.value<LibXISF::Float32>() /= 1000.0f;
+
         updateProperty(prop);
         return true;
     }
@@ -663,10 +758,14 @@ void XISFReaderPrivate::parseCompression(const pugi::xml_node &node, DataBlock &
             dataBlock.codec = DataBlock::LZ4HC;
         else if(compression[0].find("lz4") == 0)
             dataBlock.codec = DataBlock::LZ4;
+#ifdef HAVE_ZSTD
+        else if(compression[0].find("zstd") == 0)
+            dataBlock.codec = DataBlock::ZSTD;
+#endif
         else
             throw Error("Unknown compression codec");
 
-        dataBlock.uncompressedSize = std::stoul(compression[1]);
+        dataBlock.uncompressedSize = std::stoull(compression[1]);
 
         if(compression[0].find("+sh") != std::string::npos)
         {
@@ -675,6 +774,18 @@ void XISFReaderPrivate::parseCompression(const pugi::xml_node &node, DataBlock &
             else
                 throw Error("Missing byte shuffling size");
         }
+
+        if(node.attribute("subblocks"))
+        {
+            std::vector<std::string> subblocks = splitString(node.attribute("subblocks").as_string(), ':');
+            for(auto &block : subblocks)
+            {
+                size_t pos = 0;
+                size_t comp = std::stoull(block, &pos);
+                size_t deco = std::stoull(block.substr(pos+1));
+                dataBlock.subblocks.push_back({comp, deco});
+            }
+        }
     }
 }
 
@@ -696,8 +807,8 @@ DataBlock XISFReaderPrivate::parseDataBlock(const pugi::xml_node &node)
     }
     else if(location.size() >= 3 && location[0] == "attachment")
     {
-        dataBlock.attachmentPos = std::stoul(location[1]);
+        dataBlock.attachmentPos = std::stoull(location[1]);
-        dataBlock.attachmentSize = std::stoul(location[2]);
+        dataBlock.attachmentSize = std::stoull(location[2]);
     }
     else
     {
@@ -773,9 +884,9 @@ Image XISFReaderPrivate::parseImage(const pugi::xml_node &node)
 
     std::vector<std::string> geometry = splitString(node.attribute("geometry").as_string(), ':');
     if(geometry.size() != 3)throw Error("We support only 2D images");
-    image._width = std::stoul(geometry[0]);
+    image._width = std::stoull(geometry[0]);
-    image._height = std::stoul(geometry[1]);
+    image._height = std::stoull(geometry[1]);
-    image._channelCount = std::stoul(geometry[2]);
+    image._channelCount = std::stoull(geometry[2]);
     if(!image._width || !image._height || !image._channelCount)throw Error("Invalid image geometry");
 
     std::vector<std::string> bounds = splitString(node.attribute("bounds").as_string(), ':');
@@ -819,7 +930,15 @@ void XISFReaderPrivate::readAttachment(DataBlock &dataBlock)
 {
     ByteArray data(dataBlock.attachmentSize);
     _io->seekg(dataBlock.attachmentPos);
-    _io->read(data.data(), dataBlock.attachmentSize);
+    size_t size = dataBlock.attachmentSize;
+    char *ptr = data.data();
+    while(size > 0)
+    {
+        size_t s = std::min(size, GiB);
+        _io->read(ptr, s);
+        size -= s;
+        ptr += s;
+    }
     dataBlock.decompress(data);
 }
 
@@ -868,7 +987,15 @@ void XISFWriterPrivate::save(std::ostream &io)
 
     for(auto &image : _images)
     {
-        io.write(image._dataBlock.data.constData(), image._dataBlock.data.size());
+        const char *ptr = image._dataBlock.data.constData();
+        size_t size = image._dataBlock.data.size();
+        while(size > 0)
+        {
+            size_t s = std::min(size, GiB);
+            io.write(ptr, s);
+            ptr += s;
+            size -= s;
+        }
     }
 }
 
@@ -916,7 +1043,7 @@ void XISFWriterPrivate::writeHeader()
         offset += image._dataBlock.data.size();
     }
 
-    uint32_t headerSize = size - sizeof(signature);
+    uint32_t headerSize = header.size() - sizeof(signature);
     header.resize(size, 0);
     header.replace(8, sizeof(uint32_t), (const char*)&headerSize, sizeof(uint32_t));
 
@@ -957,7 +1084,8 @@ void XISFWriterPrivate::writeImageElement(pugi::xml_node &node, const Image &ima
     if(image._iccProfile.size())
     {
         ByteArray base64 = image._iccProfile;
-        base64.decodeBase64();
+        base64.encodeBase64();
+        base64.append('\0');
         pugi::xml_node icc_node = image_node.append_child("ICCProfile");
         icc_node.append_attribute("location").set_value("inline:base64");
         icc_node.append_child(pugi::node_pcdata).set_value(base64.data());
@@ -988,6 +1116,12 @@ void XISFWriterPrivate::writeDataBlockAttributes(pugi::xml_node &image_node, con
         codec = "lz4";
     else if(dataBlock.codec == DataBlock::LZ4HC)
         codec = "lz4hc";
+    else if(dataBlock.codec == DataBlock::ZSTD)
+#ifdef HAVE_ZSTD
+        codec = "zstd";
+#else
+        throw Error("ZSTD support not compiled");
+#endif
 
     if(dataBlock.byteShuffling > 1)
         codec += "+sh";
@@ -1000,6 +1134,19 @@ void XISFWriterPrivate::writeDataBlockAttributes(pugi::xml_node &image_node, con
 
         image_node.append_attribute("compression").set_value(codec.c_str());
     }
+
+    if(!dataBlock.subblocks.empty())
+    {
+        std::string subblocks;
+        for(auto i = dataBlock.subblocks.begin(); i != dataBlock.subblocks.end(); i++)
+        {
+            if(i != dataBlock.subblocks.begin())
+                subblocks += ":";
+
+            subblocks += std::to_string(i->first) + "," + std::to_string(i->second);
+        }
+        image_node.append_attribute("subblocks").set_value(subblocks.c_str());
+    }
 }
 
 void XISFWriterPrivate::writePropertyElement(pugi::xml_node &node, const Property &property)
@@ -1149,6 +1296,10 @@ struct Init
                 compressionCodecOverride = DataBlock::LZ4HC;
             else if(compression.find("lz4") == 0)
                 compressionCodecOverride = DataBlock::LZ4;
+#ifdef HAVE_ZSTD
+            else if(compression.find("zstd") == 0)
+                compressionCodecOverride = DataBlock::ZSTD;
+#endif
 
             if(compression.find("+sh") != std::string::npos)
                 byteShuffleOverride = true;

libxisf.h

@@ -172,18 +172,22 @@ struct LIBXISF_EXPORT DataBlock
         None,
         Zlib,
         LZ4,
-        LZ4HC
+        LZ4HC,
+        ZSTD
     };
     bool embedded = false;
     uint32_t byteShuffling = 0;
     uint64_t attachmentPos = 0;
     uint64_t attachmentSize = 0;
     uint64_t uncompressedSize = 0;
+    std::vector<std::pair<uint64_t, uint64_t>> subblocks;
     CompressionCodec codec = None;
     int compressLevel = -1;
     ByteArray data;
     void decompress(const ByteArray &input, const std::string &encoding = "");
     void compress(int sampleFormatSize);
+    /// ZSTD compression can be disabled at compile time
+    static bool CompressionCodecSupported(CompressionCodec codec);
 };
 
 struct LIBXISF_EXPORT Property

streambuffer.cpp

@@ -89,7 +89,7 @@ std::streamsize StreamBuffer::xsgetn(char_type *s, std::streamsize n)
     std::streamsize len = egptr() - gptr();
     if(len > 0)
     {
-        std::streamsize c = std::min(n, len);
+        std::streamsize c = n < len ? n : len;
         std::memcpy(s, gptr(), c);
         gbump(c);
         ret = c;

test/benchmark.cpp

@@ -78,6 +78,17 @@ void benchmarkType(float avg, float stdDev)
         std::cout << "LZ4HC compression   \tElapsed time: " << timer.elapsed() << " " << "ms\tSpeed: "
                   << size/1024.0/1.024/timer.elapsed() << "MiB/s\tRatio: " << baseSize/xisfImage.size() << std::endl;
     }
+    if(DataBlock::CompressionCodecSupported(DataBlock::ZSTD))
+    {
+        image.setCompression(DataBlock::ZSTD);
+        timer.start();
+        XISFWriter writer;
+        writer.writeImage(image);
+        ByteArray xisfImage;
+        writer.save(xisfImage);
+        std::cout << "ZSTD compression   \tElapsed time: " << timer.elapsed() << " " << "ms\tSpeed: "
+                  << size/1024.0/1.024/timer.elapsed() << "MiB/s\tRatio: " << baseSize/xisfImage.size() << std::endl;
+    }
     image.setByteshuffling(true);
     {
         image.setCompression(DataBlock::Zlib);
@@ -109,6 +120,17 @@ void benchmarkType(float avg, float stdDev)
         std::cout << "LZ4HC compression SH\tElapsed time: " << timer.elapsed() << " " << "ms\tSpeed: "
                   << size/1024.0/1.024/timer.elapsed() << "MiB/s\tRatio: " << baseSize/xisfImage.size() << std::endl;
     }
+    if(DataBlock::CompressionCodecSupported(DataBlock::ZSTD))
+    {
+        image.setCompression(DataBlock::ZSTD);
+        timer.start();
+        XISFWriter writer;
+        writer.writeImage(image);
+        ByteArray xisfImage;
+        writer.save(xisfImage);
+        std::cout << "ZSTD compression SH\tElapsed time: " << timer.elapsed() << " " << "ms\tSpeed: "
+                  << size/1024.0/1.024/timer.elapsed() << "MiB/s\tRatio: " << baseSize/xisfImage.size() << std::endl;
+    }
 }
 
 void benchmark()

utils.cpp

@@ -52,7 +52,7 @@ void sha1(uint8_t *data, size_t len, uint8_t *hash)
     nlen += 64 - nlen % 64;
     tmp.resize(nlen, 0);
 
-    size_t ml = len * 8;
+    uint64_t ml = len * 8;
     tmp[nlen - 1] = ml & 0xff;
     tmp[nlen - 2] = ml >>  8 & 0xff;
     tmp[nlen - 3] = ml >> 16 & 0xff;

variant.cpp

@@ -108,8 +108,8 @@ static std::map<Variant::Type, const char*> idToType = {
     {Variant::Type::UI32Matrix,    "UI32Matrix"},
     {Variant::Type::I64Matrix,     "I64Matrix"},
     {Variant::Type::UI64Matrix,    "UI64Matrix"},
-    {Variant::Type::F32Matrix,     "I8Matrix"},
+    {Variant::Type::F32Matrix,     "F32Matrix"},
-    {Variant::Type::F64Matrix,     "UI8Matrix"},
+    {Variant::Type::F64Matrix,     "F64Matrix"},
 };
 
 template<typename T>
@@ -196,9 +196,12 @@ void deserializeVariant(const pugi::xml_node &node, Variant &variant, const Byte
     std::string type = node.attribute("type").as_string();
     Variant::Type typeId = typeToId[type];
 
-    if(typeId == Variant::Type::String && !node.attribute("location"))
+    if(typeId == Variant::Type::String)
     {
-        variant.setValue(node.text().as_string());
+        if(!node.attribute("location"))
+            variant.setValue(node.text().as_string());
+        else
+            variant.setValue(String(data.constData(), data.size()));
     }
     else if(node.attribute("value"))
     {
@@ -400,6 +403,7 @@ Variant variantFromString(Variant::Type type, const String &str)
 
 Variant::Type Variant::type() const
 {
+    int idx = _value.index();
     return (Variant::Type)_value.index();
 }
 
@@ -430,6 +434,32 @@ String Variant::toString() const
         return ss.str();
     };
 
+    auto matrixString = [](auto matrix) {
+        std::stringstream ss;
+        ss << "{";
+        for(int i=0; i<matrix.rows(); i++)
+        {
+            ss << "{";
+            for(int o=0; o<matrix.cols(); o++)
+            {
+#if __GNUC__ >= 11 || __clang__
+                char str[128] = {0};
+                char *end = str + sizeof(str);
+                std::to_chars(str, end, matrix(i, o));
+                ss << str;
+#else
+                ss << std::to_string(matrix(i, o));
+#endif
+                if(o < matrix.cols() - 1)
+                    ss << ",";
+            }
+            ss << "}";
+            if(i < matrix.rows() - 1)
+                ss << ",";
+        }
+        return ss.str();
+    };
+
     switch(type())
     {
     case Variant::Type::Int8: toChars<Int8>(_value, str, end); break;
@@ -454,6 +484,8 @@ String Variant::toString() const
     case Variant::Type::UI64Vector: string = vectorString(std::get<UI64Vector>(_value)); break;
     case Variant::Type::F32Vector: string = vectorString(std::get<F32Vector>(_value)); break;
     case Variant::Type::F64Vector: string = vectorString(std::get<F64Vector>(_value)); break;
+    case Variant::Type::F32Matrix: string = matrixString(std::get<F32Matrix>(_value)); break;
+    case Variant::Type::F64Matrix: string = matrixString(std::get<F64Matrix>(_value)); break;
     case Variant::Type::String: string = std::get<String>(_value); break;
     case Variant::Type::TimePoint:
     {