3D DWT works for haar filter but not for Daubechies coefficients

Question

I am currently trying to implement a 3D Discrete Wavelet transform which works perfectly with haar coefficients but getting inaccurate values for any other non-symmetrical filter like db2 and db3.

I have tried different boundary handling techniques and working on each subband separately but kept on getting the same error. I am comparing my results with the 3D DWT feature on pywavlets while also passing my transform through the inverse transform to get my original data back but I am still getting an overall MSE of around 40K when comparing the inverse result with the original. Does anyone have an idea why ?

#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;
}

#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;
            }

        }
    }
}

First original 15 brain slices

Transformed 15 brain slices