Skip to content

/CANDLE-J

Permalink
Branch: master
Find file Copy path

CANDLE-J/ONLM_Mod.c

Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
564 lines (429 sloc) 16.2 KB
Raw Blame History
/* Pierrick Coupe - pierrick.coupe@gmail.com */
/* Jose V. Manjon - jmanjon@fis.upv.es */
/* Brain Imaging Center, Montreal Neurological Institute. */
/* Mc Gill University */
/* */
/* Copyright (C) 2010 Pierrick Coupe and Jose V. Manjon */
/* Details on ONLM filter */
/***************************************************************************
* The ONLM filter is described in: *
* *
* P. Coupe, P. Yger, S. Prima, P. Hellier, C. Kervrann, C. Barillot. *
* An Optimized Blockwise Non Local Means Denoising Filter for 3D Magnetic*
* Resonance Images. IEEE Transactions on Medical Imaging, 27(4):425-441, *
* April 2008 *
***************************************************************************/
/***************************************************************************
* This file was modified by Haider Khan - haiderriazkhan@hotmail.com *
* *
* The ONLM C program can now be called from Java via JNA *
***************************************************************************/
#include "math.h"
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
/* Multithreading stuff*/
#ifdef _WIN32
#include <windows.h>
#include <process.h>
#else
#include <pthread.h>
#endif
#include <stdbool.h>
typedef struct{
int rows;
int cols;
int slices;
float * in_image;
float * ref_image;
float * means_image;
float * var_image;
float * estimate;
unsigned short * label;
int ini;
int fin;
int radioB;
int radioS;
float *sigma;
float beta;
}myargument;
float max;
/* Function which compute the weighted average for one block */
void Average_block(float *ima,int x,int y,int z,int neighborhoodsize,float *average, float weight, int sx,int sy,int sz)
{
int x_pos,y_pos,z_pos;
bool is_outside;
int a,b,c,ns,sxy,count;
ns=2*neighborhoodsize+1;
sxy=sx*sy;
count = 0;
for (c = 0; c<ns;c++)
{
z_pos = z+c-neighborhoodsize;
for (b = 0; b<ns;b++)
{
y_pos = y+b-neighborhoodsize;
for (a = 0; a<ns;a++)
{
x_pos = x+a-neighborhoodsize;
is_outside = false;
if ((x_pos < 0) || (x_pos > sx-1)) is_outside = true;
if ((y_pos < 0) || (y_pos > sy-1)) is_outside = true;
if ((z_pos < 0) || (z_pos > sz-1)) is_outside = true;
if (is_outside)
average[count] = average[count] + ima[z*(sxy)+(y*sx)+x]*weight;
else
average[count] = average[count] + ima[z_pos*(sxy)+(y_pos*sx)+x_pos]*weight;
count++;
}
}
}
}
/* Function which computes the value assigned to each voxel */
void Value_block(float *Estimate, unsigned short *Label,int x,int y,int z,int neighborhoodsize,float *average, float global_sum, int sx,int sy,int sz)
{
int x_pos,y_pos,z_pos;
int ret;
bool is_outside;
float value = 0.0;
unsigned short label = 0;
int count=0 ;
int a,b,c,ns,sxy;
extern bool rician;
ns=2*neighborhoodsize+1;
sxy=sx*sy;
for (c = 0; c<ns;c++)
{
z_pos = z+c-neighborhoodsize;
for (b = 0; b<ns;b++)
{
y_pos = y+b-neighborhoodsize;
for (a = 0; a<ns;a++)
{
x_pos = x+a-neighborhoodsize;
is_outside = false;
if ((x_pos < 0) || (x_pos > sx-1)) is_outside = true;
if ((y_pos < 0) || (y_pos > sy-1)) is_outside = true;
if ((z_pos < 0) || (z_pos > sz-1)) is_outside = true;
if (!is_outside)
{
value = Estimate[z_pos*(sxy)+(y_pos*sx)+x_pos];
value = value + (average[count]/global_sum);
label = Label[(x_pos + y_pos*sx + z_pos*sxy)];
Estimate[z_pos*(sxy)+(y_pos*sx)+x_pos] = value;
Label[(x_pos + y_pos*sx + z_pos *sxy)] = label +1;
}
count++;
}
}
}
}
float distance(float* ima,int x,int y,int z,int nx,int ny,int nz,int f,int sx,int sy,int sz)
{
float d,acu,distancetotal;
int i,j,k,ni1,nj1,ni2,nj2,nk1,nk2;
distancetotal=0;
for(k=-f;k<=f;k++)
{
nk1=z+k;
nk2=nz+k;
if(nk1<0) nk1=-nk1;
if(nk2<0) nk2=-nk2;
if(nk1>=sz) nk1=2*sz-nk1-1;
if(nk2>=sz) nk2=2*sz-nk2-1;
for(j=-f;j<=f;j++)
{
nj1=y+j;
nj2=ny+j;
if(nj1<0) nj1=-nj1;
if(nj2<0) nj2=-nj2;
if(nj1>=sy) nj1=2*sy-nj1-1;
if(nj2>=sy) nj2=2*sy-nj2-1;
for(i=-f;i<=f;i++)
{
ni1=x+i;
ni2=nx+i;
if(ni1<0) ni1=-ni1;
if(ni2<0) ni2=-ni2;
if(ni1>=sx) ni1=2*sx-ni1-1;
if(ni2>=sx) ni2=2*sx-ni2-1;
distancetotal = distancetotal + ((ima[nk1*(sx*sy)+(nj1*sx)+ni1]-ima[nk2*(sx*sy)+(nj2*sx)+ni2])*(ima[nk1*(sx*sy)+(nj1*sx)+ni1]-ima[nk2*(sx*sy)+(nj2*sx)+ni2]));
}
}
}
acu=(2*f+1)*(2*f+1)*(2*f+1);
d=distancetotal/acu;
return d;
}
void* ThreadFunc( void* pArguments )
{
float *Estimate,*ima,*means,*variances,*sigma,*ref,*average,epsilon,totalweight,wmax,t1,d,w;
float beta;
unsigned short *Label;
int rows,cols,slices,ini,fin,v,f,i,j,k,rc,ii,jj,kk,ni,nj,nk,Ndims;
myargument arg;
arg=*(myargument *) pArguments;
rows=arg.rows;
cols=arg.cols;
slices=arg.slices;
ini=arg.ini;
fin=arg.fin;
ima=arg.in_image;
ref=arg.ref_image;
means=arg.means_image;
variances=arg.var_image;
Estimate=arg.estimate;
Label=arg.label;
v=arg.radioB;
f=arg.radioS;
sigma=arg.sigma;
beta= arg.beta;
epsilon = 0.0000001;
rc=rows*cols;
Ndims = (2*f+1)*(2*f+1)*(2*f+1);
average=(float*)malloc(Ndims*sizeof(float));
for(k=ini;k<fin;k+=2)
{
for(j=0;j<rows;j+=2)
{
for(i=0;i<cols;i+=2)
{
memset(average,0.0, sizeof(float) * Ndims );
totalweight=0.0;
wmax=0.0;
for(kk=-v;kk<=v;kk++)
{
for(jj=-v;jj<=v;jj++)
{
for(ii=-v;ii<=v;ii++)
{
ni=i+ii;
nj=j+jj;
nk=k+kk;
if(ii==0 && jj==0 && kk==0) continue;
if(ni>=0 && nj>=0 && nk>=0 && ni<cols && nj<rows && nk<slices)
{
t1 = ((2 * means[k*rc+(j*cols)+i] * means[nk*rc+(nj*cols)+ni]+ epsilon) / ( means[k*rc+(j*cols)+i]*means[k*rc+(j*cols)+i] + means[nk*rc+(nj*cols)+ni]*means[nk*rc+(nj*cols)+ni] + epsilon) ) * ((2 * sqrt(variances[k*rc+(j*cols)+i]) * sqrt(variances[nk*rc+(nj*cols)+ni]) + epsilon) / (variances[k*rc+(j*cols)+i] + variances[nk*rc+(nj*cols)+ni] + epsilon));
if(t1 > 0.9)
{
d=distance(ref,i,j,k,ni,nj,nk,f,cols,rows,slices);
if (d<=0) d = epsilon;
w = exp(-d/(beta*sigma[k*(rc)+(j*cols)+i]*sigma[k*(rc)+(j*cols)+i]));
if(w>wmax) wmax = w;
Average_block(ima,ni,nj,nk,f,average,w,cols,rows,slices);
totalweight = totalweight + w;
}
}
}
}
}
if(wmax==0.0) wmax=1.0;
Average_block(ima,i,j,k,f,average,wmax,cols,rows,slices);
totalweight = totalweight + wmax;
if(totalweight != 0.0)
Value_block(Estimate,Label,i,j,k,f,average,totalweight,cols,rows,slices);
}
}
}
free(average);
#ifdef _WIN32
_endthreadex(0);
return 0;
#else
pthread_exit(0);
return 0;
#endif
}
void cleanup(float* pVals)
{
free(pVals);
}
void ONLM(float* ima , int v , int f , float* sigma , float beta , float* ref , int rows, int cols, int slices, float** fima)
{
float *means, *variances, *Estimate;
unsigned short *Label;
float w,totalweight,wmax,d,mean,var,t1,t2,epsilon,label,estimate;
int i,j,k,ii,jj,kk,ni,nj,nk,indice,Nthreads,ini,fin;
myargument *ThreadArgs;
#ifdef _WIN32
HANDLE *ThreadList; // Handles to the worker threads
#else
pthread_t * ThreadList;
#endif
const int dims[3] = {cols, rows, slices};
int totalsize = (rows*cols*slices);
*fima = (float*)malloc(sizeof(float) * totalsize);
means = (float*)malloc(sizeof(float) * totalsize);
variances = (float*)malloc(sizeof(float) * totalsize);
Estimate = (float*)malloc(sizeof(float) * totalsize);
Label = (unsigned short*)malloc(sizeof(unsigned short) * totalsize);
for (i = 0; i < dims[2] *dims[1] * dims[0];i++)
{
Estimate[i] = 0.0;
Label[i] = 0.0;
(*fima)[i] = 0.0;
}
max=0;
for(k=0;k<dims[2];k++)
{
for(i=0;i<dims[1];i++)
{
for(j=0;j<dims[0];j++)
{
mean=0;
indice=0;
for(ii=-1;ii<=1;ii++)
{
for(jj=-1;jj<=1;jj++)
{
for(kk=-1;kk<=1;kk++)
{
ni=i+ii;
nj=j+jj;
nk=k+kk;
if(ni<0) ni=-ni;
if(nj<0) nj=-nj;
if(nk<0) nk=-nk;
if(ni>=dims[1]) ni=2*dims[1]-ni-1;
if(nj>=dims[0]) nj=2*dims[0]-nj-1;
if(nk>=dims[2]) nk=2*dims[2]-nk-1;
mean = mean + ref[nk*(dims[0]*dims[1])+(ni*dims[0])+nj];
indice=indice+1;
}
}
}
mean=mean/indice;
means[k*(dims[0]*dims[1])+(i*dims[0])+j]=mean;
if (mean > max) max =mean;
}
}
}
for(k=0;k<dims[2];k++)
{
for(i=0;i<dims[1];i++)
{
for(j=0;j<dims[0];j++)
{
var=0;
indice=0;
for(ii=-1;ii<=1;ii++)
{
for(jj=-1;jj<=1;jj++)
{
for(kk=-1;kk<=1;kk++)
{
ni=i+ii;
nj=j+jj;
nk=k+kk;
if(ni>=0 && nj>=0 && nk>0 && ni<dims[1] && nj<dims[0] && nk<dims[2])
{
var = var + (ref[nk*(dims[0]*dims[1])+(ni*dims[0])+nj]-means[k*(dims[0]*dims[1])+(i*dims[0])+j])*(ref[nk*(dims[0]*dims[1])+(ni*dims[0])+nj]-means[k*(dims[0]*dims[1])+(i*dims[0])+j]);
indice=indice+1;
}
}
}
}
var=var/(indice-1);
variances[k*(dims[0]*dims[1])+(i*dims[0])+j]=var;
}
}
}
Nthreads=8;
#ifdef _WIN32
// Reserve room for handles of threads in ThreadList
ThreadList = (HANDLE*)malloc(Nthreads* sizeof( HANDLE ));
ThreadArgs = (myargument*) malloc( Nthreads*sizeof(myargument));
for (i=0; i<Nthreads; i++)
{
/* Make Thread Structure*/
ini=(i*dims[2])/Nthreads;
fin=((i+1)*dims[2])/Nthreads;
ThreadArgs[i].cols=dims[0];
ThreadArgs[i].rows=dims[1];
ThreadArgs[i].slices=dims[2];
ThreadArgs[i].in_image=ima;
ThreadArgs[i].ref_image=ref;
ThreadArgs[i].var_image=variances;
ThreadArgs[i].means_image=means;
ThreadArgs[i].estimate=Estimate;
ThreadArgs[i].label=Label;
ThreadArgs[i].ini=ini;
ThreadArgs[i].fin=fin;
ThreadArgs[i].radioB=v;
ThreadArgs[i].radioS=f;
ThreadArgs[i].sigma=sigma;
ThreadArgs[i].beta=beta;
ThreadList[i] = (HANDLE)_beginthreadex( NULL, 0, &ThreadFunc, &ThreadArgs[i] , 0, NULL );
}
for (i=0; i<Nthreads; i++) { WaitForSingleObject(ThreadList[i], INFINITE); }
for (i=0; i<Nthreads; i++) { CloseHandle( ThreadList[i] ); }
#else
/* Reserve room for handles of threads in ThreadList*/
ThreadList = (pthread_t *) calloc(Nthreads,sizeof(pthread_t));
ThreadArgs = (myargument*) calloc( Nthreads,sizeof(myargument));
for (i=0; i<Nthreads; i++)
{
/* Make Thread Structure*/
ini=(i*dims[2])/Nthreads;
fin=((i+1)*dims[2])/Nthreads;
ThreadArgs[i].cols=dims[0];
ThreadArgs[i].rows=dims[1];
ThreadArgs[i].slices=dims[2];
ThreadArgs[i].in_image=ima;
ThreadArgs[i].ref_image=ref;
ThreadArgs[i].var_image=variances;
ThreadArgs[i].means_image=means;
ThreadArgs[i].estimate=Estimate;
ThreadArgs[i].label=Label;
ThreadArgs[i].ini=ini;
ThreadArgs[i].fin=fin;
ThreadArgs[i].radioB=v;
ThreadArgs[i].radioS=f;
ThreadArgs[i].sigma=sigma;
ThreadArgs[i].beta=beta;
}
for (i=0; i<Nthreads; i++)
{
if(pthread_create(&ThreadList[i], NULL, ThreadFunc,&ThreadArgs[i]))
{
printf("Threads cannot be created\n");
exit(1);
}
}
for (i=0; i<Nthreads; i++)
{
pthread_join(ThreadList[i],NULL);
}
#endif
free(ThreadArgs);
free(ThreadList);
free(variances);
free(means);
free(sigma);
label = 0.0;
estimate = 0.0;
/* Aggregation of the estimators (i.e. means computation) */
for (k = 0; k < dims[2]; k++ )
{
for (i = 0; i < dims[1]; i++ )
{
for (j = 0; j < dims[0]; j++ )
{
label = Label[k*(dims[0]*dims[1])+(i*dims[0])+j];
if (label == 0.0)
{
// Corrected
(*fima)[k*(dims[0]*dims[1])+(i*dims[0])+j] = ima[k*(dims[0]*dims[1])+(i*dims[0])+j];
}
else
{
estimate = Estimate[k*(dims[0]*dims[1])+(i*dims[0])+j];
estimate = (estimate/label);
(*fima)[k*(dims[0]*dims[1])+(i*dims[0])+j]=estimate;
}
}
}
}
free(Estimate);
free(Label);
}
You can’t perform that action at this time.