Свертка изображений в пространственной области

Я пытаюсь реплицировать результат этой ссылки, используя линейную свертку в пространственной области.

Сначала изображения преобразуются в 2d double массивы и затем свертываются. Изображение и ядро имеют одинаковый размер. Изображение дополняется до свертки и обрезается соответственно после свертки.

enter image description here

По сравнению со сверткой на основе FFT, выход является странным и неправильным.

Как я могу решить проблему?

Обратите внимание, что я получил следующий результат изображения Matlab, который соответствует моему выходу С# FFT:

enter image description here

,

Обновление-1: после комментария @Ben Voigt я изменил Rescale() чтобы заменить 255.0 на 1 и, следовательно, результат значительно улучшился.Но, тем не менее, выход не соответствует выходу FFT (что является правильным).
enter image description here

,

Обновление-2: После комментария @Cris Luengo, я заполнил изображение путем сшивания, а затем выполнил пространственную свертку.Результатом стало следующее:
enter image description here

Таким образом, результат хуже предыдущего.Но это имеет сходство со вторым результатом связанного ответа, что означает, что круговая свертка не является решением.

,

Обновление-3: я использовал функцию Sum() предложенную ответом @Cris Luengo.Результатом является более совершенная версия **Update-1** :
enter image description here

Но он по-прежнему не на 100% похож на версию FFT.

,

Обновление-4: После комментария @Cris Luengo, я вычитал два результата, чтобы увидеть разницу:
enter image description here, enter image description here

1. пространственная минусовая частотная область
2. частота минус пространственная область

Похоже, разница существенна, что означает, что пространственная свертка выполняется неправильно.

,

Исходный код:

(Сообщите мне, если вам нужно больше исходного кода для просмотра.)

    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 < convHeight; y++)
        {
            for (int x = 0; x < 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 < (startY + maskHeight); y++)
        {
            for (int x = startX; x < (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 < imageHeight; j++)
        {
            for (int i = 0; i < imageWidth; i++)
            {
                maxAmp = Math.Max(maxAmp, convolve[i, j]);
            }
        }

        double scale = 1.0 / maxAmp;

        for (int j = 0; j < imageHeight; j++)
        {
            for (int i = 0; i < 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;
    } 

Ответы

Ответ 1

Я нашел решение по этой ссылке. Основная идея заключалась в том, чтобы ввести offset и factor.

  • factor - сумма всех значений в ядре.
  • offset - произвольное значение для дальнейшего исправления вывода.

,

@Крис Luengo ответ также поднял действительную точку.

,

В данной ссылке указан следующий исходный код:

    private void SafeImageConvolution(Bitmap image, ConvMatrix fmat) 
    { 
        //Avoid division by 0 
        if (fmat.Factor == 0) 
            return; 

        Bitmap srcImage = (Bitmap)image.Clone(); 

        int x, y, filterx, filtery; 
        int s = fmat.Size / 2; 
        int r, g, b; 
        Color tempPix; 

        for (y = s; y < srcImage.Height - s; y++) 
        { 
            for (x = s; x < srcImage.Width - s; x++) 
            { 
                r = g = b = 0; 

                // Convolution 
                for (filtery = 0; filtery < fmat.Size; filtery++) 
                { 
                    for (filterx = 0; filterx < fmat.Size; filterx++) 
                    { 
                        tempPix = srcImage.GetPixel(x + filterx - s, y + filtery - s); 

                        r += fmat.Matrix[filtery, filterx] * tempPix.R; 
                        g += fmat.Matrix[filtery, filterx] * tempPix.G; 
                        b += fmat.Matrix[filtery, filterx] * tempPix.B; 
                    } 
                } 

                r = Math.Min(Math.Max((r / fmat.Factor) + fmat.Offset, 0), 255); 
                g = Math.Min(Math.Max((g / fmat.Factor) + fmat.Offset, 0), 255); 
                b = Math.Min(Math.Max((b / fmat.Factor) + fmat.Offset, 0), 255); 

                image.SetPixel(x, y, Color.FromArgb(r, g, b)); 
            } 
        } 
    } 

Ответ 2

Ваш текущий результат больше похож на функцию автокорреляции, чем свертка Лены. Я думаю, проблема может быть в вашей функции Sum.

Если вы посмотрите на определение суммы свертки, вы увидите, что ядро (или изображение не имеет значения) зеркалировано:

sum_m( f[n-m] g[m] )

Для одной функции m появляется со знаком плюса, а для другого - знаком минус.

Вам нужно будет изменить функцию Sum чтобы прочитать изображение mask1 в правильном порядке:

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 < (startY + maskHeight); y++)
    {
        for (int x = startX; x < (startX + maskWidth); x++)
        {
            double img = paddedImage1[x, y];
            double msk = mask1[maskWidth - x + startX - 1, maskHeight - y + startY - 1];
            sum = sum + (img * msk);
        }
    }

    return sum;
}

Другой вариант - передать зеркальную версию mask1 в эту функцию.