Здравствуйте, Sinclair, Вы писали:

S>Сходу тяжко врубиться в 20 страниц научного текста. Как выглядит код бинаризации по Sauvola? Какая часть алгоритма является "переменной" (зависит от специфики конкретной задачи), а какая — общей?

Держи

void binarize(int image_width, int image_height, unsigned char *bin_image, unsigned char *gray_image,
        int w, double k){
        int whalf = w>>1;
        
        // Calculate the integral image, and integral of the squared image
        double **integral_image,**rowsum_image;
        double **integral_sqimg,**rowsum_sqimg;

        unsigned char** gray_image_ptr = new unsigned char*[image_height];
        gray_image_ptr[0] = gray_image;
        for(int i = 1; i < image_height; i++)
            gray_image_ptr[i] = gray_image_ptr[i-1] + image_width;

        unsigned char** bin_image_ptr = new unsigned char*[image_height];
        bin_image_ptr[0] = bin_image;
        for(int i = 1; i < image_height; i++)
            bin_image_ptr[i] = bin_image_ptr[i-1] + image_width;

        integral_image = new double*[image_width];
        for(int i = 0; i < image_width; i++)
            integral_image[i] = new double[image_height];

        rowsum_image = new double*[image_width];
        for(int i = 0; i < image_width; i++)
            rowsum_image[i] = new double[image_height];

        integral_sqimg = new double*[image_width];
        for(int i = 0; i < image_width; i++)
            integral_sqimg[i] = new double[image_height];
        
        rowsum_sqimg = new double*[image_width];
        for(int i = 0; i < image_width; i++)
            rowsum_sqimg[i] = new double[image_height];
        

        int xmin,ymin,xmax,ymax;
        double diagsum,idiagsum,diff,sqdiagsum,sqidiagsum,sqdiff,area;
        double mean,std,threshold;
        
    // START 
        for(int j=0; j<image_height; j++){
            rowsum_image[0][j] = gray_image_ptr[j][0];
            rowsum_sqimg[0][j] = gray_image_ptr[j][0]*gray_image_ptr[j][0];
        }
        for(int i=1; i<image_width; i++){
            for(int j=0; j<image_height; j++){
                rowsum_image[i][j] = rowsum_image[i-1][j] + gray_image_ptr[j][i];
                rowsum_sqimg[i][j] = rowsum_sqimg[i-1][j] + gray_image_ptr[j][i]*gray_image_ptr[j][i];
            }
        }

        for(int i=0; i<image_width; i++){
            integral_image[i][0] = rowsum_image[i][0];
            integral_sqimg[i][0] = rowsum_sqimg[i][0];
        }
        for(int i=0; i<image_width; i++){
            for(int j=1; j<image_height; j++){
                integral_image[i][j] = integral_image[i][j-1] + rowsum_image[i][j];
                integral_sqimg[i][j] = integral_sqimg[i][j-1] + rowsum_sqimg[i][j];
            }
        }
    // END
        
        //Calculate the mean and standard deviation using the integral image
        
        for(int i=0; i<image_width; i++){
            for(int j=0; j<image_height; j++){
                xmin = __max(0,i-whalf);
                ymin = __max(0,j-whalf);
                xmax = __min(image_width-1,i+whalf);
                ymax = __min(image_height-1,j+whalf);
                area = (xmax-xmin+1)*(ymax-ymin+1);
                // area can't be 0 here
                // proof (assuming whalf >= 0): 
                // we'll prove that (xmax-xmin+1) > 0,
                // (ymax-ymin+1) is analogous
                // It's the same as to prove: xmax >= xmin
                // image_width - 1 >= 0         since image_width > i >= 0
                // i + whalf >= 0               since i >= 0, whalf >= 0
                // i + whalf >= i - whalf       since whalf >= 0
                // image_width - 1 >= i - whalf since image_width > i
                // --IM
  //              ASSERT(area);
                if(!xmin && !ymin){ // Point at origin
                    diff   = integral_image[xmax][ymax];
                    sqdiff = integral_sqimg[xmax][ymax];
                }
                else if(!xmin && ymin){ // first column
                    diff   = integral_image[xmax][ymax] - integral_image[xmax][ymin-1];
                    sqdiff = integral_sqimg[xmax][ymax] - integral_sqimg[xmax][ymin-1];
                }
                else if(xmin && !ymin){ // first row
                    diff   = integral_image[xmax][ymax] - integral_image[xmin-1][ymax];
                    sqdiff = integral_sqimg[xmax][ymax] - integral_sqimg[xmin-1][ymax];
                }
                else{ // rest of the image
                    diagsum    = integral_image[xmax][ymax] + integral_image[xmin-1][ymin-1];
                    idiagsum   = integral_image[xmax][ymin-1] + integral_image[xmin-1][ymax];
                    diff       = diagsum - idiagsum;
                    sqdiagsum  = integral_sqimg[xmax][ymax] + integral_sqimg[xmin-1][ymin-1];
                    sqidiagsum = integral_sqimg[xmax][ymin-1] + integral_sqimg[xmin-1][ymax];
                    sqdiff     = sqdiagsum - sqidiagsum;
                }
                
                mean = diff/area;
                std  = sqrt((sqdiff - diff*diff/area)/(area-1));
                threshold = mean*(1+k*((std/128)-1));

                if(gray_image_ptr[j][i] <= threshold) 
                    bin_image_ptr[j][i] = 0;
                else
                    bin_image_ptr[j][i] = 255;
            }
        }
        delete[] gray_image_ptr;

        delete[] bin_image_ptr;


        for(int i = 0; i < image_width; i++)
            delete[] integral_image[i];
        delete[] integral_image;

        for(int i = 0; i < image_width; i++)
            delete[] rowsum_image[i];
        delete[] rowsum_image;


        for(int i = 0; i < image_width; i++)
            delete[] integral_sqimg[i];
        delete[] integral_sqimg;
        
        for(int i = 0; i < image_width; i++)
            delete[] rowsum_sqimg[i];
        delete []rowsum_sqimg;

    }

Сама идея достаточно проста — вычисляются суммы и суммы квадратов для каждого пикселя , считая его центром окна размером в w (обычно 15-25), потом простенькая формула
Проблема в том, что если это делать в лоб, то будет O(N*M*w^2)

Поэтому вначале строятся массивы, содержащие суммы и суммы квадратов для каждого пикселя , и строятся они рекуррентно. Я пометил этот кусок // START и // END.

PD>>Нахождение линий по Хафу
PD>>https://ru.wikipedia.org/wiki/%D0%9F%D1%80%D0%B5%D0%BE%D0%B1%D1%80%D0%B0%D0%B7%D0%BE%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5_%D0%A5%D0%B0%D1%84%D0%B0
S>Опять же — интересно было бы посмотреть на то, как выглядит типичный код на С++, который ищет какие-то примитивы при помощи преобразования Хафа или

Мы его сами не писали, но вот тебе код, и даже на C#


PD>>Все это 99%, а простенький фильтр — 1%. Это тоже предложишь на LinQ переводить ? А если нет — зачем мне 1% кода переводить на LinQ ?
S>Нет. Я предложу писать код в максимально читаемой форме, где сказано "что", но не сказано "как".

Ну напишем 1% кода "в максимально читаемой форме", а 99% останутся в нечитаемой.

PD>>Эээ, не передергивай. Он как раз потратил время на оптимизацию, потому что в алгоритмах обработки изображений O(N^2) — это очень хорошо, а зачастую бывает существенно хуже. А ты этой оптимизацией и не занимался, да и не ясно, сможешь ли вообще ей заняться.
S>Я-то как раз смогу. Причём эта оптимизация пригодится не только в этой строчке этой программы, а во всех местах, где у меня вычисляется фильтр по изображению.

PD>>Ты сумел показать, что и этот фильтр можно подвести под механизм последовательной обработки, пусть понимаемой несколько расширенно. Однако когда дело дойдет до тех алгоритмов, где обработка ведется не последовательным, а более хитрым способом, ничего путного у тебя не получится
S>Ждём пример такого алгоритма. Как обычно — интересует алгоритм, в котором есть "переменная часть". В случае с фильтром у нас было задание ядра; в случае с СУБД — задание предиката для поиска.

Ну вот тебе Савола. В общем-то, тот же фильтр , так что бери его в качестве переменной части. Правда, окно там не размера 2, а обычно 15-25.
Ну и для Хафа тоже есть переменная часть. Например, можно окружности искать

PD>>Это как с последовательным доступом к файлам vs. прямым доступом. Если надо писать или читать в/из поток(а) — классы последовательного доступа замечательно работают. Если же надо файл модифицировать сложным образом — применять классы последовательного доступа, наверное, все же можно , но лучше не стоит.
S>Вот как раз "императивный for" — это алгоритм последовательного доступа.

Вообще-то, кроме for, есть еще и while

Здравствуйте, Sinclair, Вы писали:

S>Сходу тяжко врубиться в 20 страниц научного текста. Как выглядит код бинаризации по Sauvola? Какая часть алгоритма является "переменной" (зависит от специфики конкретной задачи), а какая — общей?

Держи

void binarize(int image_width, int image_height, unsigned char *bin_image, unsigned char *gray_image,
        int w, double k){
        int whalf = w>>1;
        
        // Calculate the integral image, and integral of the squared image
        double **integral_image,**rowsum_image;
        double **integral_sqimg,**rowsum_sqimg;

        unsigned char** gray_image_ptr = new unsigned char*[image_height];
        gray_image_ptr[0] = gray_image;
        for(int i = 1; i < image_height; i++)
            gray_image_ptr[i] = gray_image_ptr[i-1] + image_width;

        unsigned char** bin_image_ptr = new unsigned char*[image_height];
        bin_image_ptr[0] = bin_image;
        for(int i = 1; i < image_height; i++)
            bin_image_ptr[i] = bin_image_ptr[i-1] + image_width;

        integral_image = new double*[image_width];
        for(int i = 0; i < image_width; i++)
            integral_image[i] = new double[image_height];

        rowsum_image = new double*[image_width];
        for(int i = 0; i < image_width; i++)
            rowsum_image[i] = new double[image_height];

        integral_sqimg = new double*[image_width];
        for(int i = 0; i < image_width; i++)
            integral_sqimg[i] = new double[image_height];
        
        rowsum_sqimg = new double*[image_width];
        for(int i = 0; i < image_width; i++)
            rowsum_sqimg[i] = new double[image_height];
        

        int xmin,ymin,xmax,ymax;
        double diagsum,idiagsum,diff,sqdiagsum,sqidiagsum,sqdiff,area;
        double mean,std,threshold;
        
    // START 
        for(int j=0; j<image_height; j++){
            rowsum_image[0][j] = gray_image_ptr[j][0];
            rowsum_sqimg[0][j] = gray_image_ptr[j][0]*gray_image_ptr[j][0];
        }
        for(int i=1; i<image_width; i++){
            for(int j=0; j<image_height; j++){
                rowsum_image[i][j] = rowsum_image[i-1][j] + gray_image_ptr[j][i];
                rowsum_sqimg[i][j] = rowsum_sqimg[i-1][j] + gray_image_ptr[j][i]*gray_image_ptr[j][i];
            }
        }

        for(int i=0; i<image_width; i++){
            integral_image[i][0] = rowsum_image[i][0];
            integral_sqimg[i][0] = rowsum_sqimg[i][0];
        }
        for(int i=0; i<image_width; i++){
            for(int j=1; j<image_height; j++){
                integral_image[i][j] = integral_image[i][j-1] + rowsum_image[i][j];
                integral_sqimg[i][j] = integral_sqimg[i][j-1] + rowsum_sqimg[i][j];
            }
        }
    // END
        
        //Calculate the mean and standard deviation using the integral image
        
        for(int i=0; i<image_width; i++){
            for(int j=0; j<image_height; j++){
                xmin = __max(0,i-whalf);
                ymin = __max(0,j-whalf);
                xmax = __min(image_width-1,i+whalf);
                ymax = __min(image_height-1,j+whalf);
                area = (xmax-xmin+1)*(ymax-ymin+1);
                // area can't be 0 here
                // proof (assuming whalf >= 0): 
                // we'll prove that (xmax-xmin+1) > 0,
                // (ymax-ymin+1) is analogous
                // It's the same as to prove: xmax >= xmin
                // image_width - 1 >= 0         since image_width > i >= 0
                // i + whalf >= 0               since i >= 0, whalf >= 0
                // i + whalf >= i - whalf       since whalf >= 0
                // image_width - 1 >= i - whalf since image_width > i
                // --IM
  //              ASSERT(area);
                if(!xmin && !ymin){ // Point at origin
                    diff   = integral_image[xmax][ymax];
                    sqdiff = integral_sqimg[xmax][ymax];
                }
                else if(!xmin && ymin){ // first column
                    diff   = integral_image[xmax][ymax] - integral_image[xmax][ymin-1];
                    sqdiff = integral_sqimg[xmax][ymax] - integral_sqimg[xmax][ymin-1];
                }
                else if(xmin && !ymin){ // first row
                    diff   = integral_image[xmax][ymax] - integral_image[xmin-1][ymax];
                    sqdiff = integral_sqimg[xmax][ymax] - integral_sqimg[xmin-1][ymax];
                }
                else{ // rest of the image
                    diagsum    = integral_image[xmax][ymax] + integral_image[xmin-1][ymin-1];
                    idiagsum   = integral_image[xmax][ymin-1] + integral_image[xmin-1][ymax];
                    diff       = diagsum - idiagsum;
                    sqdiagsum  = integral_sqimg[xmax][ymax] + integral_sqimg[xmin-1][ymin-1];
                    sqidiagsum = integral_sqimg[xmax][ymin-1] + integral_sqimg[xmin-1][ymax];
                    sqdiff     = sqdiagsum - sqidiagsum;
                }
                
                mean = diff/area;
                std  = sqrt((sqdiff - diff*diff/area)/(area-1));
                threshold = mean*(1+k*((std/128)-1));

                if(gray_image_ptr[j][i] <= threshold) 
                    bin_image_ptr[j][i] = 0;
                else
                    bin_image_ptr[j][i] = 255;
            }
        }
        delete[] gray_image_ptr;

        delete[] bin_image_ptr;


        for(int i = 0; i < image_width; i++)
            delete[] integral_image[i];
        delete[] integral_image;

        for(int i = 0; i < image_width; i++)
            delete[] rowsum_image[i];
        delete[] rowsum_image;


        for(int i = 0; i < image_width; i++)
            delete[] integral_sqimg[i];
        delete[] integral_sqimg;
        
        for(int i = 0; i < image_width; i++)
            delete[] rowsum_sqimg[i];
        delete []rowsum_sqimg;

    }

Сама идея достаточно проста — вычисляются суммы и суммы квадратов для каждого пикселя , считая его центром окна размером в w (обычно 15-25), потом простенькая формула
Проблема в том, что если это делать в лоб, то будет O(N*M*w^2)

Поэтому вначале строятся массивы, содержащие суммы и суммы квадратов для каждого пикселя , и строятся они рекуррентно. Я пометил этот кусок // START и // END.

PD>>Нахождение линий по Хафу
PD>>https://ru.wikipedia.org/wiki/%D0%9F%D1%80%D0%B5%D0%BE%D0%B1%D1%80%D0%B0%D0%B7%D0%BE%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5_%D0%A5%D0%B0%D1%84%D0%B0
S>Опять же — интересно было бы посмотреть на то, как выглядит типичный код на С++, который ищет какие-то примитивы при помощи преобразования Хафа или

Мы его сами не писали, но вот тебе код, и даже на C#


PD>>Все это 99%, а простенький фильтр — 1%. Это тоже предложишь на LinQ переводить ? А если нет — зачем мне 1% кода переводить на LinQ ?
S>Нет. Я предложу писать код в максимально читаемой форме, где сказано "что", но не сказано "как".

Ну напишем 1% кода "в максимально читаемой форме", а 99% останутся в нечитаемой.

PD>>Эээ, не передергивай. Он как раз потратил время на оптимизацию, потому что в алгоритмах обработки изображений O(N^2) — это очень хорошо, а зачастую бывает существенно хуже. А ты этой оптимизацией и не занимался, да и не ясно, сможешь ли вообще ей заняться.
S>Я-то как раз смогу. Причём эта оптимизация пригодится не только в этой строчке этой программы, а во всех местах, где у меня вычисляется фильтр по изображению.

PD>>Ты сумел показать, что и этот фильтр можно подвести под механизм последовательной обработки, пусть понимаемой несколько расширенно. Однако когда дело дойдет до тех алгоритмов, где обработка ведется не последовательным, а более хитрым способом, ничего путного у тебя не получится
S>Ждём пример такого алгоритма. Как обычно — интересует алгоритм, в котором есть "переменная часть". В случае с фильтром у нас было задание ядра; в случае с СУБД — задание предиката для поиска.

Ну вот тебе Савола. В общем-то, тот же фильтр , так что бери его в качестве переменной части. Правда, окно там не размера 2, а обычно 15-25.
Ну и для Хафа тоже есть переменная часть. Например, можно окружности искать

PD>>Это как с последовательным доступом к файлам vs. прямым доступом. Если надо писать или читать в/из поток(а) — классы последовательного доступа замечательно работают. Если же надо файл модифицировать сложным образом — применять классы последовательного доступа, наверное, все же можно , но лучше не стоит.
S>Вот как раз "императивный for" — это алгоритм последовательного доступа.

Вообще-то, кроме for, есть еще и while

Кстати, как насчет binary search с помощью LinQ ? Разумеется, без написания цикла.