RasterProcessTool/Toolbox/SimulationSARTool/SimulationSAR/BPBasic0_CUDA.cu

#include <cstdio>
#include <cufft.h>
#include <cmath>
#include <cuda_runtime.h>
#include <iostream>
#include <memory>
#include <vector>
#include <cmath>
#include <complex>
#include <device_launch_parameters.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#include <cufft.h>
#include <cufftw.h>
#include <cufftXt.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#include "BaseConstVariable.h"
#include "GPUTool.cuh"
#include "GPUBPTool.cuh"
#include "BPBasic0_CUDA.cuh"



__global__ void phaseCompensationKernel(cufftComplex* phdata, const double* Freq, double r, int K, int Na) {
    int freqIdx = blockIdx.x * blockDim.x + threadIdx.x;
    int pulseIdx = blockIdx.y * blockDim.y + threadIdx.y;

    if (freqIdx >= K || pulseIdx >= Na) return;

    int idx = pulseIdx * K + freqIdx;
    double phase = 4 * PI * Freq[freqIdx] * r / c;
    double cos_phase = cosf(phase);
    double sin_phase = sinf(phase);

    cufftComplex ph = phdata[idx];
    double new_real = ph.x * cos_phase - ph.y * sin_phase;
    double new_imag = ph.x * sin_phase + ph.y * cos_phase;
    phdata[idx] = make_cuComplex(new_real, new_imag);
}

__global__ void fftshiftKernel(cufftComplex* data, int Nfft, int Np) {
    int pulse = blockIdx.x * blockDim.x + threadIdx.x;
    if (pulse >= Np) return;

    int half = Nfft / 2;
    for (int i = 0; i < half; ++i) {
        cufftComplex temp = data[pulse * Nfft + i];
        data[pulse * Nfft + i] = data[pulse * Nfft + i + half];
        data[pulse * Nfft + i + half] = temp;
    }
}

__global__ void processPulseKernel(
    long prfid,
    int nx, int ny,
    const double* x_mat, const double* y_mat, const double* z_mat,
    double AntX, double AntY, double AntZ, 
    double R0, double minF,
    const cufftComplex* rc_pulse, 
    const double r_start, const double dr, const int nR,
    cufftComplex* im_final
) {
    //

    long long idx = blockIdx.x * blockDim.x + threadIdx.x;
    long long pixelcount = nx * ny;
    if (idx >= pixelcount) return;

    //printf("processPulseKernel start!!\n");

    //if (x >= nx || y >= ny) return;
    //int idx = x * ny + y;


    double dx = AntX - x_mat[idx];
    double dy = AntY - y_mat[idx];
    double dz = AntZ - z_mat[idx];
    
    //printf("processPulseKernel xmat !!\n");

    double dR = sqrtf(dx * dx + dy * dy + dz * dz) - R0;
    
    if (dR < r_start || dR >= (r_start + dr * (nR - 1))) return;
    // Linear interpolation
    double pos = (dR - r_start) / dr;
    int index = (int)floor(pos);
    double weight = pos - index;

    if (index < 0 || index >= nR - 1) return;

    cufftComplex rc_low = rc_pulse[prfid * nR +index];
    cufftComplex rc_high = rc_pulse[prfid * nR+index + 1];
    cufftComplex rc_interp;
    rc_interp.x = rc_low.x * (1 - weight) + rc_high.x * weight;
    rc_interp.y = rc_low.y * (1 - weight) + rc_high.y * weight;

    // Phase correction
    double phase = 4 * PI * minF / c * dR;
    double cos_phase = cos(phase);
    double sin_phase = sin(phase);

    cufftComplex phCorr;
    phCorr.x = rc_interp.x * cos_phase - rc_interp.y * sin_phase;
    phCorr.y = rc_interp.x * sin_phase + rc_interp.y * cos_phase;

    // Accumulate
    im_final[idx].x += phCorr.x;
    im_final[idx].y += phCorr.y;
    //printf("r_start=%e;dr=%e;nR=%d\n", r_start, dr, nR);
	if (abs(phCorr.x) > 1e-14 || abs(phCorr.y > 1e-14)) {
        printf(
            "[DEBUG] prfid=%-4ld | idx=%-8lld\n"
            "  Ant: X=%-18.10e Y=%-18.10e Z=%-18.10e\n"
            "  Pix: X=%-18.10e Y=%-18.10e Z=%-18.10e\n"
            "  dR=%-18.10e | pos=%-8.4f[%-6d+%-8.6f]\n"
            "  RC: low=(%-18.10e,%-18.10e) high=(%-18.10e,%-18.10e)\n"
            "  => interp=(%-18.10e,%-18.10e)\n"
            "  Phase: val=%-18.10e | corr=(%-18.10e,%-18.10e)\n"
            "  Final: im=(%-18.10e,%-18.10e)\n"
            "----------------------------------------\n",
            prfid, idx,
            AntX, AntY, AntZ,
            x_mat[idx], y_mat[idx], z_mat[idx],
            dR,
            pos, index, weight,
            rc_low.x, rc_low.y,
            rc_high.x, rc_high.y,
            rc_interp.x, rc_interp.y,
            phase,
            phCorr.x, phCorr.y,
            im_final[idx].x, im_final[idx].y
        );
	}
}

void bpBasic0CUDA(GPUDATA& data, int flag,double* h_R) {
    // Phase compensation
    if (flag == 1) {
        dim3 block(16, 16);
        dim3 grid((data.K + 15) / 16, (data.Np + 15) / 16);
        phaseCompensationKernel << <grid, block >> > (data.phdata, data.Freq, data.R0, data.K, data.Np);
        cudaCheckError(cudaDeviceSynchronize());
        //data.R0 = data.r;  // <20><><EFBFBD><EFBFBD>data.r<><72><EFBFBD><EFBFBD>ȷ<EFBFBD><C8B7><EFBFBD><EFBFBD>
    }

    // FFT<46><54><EFBFBD><EFBFBD>
    cufftHandle plan;
    cufftPlan1d(&plan, data.Nfft, CUFFT_C2C, data.Np);
    cufftExecC2C(plan, data.phdata, data.phdata, CUFFT_INVERSE);
    cufftDestroy(plan);

    // FFT<46><54>λ
    dim3 blockShift(256);
    dim3 gridShift((data.Np + 255) / 256);
    fftshiftKernel << <gridShift, blockShift >> > (data.phdata, data.Nfft, data.Np);
    cudaCheckError(cudaDeviceSynchronize());

    printf("fft finished!!\n");
    // ͼ<><CDBC><EFBFBD>ؽ<EFBFBD>
    double r_start = data.r_vec[0];
    double dr = (data.r_vec[data.Nfft - 1] - r_start) / (data.Nfft - 1);
    printf("dr = %f\n",dr);
    long pixelcount = data.nx* data.ny;
    long grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE;
    printf("grid finished!!\n");



    //double* d_R = (double*)mallocCUDADevice(sizeof(double) * data.nx * data.ny);
    printf("r_start=%e;dr=%e;nR=%d\n", r_start, dr, data.Nfft);
    printf("BPimage .....\n");
    for (long ii = 0; ii < data.Np; ++ii) {

        processPulseKernel << <grid_size, BLOCK_SIZE >> > (
            ii,
            data.nx, data.ny,
            data.x_mat, data.y_mat, data.z_mat,
            data.AntX[ii], data.AntY[ii], data.AntZ[ii],
            data.R0, data.minF[ii],
            data.phdata,
            r_start, dr, data.Nfft,
            data.im_final
            //,d_R
            );
        PrintLasterError("processPulseKernel");
        if (ii % 1000==0) {
            printf("\rPRF(%f %)  %d / %d\t\t\t\t",(ii*100.0/data.Np), ii,data.Np);
        }


       // DeviceToHost(h_R, d_R, sizeof(double) * data.nx * data.ny);

       // double minR = h_R[0], maxR = h_R[0];
       // for (long i = 0; i < data.nx * data.ny; i++) {
       //     if (minR > h_R[i]) { minR = h_R[i]; }
       //     if (maxR < h_R[i]) { maxR = h_R[i]; }
       // }

       //printf("prfid=%d; R=[ %e , %e ]\n", ii,minR, maxR);
       //break;

    }
    //FreeCUDADevice(d_R);
    
   
    cudaCheckError(cudaDeviceSynchronize());
}

 


void initGPUData(GPUDATA& h_data, GPUDATA& d_data) {
	d_data.AntX =h_data.AntX; //(double*)mallocCUDADevice(sizeof(double) * h_data.Np);
    d_data.AntY = h_data.AntY;//(double*)mallocCUDADevice(sizeof(double) * h_data.Np);
    d_data.AntZ = h_data.AntZ;// (double*)mallocCUDADevice(sizeof(double) * h_data.Np);
	d_data.minF = h_data.minF;// (double*)mallocCUDADevice(sizeof(double) * h_data.Np);
	d_data.x_mat = (double*)mallocCUDADevice(sizeof(double) * h_data.nx * h_data.ny);
	d_data.y_mat = (double*)mallocCUDADevice(sizeof(double) * h_data.nx * h_data.ny);
	d_data.z_mat = (double*)mallocCUDADevice(sizeof(double) * h_data.nx * h_data.ny);
    d_data.r_vec = h_data.r_vec;// (double*)mallocCUDADevice(sizeof(double) * h_data.Nfft);
	d_data.Freq = (double*)mallocCUDADevice(sizeof(double) * h_data.Nfft);
	d_data.phdata = (cufftComplex*)mallocCUDADevice(sizeof(cufftComplex) * h_data.K * h_data.Np);
	d_data.im_final = (cufftComplex*)mallocCUDADevice(sizeof(cufftComplex) * h_data.nx * h_data.ny);

    //HostToDevice(h_data.AntX, d_data.AntX,sizeof(double) * h_data.Np);
    //HostToDevice(h_data.AntY, d_data.AntY,sizeof(double) * h_data.Np);
    //HostToDevice(h_data.AntZ, d_data.AntZ,sizeof(double) * h_data.Np);
    //HostToDevice(h_data.minF, d_data.minF,sizeof(double) * h_data.Np);
    HostToDevice(h_data.x_mat, d_data.x_mat,sizeof(double) * h_data.nx * h_data.ny);
    HostToDevice(h_data.y_mat, d_data.y_mat,sizeof(double) * h_data.nx * h_data.ny);
    HostToDevice(h_data.z_mat, d_data.z_mat,sizeof(double) * h_data.nx * h_data.ny);
    HostToDevice(h_data.Freq, d_data.Freq, sizeof(double) * h_data.Nfft);
    //HostToDevice(h_data.r_vec, d_data.r_vec, sizeof(double) * h_data.Nfft);
    HostToDevice(h_data.phdata, d_data.phdata, sizeof(cufftComplex) * h_data.K * h_data.Np);
    HostToDevice(h_data.im_final, d_data.im_final, sizeof(cufftComplex) * h_data.nx * h_data.ny);
 
    // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    d_data.Nfft = h_data.Nfft;
    d_data.K = h_data.K;
    d_data.Np = h_data.Np;
    d_data.nx = h_data.nx;
    d_data.ny = h_data.ny;
    d_data.R0 = h_data.R0;
    d_data.deltaF = h_data.deltaF;
}

void freeGPUData(GPUDATA& d_data) {

    //FreeCUDADevice((d_data.AntX));
    //FreeCUDADevice((d_data.AntY));
    //FreeCUDADevice((d_data.AntZ));
    //FreeCUDADevice((d_data.minF));
    FreeCUDADevice((d_data.x_mat));
    FreeCUDADevice((d_data.y_mat));
    FreeCUDADevice((d_data.z_mat));
    //FreeCUDADevice((d_data.r_vec));
    FreeCUDADevice((d_data.Freq));
    FreeCUDADevice((d_data.phdata));
    FreeCUDADevice((d_data.im_final));
 
}

void freeHostData(GPUDATA& h_data) {
    //FreeCUDAHost((h_data.AntX)); 
    //FreeCUDAHost((h_data.AntY));
    //FreeCUDAHost((h_data.AntZ));
    FreeCUDAHost((h_data.minF));
    //FreeCUDAHost((h_data.x_mat));
    //FreeCUDAHost((h_data.y_mat));
    //FreeCUDAHost((h_data.z_mat));
    FreeCUDAHost((h_data.r_vec));
    FreeCUDAHost((h_data.Freq));
    FreeCUDAHost((h_data.phdata));
    FreeCUDAHost((h_data.im_final));
}

void BPBasic0(GPUDATA& h_data)
{
    GPUDATA d_data;

    initGPUData(h_data, d_data);

    bpBasic0CUDA(d_data, 0);

    DeviceToHost(h_data.im_final, d_data.im_final, sizeof(cuComplex) * h_data.nx * h_data.ny);
    freeGPUData(d_data);
}

//int main() {
//    GPUDATA h_data, d_data;
//
//    // <20><>ʼ<EFBFBD><CABC><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
//    h_data.Nfft = 1024;
//    h_data.K = 512;
//    // ... <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʼ<EFBFBD><CABC>
//
//    // <20><>ʼ<EFBFBD><CABC><EFBFBD>豸<EFBFBD>ڴ<EFBFBD>
//    initGPUData(h_data, d_data);
//
//    // ִ<><D6B4><EFBFBD>㷨
//    bpBasic0CUDA(d_data, 0);
//
//    // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
//    cudaCheckError(cudaMemcpy(h_data.im_final, d_data.im_final,
//        sizeof(cufftComplex) * h_data.nx * h_data.ny, cudaMemcpyDeviceToHost));
//
//    // <20>ͷ<EFBFBD><CDB7><EFBFBD>Դ
//    freeGPUData(d_data);
//
//    return 0;
//}