Sum neighborhood of every matrix element efficiently

Question

I programmed a piece of code where I have a huge 3D matrix in C++ using boost::multi_array.

For every matrix element I want to sum up the neighborhood in a certain distance dist. Every element is weighted according to its distance to the center. The center element should not be included in the sum.

The distance dist is given by the user and can vary. My program is doing the right calculations but is slow when the matrix gets big. I have sometimes matrices with more than 100000 elements...

So my question is: Is there any way to make this computation faster? Maybe also by using another library?

The part consists mainly of two functions. In the first function I access every matrix element and calculate the sum of the neighborhood. The inputMatrix is a 3D boost multi array:

boost::multi_array<float, 3> inputMatrix = imageMatrix;
T actualElement;
int posActualElement;
for (int depth = 0; depth<inputMatrix.shape()[2]; depth++) {
    for (int row = 0; row<inputMatrix.shape()[0]; row++) {
        for (int col = 0; col<inputMatrix.shape()[1]; col++) {

            indexOfElement[0] = row;
            indexOfElement[1] = col;
            indexOfElement[2] = depth;
            //get actual Element if it is the centre of a whole neighborhood
            actualElement = inputMatrix[row][col][depth];
            if (!std::isnan(actualElement)) {
                //get the sum of the actual neighborhood
                sumOfActualNeighborhood = getNeighborhood3D(inputMatrix, indexOfElement);
      }
    }
  }
}

The function neighborhood3D looks as follows:

template <class T, size_t R>
T NGTDMFeatures3D<T, R>::getNeighborhood3D(boost::multi_array<T, R> inputMatrix, int *indexOfElement) {
    std::vector<T> neighborhood;
    T actElement;
    float weight;
    for (int k = -dist; k<dist + 1; k++) {
        for (int i = -dist; i<dist + 1; i++) {
            for (int j = -dist; j<dist + 1; j++) {

                if (i != 0 || j != 0 || k != 0) {
                    if (indexOfElement[0] + i>-1 && indexOfElement[0] + i<inputMatrix.shape()[0] && indexOfElement[1] + j>-1 && indexOfElement[1] + j<inputMatrix.shape()[1] && indexOfElement[2] + k>-1 && indexOfElement[2] + k<inputMatrix.shape()[2]) {
                        actElement = inputMatrix[indexOfElement[0] + i][indexOfElement[1] + j][indexOfElement[2] + k];
                        if (!std::isnan(actElement)) {
                            weight = calculateWeight3D(i, j, k, normNGTDM, actualSpacing);
                            neighborhood.push_back(weight*actElement);
                        }
                    }
                }
            }
        }
    }
    T sum = accumulate(neighborhood.begin(), neighborhood.end(), 0);
    sum = sum / neighborhood.size();
    return sum;
}