Image Sequence processing in CUDA

Question

My task is to implement an image reconstruction algorithm code using CUDA. I am provided with a code in C for the same. The input to the code is a DAT file which contains 360 images of size 640 x 480.The code goes something like this:

    FILE *in,*out;
    float *i_data,*o_data;
    i_data=(float *)malloc(mem_size);
    for(int projection=0;projection<360;projection++)
    {
      in=fopen("filename.dat","rb");
      fread(i_data,mem_size,1,in);
      ... some math ...
      for(int slice_no=-240;slice_no<240:slice_no++)
      {
          for (i=-320;i<320;i++)
          for (j=-320;j<320;j++)
          {
             // do some operations
             (*(o_data*slice_no)+(j+320)+(i+240))+=(*(i_data*value)+(j+240)+(i+320));
             // some more math
          }
      } 
    }

The output float pointer is written back to a dat file. If I have to parallelize these loops, how would I do that in CUDA? I tried implementing the inner two for loops in CUDA using 640 blocks each of 640 threads. How do I give the thread index to the pointer operation inside the loop. I tried giving

         int i=blockIdx.x;
         int j=threadIdx.x;

     and 

         kernel<<<640,640>>>

But this gives wrong values in the output pointer. Most are NAN. Except the line with pointers shown in the above snippet, I was able to implement the other math successfully.

Could anyone please help me doing this? What is that I am doing wrong in this code? Also is it possible to parallelize all the for loops here?

Answer 1

It's not efficient to constraint the inner loop to run within a single block. Also, you'll get much more benefit if you parallelize all three loops (or 4 loops if you can load more data from disk and send many images at once to the gpu).

Your problem looks like

for(x = minx; x < maxx; ++x)
{
  for(y = miny; y < maxy; ++y)
  {
    for(z = minz; z < maxz; ++z)
    {
      // do some math
    }
   }
}

The number of threads in your kernel should be:

num_thread = (maxx-minx) * (maxy-miny) * (maxz-minz);

You should be independent from the size of a block. Set it to a constant like blockDim.x = 256 (experiment to find a better constant).

The number of blocks will be

gridDim.x =(num_thread + blockDim.x) / blockDim.x

Note that I'm reasoning as if the three loops were unrolled into single big loop but you can make 2D and even 3D blocks of threads if your driver allows it.

In your kernel, from the built-in variables, you can compute the global index inside the big unrolled loop as

int index = blockIdx.x * blockDim.x + threadIdx.x

Make sure you're not out of range

if(index >= num_thread) return; // do nothing

Now using index you can recover x, y and z (slice_no, i, j) as

x = index / ( (maxy-miny) * (maxz-minz) ) + minx;
y = ( index % ( (maxy-miny) * (maxz-minz) ) ) / (maxz-minz) + miny;
z = ( index % ( (maxy-miny) * (maxz-minz) ) ) % (maxz-minz) + minz;