// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// png.hpp — PNG 画像ファイルの読み書き
//
// 外部ライブラリ不要の自前 deflate (動的 Huffman + LZ77) 実装。
// RFC 1951 (DEFLATE), RFC 1950 (zlib), RFC 2083 (PNG) 準拠。
//
// API:
//   savePng(path, Image<T>)   — 8bit RGB/RGBA/Gray PNG として保存
//   loadPng<T>(path)          — PNG を Image<T> として読み込み

#ifndef CALX_IMAGE_IO_PNG_HPP
#define CALX_IMAGE_IO_PNG_HPP

#include <image/Image.hpp>
#include <image/io/bmp.hpp>   // detail::toU8, detail::fromU8

#include <algorithm>
#include <array>
#include <cstdint>
#include <cstring>
#include <fstream>
#include <numeric>
#include <stdexcept>
#include <string>
#include <vector>

namespace calx::image {

// ========================================================================
// CRC-32 (PNG チャンク用)
// ========================================================================

namespace detail_png {

inline uint32_t crc32(const uint8_t* data, std::size_t len) {
    static const auto table = [] {
        std::array<uint32_t, 256> t{};
        for (uint32_t i = 0; i < 256; ++i) {
            uint32_t c = i;
            for (int j = 0; j < 8; ++j)
                c = (c & 1) ? (0xEDB88320u ^ (c >> 1)) : (c >> 1);
            t[i] = c;
        }
        return t;
    }();
    uint32_t crc = 0xFFFFFFFFu;
    for (std::size_t i = 0; i < len; ++i)
        crc = table[(crc ^ data[i]) & 0xFF] ^ (crc >> 8);
    return crc ^ 0xFFFFFFFFu;
}

inline uint32_t adler32(const uint8_t* data, std::size_t len) {
    uint32_t a = 1, b = 0;
    constexpr uint32_t MOD = 65521;
    for (std::size_t i = 0; i < len; ++i) {
        a = (a + data[i]) % MOD;
        b = (b + a) % MOD;
    }
    return (b << 16) | a;
}

// ========================================================================
// Huffman ツリー構築
// ========================================================================

struct HuffCode {
    uint16_t code;
    uint8_t  len;
};

// 頻度 → コード長 (制限付き: maxBits 以下)
inline void freqToLengths(const uint32_t* freq, int n,
                          uint8_t* lengths, int maxBits) {
    // Package-Merge アルゴリズム (長さ制限付き最適 Huffman)
    // 簡略化: まず単純な方法でコード長を計算し、maxBits に制限

    struct Node {
        uint32_t weight;
        int symbol;
        int left, right;  // -1 = leaf
    };

    // 使用シンボルを収集
    std::vector<int> symbols;
    for (int i = 0; i < n; ++i)
        if (freq[i] > 0) symbols.push_back(i);

    std::fill(lengths, lengths + n, uint8_t(0));

    if (symbols.empty()) return;
    if (symbols.size() == 1) {
        lengths[symbols[0]] = 1;
        return;
    }

    // ノードプール
    std::vector<Node> nodes;
    nodes.reserve(symbols.size() * 2);
    for (int s : symbols)
        nodes.push_back({freq[s], s, -1, -1});

    // Huffman 木構築 (優先度キュー代わりにソート)
    std::vector<int> active(symbols.size());
    std::iota(active.begin(), active.end(), 0);

    while (active.size() > 1) {
        // 最小 2 つを取り出し
        std::sort(active.begin(), active.end(), [&](int a, int b) {
            return nodes[a].weight > nodes[b].weight;
        });
        int a = active.back(); active.pop_back();
        int b = active.back(); active.pop_back();
        int idx = static_cast<int>(nodes.size());
        nodes.push_back({nodes[a].weight + nodes[b].weight, -1, a, b});
        active.push_back(idx);
    }

    // 深さ計算 (スタックで再帰を回避)
    {
        struct Frame { int idx; int depth; };
        std::vector<Frame> stack;
        stack.push_back({active[0], 0});
        while (!stack.empty()) {
            auto [idx, depth] = stack.back();
            stack.pop_back();
            if (nodes[idx].left == -1) {
                lengths[nodes[idx].symbol] =
                    static_cast<uint8_t>(std::min(depth, 255));
            } else {
                stack.push_back({nodes[idx].left, depth + 1});
                stack.push_back({nodes[idx].right, depth + 1});
            }
        }
    }

    // maxBits 制限
    // 1. オーバーフローをカウントしつつクランプ
    int bl_count[16] = {};  // maxBits ≤ 15
    int overflowCount = 0;
    for (int i = 0; i < n; ++i) {
        if (lengths[i] == 0) continue;
        if (lengths[i] > maxBits) {
            lengths[i] = static_cast<uint8_t>(maxBits);
            overflowCount++;
        }
        bl_count[lengths[i]]++;
    }

    if (overflowCount == 0) return;

    // 2. オーバーフロー修正
    //    各反復: 短いコードを 1 段深くして 2 つに分割、maxBits のコードを 1 つ除去
    //    → Kraft 和を 2^(-maxBits) ずつ削減
    //    Kraft 和が 1 以下になるまでループ (深い overflow でも確実に修正)
    auto kraftOk = [&]() {
        // Kraft 和を 2^maxBits 倍した整数で判定 (浮動小数点誤差回避)
        uint64_t sum = 0;
        for (int b = 1; b <= maxBits; ++b)
            sum += static_cast<uint64_t>(bl_count[b]) << (maxBits - b);
        return sum <= (uint64_t(1) << maxBits);
    };

    while (!kraftOk()) {
        int bits = maxBits - 1;
        while (bits > 0 && bl_count[bits] == 0) bits--;
        if (bits == 0) break;       // 安全弁
        bl_count[bits]--;           // 短いコードを 1 つ除去
        bl_count[bits + 1] += 2;    // 1 段深くして 2 つに分割
        bl_count[maxBits]--;        // maxBits のオーバーフロー 1 つ除去
    }

    // 3. bl_count に基づいてシンボルにコード長を再割当
    //    頻度昇順 (低頻度 → 長いコード)
    struct SymFreq { int idx; uint32_t f; };
    std::vector<SymFreq> sorted;
    for (int i = 0; i < n; ++i)
        if (freq[i] > 0) sorted.push_back({i, freq[i]});
    std::sort(sorted.begin(), sorted.end(), [](const SymFreq& a, const SymFreq& b) {
        return a.f < b.f || (a.f == b.f && a.idx < b.idx);
    });

    std::fill(lengths, lengths + n, uint8_t(0));
    int si = 0;
    for (int len = maxBits; len >= 1 && si < static_cast<int>(sorted.size()); --len) {
        for (int j = 0; j < bl_count[len] && si < static_cast<int>(sorted.size()); ++j) {
            lengths[sorted[si].idx] = static_cast<uint8_t>(len);
            ++si;
        }
    }
}

// コード長 → canonical Huffman コード
inline void lengthsToCodes(const uint8_t* lengths, int n,
                           HuffCode* codes, int maxBits) {
    std::vector<int> blCount(maxBits + 1, 0);
    for (int i = 0; i < n; ++i)
        if (lengths[i] > 0) blCount[lengths[i]]++;

    std::vector<uint32_t> nextCode(maxBits + 1, 0);
    uint32_t code = 0;
    for (int bits = 1; bits <= maxBits; ++bits) {
        code = (code + blCount[bits - 1]) << 1;
        nextCode[bits] = code;
    }

    for (int i = 0; i < n; ++i) {
        if (lengths[i] == 0) {
            codes[i] = {0, 0};
        } else {
            // ビット反転 (deflate はビット逆順)
            uint16_t c = static_cast<uint16_t>(nextCode[lengths[i]]++);
            uint16_t rev = 0;
            for (int b = 0; b < lengths[i]; ++b)
                rev |= ((c >> b) & 1) << (lengths[i] - 1 - b);
            codes[i] = {rev, lengths[i]};
        }
    }
}

// ========================================================================
// ビットストリームライター
// ========================================================================

class BitWriter {
public:
    void writeBits(uint32_t value, int nbits) {
        buf_ |= static_cast<uint64_t>(value) << bitPos_;
        bitPos_ += nbits;
        flush();
    }

    void writeCode(HuffCode c) {
        writeBits(c.code, c.len);
    }

    void alignByte() {
        if (bitPos_ % 8 != 0)
            bitPos_ = ((bitPos_ + 7) / 8) * 8;
        flush();
    }

    const std::vector<uint8_t>& data() const { return out_; }

    void finish() {
        while (bitPos_ > 0) {
            out_.push_back(static_cast<uint8_t>(buf_ & 0xFF));
            buf_ >>= 8;
            bitPos_ = (bitPos_ >= 8) ? bitPos_ - 8 : 0;
        }
    }

private:
    std::vector<uint8_t> out_;
    uint64_t buf_ = 0;
    int bitPos_ = 0;

    void flush() {
        while (bitPos_ >= 8) {
            out_.push_back(static_cast<uint8_t>(buf_ & 0xFF));
            buf_ >>= 8;
            bitPos_ -= 8;
        }
    }
};

// ========================================================================
// Deflate エンコーダで使う定数テーブル
// ========================================================================

// Length コード (257-285) のベース値と追加ビット数
struct LenEntry { uint16_t base; uint8_t extra; };

inline const LenEntry* lenTable() {
    static const LenEntry t[] = {
        {3,0},{4,0},{5,0},{6,0},{7,0},{8,0},{9,0},{10,0},
        {11,1},{13,1},{15,1},{17,1},{19,2},{23,2},{27,2},{31,2},
        {35,3},{43,3},{51,3},{59,3},{67,4},{83,4},{99,4},{115,4},
        {131,5},{163,5},{195,5},{227,5},{258,0}
    };
    return t;
}

// Distance コード (0-29) のベース値と追加ビット数
struct DistEntry { uint16_t base; uint8_t extra; };

inline const DistEntry* distTable() {
    static const DistEntry t[] = {
        {1,0},{2,0},{3,0},{4,0},{5,1},{7,1},{9,2},{13,2},
        {17,3},{25,3},{33,4},{49,4},{65,5},{97,5},{129,6},{193,6},
        {257,7},{385,7},{513,8},{769,8},{1025,9},{1537,9},{2049,10},{3073,10},
        {4097,11},{6145,11},{8193,12},{12289,12},{16385,13},{24577,13}
    };
    return t;
}

// 長さ → length コード (257+)
inline int encodeLenCode(int length) {
    auto* t = lenTable();
    for (int i = 28; i >= 0; --i)
        if (length >= t[i].base) return 257 + i;
    return 257;
}

// 距離 → distance コード
inline int encodeDistCode(int dist) {
    auto* t = distTable();
    for (int i = 29; i >= 0; --i)
        if (dist >= t[i].base) return i;
    return 0;
}

// ========================================================================
// LZ77 エンコーダ (ハッシュチェーン)
// ========================================================================

struct LZ77Token {
    enum Type : uint8_t { LITERAL, MATCH };
    Type type;
    union {
        uint8_t  literal;
        struct { uint16_t length; uint16_t distance; } match;
    };
};

inline std::vector<LZ77Token> lz77Encode(const uint8_t* data, std::size_t len) {
    constexpr int WINDOW = 32768;
    constexpr int MAX_MATCH = 258;
    constexpr int MIN_MATCH = 3;
    constexpr int HASH_SIZE = 1 << 15;
    constexpr int HASH_MASK = HASH_SIZE - 1;
    constexpr int MAX_CHAIN = 128;  // チェーン探索の上限

    std::vector<LZ77Token> tokens;
    tokens.reserve(len / 2);

    std::vector<int> head(HASH_SIZE, -1);
    std::vector<int> prev(len, -1);

    auto hashFunc = [](const uint8_t* p) -> int {
        return ((p[0] << 10) ^ (p[1] << 5) ^ p[2]) & HASH_MASK;
    };

    std::size_t pos = 0;
    while (pos < len) {
        int bestLen = 0;
        int bestDist = 0;

        if (pos + MIN_MATCH <= len) {
            int h = hashFunc(data + pos);
            int chain = head[h];
            int chainCount = 0;

            while (chain >= 0 && chainCount < MAX_CHAIN) {
                int dist = static_cast<int>(pos) - chain;
                if (dist > WINDOW) break;

                // マッチ長を計測
                int maxL = std::min(MAX_MATCH, static_cast<int>(len - pos));
                int l = 0;
                while (l < maxL && data[pos + l] == data[chain + l]) ++l;

                if (l > bestLen) {
                    bestLen = l;
                    bestDist = dist;
                    if (l == MAX_MATCH) break;
                }

                chain = prev[chain];
                ++chainCount;
            }

            // ハッシュチェーンに登録
            if (pos < len - 2) {
                prev[pos] = head[h];
                head[h] = static_cast<int>(pos);
            }
        }

        if (bestLen >= MIN_MATCH) {
            LZ77Token tok;
            tok.type = LZ77Token::MATCH;
            tok.match.length = static_cast<uint16_t>(bestLen);
            tok.match.distance = static_cast<uint16_t>(bestDist);
            tokens.push_back(tok);

            // マッチした範囲のハッシュも登録 (lazy matching 簡略版)
            for (int i = 1; i < bestLen && pos + i + 2 < len; ++i) {
                int hh = hashFunc(data + pos + i);
                prev[pos + i] = head[hh];
                head[hh] = static_cast<int>(pos + i);
            }
            pos += bestLen;
        } else {
            LZ77Token tok;
            tok.type = LZ77Token::LITERAL;
            tok.literal = data[pos];
            tokens.push_back(tok);
            ++pos;
        }
    }

    return tokens;
}

// ========================================================================
// Deflate 圧縮 (動的 Huffman)
// ========================================================================

inline std::vector<uint8_t> deflateCompress(const uint8_t* data, std::size_t len) {
    // LZ77
    auto tokens = lz77Encode(data, len);

    // 頻度カウント
    constexpr int LITLEN_COUNT = 286;  // 0-285
    constexpr int DIST_COUNT = 30;     // 0-29
    uint32_t litFreq[LITLEN_COUNT] = {};
    uint32_t distFreq[DIST_COUNT] = {};

    for (const auto& tok : tokens) {
        if (tok.type == LZ77Token::LITERAL) {
            litFreq[tok.literal]++;
        } else {
            int lc = encodeLenCode(tok.match.length);
            litFreq[lc]++;
            int dc = encodeDistCode(tok.match.distance);
            distFreq[dc]++;
        }
    }
    litFreq[256] = 1;  // end-of-block

    // Huffman コード長計算
    uint8_t litLengths[LITLEN_COUNT] = {};
    uint8_t distLengths[DIST_COUNT] = {};
    freqToLengths(litFreq, LITLEN_COUNT, litLengths, 15);
    freqToLengths(distFreq, DIST_COUNT, distLengths, 15);

    // distance が全て 0 の場合 (リテラルのみ)
    int maxDistCode = 0;
    for (int i = DIST_COUNT - 1; i >= 0; --i) {
        if (distLengths[i] > 0) { maxDistCode = i; break; }
    }
    int HDIST = maxDistCode + 1;
    if (HDIST < 1) {
        HDIST = 1;
        distLengths[0] = 1;  // ダミー
    }

    // HLIT: lit/len コード数
    int maxLitCode = 0;
    for (int i = LITLEN_COUNT - 1; i >= 0; --i) {
        if (litLengths[i] > 0) { maxLitCode = i; break; }
    }
    int HLIT = std::max(maxLitCode + 1, 257);

    // コード長列をランレングス符号化
    std::vector<uint8_t> allLengths(HLIT + HDIST);
    for (int i = 0; i < HLIT; ++i) allLengths[i] = litLengths[i];
    for (int i = 0; i < HDIST; ++i) allLengths[HLIT + i] = distLengths[i];

    // RLE 符号化 (16: 直前を繰り返し, 17: 0を3-10個, 18: 0を11-138個)
    struct RleEntry { uint8_t symbol; uint8_t extra; uint8_t extraBits; };
    std::vector<RleEntry> rle;

    for (int i = 0; i < static_cast<int>(allLengths.size()); ) {
        if (allLengths[i] == 0) {
            int run = 1;
            while (i + run < static_cast<int>(allLengths.size()) &&
                   allLengths[i + run] == 0 && run < 138) ++run;
            if (run >= 11) {
                rle.push_back({18, static_cast<uint8_t>(run - 11), 7});
                i += run;
            } else if (run >= 3) {
                rle.push_back({17, static_cast<uint8_t>(run - 3), 3});
                i += run;
            } else {
                for (int j = 0; j < run; ++j)
                    rle.push_back({0, 0, 0});
                i += run;
            }
        } else {
            uint8_t val = allLengths[i];
            rle.push_back({val, 0, 0});
            ++i;
            int run = 0;
            while (i + run < static_cast<int>(allLengths.size()) &&
                   allLengths[i + run] == val && run < 6) ++run;
            if (run >= 3) {
                rle.push_back({16, static_cast<uint8_t>(run - 3), 2});
                i += run;
            }
        }
    }

    // コード長の Huffman (0-18 のシンボル)
    constexpr int CL_COUNT = 19;
    uint32_t clFreq[CL_COUNT] = {};
    for (const auto& r : rle) clFreq[r.symbol]++;

    uint8_t clLengths[CL_COUNT] = {};
    freqToLengths(clFreq, CL_COUNT, clLengths, 7);

    // deflate 規定の並び順
    static const int clOrder[CL_COUNT] = {
        16, 17, 18, 0, 8, 7, 9, 6, 10, 5, 11, 4, 12, 3, 13, 2, 14, 1, 15
    };
    int HCLEN = CL_COUNT;
    while (HCLEN > 4 && clLengths[clOrder[HCLEN - 1]] == 0) --HCLEN;

    // Huffman コード生成
    HuffCode litCodes[LITLEN_COUNT];
    HuffCode distCodes[DIST_COUNT];
    HuffCode clCodes[CL_COUNT];
    lengthsToCodes(litLengths, LITLEN_COUNT, litCodes, 15);
    lengthsToCodes(distLengths, DIST_COUNT, distCodes, 15);
    lengthsToCodes(clLengths, CL_COUNT, clCodes, 7);

    // ビットストリーム出力
    BitWriter bw;

    // BFINAL=1, BTYPE=10 (動的 Huffman)
    bw.writeBits(1, 1);   // BFINAL
    bw.writeBits(2, 2);   // BTYPE = dynamic

    // HLIT, HDIST, HCLEN
    bw.writeBits(HLIT - 257, 5);
    bw.writeBits(HDIST - 1, 5);
    bw.writeBits(HCLEN - 4, 4);

    // コード長テーブル
    for (int i = 0; i < HCLEN; ++i)
        bw.writeBits(clLengths[clOrder[i]], 3);

    // RLE 符号化されたコード長
    for (const auto& r : rle) {
        bw.writeCode(clCodes[r.symbol]);
        if (r.extraBits > 0)
            bw.writeBits(r.extra, r.extraBits);
    }

    // トークンを出力
    auto* lt = lenTable();
    auto* dt = distTable();

    for (const auto& tok : tokens) {
        if (tok.type == LZ77Token::LITERAL) {
            bw.writeCode(litCodes[tok.literal]);
        } else {
            int lc = encodeLenCode(tok.match.length);
            bw.writeCode(litCodes[lc]);
            int li = lc - 257;
            if (lt[li].extra > 0)
                bw.writeBits(tok.match.length - lt[li].base, lt[li].extra);

            int dc = encodeDistCode(tok.match.distance);
            bw.writeCode(distCodes[dc]);
            if (dt[dc].extra > 0)
                bw.writeBits(tok.match.distance - dt[dc].base, dt[dc].extra);
        }
    }

    // end-of-block
    bw.writeCode(litCodes[256]);
    bw.finish();

    return bw.data();
}

// ========================================================================
// Deflate 展開 (読み込み用)
// ========================================================================

class BitReader {
public:
    BitReader(const uint8_t* data, std::size_t len)
        : data_(data), len_(len) {}

    uint32_t readBits(int n) {
        while (bitPos_ < n) {
            if (bytePos_ >= len_)
                throw std::runtime_error("BitReader: unexpected end");
            buf_ |= static_cast<uint64_t>(data_[bytePos_++]) << bitPos_;
            bitPos_ += 8;
        }
        uint32_t val = static_cast<uint32_t>(buf_) & ((1u << n) - 1);
        buf_ >>= n;
        bitPos_ -= n;
        return val;
    }

    void alignByte() {
        int skip = bitPos_ % 8;
        if (skip > 0) { buf_ >>= skip; bitPos_ -= skip; }
    }

    // 読みすぎたビットを戻す
    void unreadBits(uint32_t value, int n) {
        buf_ = (buf_ << n) | static_cast<uint64_t>(value);
        bitPos_ += n;
    }

private:
    const uint8_t* data_;
    std::size_t len_;
    std::size_t bytePos_ = 0;
    uint64_t buf_ = 0;
    int bitPos_ = 0;
};

// Huffman デコードテーブル (最大 15 ビット)
struct HuffDecoder {
    // テーブルルックアップ方式: LSB ファーストでインデックス
    std::vector<int16_t> table;   // symbol or -1
    std::vector<uint8_t> codeLens; // 各エントリの実際のコード長
    int maxBits = 0;

    void build(const uint8_t* lengths, int n) {
        maxBits = *std::max_element(lengths, lengths + n);
        if (maxBits == 0) return;
        if (maxBits > 15) maxBits = 15;

        // canonical コードを再構築
        std::vector<int> blCount(maxBits + 1, 0);
        for (int i = 0; i < n; ++i)
            if (lengths[i] > 0) blCount[lengths[i]]++;

        std::vector<int> nextCode(maxBits + 1, 0);
        int code = 0;
        for (int bits = 1; bits <= maxBits; ++bits) {
            code = (code + blCount[bits - 1]) << 1;
            nextCode[bits] = code;
        }

        // LSB ファーストのルックアップテーブルを構築
        std::size_t tableSize = std::size_t(1) << maxBits;
        table.assign(tableSize, -1);
        codeLens.assign(tableSize, 0);

        for (int i = 0; i < n; ++i) {
            if (lengths[i] == 0) continue;
            int c = nextCode[lengths[i]]++;
            int len = lengths[i];
            // canonical code を LSB ファーストに反転
            int rev = 0;
            for (int b = 0; b < len; ++b)
                rev |= ((c >> (len - 1 - b)) & 1) << b;
            // 残りビットの全パターンを埋める
            int shift = maxBits - len;
            for (int fill = 0; fill < (1 << shift); ++fill) {
                int idx = rev | (fill << len);
                table[idx] = static_cast<int16_t>(i);
                codeLens[idx] = static_cast<uint8_t>(len);
            }
        }
    }

    int decode(BitReader& br) const {
        if (maxBits == 0) return -1;
        // maxBits ビットをピーク (LSB ファースト)
        uint32_t bits = br.readBits(maxBits);
        int16_t sym = table[bits];
        if (sym < 0) throw std::runtime_error("Huffman decode error");
        // 実際に消費したのは codeLens[bits] ビットだけ
        // 余分に読んだビットを返す
        int consumed = codeLens[bits];
        int excess = maxBits - consumed;
        if (excess > 0)
            br.unreadBits(bits >> consumed, excess);
        return sym;
    }
};

inline std::vector<uint8_t> deflateDecompress(const uint8_t* data, std::size_t len) {
    BitReader br(data, len);
    std::vector<uint8_t> out;
    out.reserve(len * 4);

    auto* lt = lenTable();
    auto* dt = distTable();

    bool bfinal = false;
    while (!bfinal) {
        bfinal = br.readBits(1) != 0;
        uint32_t btype = br.readBits(2);

        if (btype == 0) {
            // 非圧縮
            br.alignByte();
            uint32_t blockLen = br.readBits(16);
            br.readBits(16);  // NLEN (無視)
            for (uint32_t i = 0; i < blockLen; ++i)
                out.push_back(static_cast<uint8_t>(br.readBits(8)));
        } else if (btype == 1 || btype == 2) {
            HuffDecoder litDec, distDec;

            if (btype == 1) {
                // 固定 Huffman
                uint8_t litLen[288] = {};
                for (int i = 0; i <= 143; ++i)   litLen[i] = 8;
                for (int i = 144; i <= 255; ++i) litLen[i] = 9;
                for (int i = 256; i <= 279; ++i) litLen[i] = 7;
                for (int i = 280; i <= 287; ++i) litLen[i] = 8;
                litDec.build(litLen, 288);

                uint8_t distLen[32] = {};
                for (int i = 0; i < 32; ++i) distLen[i] = 5;
                distDec.build(distLen, 32);
            } else {
                // 動的 Huffman
                int hlit = br.readBits(5) + 257;
                int hdist = br.readBits(5) + 1;
                int hclen = br.readBits(4) + 4;

                static const int clOrder[19] = {
                    16, 17, 18, 0, 8, 7, 9, 6, 10, 5, 11, 4, 12, 3, 13, 2, 14, 1, 15
                };
                uint8_t clLen[19] = {};
                for (int i = 0; i < hclen; ++i)
                    clLen[clOrder[i]] = static_cast<uint8_t>(br.readBits(3));

                HuffDecoder clDec;
                clDec.build(clLen, 19);

                // コード長を復号
                std::vector<uint8_t> codeLens(hlit + hdist, 0);
                for (int i = 0; i < hlit + hdist; ) {
                    int sym = clDec.decode(br);
                    if (sym < 16) {
                        codeLens[i++] = static_cast<uint8_t>(sym);
                    } else if (sym == 16) {
                        int rep = br.readBits(2) + 3;
                        uint8_t prev = (i > 0) ? codeLens[i - 1] : 0;
                        for (int j = 0; j < rep; ++j) codeLens[i++] = prev;
                    } else if (sym == 17) {
                        int rep = br.readBits(3) + 3;
                        for (int j = 0; j < rep; ++j) codeLens[i++] = 0;
                    } else if (sym == 18) {
                        int rep = br.readBits(7) + 11;
                        for (int j = 0; j < rep; ++j) codeLens[i++] = 0;
                    }
                }

                litDec.build(codeLens.data(), hlit);
                distDec.build(codeLens.data() + hlit, hdist);
            }

            // デコード
            while (true) {
                int sym = litDec.decode(br);
                if (sym < 256) {
                    out.push_back(static_cast<uint8_t>(sym));
                } else if (sym == 256) {
                    break;
                } else {
                    int li = sym - 257;
                    int length = lt[li].base;
                    if (lt[li].extra > 0)
                        length += br.readBits(lt[li].extra);

                    int dc = distDec.decode(br);
                    int distance = dt[dc].base;
                    if (dt[dc].extra > 0)
                        distance += br.readBits(dt[dc].extra);

                    for (int i = 0; i < length; ++i)
                        out.push_back(out[out.size() - distance]);
                }
            }
        } else {
            throw std::runtime_error("deflate: invalid block type");
        }
    }

    return out;
}

// ========================================================================
// PNG フィルタ
// ========================================================================

// Paeth 予測関数
inline uint8_t paethPredictor(int a, int b, int c) {
    int p = a + b - c;
    int pa = std::abs(p - a);
    int pb = std::abs(p - b);
    int pc = std::abs(p - c);
    if (pa <= pb && pa <= pc) return static_cast<uint8_t>(a);
    if (pb <= pc) return static_cast<uint8_t>(b);
    return static_cast<uint8_t>(c);
}

// 各フィルタの適用 → 絶対値和が最小のフィルタを選択
inline uint8_t chooseBestFilter(const uint8_t* row, const uint8_t* prevRow,
                                int stride, int bpp) {
    uint32_t bestSum = UINT32_MAX;
    uint8_t bestFilter = 0;

    // フィルタ 0: None
    {
        uint32_t s = 0;
        for (int i = 0; i < stride; ++i)
            s += row[i];
        if (s < bestSum) { bestSum = s; bestFilter = 0; }
    }

    // フィルタ 1: Sub
    {
        uint32_t s = 0;
        for (int i = 0; i < stride; ++i) {
            int a = (i >= bpp) ? row[i - bpp] : 0;
            s += static_cast<uint8_t>(row[i] - a);
        }
        if (s < bestSum) { bestSum = s; bestFilter = 1; }
    }

    // フィルタ 2: Up
    {
        uint32_t s = 0;
        for (int i = 0; i < stride; ++i) {
            int b = prevRow ? prevRow[i] : 0;
            s += static_cast<uint8_t>(row[i] - b);
        }
        if (s < bestSum) { bestSum = s; bestFilter = 2; }
    }

    // フィルタ 3: Average
    {
        uint32_t s = 0;
        for (int i = 0; i < stride; ++i) {
            int a = (i >= bpp) ? row[i - bpp] : 0;
            int b = prevRow ? prevRow[i] : 0;
            s += static_cast<uint8_t>(row[i] - (a + b) / 2);
        }
        if (s < bestSum) { bestSum = s; bestFilter = 3; }
    }

    // フィルタ 4: Paeth
    {
        uint32_t s = 0;
        for (int i = 0; i < stride; ++i) {
            int a = (i >= bpp) ? row[i - bpp] : 0;
            int b = prevRow ? prevRow[i] : 0;
            int c = (prevRow && i >= bpp) ? prevRow[i - bpp] : 0;
            s += static_cast<uint8_t>(row[i] - paethPredictor(a, b, c));
        }
        if (s < bestSum) { bestFilter = 4; }
    }

    return bestFilter;
}

inline void applyFilter(uint8_t filter, const uint8_t* row,
                         const uint8_t* prevRow,
                         uint8_t* out, int stride, int bpp) {
    switch (filter) {
    case 0:
        std::memcpy(out, row, stride);
        break;
    case 1:
        for (int i = 0; i < stride; ++i) {
            int a = (i >= bpp) ? row[i - bpp] : 0;
            out[i] = static_cast<uint8_t>(row[i] - a);
        }
        break;
    case 2:
        for (int i = 0; i < stride; ++i) {
            int b = prevRow ? prevRow[i] : 0;
            out[i] = static_cast<uint8_t>(row[i] - b);
        }
        break;
    case 3:
        for (int i = 0; i < stride; ++i) {
            int a = (i >= bpp) ? row[i - bpp] : 0;
            int b = prevRow ? prevRow[i] : 0;
            out[i] = static_cast<uint8_t>(row[i] - (a + b) / 2);
        }
        break;
    case 4:
        for (int i = 0; i < stride; ++i) {
            int a = (i >= bpp) ? row[i - bpp] : 0;
            int b = prevRow ? prevRow[i] : 0;
            int c = (prevRow && i >= bpp) ? prevRow[i - bpp] : 0;
            out[i] = static_cast<uint8_t>(row[i] - paethPredictor(a, b, c));
        }
        break;
    }
}

inline void unfilterRow(uint8_t filter, uint8_t* row,
                        const uint8_t* prevRow, int stride, int bpp) {
    switch (filter) {
    case 0: break;
    case 1:
        for (int i = bpp; i < stride; ++i)
            row[i] += row[i - bpp];
        break;
    case 2:
        if (prevRow)
            for (int i = 0; i < stride; ++i) row[i] += prevRow[i];
        break;
    case 3:
        for (int i = 0; i < stride; ++i) {
            int a = (i >= bpp) ? row[i - bpp] : 0;
            int b = prevRow ? prevRow[i] : 0;
            row[i] += static_cast<uint8_t>((a + b) / 2);
        }
        break;
    case 4:
        for (int i = 0; i < stride; ++i) {
            int a = (i >= bpp) ? row[i - bpp] : 0;
            int b = prevRow ? prevRow[i] : 0;
            int c = (prevRow && i >= bpp) ? prevRow[i - bpp] : 0;
            row[i] += paethPredictor(a, b, c);
        }
        break;
    }
}

// ========================================================================
// PNG チャンク書き出しヘルパー
// ========================================================================

inline void writeBE32(uint8_t* p, uint32_t v) {
    p[0] = (v >> 24) & 0xFF;
    p[1] = (v >> 16) & 0xFF;
    p[2] = (v >>  8) & 0xFF;
    p[3] =  v        & 0xFF;
}

inline uint32_t readBE32(const uint8_t* p) {
    return (uint32_t(p[0]) << 24) | (uint32_t(p[1]) << 16) |
           (uint32_t(p[2]) << 8)  |  uint32_t(p[3]);
}

inline void writeChunk(std::ofstream& ofs, const char* type,
                       const uint8_t* data, std::size_t len) {
    uint8_t hdr[4];
    writeBE32(hdr, static_cast<uint32_t>(len));
    ofs.write(reinterpret_cast<const char*>(hdr), 4);
    ofs.write(type, 4);
    if (len > 0)
        ofs.write(reinterpret_cast<const char*>(data), len);

    // CRC (type + data)
    std::vector<uint8_t> crcBuf(4 + len);
    std::memcpy(crcBuf.data(), type, 4);
    if (len > 0)
        std::memcpy(crcBuf.data() + 4, data, len);
    uint32_t crc = crc32(crcBuf.data(), crcBuf.size());
    writeBE32(hdr, crc);
    ofs.write(reinterpret_cast<const char*>(hdr), 4);
}

} // namespace detail_png

// ========================================================================
// savePng
// ========================================================================

template<typename T>
void savePng(const std::string& path, const Image<T>& img) {
    if (img.empty()) throw std::runtime_error("savePng: empty image");

    int w = static_cast<int>(img.width());
    int h = static_cast<int>(img.height());
    int ch = img.numChannels();

    // PNG カラータイプ
    uint8_t colorType;
    int bpp;  // bytes per pixel
    if (ch >= 4) {
        colorType = 6;  // RGBA
        bpp = 4;
    } else if (ch >= 3) {
        colorType = 2;  // RGB
        bpp = 3;
    } else {
        colorType = 0;  // Grayscale
        bpp = 1;
    }

    int stride = w * bpp;

    // ピクセルデータを uint8_t 行列に変換 + フィルタ適用
    std::vector<uint8_t> rawRows(h * stride);
    for (int y = 0; y < h; ++y) {
        uint8_t* row = rawRows.data() + y * stride;
        for (int x = 0; x < w; ++x) {
            if (bpp >= 3) {
                row[x * bpp + 0] = detail::toU8(img.at(y, x, 0));
                row[x * bpp + 1] = detail::toU8(img.at(y, x, 1));
                row[x * bpp + 2] = detail::toU8(img.at(y, x, 2));
                if (bpp == 4)
                    row[x * bpp + 3] = detail::toU8(img.at(y, x, 3));
            } else {
                row[x] = detail::toU8(img.at(y, x));
            }
        }
    }

    // フィルタ適用 + 結合
    std::vector<uint8_t> filtered;
    filtered.reserve(h * (stride + 1));
    std::vector<uint8_t> filteredRow(stride);

    for (int y = 0; y < h; ++y) {
        const uint8_t* row = rawRows.data() + y * stride;
        const uint8_t* prev = (y > 0) ? rawRows.data() + (y - 1) * stride : nullptr;

        uint8_t bestF = detail_png::chooseBestFilter(row, prev, stride, bpp);
        filtered.push_back(bestF);
        detail_png::applyFilter(bestF, row, prev, filteredRow.data(), stride, bpp);
        filtered.insert(filtered.end(), filteredRow.begin(), filteredRow.end());
    }

    // Deflate 圧縮
    auto deflated = detail_png::deflateCompress(filtered.data(), filtered.size());

    // zlib ラッパー (CMF + FLG + deflate + ADLER32)
    std::vector<uint8_t> zlibData;
    zlibData.reserve(deflated.size() + 6);
    zlibData.push_back(0x78);  // CMF: deflate, window=32768
    zlibData.push_back(0x9C);  // FLG: default compression
    zlibData.insert(zlibData.end(), deflated.begin(), deflated.end());

    uint32_t adler = detail_png::adler32(filtered.data(), filtered.size());
    uint8_t adlerBE[4];
    detail_png::writeBE32(adlerBE, adler);
    zlibData.insert(zlibData.end(), adlerBE, adlerBE + 4);

    // PNG ファイル書き出し
    std::ofstream ofs(path, std::ios::binary);
    if (!ofs) throw std::runtime_error("savePng: cannot open " + path);

    // シグネチャ
    static const uint8_t sig[8] = {137, 80, 78, 71, 13, 10, 26, 10};
    ofs.write(reinterpret_cast<const char*>(sig), 8);

    // IHDR
    uint8_t ihdr[13];
    detail_png::writeBE32(ihdr + 0, w);
    detail_png::writeBE32(ihdr + 4, h);
    ihdr[8] = 8;          // bit depth
    ihdr[9] = colorType;
    ihdr[10] = 0;         // compression method
    ihdr[11] = 0;         // filter method
    ihdr[12] = 0;         // interlace method
    detail_png::writeChunk(ofs, "IHDR", ihdr, 13);

    // IDAT
    detail_png::writeChunk(ofs, "IDAT", zlibData.data(), zlibData.size());

    // IEND
    detail_png::writeChunk(ofs, "IEND", nullptr, 0);
}

// ========================================================================
// loadPng
// ========================================================================

template<typename T>
Image<T> loadPng(const std::string& path) {
    std::ifstream ifs(path, std::ios::binary);
    if (!ifs) throw std::runtime_error("loadPng: cannot open " + path);

    // シグネチャ確認
    uint8_t sig[8];
    ifs.read(reinterpret_cast<char*>(sig), 8);
    static const uint8_t expected[8] = {137, 80, 78, 71, 13, 10, 26, 10};
    if (std::memcmp(sig, expected, 8) != 0)
        throw std::runtime_error("loadPng: not a PNG file");

    int w = 0, h = 0, bitDepth = 0, colorType = 0;
    std::vector<uint8_t> compressedData;

    // チャンクを読む
    while (ifs) {
        uint8_t chunkHdr[8];
        if (!ifs.read(reinterpret_cast<char*>(chunkHdr), 8)) break;

        uint32_t chunkLen = detail_png::readBE32(chunkHdr);
        char type[5] = {};
        std::memcpy(type, chunkHdr + 4, 4);

        std::vector<uint8_t> chunkData(chunkLen);
        if (chunkLen > 0)
            ifs.read(reinterpret_cast<char*>(chunkData.data()), chunkLen);

        uint8_t crcBuf[4];
        ifs.read(reinterpret_cast<char*>(crcBuf), 4);  // CRC (読み飛ばし)

        if (std::strcmp(type, "IHDR") == 0) {
            w = static_cast<int>(detail_png::readBE32(chunkData.data()));
            h = static_cast<int>(detail_png::readBE32(chunkData.data() + 4));
            bitDepth = chunkData[8];
            colorType = chunkData[9];
            if (bitDepth != 8)
                throw std::runtime_error("loadPng: only 8-bit depth supported");
        } else if (std::strcmp(type, "IDAT") == 0) {
            compressedData.insert(compressedData.end(),
                                  chunkData.begin(), chunkData.end());
        } else if (std::strcmp(type, "IEND") == 0) {
            break;
        }
    }

    if (w == 0 || h == 0)
        throw std::runtime_error("loadPng: missing IHDR");

    // zlib ヘッダをスキップ (2 バイト) して deflate 展開
    if (compressedData.size() < 6)
        throw std::runtime_error("loadPng: IDAT too short");

    auto decompressed = detail_png::deflateDecompress(
        compressedData.data() + 2,
        compressedData.size() - 6  // 2 header + 4 adler
    );

    // チャネル数
    int ch, bpp;
    switch (colorType) {
    case 0: ch = 1; bpp = 1; break;  // Grayscale
    case 2: ch = 3; bpp = 3; break;  // RGB
    case 4: ch = 2; bpp = 2; break;  // Grayscale + Alpha
    case 6: ch = 4; bpp = 4; break;  // RGBA
    default:
        throw std::runtime_error("loadPng: unsupported color type");
    }

    int stride = w * bpp;
    if (static_cast<int>(decompressed.size()) < h * (stride + 1))
        throw std::runtime_error("loadPng: decompressed data too short");

    // フィルタ解除 + Image 構築
    Image<T> img(h, w, ch);
    std::vector<uint8_t> curRow(stride), prevRow(stride, 0);

    for (int y = 0; y < h; ++y) {
        const uint8_t* src = decompressed.data() + y * (stride + 1);
        uint8_t filter = src[0];
        std::memcpy(curRow.data(), src + 1, stride);

        detail_png::unfilterRow(filter, curRow.data(),
                                y > 0 ? prevRow.data() : nullptr, stride, bpp);

        for (int x = 0; x < w; ++x) {
            for (int c = 0; c < ch; ++c)
                img.at(y, x, c) = detail::fromU8<T>(curRow[x * bpp + c]);
        }

        std::swap(curRow, prevRow);
    }

    return img;
}

} // namespace calx::image

#endif // CALX_IMAGE_IO_PNG_HPP
