Problem with addition of neighboring pixels

Question

I am trying to implement de-noising using non local means algorithm. I am having an issue in regards to addition of pixels.

// Find the patch differences by blurring the difference images
            d_y_channels[0].at<float>(i, j) = d_channels[0].at<float>(max(i - 1, 0), j) + d_channels[0].at<float>(i, j) + d_channels[0].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[1].at<float>(i, j) = d_channels[1].at<float>(max(i - 1, 0), j) + d_channels[1].at<float>(i, j) + d_channels[1].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[2].at<float>(i, j) = d_channels[2].at<float>(max(i - 1, 0), j) + d_channels[2].at<float>(i, j) + d_channels[2].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[3].at<float>(i, j) = d_channels[3].at<float>(max(i - 1, 0), j) + d_channels[3].at<float>(i, j) + d_channels[3].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[4].at<float>(i, j) = d_channels[4].at<float>(max(i - 1, 0), j) + d_channels[4].at<float>(i, j) + d_channels[4].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[5].at<float>(i, j) = d_channels[5].at<float>(max(i - 1, 0), j) + d_channels[5].at<float>(i, j) + d_channels[5].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[6].at<float>(i, j) = d_channels[6].at<float>(max(i - 1, 0), j) + d_channels[6].at<float>(i, j) + d_channels[6].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[7].at<float>(i, j) = d_channels[7].at<float>(max(i - 1, 0), j) + d_channels[7].at<float>(i, j) + d_channels[7].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[8].at<float>(i, j) = d_channels[8].at<float>(max(i - 1, 0), j) + d_channels[8].at<float>(i, j) + d_channels[8].at<float>(min(i + 1, img_rev.rows - 1), j);

When I tried to print the output using the following code,

if (i == 154 && j == 35) {
                cout << endl << "At 154 and 35";
                cout << endl << d_y_channels[0].at<float>(i, j) << " = " << d_channels[0].at<float>(max(i - 1, 0), j) << " + " << d_channels[0].at<float>(i, j) << " + " << d_channels[0].at<float>(min(i + 1, img_rev.rows - 1), j);
                cout << endl << d_y_channels[0].at<float>(i, j) << "," << d_y_channels[1].at<float>(i, j) << "," << d_y_channels[2].at<float>(i, j);
                cout << endl << d_y_channels[3].at<float>(i, j) << "," << d_y_channels[4].at<float>(i, j) << "," << d_y_channels[5].at<float>(i, j);
                cout << endl << d_y_channels[6].at<float>(i, j) << "," << d_y_channels[7].at<float>(i, j) << "," << d_y_channels[8].at<float>(i, j);

            }

it gave me this,

As you can see, the addition is not properly updated.

However, When I modify the line

d_y_channels[0].at<float>(i, j) = d_channels[0].at<float>(max(i - 1, 0), j) + d_channels[0].at<float>(i, j) + d_channels[0].at<float>(min(i + 1, img_rev.rows - 1), j);

to

d_y_channels[0].at<float>(i, j) = d_channels[0].at<float>(max(i - 1, 0), j) + d_channels[0].at<float>(i, j);// +d_channels[0].at<float>(min(i + 1, img_rev.rows - 1), j);

It is giving me an apt output.

I don't understand where it is going wrong. I checked the entire code and declarations multiple times, but i am not able to figure out why its going wrong with d_channels[0].at<float>(min(i + 1, img_rev.rows - 1), j).

I would really be thankful if someone can help me with this.

Thanks in advance.

I have attached the entire code(with declarations) below for better understanding

// image_read.cpp : Defines the entry point for the console application.
//

//#include "stdafx.h"
#include "stdafx.h"
#include <opencv2/opencv.hpp>
#include <iostream>
#include <fstream>
using namespace cv;
using namespace std;
#include <vector>
#include <string>
#include "math.h"
#include <chrono>    

int main(int argc, char** argv)
{
    Mat img_rev = imread("C:/Users/20181217/Desktop/images/imgs/denoise_input_h.png");
    //Mat ideal = imread("C:/Users/20181217/Desktop/images/imgs/den_check_2.png");
    int search_area = 3;
    int patch_size = 3;//dx
    float sigma = 0.12;//dy
    cout << img_rev.rows << "," << img_rev.cols << endl;
    cvtColor(img_rev, img_rev, COLOR_BGR2RGB);
    img_rev.convertTo(img_rev, CV_32F);
    img_rev = img_rev / 256.0;
    array<Mat, 3> channels;
    split(img_rev, channels);
    Mat red = channels[0].clone();
    Mat green = channels[1].clone();
    Mat blue = channels[2].clone();

    //define nine channels, one for each in the 3x3 search window
    Mat dc_red_four(img_rev.rows, img_rev.cols, CV_32FC(9));
    Mat dc_green_four(img_rev.rows, img_rev.cols, CV_32FC(9));
    Mat dc_blue_four(img_rev.rows, img_rev.cols, CV_32FC(9));
    cout << dc_blue_four.channels();

    array<Mat, 9> red_channels;
    split(dc_red_four, red_channels);

    array<Mat, 9> green_channels;
    split(dc_green_four, green_channels);

    array<Mat, 9> blue_channels;
    split(dc_blue_four, blue_channels);

    Mat d(img_rev.rows, img_rev.cols, CV_32FC(9));
    array<Mat, 9> d_channels;
    split(d, d_channels);

    Mat blur_d_y(img_rev.rows, img_rev.cols, CV_32FC(9));
    array<Mat, 9> d_y_channels;
    split(blur_d_y, d_y_channels);
    /*
    0 = -1,-1
    1 = -1, 0
    2 = -1, 1
    3 =  0,-1
    4 =  0, 0
    5 =  0, 1
    6 =  1,-1
    7 =  1, 0
    8 =  1, 1
    */

    for (int i = 0; i < img_rev.rows; i++)
    {
        for (int j = 0; j < img_rev.cols; j++)
        {
            //creating the 3x3 window and squaring the distance
            red_channels[0].at<float>(i, j) = (red.at<float>(i, j) - red.at<float>(max(i - 1, 0), max(j - 1, 0))) * (red.at<float>(i, j) - red.at<float>(max(i - 1, 0), max(j - 1, 0)));
            red_channels[1].at<float>(i, j) = (red.at<float>(i, j) - red.at<float>(max(i - 1, 0), j)) * (red.at<float>(i, j) - red.at<float>(max(i - 1, 0), j));
            red_channels[2].at<float>(i, j) = (red.at<float>(i, j) - red.at<float>(max(i - 1, 0), min(j + 1, img_rev.cols - 1))) * (red.at<float>(i, j) - red.at<float>(max(i - 1, 0), min(j + 1, img_rev.cols - 1)));

            red_channels[3].at<float>(i, j) = (red.at<float>(i, j) - red.at<float>(i, max(j - 1, 0))) * (red.at<float>(i, j) - red.at<float>(i, max(j - 1, 0)));
            red_channels[4].at<float>(i, j) = (red.at<float>(i, j) - red.at<float>(i, j)) * (red.at<float>(i, j) - red.at<float>(i, j));
            red_channels[5].at<float>(i, j) = (red.at<float>(i, j) - red.at<float>(i, min(j + 1, img_rev.cols - 1))) * (red.at<float>(i, j) - red.at<float>(i, min(j + 1, img_rev.cols - 1)));

            red_channels[6].at<float>(i, j) = (red.at<float>(i, j) - red.at<float>(min(i + 1, img_rev.rows - 1), max(j - 1, 0))) * (red.at<float>(i, j) - red.at<float>(min(i + 1, img_rev.rows - 1), max(j - 1, 0)));
            red_channels[7].at<float>(i, j) = (red.at<float>(i, j) - red.at<float>(min(i + 1, img_rev.rows - 1), j)) * (red.at<float>(i, j) - red.at<float>(min(i + 1, img_rev.rows - 1), j));
            red_channels[8].at<float>(i, j) = (red.at<float>(i, j) - red.at<float>(min(i + 1, img_rev.rows - 1), min(j + 1, img_rev.cols - 1))) * (red.at<float>(i, j) - red.at<float>(min(i + 1, img_rev.rows - 1), min(j + 1, img_rev.cols - 1)));

            //green
            green_channels[0].at<float>(i, j) = (green.at<float>(i, j) - green.at<float>(max(i - 1, 0), max(j - 1, 0))) * (green.at<float>(i, j) - green.at<float>(max(i - 1, 0), max(j - 1, 0)));
            green_channels[1].at<float>(i, j) = (green.at<float>(i, j) - green.at<float>(max(i - 1, 0), j)) * (green.at<float>(i, j) - green.at<float>(max(i - 1, 0), j));
            green_channels[2].at<float>(i, j) = (green.at<float>(i, j) - green.at<float>(max(i - 1, 0), min(j + 1, img_rev.cols - 1))) * (green.at<float>(i, j) - green.at<float>(max(i - 1, 0), min(j + 1, img_rev.cols - 1)));

            green_channels[3].at<float>(i, j) = (green.at<float>(i, j) - green.at<float>(i, max(j - 1, 0))) * (green.at<float>(i, j) - green.at<float>(i, max(j - 1, 0)));
            green_channels[4].at<float>(i, j) = (green.at<float>(i, j) - green.at<float>(i, j)) * (green.at<float>(i, j) - green.at<float>(i, j));
            green_channels[5].at<float>(i, j) = (green.at<float>(i, j) - green.at<float>(i, min(j + 1, img_rev.cols - 1))) * (green.at<float>(i, j) - green.at<float>(i, min(j + 1, img_rev.cols - 1)));

            green_channels[6].at<float>(i, j) = (green.at<float>(i, j) - green.at<float>(min(i + 1, img_rev.rows - 1), max(j - 1, 0))) * (green.at<float>(i, j) - green.at<float>(min(i + 1, img_rev.rows - 1), max(j - 1, 0)));
            green_channels[7].at<float>(i, j) = (green.at<float>(i, j) - green.at<float>(min(i + 1, img_rev.rows - 1), j)) * (green.at<float>(i, j) - green.at<float>(min(i + 1, img_rev.rows - 1), j));
            green_channels[8].at<float>(i, j) = (green.at<float>(i, j) - green.at<float>(min(i + 1, img_rev.rows - 1), min(j + 1, img_rev.cols - 1))) * (green.at<float>(i, j) - green.at<float>(min(i + 1, img_rev.rows - 1), min(j + 1, img_rev.cols - 1)));

            //blue
            blue_channels[0].at<float>(i, j) = (blue.at<float>(i, j) - blue.at<float>(max(i - 1, 0), max(j - 1, 0))) * (blue.at<float>(i, j) - blue.at<float>(max(i - 1, 0), max(j - 1, 0)));
            blue_channels[1].at<float>(i, j) = (blue.at<float>(i, j) - blue.at<float>(max(i - 1, 0), j)) * (blue.at<float>(i, j) - blue.at<float>(max(i - 1, 0), j));
            blue_channels[2].at<float>(i, j) = (blue.at<float>(i, j) - blue.at<float>(max(i - 1, 0), min(j + 1, img_rev.cols - 1))) * (blue.at<float>(i, j) - blue.at<float>(max(i - 1, 0), min(j + 1, img_rev.cols - 1)));

            blue_channels[3].at<float>(i, j) = (blue.at<float>(i, j) - blue.at<float>(i, max(j - 1, 0))) * (blue.at<float>(i, j) - blue.at<float>(i, max(j - 1, 0)));
            blue_channels[4].at<float>(i, j) = (blue.at<float>(i, j) - blue.at<float>(i, j)) * (blue.at<float>(i, j) - blue.at<float>(i, j));
            blue_channels[5].at<float>(i, j) = (blue.at<float>(i, j) - blue.at<float>(i, min(j + 1, img_rev.cols - 1))) * (blue.at<float>(i, j) - blue.at<float>(i, min(j + 1, img_rev.cols - 1)));

            blue_channels[6].at<float>(i, j) = (blue.at<float>(i, j) - blue.at<float>(min(i + 1, img_rev.rows - 1), max(j - 1, 0))) * (blue.at<float>(i, j) - blue.at<float>(min(i + 1, img_rev.rows - 1), max(j - 1, 0)));
            blue_channels[7].at<float>(i, j) = (blue.at<float>(i, j) - blue.at<float>(min(i + 1, img_rev.rows - 1), j)) * (blue.at<float>(i, j) - blue.at<float>(min(i + 1, img_rev.rows - 1), j));
            blue_channels[8].at<float>(i, j) = (blue.at<float>(i, j) - blue.at<float>(min(i + 1, img_rev.rows - 1), min(j + 1, img_rev.cols - 1))) * (blue.at<float>(i, j) - blue.at<float>(min(i + 1, img_rev.rows - 1), min(j + 1, img_rev.cols - 1)));

            //sum across colour channels
            d_channels[0].at<float>(i, j) = red_channels[0].at<float>(i, j) + green_channels[0].at<float>(i, j) + blue_channels[0].at<float>(i, j);
            d_channels[1].at<float>(i, j) = red_channels[1].at<float>(i, j) + green_channels[1].at<float>(i, j) + blue_channels[1].at<float>(i, j);
            d_channels[2].at<float>(i, j) = red_channels[2].at<float>(i, j) + green_channels[2].at<float>(i, j) + blue_channels[2].at<float>(i, j);
            d_channels[3].at<float>(i, j) = red_channels[3].at<float>(i, j) + green_channels[3].at<float>(i, j) + blue_channels[3].at<float>(i, j);
            d_channels[4].at<float>(i, j) = red_channels[4].at<float>(i, j) + green_channels[4].at<float>(i, j) + blue_channels[4].at<float>(i, j);
            d_channels[5].at<float>(i, j) = red_channels[5].at<float>(i, j) + green_channels[5].at<float>(i, j) + blue_channels[5].at<float>(i, j);
            d_channels[6].at<float>(i, j) = red_channels[6].at<float>(i, j) + green_channels[6].at<float>(i, j) + blue_channels[6].at<float>(i, j);
            d_channels[7].at<float>(i, j) = red_channels[7].at<float>(i, j) + green_channels[7].at<float>(i, j) + blue_channels[7].at<float>(i, j);
            d_channels[8].at<float>(i, j) = red_channels[8].at<float>(i, j) + green_channels[8].at<float>(i, j) + blue_channels[8].at<float>(i, j);

            // Find the patch differences by blurring the difference images
            d_y_channels[0].at<float>(i, j) = d_channels[0].at<float>(max(i - 1, 0), j) + d_channels[0].at<float>(i, j);// +d_channels[0].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[1].at<float>(i, j) = d_channels[1].at<float>(max(i - 1, 0), j) + d_channels[1].at<float>(i, j) + d_channels[1].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[2].at<float>(i, j) = d_channels[2].at<float>(max(i - 1, 0), j) + d_channels[2].at<float>(i, j) + d_channels[2].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[3].at<float>(i, j) = d_channels[3].at<float>(max(i - 1, 0), j) + d_channels[3].at<float>(i, j) + d_channels[3].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[4].at<float>(i, j) = d_channels[4].at<float>(max(i - 1, 0), j) + d_channels[4].at<float>(i, j) + d_channels[4].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[5].at<float>(i, j) = d_channels[5].at<float>(max(i - 1, 0), j) + d_channels[5].at<float>(i, j) + d_channels[5].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[6].at<float>(i, j) = d_channels[6].at<float>(max(i - 1, 0), j) + d_channels[6].at<float>(i, j) + d_channels[6].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[7].at<float>(i, j) = d_channels[7].at<float>(max(i - 1, 0), j) + d_channels[7].at<float>(i, j) + d_channels[7].at<float>(min(i + 1, img_rev.rows - 1), j);
            d_y_channels[8].at<float>(i, j) = d_channels[8].at<float>(max(i - 1, 0), j) + d_channels[8].at<float>(i, j) + d_channels[8].at<float>(min(i + 1, img_rev.rows - 1), j);

            

            if (i == 154 && j == 35) {
                cout << endl << "At 154 and 35";
                cout << endl << d_y_channels[0].at<float>(i, j) << " = " << d_channels[0].at<float>(max(i - 1, 0), j) << " + " << d_channels[0].at<float>(i, j);
                cout << endl << d_y_channels[0].at<float>(i, j) << "," << d_y_channels[1].at<float>(i, j) << "," << d_y_channels[2].at<float>(i, j);
                cout << endl << d_y_channels[3].at<float>(i, j) << "," << d_y_channels[4].at<float>(i, j) << "," << d_y_channels[5].at<float>(i, j);
                cout << endl << d_y_channels[6].at<float>(i, j) << "," << d_y_channels[7].at<float>(i, j) << "," << d_y_channels[8].at<float>(i, j);

            }


        }
    }




    //cvtColor(img_rev, img_rev, COLOR_BGR2RGB);
    //cout << img_rev;

    //cout << ideal;
    cout << endl << "finished";
    waitKey(100000000);
    return 0;
}

Answer 1

Do you know what -4.316e+08 means? It is 400 million. Adding 400 million + 0.04 is 400 million. This addition is totally correct.

The problem is that you are computing d_channels sequentially, and within the same loop you want to compute d_y_channels that needs a value of d_channels that you haven't computed yet. Thus, you read an uninitialized value, which in this case is -4.316e+08, but could be anything.

You should do the blurring in a separate loop. This way, you first compute all of d_channels, and when you do the blurring, all the values you need are already computed.