3D离散小波变换Haar可用但Daubechies系数结果异常求助

我正在实现3D离散小波变换（3D DWT），使用Haar系数时功能正常，但使用db2、db3等非对称滤波器时结果不准确。尝试了多种边界处理方式，也单独处理各子带，但问题仍未解决。将实现结果与pywavelets的3D DWT功能对比，同时对变换结果执行逆变换以还原原始数据，但逆变换结果与原始数据的*均方误差（MSE）*约为40K。请问该问题的可能原因是什么？

核心实现代码 #

DWT类实现 #

#include "DWT.h"

// Constructor for the DWT class to be used for convolving the filters
DWT::DWT(const float* lpf, const float* hpf, size_t filter_size)
    : convolve(lpf, hpf, filter_size) {}

// Perform the 3D Discrete Wavelet Transform
Array3D<float> DWT::dwt_3d(const Array3D<float>& data, int levels) const {
    if (levels < 1) {
        throw invalid_argument("Levels must be greater than or equal to 1");
    }

    Array3D<float> result = data;

    size_t depth = data.get_depth();
    size_t rows = data.get_rows();
    size_t cols = data.get_cols();

    for (int level = 0; level < levels; ++level) {
        // Convolve and subsample ONLY within the bounds of the current level
        convolve.dim0(result, depth, rows, cols);
        convolve.dim1(result, depth, rows, cols);
        convolve.dim2(result, depth, rows, cols);

        // Calculate new bounds for the next level's LLL subband
        depth = (depth + 1) / 2;
        rows = (rows + 1) / 2;
        cols = (cols + 1) / 2;
    }

    return result;
}

Convolve类实现 #

#include "convolve.h"

// Constructor for the Convolve class
Convolve::Convolve(const float* lpf, const float* hpf, size_t filter_size)
    : lpf(lpf), hpf(hpf), filter_size(filter_size) {}

// Convolution along the first dimension (rows)
void Convolve::dim0(Array3D<float>& data, size_t depth_limit, size_t row_limit, size_t col_limit) const {
    Array3D<float> temp(data); // Create a temporary copy to avoid overrding the original data

    for (size_t d = 0; d < depth_limit; ++d) {
        for (size_t c = 0; c < col_limit; ++c) {

            for (size_t i = 0; i < row_limit / 2; ++i) {
                float sum_low = 0.0f;
                float sum_high = 0.0f;

                for (size_t j = 0; j < filter_size; ++j) {
                    size_t index = (2 * i + j) % row_limit;
                    float input_val = temp(d, index, c);
                    sum_low += lpf[j] * input_val;
                    sum_high += hpf[j] * input_val;
                }

                data(d, i, c) = sum_low;
                data(d, i + row_limit / 2, c) = sum_high;
            }

         }
    }
}

// Convolution along the second dimension (columns)
void Convolve::dim1(Array3D<float>& data, size_t depth_limit, size_t row_limit, size_t col_limit) const {
    Array3D<float> temp(data); // Create a temporary copy to avoid overrding the original data

    for (size_t d = 0; d < depth_limit; ++d) {
        for (size_t r = 0; r < row_limit; ++r) {

            for (size_t i = 0; i < col_limit / 2; ++i) {
                float sum_low = 0.0f;
                float sum_high = 0.0f;

                for (size_t j = 0; j < filter_size; ++j) {
                    size_t index = (2 * i + j) % col_limit;
                    float input_val = temp(d, r, index);
                    sum_low += lpf[j] * input_val;
                    sum_high += hpf[j] * input_val;
                }

                data(d, r, i) = sum_low;
                data(d, r, i + col_limit / 2) = sum_high;
            }

        }
    }
}

// Convolution along the third dimension (depths)
void Convolve::dim2(Array3D<float>& data, size_t depth_limit, size_t row_limit, size_t col_limit) const {
    Array3D<float> temp(data); // Create a temporary copy to avoid overrding the original data

    for (size_t r = 0; r < row_limit; ++r) {
        for (size_t c = 0; c < col_limit; ++c) { 

            for (size_t i = 0; i < depth_limit / 2; ++i) {
                float sum_low = 0.0f;
                float sum_high = 0.0f;

                for (size_t j = 0; j < filter_size; ++j) {
                    size_t index = (2 * i + j) % depth_limit;
                    float input_val = temp(index, r, c);
                    sum_low += lpf[j] * input_val;
                    sum_high += hpf[j] * input_val;
                }

                data(i, r, c) = sum_low;
                data(i + depth_limit / 2, r, c) = sum_high;
            }

        }
    }
}

数据示例 #

• 原始15层脑切片：
• 变换后的15层脑切片：

内容的提问来源于stack exchange，提问作者Antonio Galdes