Geometry & Recognition

void retinex_scale_distribution(const int nscales/*=3*/, const int s/*=240*/,
                                double scales []) 
{    
    //  ASSERT(nscales>=3);
    const double size_step = double(s) / nscales;
    for (int i = 0; i < nscales; ++i)
        scales[i] = 2. + i * size_step;
}

// recursive Gaussian convolution;

void gauss_smooth (double *in, int size, int rowstride, double *out, double b[5]) {    
    /* forward pass */
    const int bufsize = size + 3;
    std::vector<double> w1(bufsize, 0), w2(bufsize, 0);
    w1[0] = in[0]; w1[1] = in[0]; w1[2] = in[0];
    for (int i = 0, n = 3; i < size ; i++, n++)
        w1[n] = b[4] * in[i*rowstride] + \
               (b[1] * w1[n-1] + b[2] * w1[n-2] + b[3] * w1[n-3]) / b[0];
    /* backward pass */
    w2[size+0] = w1[size+2];
    w2[size+1] = w1[size+2];
    w2[size+2] = w1[size+2];
    for (int i = size - 1, n = i; i >= 0; i--, n--)
        w2[n] = out[i*rowstride] = b[4] * w1[n] + \
                                  (b[1] * w2[n+1] + b[2] * w2[n+2] + b[3] * w2[n+3]) / b[0];
}

// 이미지의 평균과 표준편차 계산;

void image_statistics(double *img, int size /*=width*height*/,
                      double &mean, double &std) {    
    double s = 0, ss = 0;
    for (int i = size; i-- > 0;) {
        double a = img[i];
        s += a;
        ss += SQR(a) ;
    } ;
    mean = s / size ;
    std = sqrt((ss - s * s / size) / size);
}
//
void rescale_range(double *data, int size) {
    double mean, sig;
    image_statistics(data, size, mean, sig);
    const double max_val = mean + 1.2 * sig;
    const double min_val = mean - 1.2 * sig;    
    double range = max_val - min_val;
    if (!range) range = 1;
    range = 255 / range ;
    // change the range;
    for (int i = size; i-->0;)
        data[i] = (data[i] - min_val) * range;
}

// 메인;

void retinex_process(double *src, int width, int height, double *dst) {
    const int nfilter = 3;
    const double weight = 1.0 / double(nfilter); //equal-weight for each scale; 
    const int default_scale = 240;
    double sigma[nfilter];
    retinex_scale_distribution(nfilter, default_scale, sigma);
    std::vector<double> in(width*height);
    std::vector<double> out(width*height, 0);
    // scale-space gauss_smooth;
    double c[5];
    for (int i = 0; i < nfilter; i++) {
        compute_coefs3(c, sigma[i]);
        // copy src to temporay buffer(=in);
        for (int pos = width*height; pos-->0; ) 
            in[pos] = double(src[pos]) + 1;
        // (1) horizontal convolution(stride = 1 for grey);
        for (int y = 0; y < height; y++) 
            gauss_smooth(&in[y * width], width, 1, &out[y * width], c);
        // (2) vertical convolution(stride = width for grey);
        in.swap(out);  //copy out to in;
        for (int x = 0; x < width; x++) 
            gauss_smooth(&in[x], height, width, &out[x], c);
        // 각 스케일에서 반사 성분을 누적(equal weight);
        for (int pos = width*height; pos-->0;) 
            dst[pos] += weight * (log(src[pos] + 1) - log(out[pos]));
    }
    const double alpha = 128;
    const double gain = 1;
    const double offset = 0;
    for (int pos = width*height; pos-->0;) {
        double logI = log(double(src[pos]) + 1);
    	dst[pos] = gain * ((log(alpha* (src[pos] + 1.)) - logI) * dst[pos]) + offset;
    }
    // scale to [0,255];
    rescale_range(dst, width*height); 
};

/** 2020.11.24.일 수정됨;**/

원본 이미지;

가우시안 Convolution 결과(sigma[0]=2);

가우시안 Convolution 결과(sigma[1]=82);

가우시안 Convolution 결과(sigma[2]=162);

결과 이미지(scale=240)