title: Image convolution in spatial domain
body: <p>I am trying to replicate the outcome of <a href="https://stackoverflow.com/q/38709810/159072"><strong>this link</strong></a> using <strong>linear</strong> convolution in <strong>spatial-domain</strong>. </p>
  
  <p>Images are first converted to 2d <code>double</code> arrays and then convolved. Image and kernel are of the same size. The image is padded before convolution and cropped accordingly after the convolution. </p>
  
  <p><a href="https://i.stack.imgur.com/Nn5I8.png" rel="noreferrer"><img src="https://i.stack.imgur.com/Nn5I8.png" alt="enter image description here"></a></p>
  
  <p>As compared to the FFT-based convolution, the output is weird and <strong>incorrect</strong>. </p>
  
  <p><strong>How can I solve the issue?</strong></p>
  
  <p>Note that I obtained the following image output from Matlab which matches my C# FFT output:   </p>
  
  <p><a href="https://i.stack.imgur.com/7LltS.jpg" rel="noreferrer"><img src="https://i.stack.imgur.com/7LltS.jpg" alt="enter image description here"></a> </p>
  
  <p>.    </p>
  
  <p><sup><strong>Update-1:</strong> Following @<strong>Ben Voigt</strong>'s comment, I changed the <code>Rescale()</code> function to replace <code>255.0</code> with <code>1</code> and thus the output is improved substantially. But, still, the output doesn't match the FFT output (which is the correct one).</sup><br>
  <a href="https://i.stack.imgur.com/zRMTf.png" rel="noreferrer"><img src="https://i.stack.imgur.com/zRMTf.png" alt="enter image description here"></a>  </p>
  
  <p>.  </p>
  
  <p><sup><strong>Update-2:</strong> Following @<strong>Cris Luengo</strong>'s comment, I have padded the image by stitching and then performed spatial convolution. The outcome has been as follows:</sup><br>
  <a href="https://i.stack.imgur.com/wMGnQ.png" rel="noreferrer"><img src="https://i.stack.imgur.com/wMGnQ.png" alt="enter image description here"></a></p>
  
  <p><sup>So, the output is worse than the previous one. But, this has a similarity with the 2nd output of the <a href="https://stackoverflow.com/a/38724731/159072">linked answer</a> which means a circular convolution is not the solution.</sup>    </p>
  
  <p>.  </p>
  
  <p><sup><strong>Update-3</strong>: I have used the <code>Sum()</code> function proposed by @<strong>Cris Luengo</strong>'s answer. The result is a more improved version of <code>**Update-1**</code>:</sup><br>
  <a href="https://i.stack.imgur.com/aqneV.png" rel="noreferrer"><img src="https://i.stack.imgur.com/aqneV.png" alt="enter image description here"></a>   </p>
  
  <p><sup>But, it is still not 100% similar to the FFT version.</sup>    </p>
  
  <p>.  </p>
  
  <p><sup><strong>Update-4</strong>: Following @<strong>Cris Luengo</strong>'s comment, I have subtracted the two outcomes to see the difference:</sup><br>
  <a href="https://i.stack.imgur.com/o2zSz.png" rel="noreferrer"><img src="https://i.stack.imgur.com/o2zSz.png" alt="enter image description here"></a>, <a href="https://i.stack.imgur.com/u7XkX.png" rel="noreferrer"><img src="https://i.stack.imgur.com/u7XkX.png" alt="enter image description here"></a></p>
  
  <p><sup>
   1. spatial minus frequency domain<br>
   2. frequency minus spatial domain
  </sup>  </p>
  
  <p><sub>Looks like, the difference is substantial which means, spatial convolution is not being done correctly.</sub>    </p>
  
  <p>.  </p>
  
  <p><strong>Source Code:</strong>   </p>
  
  <p><sup>(Notify me if you need more source code to see.)</sup>      </p>
  
  <pre><code>    public static double[,] LinearConvolutionSpatial(double[,] image, double[,] mask)
      {
          int maskWidth = mask.GetLength(0);
          int maskHeight = mask.GetLength(1);
  
          double[,] paddedImage = ImagePadder.Pad(image, maskWidth);
  
          double[,] conv = Convolution.ConvolutionSpatial(paddedImage, mask);
  
          int cropSize = (maskWidth/2);
  
          double[,] cropped = ImageCropper.Crop(conv, cropSize);
  
          return conv;
      } 
      static double[,] ConvolutionSpatial(double[,] paddedImage1, double[,] mask1)
      {
          int imageWidth = paddedImage1.GetLength(0);
          int imageHeight = paddedImage1.GetLength(1);
  
          int maskWidth = mask1.GetLength(0);
          int maskHeight = mask1.GetLength(1);
  
          int convWidth = imageWidth - ((maskWidth / 2) * 2);
          int convHeight = imageHeight - ((maskHeight / 2) * 2);
  
          double[,] convolve = new double[convWidth, convHeight];
  
          for (int y = 0; y &lt; convHeight; y++)
          {
              for (int x = 0; x &lt; convWidth; x++)
              {
                  int startX = x;
                  int startY = y;
  
                  convolve[x, y] = Sum(paddedImage1, mask1, startX, startY);
              }
          }
  
          Rescale(convolve);
  
          return convolve;
      } 
  
      static double Sum(double[,] paddedImage1, double[,] mask1, int startX, int startY)
      {
          double sum = 0;
  
          int maskWidth = mask1.GetLength(0);
          int maskHeight = mask1.GetLength(1);
  
          for (int y = startY; y &lt; (startY + maskHeight); y++)
          {
              for (int x = startX; x &lt; (startX + maskWidth); x++)
              {
                  double img = paddedImage1[x, y];
                  double msk = mask1[x - startX, y - startY];
                  sum = sum + (img * msk);
              }
          }
  
          return sum;
      }
  
      static void Rescale(double[,] convolve)
      {
          int imageWidth = convolve.GetLength(0);
          int imageHeight = convolve.GetLength(1);
  
          double maxAmp = 0.0;
  
          for (int j = 0; j &lt; imageHeight; j++)
          {
              for (int i = 0; i &lt; imageWidth; i++)
              {
                  maxAmp = Math.Max(maxAmp, convolve[i, j]);
              }
          }
  
          double scale = 1.0 / maxAmp;
  
          for (int j = 0; j &lt; imageHeight; j++)
          {
              for (int i = 0; i &lt; imageWidth; i++)
              {
                  double d = convolve[i, j] * scale;
                  convolve[i, j] = d;
              }
          }
      } 
  
      public static Bitmap ConvolveInFrequencyDomain(Bitmap image1, Bitmap kernel1)
      {
          Bitmap outcome = null;
  
          Bitmap image = (Bitmap)image1.Clone();
          Bitmap kernel = (Bitmap)kernel1.Clone();
  
          //linear convolution: sum. 
          //circular convolution: max
          uint paddedWidth = Tools.ToNextPow2((uint)(image.Width + kernel.Width));
          uint paddedHeight = Tools.ToNextPow2((uint)(image.Height + kernel.Height));
  
          Bitmap paddedImage = ImagePadder.Pad(image, (int)paddedWidth, (int)paddedHeight);
          Bitmap paddedKernel = ImagePadder.Pad(kernel, (int)paddedWidth, (int)paddedHeight);
  
          Complex[,] cpxImage = ImageDataConverter.ToComplex(paddedImage);
          Complex[,] cpxKernel = ImageDataConverter.ToComplex(paddedKernel);
  
          // call the complex function
          Complex[,] convolve = Convolve(cpxImage, cpxKernel);
  
          outcome = ImageDataConverter.ToBitmap(convolve);
  
          outcome = ImageCropper.Crop(outcome, (kernel.Width/2)+1);
  
          return outcome;
      } 
  </code></pre>