Homec4science
Diffusion iCAPs public (master) Browse

iCAPs/kernel.cu
31022b3a13e6master

Last Change
Raw File

kernel.cu
View Options

/*
Author : E. Orliac, SCITAS, EPFL
Date : 08.11.2017
Purpose:
Remarks:
*/
#include "cuda_runtime.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <device_launch_parameters.h>
#include <device_functions.h>
#include <math.h>
#include <thrust/sort.h>
#include <thrust/device_ptr.h>
#include <thrust/execution_policy.h>
#define DB6 6
#define DB3 3
__device__ __constant__ double d_g6[DB6];
__device__ __constant__ double d_h6[DB6];
__device__ __constant__ double d_unity[1];
/*
__inline__ __device__ void reverse_vector(double *data, const int size) {
unsigned j;
double tmp;
for(unsigned i=0; i<(size/2); i++) {
j = size - i - 1;
tmp = data[j];
data[j] = data[i];
data[i] = tmp;
}
}
*/
__device__ void filter(double *out,
double *in, const int n,
double *num, const int lnum,
double *den, const int lden,
const bool reverse) {
int i, j;
double acc;
if (lden == 1 && den[0] == 1.0) {
if (reverse) {
for(i=0; i<n; i++) {
out[n-1-i] = num[0]*in[n-1-i];
for(j=1; j<lnum; j++) {
if (j<=i) {
out[n-1-i] += num[j]*in[n-1-i+j];
}
}
}
} else {
/*
for(i=0; i<n; i++) {
out[i] = 0.0;
for(j=0; j<lnum; j++) {
if (j<=i) {
out[i] += num[j]*in[i-j];
}
}
}
*/
for(i=0; i<n; i+=2) {
out[i] = num[0]*in[i];
out[i+1] = num[0]*in[i+1];
for(j=1; j<lnum; j++) {
if (j<=i) {
out[i] += num[j]*in[i-j];
}
if (j<=i+1) {
out[i+1] += num[j]*in[i-j+1];
}
}
}
}
} else if (lnum == 1 && num[0] == 1.0) {
for(i=0; i<n; i++)
out[i] = 0.0;
if (reverse) {
for(i=0; i<n; i++) {
acc = in[n-1-i];
for(j=1; j<lden; j++) {
if (j<=i) {
acc -= den[j]*out[i-j];
}
}
out[i] = acc;
}
} else {
for(i=0; i<n; i++) {
acc = in[i];
for(j=1; j<lden; j++) {
if (j<=i) {
acc -= den[j]*out[i-j];
}
}
out[i] = acc;
}
}
} else {
printf("Fatal. Unknown case in kernel:filter.\n");
assert(false);
}
}
__device__ void filter_boundary(double *tmp,
double *out,
double *data, const int n,
double *num, const int lnum,
const int lden,
double *den1, const int lden1,
double *den2, const int lden2,
const int condition) {
unsigned j = 0;
double dtmp = 0.0;
if (condition == 1) {
filter(out, data, n, num, lnum, d_unity, 1, true);
} else if (condition == 0) {
filter(out, data, n, num, lnum, d_unity, 1, false);
} else {
printf("Fatal. Unknown case in kernel filter_boundary.\n");
assert(false);
}
if (lden == 2) {
/* den1 is causal; den2 is non-causal */
if (lden1+lden2 > 2) {
// shiftnc = lden2 - 1
if (lden2 - 1 != 0) {
printf("FATAL. Non-zero shiftnc! Need to add missing implementation.\n");
assert(false);
}
if (condition == 0) {
filter(tmp, out, n, d_unity, 1, den1, lden1, false);
filter(out, tmp, n, d_unity, 1, den2, lden2, true);
for (unsigned i=0; i<n; i+=2) {
out[i] *= den2[lden2-1];
out[i+1] *= den2[lden2-1];
}
} else {
filter(tmp, out, n, d_unity, 1, den1, lden1, true);
filter(out, tmp, n, d_unity, 1, den2, lden2, true);
//EO: combine revert and scaling
for(unsigned i=0; i<(n/2); i++) {
j = n - i - 1;
dtmp = out[j];
out[j] = out[i] * den2[lden2-1];
out[i] = dtmp * den2[lden2-1];
}
}
} else {
printf("fatal: lden1+lden2 <= 2.");
assert(false);
}
} else {
printf("fatal: lden not equal to 2.");
assert(false);
}
}
__global__ void my_kernel(double* __restrict__ tcIn,
const int tcLength,
const int voxelNb,
double *num,
const int lnum,
const int lden,
double *den1, const int lden1,
double *den2, const int lden2,
const double lambdaTemp,
const double maxeig,
const int cost_save,
const int nit,
const double *noiseFinIn,
double *tcOut,
double *noiseFinOut) {
double acc, nv, noise_estimate, precision, t_l, lambda_t;
double lambda = 0.0;
double t = 1.0;
double maxeig_ = 1.0 / maxeig;
double tclsqrt_ = sqrt(1.0 / tcLength);
double *z = (double*) malloc(tcLength*sizeof(double));
assert(z != NULL);
double *z_l = (double*) malloc(tcLength*sizeof(double));
assert(z_l != NULL);
double *z3 = (double*) malloc(tcLength*sizeof(double));
assert(z3 != NULL);
double *s = (double*) malloc(tcLength*sizeof(double));
assert(s != NULL);
double *tmp = (double*) malloc(tcLength*sizeof(double));
assert(tmp != NULL);
unsigned tid = blockIdx.x * blockDim.x + threadIdx.x;
if (tid >= voxelNb)
return;
unsigned tcinid = tid * tcLength;
/*
Compute Daubechies 6 level 1 detailed coefficients
Replicates Matlab sequence (see MyTemporal.m):
[coef, len] = wavedec(s,1,'db3');
coef(1:len(1)) = [];
! Assumes default dwtmode('sym') is set.
! Requires even length input signal.
*/
// Requires even length input signal
assert(tcLength%2==0);
// Number of detailed coefficients to be computed
unsigned nc = floor((tcLength-1)*0.5) + DB3;
double *coeff = (double*) malloc(nc * sizeof(double));
assert(coeff != NULL);
// Number of coefficients that do not depend on padded values
unsigned n1 = nc - 4;
coeff[0] = 0.0;
coeff[0] += d_g6[0] * tcIn[tcinid + 3];
coeff[0] += d_g6[1] * tcIn[tcinid + 2];
coeff[0] += d_g6[2] * tcIn[tcinid + 1];
coeff[0] += d_g6[3] * tcIn[tcinid];
coeff[0] += d_g6[4] * tcIn[tcinid];
coeff[0] += d_g6[5] * tcIn[tcinid + 1];
coeff[1] = 0.0;
coeff[1] += d_g6[0] * tcIn[tcinid + 1];
coeff[1] += d_g6[1] * tcIn[tcinid];
coeff[1] += d_g6[2] * tcIn[tcinid];
coeff[1] += d_g6[3] * tcIn[tcinid + 1];
coeff[1] += d_g6[4] * tcIn[tcinid + 2];
coeff[1] += d_g6[5] * tcIn[tcinid + 3];
for (unsigned i=0; i < n1; i++) {
coeff[i+2] = 0.0;
for (unsigned j=0; j<DB6; j++) {
coeff[i+2] += d_g6[j] * tcIn[tcinid + 2*i + j];
}
}
coeff[n1+2] = 0.0;
coeff[n1+2] += d_g6[0] * tcIn[tcinid + tcLength - 4];
coeff[n1+2] += d_g6[1] * tcIn[tcinid + tcLength - 3];
coeff[n1+2] += d_g6[2] * tcIn[tcinid + tcLength - 2];
coeff[n1+2] += d_g6[3] * tcIn[tcinid + tcLength - 1];
coeff[n1+2] += d_g6[4] * tcIn[tcinid + tcLength - 1];
coeff[n1+2] += d_g6[5] * tcIn[tcinid + tcLength - 2];
coeff[n1+3] = 0.0;
coeff[n1+3] += d_g6[0] * tcIn[tcinid + tcLength - 2];
coeff[n1+3] += d_g6[1] * tcIn[tcinid + tcLength - 1];
coeff[n1+3] += d_g6[2] * tcIn[tcinid + tcLength - 1];
coeff[n1+3] += d_g6[3] * tcIn[tcinid + tcLength - 2];
coeff[n1+3] += d_g6[4] * tcIn[tcinid + tcLength - 3];
coeff[n1+3] += d_g6[5] * tcIn[tcinid + tcLength - 4];
thrust::sort(thrust::seq, coeff, coeff+nc);
double median = 0.0;
const size_t lhs = (nc - 1) * 0.5;
const size_t rhs = nc * 0.5;
if (lhs == rhs) {
median = coeff[lhs];
} else {
median = (coeff[lhs] + coeff[rhs]) * 0.5;
}
for(unsigned i=0; i<nc; i++) {
coeff[i] = fabs(coeff[i] - median);
}
thrust::sort(thrust::seq, coeff, coeff+nc);
if (lhs == rhs) {
median = coeff[lhs];
} else {
median = (coeff[lhs] + coeff[rhs]) * 0.5;
}
free(coeff);
__syncthreads();
lambda_t = median * lambdaTemp;
noise_estimate = lambda_t;
precision = noise_estimate * 1E-5;
if (noiseFinIn[tid] == 0.0) {
//printf("case 0.0 in cuda: lambda = %15.10f\n", lambda_t);
lambda = lambda_t;
} else {
//printf("case != 0.0 in cuda: lambda = %15.10f (vox %i)\n", noiseFinIn[tid], tid);
lambda = noiseFinIn[tid];
}
for (unsigned i=0; i<tcLength; i+=2) {
s[i] = 0.0;
s[i+1] = 0.0;
z[i] = tcIn[tcinid + i];
z[i+1] = tcIn[tcinid + i + 1];
}
filter_boundary(tmp, &tcOut[tcinid], &tcIn[tcinid], tcLength, num, lnum, lden, den1, lden1, den2, lden2, 0);
nv = 0.0;
for (unsigned k=0; k<nit; k++) {
for(unsigned i=0; i<tcLength; i++)
z_l[i] = z[i];
filter_boundary(tmp, z, s, tcLength, num, lnum, lden, den1, lden1, den2, lden2, 1);
filter_boundary(tmp, z3, z, tcLength, num, lnum, lden, den1, lden1, den2, lden2, 0);
for (unsigned i=0; i<tcLength; i++) {
z[i] = maxeig_ * (tcOut[tcinid+i]/lambda - z3[i]) + s[i];
if (z[i] > 1.0) z[i] = 1.0;
if (z[i] < -1.0) z[i] = -1.0;
}
t_l = t;
t = 0.5 + sqrt(0.25 + t*t);
double ct = (t_l - 1.0)/t;
for (unsigned i=0; i<tcLength; i++)
s[i] = z[i] + ct * (z[i]-z_l[i]);
nv = 0.0;
if (cost_save) {
printf("Fatal. Missing implementation if cost_save == 1 in Temporal_TA_MEX.\n");
assert(false);
} else {
filter_boundary(tmp, z3, z, tcLength, num, lnum, lden, den1, lden1, den2, lden2, 1);
acc = 0.0;
for(unsigned i=0; i<tcLength; i++)
acc += z3[i]*z3[i];
nv = lambda * sqrt(acc) * tclsqrt_;
}
if (fabs(nv-noise_estimate) > precision) {
lambda *= noise_estimate/nv;
}
}
for(unsigned i=0; i<tcLength; i++)
tcOut[tcinid + i] = tcIn[tcinid + i] - lambda*z3[i];
noiseFinOut[tid] = nv;
free(s);
free(z3);
free(z_l);
free(z);
free(tmp);
}
// Function to set constant memory for the Daubechies db6 coefficients
void set_constant_memory(const double *g6, const double *h6, const double *one) {
cudaMemcpyToSymbol(d_g6, g6, DB6 * sizeof(double));
cudaMemcpyToSymbol(d_h6, h6, DB6 * sizeof(double));
cudaMemcpyToSymbol(d_unity, one, sizeof(double));
}
c4science · Help