Cuda实现Radon变换

冷不防 2023-07-13 03:35 9阅读 0赞

有关radon变换的算法的详情，写在这里：[http://blog.csdn.net/celte/article/details/9826505][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;
    				else
    					I[p*100+k]=0;
    			}
    			else
    				I[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;
    
    }

[http_blog.csdn.net_celte_article_details_9826505]: http://blog.csdn.net/celte/article/details/9826505