hoNDArray Incorrect Mean Calculation

gadgetron / gadgetron

Gadgetron - Medical Image Reconstruction Framework

Other

237 stars 164 forks source link

Describe the bug As array size and stored value increase, the complex hoNDArray functions for mean, stddev, etc return incorrect values

To Reproduce hoNDArray<std::complex> tens(10,10,10,10,10,10,10); tens.fill(std::complex(10,10)); auto mean = Gadgetron::mean(tens); GDEBUG("Mean: %f +%fi", mean.real(), mean.imag()); Output: "Mean: 9.342177 +9.342177i"

hoNDArray<std::complex> tens(10,10,10,10,10,10,10); tens.fill(std::complex(10,10)); auto stddev = Gadgetron::stddev(tens); GDEBUG("StdDev: %f +%fi", stddev.real(), stddev.imag()); Output: "StdDev: 0.369728 +0.000000i"

Expected behavior Array operators should return correct values regardless of array size or contained values. The following returns correctly: hoNDArray<std::complex> ones(10,10,10,10,10,10,10); ones.fill(std::complex(1,1)); auto mean = Gadgetron::mean(ones); GDEBUG("Mean: %f +%fi", mean.real(), mean.imag()); Output: "Mean: 1.000000 + 1.000000i"

hoNDArray<std::complex> tens(5,5,5,5,5,5,5); tens.fill(std::complex(10,10)); auto mean = Gadgetron::mean(tens); GDEBUG("Mean: %f +%fi", mean.real(), mean.imag()); Output: "Mean: 10.00000 + 10.00000i"

System Information:

OS: Ubuntu 20.04
Install/deployment method: Docker dev container
Version: 555e98d

Additional context I think there's a buffer overrun or we exceed a max value for a fixed-point calculation somewhere, but the armadillo methods underlying hoNDArray math operations are a bit difficult to tease apart.

Adding a few examples to show the error: Code:

    GDEBUG("%d,%d,%d,%d,%d,%d,%d,%d \n", REP, E0, E1, E2, CHA, N, S, LOC);
    hoNDArray<std::complex<float>> standardDeviations(E0, E1, E2, CHA, N, S, LOC);
    int numStdDevs = 0;
    std::complex<float> sumStdDevs = 0;

    for (size_t e0 = 0; e0 < E0; e0++) {
        for (size_t e1 = 0; e1 < E1; e1++) {
            for (size_t e2 = 0; e2 < E2; e2++) {
                for (size_t cha = 0; cha < CHA; cha++) {
                    for (size_t n = 0; n < N; n++) {
                        for (size_t s = 0; s < S; s++) {
                            for (size_t loc = 0; loc < LOC; loc++) {
                                hoNDArray<std::complex<float>> replicaData(REP);
                                for (size_t rep = 0; rep < REP; rep++) {
                                    replicaData[rep] = pseudoreplicaList[rep].rbit_[0].data_.data_[e0, e1, e2, cha, n, s, loc];
                                }
                                //standardDeviations[e0, e1, e2, cha, n, s, loc] = stddev(replicaData);
                                standardDeviations[e0, e1, e2, cha, n, s, loc] = std::complex<float>(1.0,1.0);
                                sumStdDevs += std::complex<float>(1.0,1.0);
                                numStdDevs++;
                            }
                        }
                    }
                }
            }
        }
    }
    GDEBUG("Count: %d\n", numStdDevs);
    GDEBUG("Sum: %f + i%f\n", sumStdDevs.real(), sumStdDevs.imag());
    GDEBUG("Manual Mean: %f + i%f\n", sumStdDevs.real()/float(numStdDevs), sumStdDevs.imag()/float(numStdDevs));

    auto meanOfStdDevs = mean(standardDeviations);
    GDEBUG("Built-in Mean: %f + i%f\n", meanOfStdDevs.real(), meanOfStdDevs.imag());

Outputs:

01-20 14:25:11.530 DEBUG [StandardDeviationGadget.cpp:18] ReconDataList Length: 6
01-20 14:25:11.530 DEBUG [StandardDeviationGadget.cpp:28] 6,128,128,1,16,10,1,1 
01-20 14:25:12.093 DEBUG [StandardDeviationGadget.cpp:56] Count: 2621440
01-20 14:25:12.093 DEBUG [StandardDeviationGadget.cpp:57] Sum: 2621440.000000 + i2621440.000000
01-20 14:25:12.093 DEBUG [StandardDeviationGadget.cpp:58] Manual Mean: 1.000000 + i1.000000
01-20 14:25:12.096 DEBUG [StandardDeviationGadget.cpp:61] Built-in Mean: 0.000000 + i0.000000

gadgetron / gadgetron

hoNDArray Incorrect Mean Calculation #1189