C matrix libreries

Hello people
I found a function to binarize a gray image in matlab (Niblack algorithm) here.

I make a implementation in C and get the same results but the excution time is very huge.
in matlab. 31 seconds and
my C function 1 hour.
i also use and implementation fourier lib and is slow too.

Code:

---

unsigned char **im2bw_niblack(        double **mat,
                                int nrows, int ncols,
                                int sizeWindow,
                                double weight)
{

    int halfSizeWindow = sizeWindow / 2;
    int
        minX = halfSizeWindow, minY = halfSizeWindow,  // Esquina superior Izquierda
        maxX = nrows-1-halfSizeWindow, maxY = ncols-1- halfSizeWindow;  // Esquina inferior Derecha
    unsigned char **ucmat;                  //Matriz Binarizada de salida.
    ucmat=ucmatrix(0,nrows-1,0,ncols-1);    // Se inicializa la matriz.

    // Se calcula la mediana
    double sumPixelsWindow=0;        // Suma de los pixeles en la vecindad.
    double sum2PixelsWindow=0;        // Suma de los cuadrados de la vecindad.
    int i=0,j=0,k=0,l=0;
    double localMean=0, localVar=0, localStd=0;
    double localValue=0;
    for(j=minY; j<=maxY; j++){
        for(i=minX; i<=maxX; i++){
            // Se calcula la media para cada pixel
            sumPixelsWindow =0;
            sum2PixelsWindow =0;

            for(l=j-halfSizeWindow;l<=j+halfSizeWindow;l++){
                for(k=i-halfSizeWindow;k<=i+halfSizeWindow;k++){
                    sumPixelsWindow += mat[k][l];      // Se acumulan las vecindades
                    sum2PixelsWindow += mat[k][l]*mat[k][l];

                }
            }
           
            localMean  = sumPixelsWindow/(sizeWindow*sizeWindow);  //Se halla la media.
            localVar    = sum2PixelsWindow/(sizeWindow*sizeWindow) - localMean*localMean;  // Se calcula la varianza.
            localStd    = sqrt(fabs(localVar));          // Se calcula la desviacion estandard.

            localValue = localMean+weight*localStd;

            // Se binariza
            if(mat[i][j]< localValue){
                ucmat[i][j] = 0;           
            }else{
                ucmat[i][j] = 1;
            }

            printf("ucmat[i][j] = %c\n",ucmat[i][j]);
           
        }
    }

    return (ucmat);

} // end sub im2bw_niblack()

---

thanks.

This is very inefficient.

Code:

---

            for(l=j-halfSizeWindow;l<=j+halfSizeWindow;l++){
                for(k=i-halfSizeWindow;k<=i+halfSizeWindow;k++){
                    sumPixelsWindow += mat[k][l];      // Se acumulan las vecindades
                    sum2PixelsWindow += mat[k][l]*mat[k][l];

                }
            }

---

You've already calculated most of that summation for the previous iteration of the enclosing loop.

You could optimize that in several ways: since the window will just shift one to the right for each iteration of the inner loop (after the first iteration, of course), you could subtract the sum of the first column, add the sum of the last column, and you'll have the right value. Since you use the square of the elements frequently, it might be worth it to pre-calculate the square of each element in the matrix that you'll be concerned with.

Also, be careful when sizeWindow is odd. I don't think your window summation loop will work properly then.

thanks dwks thats a really good idea.
I will try and then post again.

This is my implementation

Code:

---


unsigned char **im2bw_niblack(        double **mat,
                                int nrows, int ncols,
                                int sizeWindow,
                                double weight)
{
    int halfSizeWindow = sizeWindow / 2;
    int
        minX = halfSizeWindow, minY = halfSizeWindow,  // Esquina superior Izquierda
        maxX = ncols-1-halfSizeWindow, maxY = nrows-1- halfSizeWindow;  // Esquina inferior Derecha
    unsigned char **ucmat;                      //Matriz Binarizada de salida.
    ucmat=ucmatrix(0,nrows-1,0,ncols-1);        // Se inicializa la matriz.

    // Se calcula la mediana
    double  sumPixelsWindow=0;                  // Suma de los pixeles en la vecindad.
    double  sum2PixelsWindow=0;                  // Suma de los cuadrados de la vecindad.
    double  localMean=0, localVar=0, localStd=0; // Valor Medio, Varianza, desviacion Standard
    double  localValue=0;
    double  mainSumPixelWindow=0,mainSum2PixelWindow=0;
    int    i=0,j=0,k=0,l=0;

    double  sumCols[ncols],sum2Cols[ncols];
   
    for(j=0;j<ncols;j++){
        sumCols[j]  = 0;
        sum2Cols[j] = 0;
    }
    //  Se acumulan los valores las columnas.
    for (i=0;i<sizeWindow;i++){
            for (j=0;j<ncols;j++){
                    sumCols[j]  += mat[i][j];
                    sum2Cols[j] += mat[i][j]*mat[i][j];
            }
    }

    //        Calcular el mainSumPixelWindow
    for(j=0;j<sizeWindow;j++){
            mainSumPixelWindow        += sumCols[j];
            mainSum2PixelWindow        += sum2Cols[j];
    }
    for(i=minY; i<=maxY; i++){
        for(j=minX; j<=maxX; j++){

            // Se calcula el sumFirstCol y sumLastCol
            if(i==minX){
                if(j==minY){        //******* i==minX & j==minY

                    sumPixelsWindow        =  mainSumPixelWindow;
                    sum2PixelsWindow        =  mainSum2PixelWindow;
                }else{              //******* i==minX & j=!minY
                   
                    sumPixelsWindow        += (sumCols[j+halfSizeWindow]  - sumCols[j-halfSizeWindow-1]);
                    sum2PixelsWindow        += (sum2Cols[j+halfSizeWindow]  - sum2Cols[j-halfSizeWindow-1]);
                }
            }else{

                if(j==minY){        //******* i!=minX & j==minY
                    // Actualizar el valor de sumPixelsCols
                    for (k=0;k<ncols;k++){
                   
                      sumCols[k]  += (mat[i+halfSizeWindow][k] - mat[i-halfSizeWindow-1][k]);
                      sum2Cols[k] += (mat[i+halfSizeWindow][k]*mat[i+halfSizeWindow][k] - mat[i-halfSizeWindow-1][k]*mat[i-halfSizeWindow-1][k]);
                    }
                    //        Calcular el mainSumPixelWindow
                    mainSumPixelWindow        = 0;
                    mainSum2PixelWindow = 0;
                    for(k=0;k<sizeWindow;k++){
                            mainSumPixelWindow        += sumCols[k];
                            mainSum2PixelWindow        += sum2Cols[k];
                    }

                    sumPixelsWindow        = mainSumPixelWindow;
                    sum2PixelsWindow        = mainSum2PixelWindow;
                }else{              //******* i!=minX & j!=minY
                    sumPixelsWindow  += (sumCols[j+halfSizeWindow] - sumCols[j-halfSizeWindow-1]);
                    sum2PixelsWindow += (sum2Cols[j+halfSizeWindow]  - sum2Cols[j-halfSizeWindow-1]);
                }
            }
       
            localMean  = sumPixelsWindow/(sizeWindow*sizeWindow);  //Se halla la media.
            localVar    = sum2PixelsWindow/(sizeWindow*sizeWindow) - localMean*localMean;  // Se calcula la varianza.
            localStd    = sqrt(fabs(localVar));          // Se calcula la desviacion estandard.

            localValue = localMean+weight*localStd;
            // printf("mat[%d][%d] = %f\tlocalValue = %f\tsumPixelsWindow= %f\tlocalMean =  %f\tlocalVar = %f\tlocalStd = %f\n",i,j,mat[i][j],localValue,sumPixelsWindow,localMean,localVar,localStd);
            // printf("i=%d\tj=%d\t%f\t%f\t%f\n",i,j,mat[i][j],localValue,sumPixelsWindow);


            if(mat[i][j]< localValue){
                ucmat[i][j] = 0;
            }else{
                ucmat[i][j] = 1;
            }
        }
    }

    for(i=0;i<nrows;i++){
        for(j=0;j<ncols;j++){
            if(i<minY) ucmat[i][j] = 1;
            if(i>maxY) ucmat[i][j] = 1;
            if(j<minX) ucmat[i][j] = 1;
            if(j>maxX) ucmat[i][j] = 1;
        }

    }

   

    printf("Binarizacion culminada\n");
    return (ucmat);

} // end sub im2bw_niblack()

---

Thanks for help me
bests regards

Looks good, though I would add a few more comments in to explain what you're doing. Anyway, is it sufficiently fast now? :)

You should realize too that a program like matlab will have so many optimizations in it that you probably won't be able to achieve its performance . . . but why should that stop you from trying?

Yes is really fast, for my images takes 0.12 seconds. Oops sorry for the comments but my native language is spanish. ;)

Thanks for the help and comments.
bye.