function [compressed_data, reconstructed_img, compression_info] = dct_color_compress(img, quality_factor)
    % 基於DCT的彩色圖像壓縮
    % 輸入:
    %   img - 輸入彩色圖像
    %   quality_factor - 質量因子 (1-100)
    % 輸出:
    %   compressed_data - 壓縮後的數據
    %   reconstructed_img - 重建圖像
    %   compression_info - 壓縮信息
    
    % 參數驗證
    if nargin < 2
        quality_factor = 50;
    end
    quality_factor = max(1, min(100, quality_factor));
    
    % 轉換為double類型
    original_img = im2double(img);
    
    % 顏色空間轉換: RGB -> YCbCr
    ycbcr_img = rgb2ycbcr(original_img);
    
    % 提取Y, Cb, Cr分量
    Y = ycbcr_img(:, :, 1);
    Cb = ycbcr_img(:, :, 2);
    Cr = ycbcr_img(:, :, 3);
    
    % 色度下采樣 (4:2:0)
    [Cb_subsampled, Cr_subsampled] = chroma_subsample(Cb, Cr);
    
    % 對每個分量進行DCT壓縮
    fprintf('正在壓縮Y分量...\n');
    [compressed_Y, info_Y] = dct_compress_channel(Y, quality_factor, 'luminance');
    
    fprintf('正在壓縮Cb分量...\n');
    [compressed_Cb, info_Cb] = dct_compress_channel(Cb_subsampled, quality_factor, 'chrominance');
    
    fprintf('正在壓縮Cr分量...\n');
    [compressed_Cr, info_Cr] = dct_compress_channel(Cr_subsampled, quality_factor, 'chrominance');
    
    % 組合壓縮數據
    compressed_data.Y = compressed_Y;
    compressed_data.Cb = compressed_Cb;
    compressed_data.Cr = compressed_Cr;
    compressed_data.original_size = size(img);
    compressed_data.quality_factor = quality_factor;
    
    % 重建圖像
    fprintf('重建圖像...\n');
    reconstructed_img = dct_color_decompress(compressed_data);
    
    % 計算壓縮信息
    compression_info = calculate_compression_info(original_img, compressed_data, info_Y, info_Cb, info_Cr);
    
    fprintf('壓縮完成!\n');
end

function [Cb_subsampled, Cr_subsampled] = chroma_subsample(Cb, Cr)
    % 色度下采樣 (4:2:0)
    % 在水平和垂直方向都進行2:1下采樣
    
    % 使用抗混疊濾波
    h = fspecial('gaussian', [3 3], 0.5);
    Cb_filtered = imfilter(Cb, h, 'replicate');
    Cr_filtered = imfilter(Cr, h, 'replicate');
    
    % 下采樣
    Cb_subsampled = Cb_filtered(1:2:end, 1:2:end);
    Cr_subsampled = Cr_filtered(1:2:end, 1:2:end);
end

function [Cb_upsampled, Cr_upsampled] = chroma_upsample(Cb_subsampled, Cr_subsampled, target_size)
    % 色度上採樣
    % 使用雙線性插值恢復原始尺寸
    
    Cb_upsampled = imresize(Cb_subsampled, target_size, 'bilinear');
    Cr_upsampled = imresize(Cr_subsampled, target_size, 'bilinear');
end

function [compressed_channel, info] = dct_compress_channel(channel, quality_factor, channel_type)
    % 對單個通道進行DCT壓縮
    % 輸入:
    %   channel - 輸入通道數據
    %   quality_factor - 質量因子
    %   channel_type - 'luminance' 或 'chrominance'
    
    [height, width] = size(channel);
    block_size = 8;
    
    % 確保圖像尺寸是8的倍數
    padded_height = ceil(height / block_size) * block_size;
    padded_width = ceil(width / block_size) * block_size;
    channel_padded = padarray(channel, [padded_height-height, padded_width-width], 'replicate', 'post');
    
    % 獲取量化表
    quant_table = get_quantization_table(quality_factor, channel_type);
    
    % 分塊處理
    compressed_blocks = cell(padded_height/block_size, padded_width/block_size);
    dc_coefficients = [];
    ac_coefficients = [];
    
    for i = 1:block_size:padded_height
        for j = 1:block_size:padded_width
            % 提取8x8塊
            block = channel_padded(i:i+block_size-1, j:j+block_size-1);
            
            % 電平偏移 (0-255 -> -128 to 127)
            block_shifted = block - 128;
            
            % DCT變換
            dct_coeffs = dct2(block_shifted);
            
            % 量化
            quantized_coeffs = round(dct_coeffs ./ quant_table);
            
            % 保存壓縮後的塊
            compressed_blocks{(i-1)/block_size+1, (j-1)/block_size+1} = quantized_coeffs;
            
            % 提取DC和AC係數用於熵編碼分析
            dc_coefficients = [dc_coefficients; quantized_coeffs(1,1)];
            ac_block = quantized_coeffs(:);
            ac_block(1) = []; % 移除DC係數
            ac_coefficients = [ac_coefficients; ac_block];
        end
    end
    
    % 組織壓縮數據
    compressed_channel.blocks = compressed_blocks;
    compressed_channel.quant_table = quant_table;
    compressed_channel.original_size = size(channel);
    compressed_channel.padded_size = [padded_height, padded_width];
    compressed_channel.block_size = block_size;
    
    % 統計信息
    info.dc_coefficients = dc_coefficients;
    info.ac_coefficients = ac_coefficients;
    info.non_zero_ac = sum(ac_coefficients ~= 0);
    info.total_ac = length(ac_coefficients);
end

function quant_table = get_quantization_table(quality_factor, channel_type)
    % 生成量化表
    % 基於JPEG標準量化表
    
    % 標準亮度量化表
    luminance_quant = [16  11  10  16  24  40  51  61;
                       12  12  14  19  26  58  60  55;
                       14  13  16  24  40  57  69  56;
                       14  17  22  29  51  87  80  62;
                       18  22  37  56  68  109 103 77;
                       24  35  55  64  81  104 113 92;
                       49  64  78  87  103 121 120 101;
                       72  92  95  98  112 100 103 99];
    
    % 標準色度量化表
    chrominance_quant = [17  18  24  47  99  99  99  99;
                         18  21  26  66  99  99  99  99;
                         24  26  56  99  99  99  99  99;
                         47  66  99  99  99  99  99  99;
                         99  99  99  99  99  99  99  99;
                         99  99  99  99  99  99  99  99;
                         99  99  99  99  99  99  99  99;
                         99  99  99  99  99  99  99  99];
    
    % 選擇量化表
    if strcmpi(channel_type, 'luminance')
        base_quant = luminance_quant;
    else
        base_quant = chrominance_quant;
    end
    
    % 根據質量因子調整量化表
    if quality_factor < 50
        scale_factor = 50 / quality_factor;
    else
        scale_factor = (100 - quality_factor) / 50;
    end
    
    quant_table = max(1, min(255, round(base_quant * scale_factor)));
end

function reconstructed_img = dct_color_decompress(compressed_data)
    % DCT彩色圖像解壓縮
    
    % 提取壓縮數據
    compressed_Y = compressed_data.Y;
    compressed_Cb = compressed_data.Cb;
    compressed_Cr = compressed_data.Cr;
    original_size = compressed_data.original_size;
    
    % 解壓縮各分量
    fprintf('解壓縮Y分量...\n');
    Y_reconstructed = dct_decompress_channel(compressed_Y);
    
    fprintf('解壓縮Cb分量...\n');
    Cb_subsampled = dct_decompress_channel(compressed_Cb);
    
    fprintf('解壓縮Cr分量...\n');
    Cr_subsampled = dct_decompress_channel(compressed_Cr);
    
    % 色度上採樣
    [Cb_upsampled, Cr_upsampled] = chroma_upsample(Cb_subsampled, Cr_subsampled, original_size(1:2));
    
    % 組合YCbCr圖像
    ycbcr_reconstructed = cat(3, Y_reconstructed, Cb_upsampled, Cr_upsampled);
    
    % 轉換回RGB
    reconstructed_img = ycbcr2rgb(ycbcr_reconstructed);
    
    % 裁剪到原始尺寸
    reconstructed_img = reconstructed_img(1:original_size(1), 1:original_size(2), :);
end

function channel_reconstructed = dct_decompress_channel(compressed_channel)
    % 單通道DCT解壓縮
    
    blocks = compressed_channel.blocks;
    quant_table = compressed_channel.quant_table;
    block_size = compressed_channel.block_size;
    padded_size = compressed_channel.padded_size;
    
    % 重建圖像
    channel_reconstructed = zeros(padded_size);
    
    for i = 1:block_size:padded_size(1)
        for j = 1:block_size:padded_size(2)
            % 獲取量化後的係數
            quantized_coeffs = blocks{(i-1)/block_size+1, (j-1)/block_size+1};
            
            % 反量化
            dct_coeffs = quantized_coeffs .* quant_table;
            
            % 逆DCT變換
            block_reconstructed = idct2(dct_coeffs);
            
            % 電平偏移恢復
            block_reconstructed = block_reconstructed + 128;
            
            % 放回重建圖像
            channel_reconstructed(i:i+block_size-1, j:j+block_size-1) = block_reconstructed;
        end
    end
    
    % 裁剪到原始尺寸
    original_size = compressed_channel.original_size;
    channel_reconstructed = channel_reconstructed(1:original_size(1), 1:original_size(2));
end

function compression_info = calculate_compression_info(original_img, compressed_data, info_Y, info_Cb, info_Cr)
    % 計算壓縮性能指標
    
    % 原始數據量 (字節)
    original_size_bytes = numel(original_img) * 8; % 假設8位/像素
    
    % 估算壓縮後數據量
    % 這裏簡化計算，實際應該包括量化表、霍夫曼表等開銷
    compressed_size_estimate = 0;
    
    % Y分量
    num_blocks_Y = numel(compressed_data.Y.blocks);
    non_zero_coeffs_Y = info_Y.non_zero_ac + num_blocks_Y; % 包括DC係數
    
    % Cb分量
    num_blocks_Cb = numel(compressed_data.Cb.blocks);
    non_zero_coeffs_Cb = info_Cb.non_zero_ac + num_blocks_Cb;
    
    % Cr分量
    num_blocks_Cr = numel(compressed_data.Cr.blocks);
    non_zero_coeffs_Cr = info_Cr.non_zero_ac + num_blocks_Cr;
    
    % 假設每個非零係數平均需要12位 (包括遊程和值)
    compressed_size_estimate = (non_zero_coeffs_Y + non_zero_coeffs_Cb + non_zero_coeffs_Cr) * 1.5;
    
    % 壓縮比
    compression_ratio = original_size_bytes / compressed_size_estimate;
    
    % 保存壓縮信息
    compression_info.original_size_bytes = original_size_bytes;
    compression_info.compressed_size_estimate = compressed_size_estimate;
    compression_info.compression_ratio = compression_ratio;
    compression_info.bit_rate = compressed_size_estimate * 8 / (numel(original_img) / 3); % bpp
    
    % 係數統計
    compression_info.non_zero_coeffs.total = non_zero_coeffs_Y + non_zero_coeffs_Cb + non_zero_coeffs_Cr;
    compression_info.non_zero_coeffs.Y = non_zero_coeffs_Y;
    compression_info.non_zero_coeffs.Cb = non_zero_coeffs_Cb;
    compression_info.non_zero_coeffs.Cr = non_zero_coeffs_Cr;
    
    fprintf('壓縮性能分析:\n');
    fprintf('原始數據量: %.2f KB\n', original_size_bytes/1024);
    fprintf('壓縮後估算: %.2f KB\n', compressed_size_estimate/1024);
    fprintf('壓縮比: %.2f:1\n', compression_ratio);
    fprintf('比特率: %.2f bpp\n', compression_info.bit_rate);
end

function analyze_dct_coefficients(compressed_data)
    % 分析DCT係數分佈
    
    figure('Position', [100, 100, 1200, 800]);
    
    % Y分量係數分析
    subplot(2,3,1);
    dc_y = [];
    for i = 1:numel(compressed_data.Y.blocks)
        block = compressed_data.Y.blocks{i};
        dc_y = [dc_y; block(1,1)];
    end
    hist(dc_y, 50);
    title('Y分量DC係數分佈');
    xlabel('DC係數值');
    ylabel('頻數');
    
    subplot(2,3,2);
    ac_y = [];
    for i = 1:numel(compressed_data.Y.blocks)
        block = compressed_data.Y.blocks{i};
        ac_block = block(:);
        ac_block(1) = []; % 移除DC
        ac_y = [ac_y; ac_block];
    end
    hist(ac_y, 100);
    title('Y分量AC係數分佈');
    xlabel('AC係數值');
    ylabel('頻數');
    
    % 顯示量化表
    subplot(2,3,3);
    imagesc(compressed_data.Y.quant_table);
    colorbar;
    title('Y分量量化表');
    axis image;
    
    % 顯示一個示例塊的DCT係數
    subplot(2,3,4);
    example_block = compressed_data.Y.blocks{1,1};
    imagesc(log(abs(example_block) + 1));
    colorbar;
    title('示例塊DCT係數 (對數尺度)');
    axis image;
    
    % 零係數比例
    subplot(2,3,5);
    zero_ratio_y = sum(ac_y == 0) / length(ac_y);
    zero_ratio_cb = sum([compressed_data.Cb.blocks{:}](:) == 0) / numel([compressed_data.Cb.blocks{:}]);
    zero_ratio_cr = sum([compressed_data.Cr.blocks{:}](:) == 0) / numel([compressed_data.Cr.blocks{:}]);
    
    bar([zero_ratio_y, zero_ratio_cb, zero_ratio_cr]);
    set(gca, 'XTickLabel', {'Y', 'Cb', 'Cr'});
    title('零係數比例');
    ylabel('比例');
    
    % 係數能量分佈
    subplot(2,3,6);
    energy_y = accum_dct_energy(compressed_data.Y.blocks);
    plot(1:64, energy_y, 'b-', 'LineWidth', 2);
    title('DCT係數累積能量分佈');
    xlabel('Zig-zag掃描順序');
    ylabel('累積能量比例');
    grid on;
end

function energy_accum = accum_dct_energy(blocks)
    % 計算DCT係數累積能量分佈
    
    total_energy = 0;
    energy_by_position = zeros(64, 1);
    
    for i = 1:numel(blocks)
        block = blocks{i};
        zigzag_coeffs = zigzag_scan(block);
        energy = zigzag_coeffs.^2;
        total_energy = total_energy + sum(energy);
        energy_by_position = energy_by_position + energy;
    end
    
    % 按能量排序
    [~, sort_idx] = sort(energy_by_position, 'descend');
    
    % 計算累積能量
    energy_accum = zeros(64, 1);
    for i = 1:64
        energy_accum(i) = sum(energy_by_position(sort_idx(1:i))) / total_energy;
    end
end

function zigzag = zigzag_scan(block)
    % Zig-zag掃描
    zigzag = zeros(64, 1);
    index = 1;
    
    for sum = 0:2:14
        if rem(sum, 2) == 0
            for i = max(1, sum-7):min(8, sum+1)
                j = sum - i + 2;
                if j >= 1 && j <= 8
                    zigzag(index) = block(i, j);
                    index = index + 1;
                end
            end
        else
            for i = max(1, sum-7):min(8, sum+1)
                j = sum - i + 2;
                if j >= 1 && j <= 8
                    zigzag(index) = block(j, i);
                    index = index + 1;
                end
            end
        end
    end
end

function dct_compression_demo()
    % DCT壓縮演示函數
    
    % 讀取測試圖像
    try
        img = imread('peppers.png');
    catch
        % 如果找不到圖像，創建一個測試圖像
        fprintf('創建測試圖像...\n');
        img = uint8(255 * rand(256, 256, 3));
    end
    
    % 測試不同質量因子
    quality_factors = [10, 25, 50, 75, 90];
    
    figure('Position', [50, 50, 1400, 900]);
    
    for i = 1:length(quality_factors)
        qf = quality_factors(i);
        
        % 壓縮和解壓縮
        [compressed_data, reconstructed_img, compression_info] = dct_color_compress(img, qf);
        
        % 計算PSNR
        psnr_val = psnr(im2double(reconstructed_img), im2double(img));
        
        % 顯示結果
        subplot(2, 3, i);
        imshow(reconstructed_img);
        title(sprintf('質量因子: %d\nPSNR: %.2f dB, 壓縮比: %.2f:1', ...
            qf, psnr_val, compression_info.compression_ratio));
        
        fprintf('質量因子 %d: PSNR = %.2f dB, 壓縮比 = %.2f:1\n', ...
            qf, psnr_val, compression_info.compression_ratio);
    end
    
    % 顯示原始圖像
    subplot(2, 3, 6);
    imshow(img);
    title('原始圖像');
    
    % 詳細分析中等質量的情況
    fprintf('\n進行詳細分析...\n');
    [compressed_data, reconstructed_img, compression_info] = dct_color_compress(img, 50);
    
    % 分析DCT係數
    figure('Position', [100, 100, 1200, 800]);
    analyze_dct_coefficients(compressed_data);
    
    % 顯示壓縮前後對比
    figure('Position', [150, 150, 1000, 400]);
    subplot(1,2,1);
    imshow(img);
    title('原始圖像');
    
    subplot(1,2,2);
    imshow(reconstructed_img);
    title('壓縮後圖像 (質量因子50)');
    
    % 計算並顯示差異
    diff_img = im2double(img) - reconstructed_img;
    diff_img = (diff_img - min(diff_img(:))) / (max(diff_img(:)) - min(diff_img(:)));
    
    figure('Position', [200, 200, 600, 400]);
    imshow(diff_img);
    title('壓縮誤差 (歸一化顯示)');
    colorbar;
end

function adaptive_quant_table = get_adaptive_quantization(block, base_quant_table, activity)
    % 基於塊活動性的自適應量化
    % activity: 塊的活動性度量 (如方差)
    
    if activity > 0.1
        % 高活動性區域，使用較細的量化
        scale_factor = 0.7;
    else
        % 平坦區域，使用較粗的量化
        scale_factor = 1.3;
    end
    
    adaptive_quant_table = max(1, round(base_quant_table * scale_factor));
end

function progressive_data = progressive_compress(img, quality_levels)
    % 漸進式DCT壓縮
    % quality_levels: 質量級別數組，如 [10, 30, 50, 70, 90]
    
    progressive_data = cell(length(quality_levels), 1);
    
    for i = 1:length(quality_levels)
        fprintf('壓縮級別 %d/%d (質量因子: %d)...\n', i, length(quality_levels), quality_levels(i));
        [compressed_data, ~, compression_info] = dct_color_compress(img, quality_levels(i));
        progressive_data{i} = compressed_data;
    end
end