# Image processing routines for an image de-noising program

I have written a multi-language software to implement an image denoising algorithm (composed of many sub algorithms) called CANDLE. I have a C program that is part of the software and would like to make sure that it is adequately optimized and safe.

It is a simple C program (no multithreading). I have ran it and everything is fine. But I am a Java programmer primarily and I want to make sure I am not doing anything wrong that just simply has not cropped up yet.

I had to trade speed of accessing array elements (by stack allocating) with space (by heap allocating) since the images can be very large.

Below is the program. The parent method is at the bottom ('estimate'). Any inputs?

#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){

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

}else{

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);

img2d = (float*)malloc(sizeof(float) * ( x_half*y_half ) );

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];

}

}

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;
}
}

}

• A small tip/word of advise: when passing the length of the array, or any other object to a function, that argument's type ought to be size_t, – Elias Van Ootegem Jun 20 '15 at 15:20

It's always cool to see image processing stuff here! Looks like some interesting stuff you're doing. There is a lot of code here, so I'll just give you some general ideas of what I think could be improved.

Give Things Useful Names

A lot of your functions and variables have meaningless or hard to understand names. For example lp2() doesn't tell me much. I can see from reading it for a few minutes that it's an integer log base 2 function. So maybe call it intLogBase2()?

Also, you have medfilt2() function that I presume does a median filter. Why "2"? Is it because it's 2D? Is there a medfilt() function that it replaced? If so, I don't see it, so why not just call it medfilt()?

Use const For Things That Don't Change

Most of your functions arguments are read-only. The function doesn't change them at all. When that's the case, you should declare them as const. It lets the compiler and anyone reading the code know that you don't intend for that variable to change in the function.

Don't Reinvent The Wheel

You wrote your own quick_sort() method. Experience tells me that it's likely to have bugs in it. (I haven't spent time analyzing it to see for sure, but I know that every time I try to implement it myself, I end up screwing up some obscure case that only comes up after shipping it!) Use the standard library function qsort() instead.

But it turns out you don't need to sort anything anyway...

There's A Faster Way To Calculate The Median

It turns out that floats were designed so that if you treat their bit patterns as integers and compare them, you'll get the correct answer. You can use this to your advantage. (I'd document in the comments of the function that you're doing it since it's non-obvious.)

You can calculate a histogram for an array of floats just like you would for integers. You can then use the histogram to find the median.

One trick is to do multiple passes. For example, cast each float to a 32-bit int. Then find the high byte of the median. Next, run the histogram again. But this time, if the high byte is less than the high byte of the median, simply count it as "less than" the median. If the high byte is greater than the median, just ignore it. If the high byte exactly equals the median, then put it into one of your 256 bins. At the end, you'll know the highest and second-highest bytes of the median by adding the number of "less than" values plus the sum of all bins up to the median. Do the same for the second-lowest byte and lowest byte.

Or, if you're OK with an approximation, bin the values into some set number of bins (like 1000, or whatever), and then just sort or run a histogram of the ones in the bin containing the 50th percentile.

I assume your convn() function is a convolution function? Is your convolution kernel separable? If so, you could do 3 1D passes over it and save quite a bit of time. If not, could you use an FFT to convert to frequency space and do the convolution there? (A convolution in frequency space is just a multiply, so it's faster.)

Do You Need To Allocate So Much Memory?

Your functions allocate and free a lot of memory. Can any of the tasks be done in place? Alternately, if the buffers are frequently the same size (either the size of the entire image or the size of a tile of the image) perhaps you could use a pool of pre-allocated buffers of the right size?

• Thanks !!! Yes, lp2 is a bad method name. The 2 in medfilt2 comes from the fact it is for 2D images. I probably should have also put a D after the two. I will also put const in the method arguments (ones where they don't change) for clarity as you suggested. Will also look into the alternative way to compute the median . Sounds interesting. – Haider Jun 21 '15 at 12:35
• Also, in terms of not reinventing the wheel, you could look into OpenCV which has a lot of image processing routines built-in and well tested by the community of people using the library. In addition, if you want to improve performance, you can look into doing the work on the GPU via a variety of technologies (OpenGL, OpenCL, CUDA, CoreImage, etc.). – user1118321 Jun 21 '15 at 17:10

## Stack vs heap

I just wanted to comment on one thing. You have two functions, medianSmall() and medianLarge(), which are identical except that first allocates the array on the stack and the second allocates the array on the heap. There's no need for two functions. You should just use the one that allocates from the heap. There's no speed advantage to putting your array on the stack. The heap version is preferable because you don't have stack overflow problems that way.

• I always thought accessing memory from the stack was faster. Which is why I made two functions. One for small arrays and one for large. Are you saying the difference is negligible (or not present at all) and there is no need for two methods in this case? – Haider Jun 21 '15 at 12:39
• @Haider Yes that's what I am saying. To prove it to yourself, you could do a timing test to see whether one way is faster than the other. – JS1 Jun 21 '15 at 16:28
• We should note that allocating the memory on the heap does have a cost over allocating it on the stack. A stack allocation is merely adding some constant amount to a pointer, whereas a heap allocation has to find a block large enough in a list of free blocks. But it's unlikely to make much difference in this case. As JS1 says, measure it and see! – user1118321 Jun 21 '15 at 17:11