Cuda实现Radon变换

有关radon变换的算法的详情，写在这里：http://blog.csdn.net/celte/article/details/9826505

用Cuda实现radon变换，可以类似C语言程序操作。。。

这里可以简单的使用一个block，每个block包含numAngles个thread（在下面的代码中是181个thread）

#include <stdio.h>#include "math.h"#include "cuda_runtime.h"#include "device_launch_parameters.h"   __global__ static void radon_cuda_core(float *gpuOutput,float *gpuInput,float *gpuAngles,int M,int N,int xOrgin,int yOrgin,int numAngles,int rFirst,int rSize){const int tid=threadIdx.x;float angle=gpuAngles[tid];float *pOutput=gpuOutput+tid*rSize;float sine=sin(angle);float cosine=cos(angle);int m,n;for(m=0;m<rSize;m++)pOutput[m]=0.0;float *pInput=gpuInput;for(n=0;n<N;n++){for(m=0;m<M;m++){float value=*pInput++;if(value!=0.0){value*=0.25;float x=n-xOrgin;float y=yOrgin-m;int r1;float delta;float r=(x-0.25)*cosine+(y-0.25)*sine-rFirst;r1=(int)r;delta=r-r1;pOutput[r1]+=value*(1.0-delta);pOutput[r1+1]+=value*delta;r=(x-0.25)*cosine+(y+0.25)*sine-rFirst;r1=(int)r;delta=r-r1;pOutput[r1]+=value*(1.0-delta);pOutput[r1+1]+=value*delta;r=(x+0.25)*cosine+(y+0.25)*sine-rFirst;r1=(int)r;delta=r-r1;pOutput[r1]+=value*(1.0-delta);pOutput[r1+1]+=value*delta;r=(x+0.25)*cosine+(y-0.25)*sine-rFirst;r1=(int)r;delta=r-r1;pOutput[r1]+=value*(1.0-delta);pOutput[r1+1]+=value*delta;}}}}static void radon_cuda(float *pPtr, float *iPtr, float *thetaPtr, int M, int N,    int xOrigin, int yOrigin, int numAngles, int rFirst, int rSize){float *gpuInput;float *gpuOutput;float *gpuAngles;cudaMalloc((void **)&gpuInput,sizeof(float)*M*N);cudaMalloc((void **)&gpuOutput,sizeof(float)*numAngles*rSize);cudaMalloc((void **)&gpuAngles,sizeof(float)*numAngles);cudaMemcpy(gpuInput,iPtr,sizeof(float)*M*N,cudaMemcpyHostToDevice);cudaMemset(gpuOutput,0,numAngles*rSize);cudaMemcpy(gpuAngles,thetaPtr,sizeof(float)*numAngles,cudaMemcpyHostToDevice);radon_cuda_core<<<1,numAngles,0>>>(gpuOutput,gpuInput,gpuAngles,M,N,xOrigin,yOrigin,numAngles,rFirst,rSize);cudaMemcpy(pPtr,gpuOutput,sizeof(float)*numAngles*rSize,cudaMemcpyDeviceToHost);cudaFree(gpuInput);cudaFree(gpuOutput);cudaFree(gpuAngles);}int main(){int M=100;int N=100;int xOrigin=((N-1)/2>0)?((N-1)/2):0;int yOrigin=((M-1)/2>0)?((M-1)/2):0;int temp1=M-1-yOrigin;int temp2=N-1-xOrigin;int rLast=(int) ceil(sqrt((float) (temp1*temp1+temp2*temp2))) + 1;int rFirst=-rLast;int rSize=rLast-rFirst+1;int numAngles=181;float *thetaPtr=(float*)calloc(numAngles,sizeof(float));float *ptr=thetaPtr;float deg2rad = 3.141592 / 180.0;int k=0;for(k=0;k<numAngles;k++)*(ptr++)=numAngles*deg2rad;float* I=(float *)calloc(100*100,sizeof(float));int p=0;for(p=0;p<100;p++){for(k=0;k<100;k++){if(k>23 && k<75){if(p>23 && p<75)I[p*100+k]=1;elseI[p*100+k]=0;}elseI[p*100+k]=0;}}float *gpu_result;gpu_result=(float *)calloc(numAngles*rSize,sizeof(float));memset(gpu_result,0,numAngles*rSize);radon_cuda(gpu_result, I, thetaPtr, M, N, xOrigin, yOrigin, numAngles, rFirst, rSize);printf("gpu计算结束。。。\n");free(I);free(thetaPtr);free(cpu_result);free(gpu_result);getchar();return 0;}