c++ - Replace a chain of image blurs with one blur -

in this question asked how implement chain of blurs in 1 single step.

then found out gaussian blur page of wikipedia that:

applying multiple, successive gaussian blurs image has same effect applying single, larger gaussian blur, radius square root of sum of squares of blur radii applied. example, applying successive gaussian blurs radii of 6 , 8 gives same results applying single gaussian blur of radius 10, since sqrt {6^{2}+8^{2}}=10.

so thought blur , singleblur same in following code:

cv::mat firstlevel; float sigma1, sigma2; //intialize firstlevel, sigma1 , sigma2 cv::mat blur = gaussianblur(firstlevel, sigma1);         blur = gaussianblur(blur, sigma2); float singlesigma = std::sqrt(std::pow(sigma1,2)+std::pow(sigma2,2)); cv::mat singleblur = gaussianblur(firstlevel, singlesigma); cv::mat diff = blur != singleblur; // equal if no elements disagree assert( cv::countnonzero(diff) == 0); 

but assert fails (and actually, example, first row of blur different first 1 of singleblur).

why?

update:

after different comments asking more information, i'll update answer.

what i'm trying parallelize this code. in particular, i'm focusing on computing blurs @ levels in advance. serial code (which works correctly) following:

   vector<mat> blurs ((par.numberofscales+3)*levels, mat());    cv::mat octavelayer = firstlevel;    int scalecycles = par.numberofscales+2;     //compute blurs @ layers (not parallelizable)    for(int i=0; i<levels; i++){        blurs[i*scalecycles+1] = octavelayer.clone();        (int j = 1; j < scalecycles; j++){            float sigma = par.sigmas[j]* sqrt(sigmastep * sigmastep - 1.0f);            blurs[j+1+i*scalecycles] = gaussianblur(blurs[j+i*scalecycles], sigma);            if(j == par.numberofscales)                octavelayer = halfimage(blurs[j+1+i*scalecycles]);        }    } 

where:

mat halfimage(const mat &input) {    mat n(input.rows/2, input.cols/2, input.type());    float *out = n.ptr<float>(0);    (int r = 0, ri = 0; r < n.rows; r++, ri += 2)       (int c = 0, ci = 0; c < n.cols; c++, ci += 2)          *out++ = input.at<float>(ri,ci);    return n; }  mat gaussianblur(const mat input, const float sigma) {    mat ret(input.rows, input.cols, input.type());    int size = (int)(2.0 * 3.0 * sigma + 1.0); if (size % 2 == 0) size++;          gaussianblur(input, ret, size(size, size), sigma, sigma, border_replicate);    return ret; } 

i'm sorry horrible indexes above, tried respect original code system (which horrible, starting counting 1 instead of 0). code above has scalecycles=5 , levels=6, 30 blurs generated in total.

this "single blur" version, first compute sigmas each blur has computed (following wikipedia's formula) , apply blur (notice still serial , not parallelizable):

   vector<mat> singleblurs ((par.numberofscales+3)*levels, mat());    vector<float> singlesigmas(scalecycles);    float acc = 0;    (int j = 1; j < scalecycles; j++){        float sigma = par.sigmas[j]* sqrt(sigmastep * sigmastep - 1.0f);        acc += pow(sigma, 2);        singlesigmas[j] = sqrt(acc);    }     octavelayer = firstlevel;    for(int i=0; i<levels; i++){        singleblurs[i*scalecycles+1] = octavelayer.clone();        (int j = 1; j < scalecycles; j++){            float sigma = singlesigmas[j];            std::cout<<"j="<<j<<" sigma="<<sigma<<std::endl;            singleblurs[j+1+i*scalecycles] = gaussianblur(singleblurs[j+i*scalecycles], sigma);            if(j == par.numberofscales)                octavelayer = halfimage(singleblurs[j+1+i*scalecycles]);        }    } 

of course code above generates 30 blurs same parameters of previous version.

and code see difference between each signgleblurs , blurs:

   assert(blurs.size() == singleblurs.size());    vector<mat> blurdiffs(blurs.size());    for(int i=1; i<levels*scalecycles; i++){         cv::mat diff;         absdiff(blurs[i], singleblurs[i], diff);         std::cout<<"i="<<i<<"diff rows="<<diff.rows<<" cols="<<diff.cols<<std::endl;         blurdiffs[i] = diff;         std::cout<<"blurs rows="<<blurs[i].rows<<" cols="<<blurs[i].cols<<std::endl;         std::cout<<"singleblurs rows="<<singleblurs[i].rows<<" cols="<<singleblurs[i].cols<<std::endl;         std::cout<<"blurdiffs rows="<<blurdiffs[i].rows<<" cols="<<blurdiffs[i].cols<<std::endl;         namedwindow( "bluediffs["+std::to_string(i)+"]", window_autosize );// create window display.         //imshow( "bluediffs["+std::to_string(i)+"]", blurdiffs[i] );                   // show our image inside it.         //waitkey(0);                                          // wait keystroke in window         mat imagef_8uc3;         std::cout<<"converting..."<<std::endl;         blurdiffs[i].convertto(imagef_8uc3, cv_8u, 255);         std::cout<<"converted"<<std::endl;         imwrite( "blurdiffs_"+std::to_string(i)+".jpg", imagef_8uc3);    } 

now, i've seen blurdiffs_1.jpg , blurdiffs_2.jpg black, suddendly blurdiffs_3.jpg until blurdiffs_29.jpg becomes whiter , whiter. reason, blurdiffs_30.jpg black.

the first (correct) version generates 1761 descriptors. second (uncorrect) version generates >2.3k descriptors.

i can't post blurdiffs matrices because (especially first ones) big , post's space limited. i'll post samples. i'll not post blurdiffs_1.jpg , blurdiffs_2.jpg because they're totally blacks. notice because of halfimage images become smaller , smaller (as expected).

blurdiffs_3.jpg:

blurdiffs_6.jpg:

blurdiffs_15.jpg:

blurdiffs_29.jpg:

how image read:

  mat tmp = imread(argv[1]);   mat image(tmp.rows, tmp.cols, cv_32fc1, scalar(0));    float *out = image.ptr<float>(0);   unsigned char *in  = tmp.ptr<unsigned char>(0);     (size_t i=tmp.rows*tmp.cols; > 0; i--)   {      *out = (float(in[0]) + in[1] + in[2])/3.0f;      out++;      in+=3;   } 

someone here suggested divide diff 255 see real difference, don't understand why of understood him correctly.

if need more details, please let me know.

a warning front: have no experience opencv, following relevant computing gaussian blur in general. , applicable threw glance @ opencv documentation wrt border treatment , use of finite kernels (fir filtering).

1. as aside: initial test sensitive round-off, have cleared issue , shown errors larger.

2. beware image border effects. pixels near edge, image virtually extended using 1 of offered methods (border_default, border_replicate, etc...). if image |abcd| , use border_replicate filtering extended image aaaa|abcd|dddd. result klmn|opqr|stuv. there new pixel values (k,l,m,n,s,t,u,v) discarded yield output image |opqr|. if apply gaussian blur, blur operate on newly extended image oooo|opqr|rrrr - different "true" intermediate result , giving result different obtained single gaussian blur larger sigma. these extension methods safe though: reflect, reflect_101, wrap.

3. using finite kernel size g(s1)*g(s2)=g(sqrt(s1^2+s2^2)) rule not hold in general due tails of kernel being cut off. can reduce error introduced increasing kernel size relative sigma, e.g.:

   int size = (int)(2.0 * 10.0 * sigma + 1.0); if (size % 2 == 0) size++; 

point 3 seems issue "biting" you. care whether g(s1)*g(s2) property preserved. both results wrong in way. affect methodology works on result in major way? note example of using 10x sigma have given may resolve difference, slower.

update: forgot add might practical resolution. compute gaussian blur using fourier transform. scheme be:

• compute fourier transform (fft) of input image
• multiply fourier transform of gaussian kernel , compute inverse fourier transform. ignore imaginary part of complex output.

you can find equation gaussian in frequency domain on wikipedia

you can perform second step separately (i.e. in parallel) each scale (sigma). border condition implied computing blur way border_wrap. if prefer can achieve same border_reflect if use discrete cosine transform (dct) instead. not know if opencv provides one. after dct-ii