# 3D Direct Convolution Implementation in C

For my project, I've written a naive C implementation of Direct3D convolution with periodic padding on the input. Unfortunately, since I'm new to C, the performance isn't so good.

By convention, all the matrices (image, kernel, result) are stored in column-major fashion. This is why I loop through them in such way so they are closer in memory. As I heard this would help.

I know the implementation is very naive. But since it's written in C, I was hoping the performance would be good, but instead it's a little disappointing. I tested it with image of size $$\100^3\$$ and kernel of size $$\10^3\$$ which totals ~1GFLOPS if you only count the multiplication and addition. This took ~7s which I believe is way below the capability of a typical CPU.

Could the performance be optimized in this routine? I'm open to anything that could help, with just a few things if you could consider:

1. The problem I'm working on can be big. Image's size can be 200 by 200 by 200 whilst the kernel's size can be 50 by 50 by 50 or even larger. I understand that one way of optimizing this is by converting this problem into a matrix multiplication problem and use the blas GEMM routine, but I'm afraid memory could not hold such a big matrix

2. Due to the nature of the problem I would prefer direct convolution instead of FFTConvolve, since my model is developed with direct convolution in mind. My impression of FFTconvolve is that it gives slightly different result than direct convolve, especially for rapidly changing image. A discrepancy I'm trying to avoid. That said, I'm in no way an expert in this. so if you have a great implementation based on FFTconvolve and / or my impression on FFTconvolve is totally biased, I would really appreciate if you could help me out.

3. The input images are assumed to be periodic, so periodic padding is necessary

4. I understand that utilizing blas / SIMD or other lower level ways would definitely help a lot here. but since I'm a newbie here I don't really know where to start. I would really appreciate if you help pointing me to the right direction if you have experience in these libraries,

int mod(int a, int b)
{
// calculate mod to get the correct index with periodic padding
int r = a % b;
return r < 0 ? r + b : r;
}
void convolve3D(const double *image, const double *kernel, const int imageDimX, const int imageDimY, const int imageDimZ, const int kernelDimX, const int kernelDimY, const int kernelDimZ, double *result)
{
int imageSize = imageDimX * imageDimY * imageDimZ;
int kernelSize = kernelDimX * kernelDimY * kernelDimZ;

int i, j, k, l, m, n;
int kernelCenterX = (kernelDimX - 1) / 2;
int kernelCenterY = (kernelDimY - 1) / 2;
int kernelCenterZ = (kernelDimZ - 1) / 2;
int xShift,yShift,zShift;
int outIndex, outI, outJ, outK;
int imageIndex = 0, kernelIndex = 0;

// Loop through each voxel
for (k = 0; k < imageDimZ; k++){
for ( j = 0; j < imageDimY; j++) {
for ( i = 0; i < imageDimX; i++) {
kernelIndex = 0;
// for each voxel, loop through each kernel coefficient
for (n = 0; n < kernelDimZ; n++){
for ( m = 0; m < kernelDimY; m++) {
for ( l = 0; l < kernelDimX; l++) {
// find the index of the corresponding voxel in the output image
xShift = l - kernelCenterX;
yShift = m - kernelCenterY;
zShift = n - kernelCenterZ;

outI = mod ((i - xShift), imageDimX);
outJ = mod ((j - yShift), imageDimY);
outK = mod ((k - zShift), imageDimZ);

outIndex = outK * imageDimX * imageDimY + outJ * imageDimX + outI;

result[outIndex] += kernel[kernelIndex]* image[imageIndex];
kernelIndex++;
}
}
}
imageIndex ++;
}
}
}
}

• as discussed on the Stack Overflow version of this question, % is very slow with divisors that aren't compile-time constants. Replacing that with a conditional add (assuming it can't wrap more than once) sped up the whole thing by 20%. (Compiling with -O3 with an unknown compiler on unknown hardware, presumably not -march=native or -ffast-math.) – Peter Cordes Jun 28 '20 at 0:56
• Does this compile? kernelIndex is declared, and later incremented, but never used anywhere. stencil and stencilIndex is only referenced, never declared or defined. Are they supposed to be kernel and kernelIndex? Requests for reviews of broken code are off-topic here. – scottbb Jun 28 '20 at 1:14
• @scottbb ah, yes, I changed it from stencil to kernel when copying my code, and missed some... I have corrected the typo – lxiangyun93 Jun 29 '20 at 4:07
• Some "little"... for (n = 0, zShift = - kernelCenterZ; n < kernelDimZ; n++, zShift++){ and the next 2 loops also. You can remove the 3 substractions in the inner body. They were calculated too often. – Holger Jun 29 '20 at 9:26
• outK and outJ don't need to be allways recalculated in the inner loop. – Holger Jun 29 '20 at 9:33