czi_file.cpp

#include "czi_file.h"

bool CziFile::is_master_file() const
{
    return ((header.filePart == 0) && (header.fileGuid == header.masterFileGuid));
}

void CziFile::print_segment_header(const SegmentHeader &segmentHeader) const
{
    printf("---------SegmentHeader---------\n");
    printf("%-25s %15s\n", "SID", segmentHeader.sId.c_str());
    printf("%-25s %15li\n", "AllocatedSize", segmentHeader.allocatedSize);
    printf("%-25s %15li\n", "UsedSize", segmentHeader.usedSize);
    printf("-------------------------------\n");
}

const char *CziFile::pixel_type_str(const PixelType px) const
{

    switch (px)
    {
    case PixelType_Gray8:
        return "Gray8";
    case PixelType_Gray16:
        return "Gray16";
    case PixelType_Gray32Float:
        return "Gray32Float";
    case PixelType_Bgr24:
        return "Bgr24";
    case PixelType_Bgr48:
        return "Bgr48";
    case PixelType_Bgr96Float:
        return "Bgr96Float";
    case PixelType_Bgra32:
        return "Bgra32";
    case PixelType_Gray64ComplexFloat:
        return "Gray64ComplexFloat";
    case PixelType_Bgr192ComplexFloat:
        return "Bgr192ComplexFloat";
    case PixelType_Gray32:
        return "Gray32";
    case PixelType_Gray64:
        return "Gray64";
    default:
        always_assert("Bad pixel type." && false);
        break;
    }

    return nullptr;
}

std::string CziFile::dimension_stack_str(const std::vector<DimensionEntryDV1> &dims, bool includeSize) const
{

    std::string result;
    for (const DimensionEntryDV1 &entry : dims)
    {
        result.append(dimension_char(entry.dimension));
        // + std::string(includeSize ? () : ""));
        if (includeSize)
            result.append("(" + std::to_string(entry.size) + ")-");
    }

    if (includeSize)
    {
        result = result.substr(0, result.length() - 1);
    }
    return result;

    // return (includeSize ?: result).c_str();
}

void CziFile::report_verbose() const
{
    report();

    printf("==========SubBlock directory==========\n");
    print_segment_header(subBlockDirectory.header);
    printf("%-25s %15i\n", "EntryCount", subBlockDirectory.entryCount);

    for (size_t entry = 0; entry < subBlockDirectory.entryCount; entry++)
    {
        printf("SubBlockId: %i\tPixelType: %s\tCompression: %s\tDimensionStack: %s\tPyramidType: %s\n",
               (int)entry,
               pixel_type_str(subBlockDirectory.entries[entry].pixelType),
               compression_type_str(subBlockDirectory.entries[entry].compression),
               dimension_stack_str(subBlockDirectory.entries[entry].dimensions, true).c_str(),
               pyramid_type_str(subBlockDirectory.entries[entry].pyramidType));
    }

    printf("======================================\n");
}

std::string byte_array_str(const ByteArray &bytes)
{
    std::string result = "";

    for (size_t i = 0; i < bytes.size(); i++)
    {
        result.append(std::to_string(bytes[i]));
        if (i % 4 == 0 && i != bytes.size() - 1)
            result.append("-");
    }
    return result;
}

void CziFile::report() const
{
    printf("Loaded CZI FILE: %s\n", fileName.c_str());

    printf("==========Summary report==========\n");
    print_segment_header(header.header);
    printf("%-25s %13i.%i\n", "FileVersion", header.fileVersion.major, header.fileVersion.minor);

    printf("%-25s %15i\n", "FilePart", header.filePart);

    printf("%-25s %30s\n", "MasterGUID", byte_array_str(header.masterFileGuid).c_str());
    printf("%-25s %30s\n", "FileGUID", byte_array_str(header.masterFileGuid).c_str());

    //TODO: Replace with GUID print.
    printf("%-25s %15i\n", "GUIDs are matching", vecUtil::vector_eq(header.masterFileGuid, header.fileGuid));
    printf("%-25s %15i\n", "UpdatePending", header.updatePending);
    printf("%-25s %15li\n", "SubblockDirPosition", header.subBlockDirectoryPosition);
    printf("%-25s %15li\n", "MetadataPosition", header.metadataPosition);
    printf("%-25s %15li\n", "AttachmentDirPosition", header.attachmentDirectoryPosition);
    printf("%-25s %15i\n", "SubBlockCount", subBlockDirectory.entryCount);
    printf("%-25s %15i\n", "ChannelCount", subBlockDirectory.channelCount);

    printf("%-25s %15i\n", "AttachmentCount", attachmentDirectory.entryCount);

    if (subBlockDirectory.entryCount > 0)
    {
        {
            std::map<PixelType, int> pixelTypeMap;
            std::map<CompressionType, int> compressionMap;
            std::map<std::string, int> dimensionMap;

            for (size_t sbId = 0; sbId < subBlockDirectory.entryCount; sbId++)
            {
                const DirectoryEntryDV entry = subBlockDirectory.entries[sbId];

                if (pixelTypeMap.count(entry.pixelType))
                    pixelTypeMap[entry.pixelType]++;
                else
                    pixelTypeMap.emplace(std::make_pair(entry.pixelType, 1));

                if (compressionMap.count(entry.compression))
                    compressionMap[entry.compression]++;
                else
                    compressionMap.emplace(std::make_pair(entry.compression, 1));

                std::string dimStackStr = dimension_stack_str(entry.dimensions, false);
                if (dimensionMap.count(dimStackStr))
                    dimensionMap[dimStackStr]++;
                else
                    dimensionMap.emplace(std::make_pair(dimStackStr, 1));
            }

            printf("Pixel type distribution:\n");
            for (auto &&ptPair : pixelTypeMap)
            {
                printf("\t%s %ix\n", pixel_type_str(ptPair.first), ptPair.second);
            }
            printf("Compression type distribution:\n");
            for (auto &&ctPair : compressionMap)
            {
                printf("\t%s %ix\n", compression_type_str(ctPair.first), ctPair.second);
            }
            printf("Dimension stack  distribution:\n");
            for (auto &&dimPair : dimensionMap)
            {
                printf("\t%s %ix\n", dimPair.first.c_str(), dimPair.second);
            }
        }

        //printf("%-25s %15i\n", "DimensionCount", subBlockDirectory.entries[0].dimensionCount);
        //printf("%-25s %15s\n", "PixelType", pixel_type_str(subBlockDirectory.entries[0].pixelType));
    }

    printf("==================================\n");
}

ImageMatrix CziFile::get_image(const uint subblockId) const
{
    auto imageBytes = get_image_data(subblockId, false);
    auto entry = subBlockDirectory.entries[subblockId];

    ImageMatrix image = ImageMatrix((uint)entry.width, (uint)entry.height, entry.pixelType, imageBytes);
    return image;
}

ByteArray CziFile::get_image_data(const uint subblockId, const bool ZCurveOrdered) const
{
    // subblockId is in bounds.
    always_assert((uint)subblockId < subBlockDirectory.entries.size());
    auto entry = subBlockDirectory.entries[subblockId];
    // Image is not empty.
    always_assert(entry.width > 0 && entry.height > 0);

    // Read image bytes.
    ByteArray imageBytes;
    {
        BinaryFileStream cziStream(fileName);
        imageBytes = cziStream.move_and_consume_bytes(entry.subBlock.dataLocation, entry.subBlock.dataSize);
    }
    // Check if read data are correct and not corrupted.
    int bytesPerPixel = get_bytes_per_pixel_type(entry.pixelType);
    ulong expectedSize = entry.width * entry.height * bytesPerPixel;
    always_assert(expectedSize == imageBytes.size());

    if (ZCurveOrdered)
    {
        always_assert((uint)entry.width <= UINT_MAX && (uint)entry.height <= UINT_MAX);
        uint componentCount = get_count_of_components_of_pixel(entry.pixelType);
        uint componentSize = get_component_size_of_pixel(entry.pixelType);

        auto zIndices = generate_ordered_z_order_indices(entry.width, entry.height, componentCount);
        always_assert(zIndices.size() == (ulong)(componentCount * entry.width * entry.height));

        ByteArray zOrderedBytes = reorder_bytes_to_z_order(imageBytes, zIndices, componentSize);

        return zOrderedBytes;
    }
    else
    {
        return imageBytes;
    }
}

ByteArray CziFile::get_continuous_image_data(const bool ZCurveOrdered) const
{
    // Find required size.
    size_t totalSize = 0;
    for (auto &&entry : subBlockDirectory.entries)
        totalSize += entry.subBlock.dataSize;

    // Allocate size.
    ByteArray data;
    data.resize(totalSize);
    size_t bufferEnd = 0;

    // Rather than simply putting images one after the other, order them by channels.
    // Compression ratios aren't better, so we won't waste time with sorting right now.
#if 0
    int loaded = 0;
    for (size_t channel = 0; channel < subBlockDirectory.channelCount; channel++)
    {
        if (ZCurveOrdered)
        {
            for (size_t subblockId = 0; subblockId < subBlockDirectory.entries.size(); subblockId++)
            {
                if (subBlockDirectory.entries[subblockId].channel != (int)channel)
                    continue;

                auto zOrderSubblockData = get_image_data(subblockId, true);
                vecUtil::vector_insert_at(data, zOrderSubblockData, bufferEnd, 0, zOrderSubblockData.size());
                bufferEnd += zOrderSubblockData.size();

                printf("\rLoaded Z-ordered subblock %i/%i", ++loaded, (int)subBlockDirectory.entries.size());
                fflush(stdout);
            }
        }
        else
        {
            BinaryFileStream stream(fileName);
            for (size_t subblockId = 0; subblockId < subBlockDirectory.entries.size(); subblockId++)
            {
                auto entry = subBlockDirectory.entries[subblockId];
                if (entry.channel != (int)channel)
                    continue;

                stream.move_to(entry.subBlock.dataLocation);
                stream.consume_into(data, bufferEnd, entry.subBlock.dataSize);
                bufferEnd += entry.subBlock.dataSize;

                printf("\rLoaded subblock %i/%i", ++loaded, (int)subBlockDirectory.entries.size());
                fflush(stdout);
            }
        }
    }
    printf("\n");

#else
    if (ZCurveOrdered)
    {
        for (size_t subblockId = 0; subblockId < subBlockDirectory.entries.size(); subblockId++)
        {
            auto zOrderSubblockData = get_image_data(subblockId, true);
            vecUtil::vector_insert_at(data, zOrderSubblockData, bufferEnd, 0, zOrderSubblockData.size());
            bufferEnd += zOrderSubblockData.size();

            printf("\rLoaded subblock %i/%i", (int)subblockId + 1, (int)subBlockDirectory.entries.size());
            fflush(stdout);
        }
        printf("\n");
    }
    else
    {
        BinaryFileStream stream(fileName);
        for (size_t subblockId = 0; subblockId < subBlockDirectory.entries.size(); subblockId++)
        {
            auto entry = subBlockDirectory.entries[subblockId];

            stream.move_to(entry.subBlock.dataLocation);
            stream.consume_into(data, bufferEnd, entry.subBlock.dataSize);
            bufferEnd += entry.subBlock.dataSize;

            printf("\rLoaded subblock %i/%i", (int)subblockId + 1, (int)subBlockDirectory.entries.size());
            fflush(stdout);
        }
        printf("\n");
    }
#endif
    return data;
}

void CziFile::dump_continuous_data(const std::string &dir) const
{
    always_assert(fs_wrapper::is_directory(dir));
    auto fName1 = dir + "/data.bin";
    auto fName2 = dir + "/z_order_data.bin";

    size_t totalSize = 0;
    for (auto &&entry : subBlockDirectory.entries)
        totalSize += entry.subBlock.dataSize;

    // Normal order
    {
        ByteArray data;
        data.resize(totalSize);

        size_t bufferEnd = 0;

        BinaryFileStream stream(fileName);
        for (size_t subblockId = 0; subblockId < subBlockDirectory.entries.size(); subblockId++)
        {
            auto entry = subBlockDirectory.entries[subblockId];

            stream.move_to(entry.subBlock.dataLocation);
            stream.consume_into(data, bufferEnd, entry.subBlock.dataSize);
            bufferEnd += entry.subBlock.dataSize;

            printf("\rLoaded subblock %i/%i", (int)subblockId + 1, (int)subBlockDirectory.entries.size());
            fflush(stdout);
        }
        printf("\n");

        // Write to file.
        write_bytes_to_file(data, fName1.c_str());
    }

    //Z-Order.
    {
        ByteArray data;
        data.resize(totalSize);

        size_t bufferEnd = 0;
        for (size_t subblockId = 0; subblockId < subBlockDirectory.entries.size(); subblockId++)
        {
            auto zOrderSubblockData = get_image_data(subblockId, true);
            vecUtil::vector_insert_at(data, zOrderSubblockData, bufferEnd, 0, zOrderSubblockData.size());
            bufferEnd += zOrderSubblockData.size();

            printf("\rLoaded subblock %i/%i", (int)subblockId + 1, (int)subBlockDirectory.entries.size());
            fflush(stdout);
        }
        printf("\n");
        // Write to file.
        write_bytes_to_file(data, fName2.c_str());
    }
}

void CziFile::absolutu_difference_between_frames(const std::string baseName) const
{
    if (subBlockDirectory.entries.size() <= 1)
    {
        printf("There aren't 2 different images in this CZI file.\n");
        return;
    }

    std::string fName;
    ImageMatrix ref = get_image(0);
    for (size_t i = 1; i < subBlockDirectory.entries.size(); i++)
    {
        fName = baseName + "diff_" + std::to_string(i - 1) + "_" + std::to_string(i) + ".ppm";

        ImageMatrix next = get_image(i);
        ImageMatrix diff_ref_next = ref.create_absolute_difference_matrix(next);

        diff_ref_next.save_as_ppm(fName.c_str());

        ref = next;
    }
}

void CziFile::dump_images(const std::string &baseName) const
{
    //#pragma omp parallel for
    for (size_t i = 0; i < subBlockDirectory.entries.size(); i++)
    {
        auto subblock = get_image(i);
        std::string fName = baseName + "subblock_" + std::to_string(i) + ".ppm";
        subblock.save_as_ppm(fName.c_str());
    }
}

void CziFile::dump_data(const std::string &baseName) const
{
    BinaryFileStream cziStream(fileName);
    int imageIndex = 0;
    for (const DirectoryEntryDV &entry : subBlockDirectory.entries)
    {
        always_assert(entry.width > 0 && entry.height > 0);

        std::string binaryFileName = baseName + std::to_string(imageIndex++) + ".bin";

        auto imageBytes = cziStream.move_and_consume_bytes(entry.subBlock.dataLocation, entry.subBlock.dataSize);

        write_bytes_to_file(imageBytes, binaryFileName.c_str());

        printf("Wrote %s\n", binaryFileName.c_str());
    }
}

void CziFile::test_compression(CompressionMethod method, bool verbose, int compressionLevel) const
{

    switch (method)
    {
    case CompressionMethod_GZIP:
        printf("Selected compression: GZIP (zlib)\n");
        break;
    case CompressionMethod_LZMA:
        printf("Selected compression: LZMA (2?)\n");
        break;
    case CompressionMethod_BZIP2:
        printf("Selected compression: BZIP2 \n");
        break;
    default:
        INVALID_CASE;
        return;
    }
    printf("\tCompression level: %i\n", compressionLevel);

    CompressionResult overallN;
    CompressionResult overallZ;
    Stopwatch stopwatch;

    size_t done = 0;

    // #pragma omp parallel for
    for (size_t i = 0; i < subBlockDirectory.entries.size(); i++)
    {
        auto data = get_image_data(i, false);
        auto dataZ = get_image_data(i, true);

        CompressionResult nResult = test_compression_method(data, method, compressionLevel);
        CompressionResult zResult = test_compression_method(dataZ, method, compressionLevel);
        // #pragma omp critical
        //         {
        done += 1;
        float perc = (float)done / (float)subBlockDirectory.entries.size();

        printf("Done: %f %%\n", perc * 100.0f);

        overallN.compressionRatio += nResult.compressionRatio;
        overallZ.compressionRatio += zResult.compressionRatio;
        overallN.originalSize += data.size();
        overallN.compressedSize += nResult.compressedSize;
        overallZ.compressedSize += zResult.compressedSize;
        overallN.percentageOfOriginalSize += nResult.percentageOfOriginalSize;
        overallZ.percentageOfOriginalSize += zResult.percentageOfOriginalSize;

        overallN.compressionTimeMS += nResult.compressionTimeMS;
        overallZ.compressionTimeMS += zResult.compressionTimeMS;

        if (verbose)
        {
            printf("Subblock %-3i Compression ratios: Normal: %8f Z-Ordered: %8f; Size(N): %5.3f%%; Size(Z): %5.3f%%; Time: %8f ms\n",
                   (int)i, nResult.compressionRatio, zResult.compressionRatio, nResult.percentageOfOriginalSize,
                   zResult.percentageOfOriginalSize, nResult.compressionTimeMS);
        }
        // }
    }
    float dataCount = (float)subBlockDirectory.entries.size();

    overallN.divide(dataCount);
    overallZ.divide(dataCount);

    printf("Overall compression ratios: Normal %8f Z-Ordered: %8f\n", overallN.compressionRatio, overallZ.compressionRatio);
    printf("Original size: %8lu B Compressed size: %8lu B Compressed Z-Order size: %8lu B\n", overallN.originalSize, overallN.compressedSize, overallZ.compressedSize);
    printf("Original size (%%): 100%% Compressed: %5.5f%% Compressed Z-Order: %5.5f%%\n", overallN.percentageOfOriginalSize, overallZ.percentageOfOriginalSize);
    printf("Average compression time: %f ms\n", overallN.compressionTimeMS);
}

std::vector<std::pair<uint, std::vector<uint>>> CziFile::get_subblocks_grouped_by_channels() const
{
    std::vector<std::pair<uint, std::vector<uint>>> result;
    result.reserve(subBlockDirectory.channelCount);

    for (uint channel = 0; channel < subBlockDirectory.channelCount; channel++)
    {
        std::vector<uint> subblocksInChannel;
        for (uint subblock = 0; subblock < subBlockDirectory.entryCount; subblock++)
        {
            int subblockChannel = subBlockDirectory.entries[subblock].channel;
            always_assert(subblockChannel >= 0);
            if ((uint)subblockChannel == channel)
            {
                subblocksInChannel.push_back(subblock);
            }
        }
        result.push_back(std::make_pair(channel, subblocksInChannel));
    }

    return result;
}

void CziFile::frames_difference() const
{
    // TODO: This will now handle only Gray16 pixels!

    auto framesByChannels = get_subblocks_grouped_by_channels();

    std::pair<uint, uint> canBeMappedCounts = std::make_pair(0, 0);

    for (const std::pair<uint, std::vector<uint>> &channelGroup : framesByChannels)
    {
        printf("Starting channel %u\n", channelGroup.first);

        //#pragma omp parallel for

        for (size_t i = 1; i < channelGroup.second.size(); i++)
        {
            // This assertion will fail if pixel type isn't Gray16
            always_assert(subBlockDirectory.entries[channelGroup.second[i - 1]].pixelType == PixelType_Gray16);
            always_assert(subBlockDirectory.entries[channelGroup.second[i]].pixelType == PixelType_Gray16);

            auto prevFrameData = get_image_data(channelGroup.second[i - 1], false);
            auto currentFrameData = get_image_data(channelGroup.second[i], false);
            always_assert(prevFrameData.size() == currentFrameData.size());

            std::vector<ushort> prevFrameValues = bytes_to_ushort_array(prevFrameData);
            std::vector<ushort> currentFrameValues = bytes_to_ushort_array(currentFrameData);
            always_assert(prevFrameValues.size() == currentFrameValues.size());

            std::vector<int> differenceArray = vecUtil::diff_vectors<ushort, int>(prevFrameValues, currentFrameValues);

            int min = vecUtil::find_min(differenceArray);
            int max = vecUtil::find_max(differenceArray);
            long maxMappedValue = (min < 0) ? (abs(min) + max) : (max);
            bool canBeMappedToUShort = maxMappedValue < USHORT_MAX;

#pragma omp critical
            {
                if (canBeMappedToUShort)
                    canBeMappedCounts.first++;
                else
                    canBeMappedCounts.second++;
            }

            ByteArray intBytes = int_array_to_bytes(differenceArray);

            auto frameCompResult = test_compression_method(currentFrameData, CompressionMethod_BZIP2, 6);
            auto diffCompResult = test_compression_method(intBytes, CompressionMethod_BZIP2, 6);

            float diffSizeDif = ((float)diffCompResult.compressedSize / (float)frameCompResult.compressedSize);

            printf("========================\nFrame [%u vs %u]\n\tCR frame: %10.5f\n\tCR  diff: %10.5f\n\tSizeDiff: %10.5f\n\tMin: %10i\n\tMax: %10i\n\tMax mapped value: %10li\n\tCan be mapped to ushort: %s\n",
                   channelGroup.second[i - 1], channelGroup.second[i], frameCompResult.compressionRatio, diffCompResult.compressionRatio, diffSizeDif,
                   min, max, maxMappedValue, canBeMappedToUShort ? "YES" : "NO");

            //printf("Max value: %7i Min value: %7i #of values: %7li\n", max, min, rangeSize);
        }
    }

    printf("%u can be mapped to short difference value.\n%u can not be mapped.\n", canBeMappedCounts.first, canBeMappedCounts.second);
}