图像去噪程序的图像处理例程

Question

我编写了一个多语言的软件来实现一种图像去噪算法(由许多子算法组成)，称为蜡烛。我有一个C程序，这是软件的一部分，并希望确保它是充分优化和安全。

这是一个简单的C程序(没有多线程)。我查过了，一切都很好。但是我主要是一个Java程序员，我想确保我没有做任何只是还没有出现的错误。

我不得不将访问数组元素的速度(通过堆栈分配)与空间(通过堆分配)进行交换，因为图像可能非常大。

下面是节目。父方法位于底部(“估计”)。有投入吗？

#include "math.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>


int lp2 (unsigned int x) {
    if (x == 0) return 0;
    unsigned int result = 1;
    while ((result < x) && (result != 0))
        result <<= 1;
    return (int) result;
}





float minimum(float *A , int size){

    float min = A[0];

    for (int i =1; i < size; i++) {

        if (A[i] < min) {

            min = A[i];

        }

    }

    return min;


}



void quick_sort (float *a, int n) {
    int i, j;
    float p, t;
    if (n < 2)
        return;
    p = a[n / 2];
    for (i = 0, j = n - 1;; i++, j--) {
        while (a[i] < p)
            i++;
        while (p < a[j])
            j--;
        if (i >= j)
            break;
        t = a[i];
        a[i] = a[j];
        a[j] = t;
    }
    quick_sort(a, i);
    quick_sort(a + i, n - i);
}





float medianSmall(float *A, int size){

    float med;
    float Sorted[size];

    memcpy(Sorted , A, sizeof(float) * size ) ;

    quick_sort(Sorted , size);

    if (size % 2) {

        med = Sorted[size/2] ;

    }else{

        med = ( Sorted[size/2] + Sorted[(size/2) - 1] ) / 2 ;
    }
    return med;
}


float medianLarge(float *A, int size){

    float *Sorted, med;

    Sorted = (float*)malloc(sizeof(float) * size );

    memcpy(Sorted , A, sizeof(float) * size ) ;

    quick_sort(Sorted , size);

    if (size % 2) {

        med = Sorted[size/2] ;

    }else{

        med = ( Sorted[size/2] + Sorted[(size/2) - 1] ) / 2 ;

    }

    free(Sorted);
    return med;

}


float * cshift3D(float*x , int N1 , int N2 , int N3 ){
    int i, j , k, counter;
    float *y;

    int n[N1];

    y = (float*)malloc(sizeof(float) * ( N1 * N2 * N3 ) );

    for (i = 0; i < N1; i++) {

        n[i] = (i + 5 ) % N1;
    }

    counter = 0;
    for (k=0; k < N3; k++){
        for (i=0; i < N1; i++){
            for (j=0; j < N2; j++){
                y[counter] = x[(k*N1*N2) + (n[i]*N2) + j];
                counter++;

            }

        }

    }

    free(x);
    return y;

}


float * permute(float*A , int rows ,int cols ,int slices ,int*p ,int per_or_iper){

    int ii, jj , kk, rows_final , cols_final , slices_final;
    float *B;

    int A_dim[3] = {rows, cols , slices};

    int B_dim[3];

    B = (float*)malloc(sizeof(float) * ( rows*cols*slices ) );

    if (per_or_iper) {

        B_dim[0] = A_dim[p[0]] ;
        B_dim[1] = A_dim[p[1]] ;
        B_dim[2] = A_dim[p[2]] ;

        rows_final = rows;
        cols_final = cols;
        slices_final = slices;



    }else{

        B_dim[p[0]] = A_dim[0];
        B_dim[p[1]] = A_dim[1];
        B_dim[p[2]] = A_dim[2];

        rows_final = B_dim[0];
        cols_final = B_dim[1];
        slices_final = B_dim[2];

    }

    int ind_val[3];
    int ind[3];
    for (kk=0; kk<slices_final; kk++){
        ind[2] = kk;
        for (ii=0; ii<rows_final; ii++){
            ind[0] = ii;
            for (jj=0; jj<cols_final; jj++){
                ind[1] = jj;
                ind_val[0] = ind[p[0]];
                ind_val[1] = ind[p[1]];
                ind_val[2] = ind[p[2]];

                if (per_or_iper) {

                    B[ind_val[2]*B_dim[0]*B_dim[1] + ind_val[0]*B_dim[1] + ind_val[1]] = A[kk*rows*cols + ii*cols + jj];

                }else{

                    B[kk*B_dim[0]*B_dim[1] + ii*B_dim[1] + jj] = A[ind_val[2]*rows*cols + ind_val[0]*cols + ind_val[1]];

                }

            }

        }

    }
    free(A);
    return B;
}


float * convn(float* xin , int rows, int cols){
    int outrows, temprows, count, counter, x , y , z , i , j;
    float s;
    float *out , *temporary;

    float hpf[10] = {0 , 0 , -0.08838834764832 , -0.08838834764832 , 0.69587998903400 , -0.69587998903400, 0.08838834764832, 0.08838834764832, 0.01122679215254 , -0.01122679215254} ;

    temprows = rows + 9;

    temporary = (float*)malloc(sizeof(float) * ( temprows * cols ) );

    count = 0;

    for (y = -9 ; y < rows; y++) {
        for (x = 0; x < cols; x++) {

            s = 0;
            counter = 0;

            for (z = y; z < (y+10); z++) {

                if (z < 0) {
                    counter++;
                    continue;
                }
                else if(z >= rows){
                    counter++;
                    continue;
                }else{

                    s += ( hpf[counter++] * xin[z*cols+x] ) ;

                }

            }

            temporary[count++] = s;
        }

    }

    count = 0;

    outrows = (temprows/2) + 1 ;

    out = (float*)malloc(sizeof(float) * ( outrows * cols ) );

    for (i = 0; i < temprows; i += 2) {
        for (j = 0; j < cols; j++) {
            out[count++] = temporary[(i*cols) + j];
        }
    }

    free(temporary);
    return out;  
}


float * afb3D_A(float*x , int d, int xx, int yy, int zz){

    int L, i, N1, N2, N3 , zloc, xloc, yloc, counter;

    float *perm,  *cshif, *iperm, *hi , *xTemp;


    int p[3];
    for(i = 0; i < 3; i++ ){

        p[i] = ((d-1)+i) % 3 ; 
    }

    perm = permute(x , xx , yy , zz , p , 1);

    L = 5;

    if (d == 1) {
        N1 = xx;
        N2 = yy;
        N3 = zz;

    }
    else if(d == 2){
        N1 = yy;
        N2 = zz;
        N3 = xx;
    }else{

        N1 = zz;
        N2 = xx;
        N3 = yy;
    }

    cshif = cshift3D(perm , N1, N2 , N3);

    hi = (float*)malloc(sizeof(float) * ( (L+(N1/2) )* N2 * N3 ) );
    xTemp = (float*)malloc(sizeof(float) * ( N1 * N2 ) );
    float *hiTemp;

    for (zloc = 0; zloc < N3; zloc++) {
        counter = 0;
        for (xloc = 0; xloc <  N1; xloc++) {
            for (yloc = 0; yloc < N2; yloc++) {
                xTemp[counter++] = cshif[zloc*N1*N2 + xloc*N2 + yloc];
            }
        }
        hiTemp = convn(xTemp, N1, N2);
        for (xloc = 0; xloc <  (L+(N1/2)); xloc++) {
            for (yloc = 0; yloc < N2; yloc++) {
                hi[zloc*(L+N1/2)*N2 + xloc*N2 + yloc] = hiTemp[xloc*N2 + yloc];
            }
        }

        free(hiTemp);    
    }

    free(xTemp);
    free(cshif);


    hiTemp = (float*)malloc(sizeof(float) * ( (L+(N1/2) )* N2 * N3) );
    memcpy(hiTemp , hi , sizeof(float) * ( (L+(N1/2) )* N2 * N3 ) );


    for (zloc = 0; zloc < N3; zloc++){
        for (xloc = 0; xloc <  L; xloc++){
            for (yloc = 0; yloc < N2; yloc++){
                hiTemp[zloc*(L+N1/2)*N2 + xloc*N2 + yloc] += hi[zloc*(L+N1/2)*N2 + (xloc + (N1/2))*N2 + yloc];
           }
        }
    }


    hi = (float*)realloc(hi , sizeof(float) * ( (N1/2) * N2 * N3) );
    for (zloc = 0; zloc < N3; zloc++){
        for (xloc = 0; xloc <  (N1/2); xloc++){
            for (yloc = 0; yloc < N2; yloc++){
                hi[zloc*(N1/2)*N2 + xloc*N2 + yloc] = hiTemp[zloc*(L+N1/2)*N2 + xloc*N2 + yloc];
            }
        }
    }

    free(hiTemp);

    iperm = permute(hi, (N1/2), N2 , N3 , p , 0);

    return iperm;

}



void pad2d(float *arr , float *newarr , int x , int y , int axis){

    int xx, yy , rowcoun, colcoun, i , row , j , col;


    xx = x + 2*axis;
    yy = y + 2*axis;

    rowcoun = 0;

    for(i=0; i < xx; i++){
        colcoun = 0;
        if(i <= (axis - 1)){

            row = axis - rowcoun;
            row = row % x;

            if (row){

                row -= 1;
                row = (x-1) - row;

            }
            rowcoun++;


        }
        else if(i <= (x+axis-1) ){

            row = i - axis;

        }else{

            row = rowcoun  - axis + 1 ;
            row = row % x ;

            if(row){

                row -= 1;

            }else{

                row = x - 1;
            }

            rowcoun++;
        }

        for(j=0; j < yy; j++){



            if(j <= (axis-1) ){

                col = axis - colcoun ;
                col = col % y ;


                if(col){

                    col -= 1;
                    col = (y-1) - col;

                }

                colcoun++;
            }
            else if(j <= (y+axis-1)){

                col = j - axis;


            }else{

                col = colcoun  - axis + 1;
                col = col % y ;
                if(col){

                    col -= 1;

                }else{

                    col = y - 1;
                }

                colcoun++;
            }


            newarr[yy*i + j] = arr[row*y + col];


        }


    }


}


void medfilt2(float *A, float *B, int rows, int cols, int axis){

    int winelem , ii , jj , xx , yy , inc, ind;



    winelem =  ( (2*axis) + 1) * ( (2*axis) + 1) ;

    float window[winelem];


    for(ii = 0; ii < rows; ii++){

        for (jj=0; jj < cols; jj++) {

            inc = 0;

            for (xx = 0; xx < ( (2*axis) + 1); xx++) {
                for (yy = 0; yy < ( (2*axis) + 1); yy++) {

                    ind = ( (ii+xx)* ((2*axis) + cols ) + (jj + yy)) ;
                    window[inc++] = A[ind];
                }
            }

            B[ii*cols + jj] = medianSmall(window , winelem);

        }


    }


}


void estimate(float*ima , int x , int y , int z , int ps, float **HHH){

    int size , xx , yy , zz , i, j, k, p1 , p2 , p3, counter , z_half , x_half, y_half , padx , pady, indexx, zi, xi, yi, val ;
    float minim, Sig;

    float *filt1, *filt2, *filt3, *padarray, *temp, *NsigMAP , *img2d , *img2dp;

    size = x*y*z;

    minim = minimum(ima , size);

    if (minim < 0) {

        for(i=0; i<size; i++){

            ima[i] = ima[i] - minim;
        }

    }


    p1 = lp2( (unsigned int)x ) ;

    p2 = lp2( (unsigned int)y ) ;

    p3 = lp2( (unsigned int)z ) ;

    // Make the image dimesnions powers of 2
    if(p1 == x & p2==y & p3 == z){

        padarray = (float*)malloc(sizeof(float) * (p1*p2*p3));
        memcpy(padarray , ima , sizeof(float) * (p1*p2*p3));


    }else{

        padarray = (float*)calloc((p1*p2*p3) , sizeof(float));

        for(k = 0; k < z; k++ ){

            for(i=0; i < x; i++){

                for(j=0; j < y; j++){

                    padarray[(k*p1*p2)+(i*p2)+ j] = ima[(k*x*y) + (i*y) + j];
                }
            }
        }
    }

    // Filter along dimension 1

    filt1 = afb3D_A(padarray, 1 , p1, p2, p3);

    p1 = p1/2;

    // Filter along dimension 2

    filt2 = afb3D_A(filt1, 2 , p1, p2, p3);

    p2 = p2/2 ;

    // Filter along dimension 3

    filt3 = afb3D_A(filt2, 3 , p1, p2, p3);

    p3 = p3/2;

    z_half = z/2 + (z % 2 != 0);
    x_half = x/2 + (x % 2 != 0);
    y_half = y/2 + (y % 2 != 0);

    // Remove Regions corresponding to zero padding
    temp = (float*)malloc(sizeof(float) * ( z_half*x_half*y_half ) );

    counter = 0;
    for(k=0; k < z_half; k++){
        for (i = 0; i < x_half; i++){
            for (j=0; j < y_half; j++){

                temp[counter++] = fabsf( filt3[k*p1*p2 + i*p2 + j] ) / 0.6745;

            }

        }
    }


    free(filt3);
    Sig = medianLarge(temp , (z_half*x_half*y_half) );
    printf("Sig:\n");
    printf("%.6f", Sig);
    padx = x_half + 2*ps;
    pady = y_half + 2*ps;


    img2d = (float*)malloc(sizeof(float) * ( x_half*y_half ) );
    img2dp = (float*)malloc(sizeof(float) * ( padx * pady ) );


    NsigMAP = (float*)malloc(sizeof(float) * ( z_half * x_half * y_half ) ) ;


    for(k=0; k < z_half; k++){
        counter = 0;
        for(i=0; i < x_half; i++ ){
            for(j=0; j < y_half; j++){

                // Get the kth Slice
                img2d[counter++] = temp[k*x_half*y_half + i*y_half + j];

            }

        }


        memset(img2dp,0, sizeof(float) * (padx*pady) );


        // Pad the 2d image
        pad2d(img2d , img2dp , x_half , y_half, ps);

        // Apply 2d median filter
        medfilt2(img2dp , img2d,  x_half , y_half , ps);

        indexx = k*x_half*y_half ;

        memcpy( &NsigMAP[indexx] , img2d, sizeof(float) * (x_half*y_half) );

    }

    free(temp);
    free(img2d);
    free(img2dp);

    *HHH = (float*)malloc(sizeof(float) * ( z * x * y ) );


    // 3-d interpolation

    for (k=0; k < z; k++) {

        zi = k/2;

        for (i=0; i < x; i++) {

            xi = i/2;

            for (j=0; j < y; j++) {

                yi = j/2;

                if (k && i && j) {

                    (*HHH)[k*x*y + y*i + j] = NsigMAP[zi*x_half*y_half + xi*y_half + yi] ;

                }else{

                    (*HHH)[k*x*y + y*i + j] = NAN;

                }

            }
        }
    }


    free(NsigMAP);

    for (val = 0; val < (x*y*z); val++) {

        if ((*HHH)[val] < Sig) {
            (*HHH)[val] = Sig;
        }
    }


}

Answer 1

在这里看到图像处理的东西总是很酷的！看上去你在做一些有趣的事情。这里有很多代码，所以我只想给大家一些我认为可以改进的东西的一般想法。

给出有用的名称

您的许多函数和变量都没有意义，或者很难理解名称。例如，lp2()没有告诉我太多。从几分钟的阅读中我可以看出，它是一个整数日志基2函数。所以也许叫它intLogBase2()？

另外，您有medfilt2()函数，我猜想它可以做一个中值滤波。为什么是"2"？是因为它是二维的吗？它是否替换了一个medfilt()函数？如果是这样的话，我没有看到它，所以为什么不叫它medfilt()呢？

将const用于不改变

的事物

大多数函数参数都是只读的。这个函数根本不改变它们。在这种情况下，您应该将它们声明为const。它让编译器和任何阅读代码的人知道，您不打算让该变量更改函数。

--不要重新发明车轮

您编写了自己的quick_sort()方法。经验告诉我，它很可能会有错误。(我还没有花时间对其进行分析以确定，但我知道每次我尝试亲自实现它时，我最终都会搞砸一些只在发货后才出现的模糊情况！)改用标准库函数qsort()。

但事实证明你不需要整理任何东西..。

有一种更快的计算中值

的方法

事实证明，floats的设计是为了让你把它们的位模式当作整数，然后比较它们，你就会得到正确的答案。你可以利用这个对你有利。(我会在你正在做的函数的注释中记录它，因为它不是显而易见的。)

您可以计算floats数组的直方图，就像计算整数一样。然后你就可以使用直方图来查找中间值。了。

其中一个诀窍是进行多次传球。例如，将每个float转换为32位int。然后找到中间值的高字节。接下来，再次运行直方图。但这一次，如果高字节小于中位数的高字节，只需将其计算为“小于”中位数。如果高字节大于中位数，只需忽略它。如果高字节正好等于中位数，那么把它放在你的256个垃圾箱中的一个。最后，通过添加“小于”值的数目加上所有回收箱的总和，您将知道中位数的最高字节和第二大字节。对第二个最低字节和最低字节也这样做。

或者，如果你对一个近似没有意见，把这些值放入一个固定数量的垃圾箱中(比如1000，或者别的什么)，然后对含有第50个百分位数的回收箱中的值进行排序或运行直方图。

我假设你的convn()函数是卷积函数？卷积核是可分离的吗？如果是这样的话，你可以做3D通过它，并节省了相当多的时间。如果没有，你能用FFT转换成频率空间并在那里进行卷积吗？(频率空间中的卷积只是一个乘积，所以速度更快。)

你需要分配这么多内存吗？

您的函数分配并释放了大量内存。这些任务中的任何一项都能到位吗？或者，如果缓冲区的大小经常是相同的(或者是整个图像的大小，或者是图像块的大小)，那么您可以使用大小合适的预分配缓冲区池吗？

Answer 2

堆栈与堆

我只想就一件事发表评论。您有两个函数，medianSmall()和medianLarge()，它们是相同的，只不过第一个函数在堆栈上分配数组，第二个函数在堆上分配数组。不需要两个函数。您应该只使用从堆中分配的那个。在堆栈上放置数组没有速度优势。堆版本更可取，因为这样不会出现堆栈溢出问题。