__global__ void ForwardprojectHorizontal( float *image, float *lines, int linesN ){ // current slice __shared__ float slice[Ny * Nz]; int slsz = Ny * Nz; float max_dist = TOR_WIDTH * TOR_WIDTH / 2.0f; for (int current_slice = 0; current_slice < Nx; ++current_slice){ // Load slice into shared memory int offset = current_slice * slsz; for (int vi = threadIdx.x; vi < slice_size; vi += blockDim.x){ slice[vi] = image[vi + offset]; } __syncthreads(); for (int line = threadIdx.x + blockIdx.x * blockDim.x; line < linesN; line += blockDim.x * gridDim.x) { // line structure CUDALor *the_line = (CUDALor*)lines + line; // line direction float l1 = the_line->dy, l2 = the_line->dz; float t = (current_slice - the_line->x0) / the_line->dx; float y = the_line->y0 + t * l1, z = the_line->z0 + t * l2; int centerY = floor(y), centerZ = floor(z); float sum = 0; for (int yy = centerY - TOR_WIDTH; yy <= centerY + TOR_WIDTH; ++yy) for (int zz = centerZ - TOR_WIDTH; zz <= centerZ + TOR_WIDTH; ++zz) { if ( yy >= 0 && yy < Ny && zz >= 0 && zz < Nz) { float dy = yy – y, dz = zz - z; float inner = dy * l1 + dz * l2; // Distance to the line, squared float d2 = dy * dy + dz * dz - inner * inner; float kern = (d2 < max_dist) ? exp(-d2 * ISIGMA) : 0; sum += slice[yy + Ny * zz] * kern; } } // Write the value back to global memory the_line->value += sum; } __syncthreads(); } }