diff --git a/BaseCommonLibrary/BaseTool/BaseTool.h b/BaseCommonLibrary/BaseTool/BaseTool.h index 7990656..1122331 100644 --- a/BaseCommonLibrary/BaseTool/BaseTool.h +++ b/BaseCommonLibrary/BaseTool/BaseTool.h @@ -132,22 +132,44 @@ void initializeMatrixWithSSE2(Eigen::MatrixXf& mat, const float* data, long rowc /** 模板函数类 ***********************************************************************************************************/ template -inline void BASECONSTVARIABLEAPI memsetInitArray(std::shared_ptr ptr, long arrcount,T ti) { +inline void memsetInitArray(std::shared_ptr ptr, long arrcount,T ti) { for (long i = 0; i < arrcount; i++) { ptr.get()[i] = ti; } } template -inline void BASECONSTVARIABLEAPI memcpyArray(std::shared_ptr srct, std::shared_ptr dest, long arrcount) { +inline void memcpyArray(std::shared_ptr srct, std::shared_ptr dest, long arrcount) { for (long i = 0; i < arrcount; i++) { dest.get()[i] = srct.get()[i]; } } +template +inline void minValueInArr(T* ptr, long arrcount, T& minvalue) { + + if (arrcount == 0)return; + minvalue = ptr[0]; + for (long i = 0; i < arrcount; i++) { + if (minvalue > ptr[i]) { + minvalue = ptr[i]; + } + } +} +template +inline void maxValueInArr(T* ptr, long arrcount, T& maxvalue) { + if (arrcount == 0)return; + + maxvalue = ptr[0]; + for (long i = 0; i < arrcount; i++) { + if (maxvalue < ptr[i]) { + maxvalue = ptr[i]; + } + } +} #endif \ No newline at end of file diff --git a/BaseCommonLibrary/BaseTool/EchoDataFormat.cpp b/BaseCommonLibrary/BaseTool/EchoDataFormat.cpp index cb5e839..aa65d70 100644 --- a/BaseCommonLibrary/BaseTool/EchoDataFormat.cpp +++ b/BaseCommonLibrary/BaseTool/EchoDataFormat.cpp @@ -531,7 +531,7 @@ std::shared_ptr EchoL0Dataset::getAntPos() return temp; } -std::shared_ptr> EchoL0Dataset::getEchoArr(long startPRF, long PRFLen) +std::shared_ptr> EchoL0Dataset::getEchoArr(long startPRF, long& PRFLen) { if (!(startPRF < this->PluseCount)) { qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_PRFIDXOUTRANGE))<PluseCount; @@ -558,6 +558,9 @@ std::shared_ptr> EchoL0Dataset::getEchoArr(long startPRF, l long width = rasterDataset->GetRasterXSize(); long height = rasterDataset->GetRasterYSize(); long band_num = rasterDataset->GetRasterCount(); + + PRFLen = (PRFLen + startPRF) < height ? PRFLen : height - startPRF; + std::shared_ptr> temp = nullptr; if (height != this->PluseCount || width != this->PlusePoints) { qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR)); diff --git a/BaseCommonLibrary/BaseTool/EchoDataFormat.h b/BaseCommonLibrary/BaseTool/EchoDataFormat.h index 843fb9a..1c7c486 100644 --- a/BaseCommonLibrary/BaseTool/EchoDataFormat.h +++ b/BaseCommonLibrary/BaseTool/EchoDataFormat.h @@ -184,7 +184,7 @@ public: // public: // ȡļ std::shared_ptr getAntPos(); - std::shared_ptr> getEchoArr(long startPRF, long PRFLen); + std::shared_ptr> getEchoArr(long startPRF, long& PRFLen); std::shared_ptr> getEchoArr(); //ļ ErrorCode saveAntPos(std::shared_ptr ptr); // עΣգдǰǷȷ diff --git a/RasterMainWidgetGUI/RasterMainWidget/RasterMainWidget.ui b/RasterMainWidgetGUI/RasterMainWidget/RasterMainWidget.ui index 70ea495..af85019 100644 --- a/RasterMainWidgetGUI/RasterMainWidget/RasterMainWidget.ui +++ b/RasterMainWidgetGUI/RasterMainWidget/RasterMainWidget.ui @@ -85,7 +85,7 @@ 0 0 906 - 22 + 23 @@ -367,7 +367,7 @@ p, li { white-space: pre-wrap; } 信息窗口 - 1 + 2 diff --git a/Toolbox/SimulationSARTool/SimulationSAR/BPBasic0_CUDA.cu b/Toolbox/SimulationSARTool/SimulationSAR/BPBasic0_CUDA.cu index baf42e3..95b4778 100644 --- a/Toolbox/SimulationSARTool/SimulationSAR/BPBasic0_CUDA.cu +++ b/Toolbox/SimulationSARTool/SimulationSAR/BPBasic0_CUDA.cu @@ -54,25 +54,19 @@ __global__ void fftshiftKernel(cufftComplex* data, int Nfft, int Np) { __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, double r_start, double dr, int nR, - cufftComplex* im_final) { - //int x = blockIdx.x * blockDim.x + threadIdx.x; - //int y = blockIdx.y * blockDim.y + threadIdx.y; + 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 +) { + // - - - int x = blockIdx.x * blockDim.x + threadIdx.x; - int y = blockIdx.y * blockDim.y + threadIdx.y; - - if (x >= nx || y >= ny) return; - - long long idx = x * ny + y; - - //long long idx = blockIdx.x * blockDim.x + threadIdx.x; - //long long pixelcount = nx * ny; - //if (idx >= pixelcount) return; + long long idx = blockIdx.x * blockDim.x + threadIdx.x; + long long pixelcount = nx * ny; + if (idx >= pixelcount) return; //printf("processPulseKernel start!!\n"); @@ -83,22 +77,19 @@ __global__ void processPulseKernel( 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; - - // Range check + if (dR < r_start || dR >= (r_start + dr * (nR - 1))) return; - // Linear interpolation double pos = (dR - r_start) / dr; - int index = (int)floorf(pos); + 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; @@ -106,9 +97,9 @@ __global__ void processPulseKernel( rc_interp.y = rc_low.y * (1 - weight) + rc_high.y * weight; // Phase correction - double phase = 4 * PI * minF * dR / c; - double cos_phase = cosf(phase); - double sin_phase = sinf(phase); + 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; @@ -117,11 +108,34 @@ __global__ void processPulseKernel( // 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) { +void bpBasic0CUDA(GPUDATA& data, int flag,double* h_R) { // Phase compensation if (flag == 1) { dim3 block(16, 16); @@ -147,11 +161,16 @@ void bpBasic0CUDA(GPUDATA& data, int flag) { // ͼؽ double r_start = data.r_vec[0]; double dr = (data.r_vec[data.Nfft - 1] - r_start) / (data.Nfft - 1); - printf("dr finished!!\n"); - size_t pixelcount = data.nx* data.ny; - size_t grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE; + 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 << > > ( @@ -163,36 +182,33 @@ void bpBasic0CUDA(GPUDATA& data, int flag) { 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 computeRvec(GPUDATA& data) { - // maxWrҪȼdeltaF - double deltaF = data.deltaF; // ȡ - double maxWr = 299792458.0f / (2.0f * deltaF); - - // r_vecˣ - double* r_vec_host = (double*)mallocCUDAHost(sizeof(double) * data.Nfft);// new double[data.Nfft]; - const double step = maxWr / data.Nfft; - const double start = -data.Nfft / 2.0f * step; - - for (int i = 0; i < data.Nfft; ++i) { - r_vec_host[i] = start + i * step; - } - - data.r_vec = r_vec_host; - -} - void initGPUData(GPUDATA& h_data, GPUDATA& d_data) { diff --git a/Toolbox/SimulationSARTool/SimulationSAR/BPBasic0_CUDA.cuh b/Toolbox/SimulationSARTool/SimulationSAR/BPBasic0_CUDA.cuh index 944d4cd..80c1f14 100644 --- a/Toolbox/SimulationSARTool/SimulationSAR/BPBasic0_CUDA.cuh +++ b/Toolbox/SimulationSARTool/SimulationSAR/BPBasic0_CUDA.cuh @@ -16,7 +16,7 @@ inline void gpuAssert(cudaError_t code, const char* file, int line) { struct GPUDATA { - size_t Nfft, K, Np, nx, ny; // ҶƵͼСͼ + long Nfft, K, Np, nx, ny; // ҶƵͼСͼ double* AntX, * AntY, * AntZ, * minF; // ꡢʼƵ double* x_mat, * y_mat, * z_mat;// double* r_vec; // 귶Χ @@ -28,8 +28,7 @@ struct GPUDATA { }; extern "C" { - void bpBasic0CUDA(GPUDATA& data, int flag); - void computeRvec(GPUDATA& data); + void bpBasic0CUDA(GPUDATA& data, int flag,double* h_R=nullptr); void initGPUData(GPUDATA& h_data, GPUDATA& d_data); void freeGPUData(GPUDATA& d_data); void freeHostData(GPUDATA& d_data); diff --git a/Toolbox/SimulationSARTool/SimulationSAR/TBPImageAlgCls.cpp b/Toolbox/SimulationSARTool/SimulationSAR/TBPImageAlgCls.cpp index cb5308a..284174c 100644 --- a/Toolbox/SimulationSARTool/SimulationSAR/TBPImageAlgCls.cpp +++ b/Toolbox/SimulationSARTool/SimulationSAR/TBPImageAlgCls.cpp @@ -10,6 +10,7 @@ #include "GPUTBPImage.cuh" #include "ImageOperatorBase.h" #include "BPBasic0_CUDA.cuh" +#include "BaseTool.h" void CreatePixelXYZ(std::shared_ptr echoL0ds, QString outPixelXYZPath) @@ -337,6 +338,21 @@ ErrorCode TBPImageAlgCls::ProcessGPU() std::shared_ptr> imgArr = this->L1ds->getImageRaster(img_rid, imrowcount); // 回波值 + double img_x_min=0, img_x_max=0; + double img_y_min=0, img_y_max=0; + double img_z_min=0, img_z_max=0; + minValueInArr(img_x.get(), imrowcount * imcolcount, img_x_min); + maxValueInArr(img_x.get(), imrowcount * imcolcount, img_x_max); + minValueInArr(img_y.get(), imrowcount * imcolcount, img_y_min); + maxValueInArr(img_y.get(), imrowcount * imcolcount, img_y_max); + minValueInArr(img_z.get(), imrowcount * imcolcount, img_z_min); + maxValueInArr(img_z.get(), imrowcount * imcolcount, img_z_max); + qDebug() << "imgX:\t" << img_x_min << " ~ " << img_x_max; + qDebug() << "imgY:\t" << img_y_min << " ~ " << img_y_max; + qDebug() << "imgZ:\t" << img_z_min << " ~ " << img_z_max; + + + //shared_complexPtrToHostCuComplex(std::complex*src, cuComplex * dst, long len) // 处理 @@ -346,13 +362,41 @@ ErrorCode TBPImageAlgCls::ProcessGPU() h_data.Nfft = freqpoint; h_data.K = freqpoint; h_data.deltaF = this->L0ds->getBandwidth() / (freqpoint - 1); - computeRvec(h_data); + + // 计算maxWr（需要先计算deltaF） + double deltaF = h_data.deltaF; // 从输入参数获取 + double maxWr = 299792458.0f / (2.0f * deltaF); + + // 生成r_vec（主机端） + double* r_vec_host = (double*)mallocCUDAHost(sizeof(double) * h_data.Nfft);// new double[data.Nfft]; + const double step = maxWr / h_data.Nfft; + const double start = -1*h_data.Nfft / 2.0f * step; + printf("nfft=%d\n", h_data.Nfft); + for (int i = 0; i < h_data.Nfft; ++i) { + r_vec_host[i] = start + i * step; + } + + h_data.r_vec = r_vec_host; h_data.x_mat = img_x.get();//地面 h_data.y_mat = img_y.get(); h_data.z_mat = img_z.get(); h_data.nx = imcolcount; h_data.ny = imrowcount; + + minValueInArr(h_data.x_mat, h_data.nx * h_data.ny, img_x_min); + maxValueInArr(h_data.x_mat, h_data.nx * h_data.ny, img_x_max); + minValueInArr(h_data.y_mat, h_data.nx * h_data.ny, img_y_min); + maxValueInArr(h_data.y_mat, h_data.nx * h_data.ny, img_y_max); + minValueInArr(h_data.z_mat, h_data.nx * h_data.ny, img_z_min); + maxValueInArr(h_data.z_mat, h_data.nx * h_data.ny, img_z_max); + qDebug() << "imgX:\t" << img_x_min << " ~ " << img_x_max; + qDebug() << "imgY:\t" << img_y_min << " ~ " << img_y_max; + qDebug() << "imgZ:\t" << img_z_min << " ~ " << img_z_max; + qDebug() << "imgXYZ" << h_data.x_mat[5016600] << " , " << h_data.y_mat[5016600] << " , " << h_data.z_mat[5016600] << " , "; + + + h_data.R0 = refRange;// 参考斜距 h_data.im_final = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * imrowcount * imcolcount); @@ -394,6 +438,16 @@ ErrorCode TBPImageAlgCls::ProcessGPU() h_data.AntX = antpx.get(); // 天线 h_data.AntY = antpy.get(); h_data.AntZ = antpz.get(); + + // checkout + if (!this->checkZeros(h_data.AntX, ehcoprfcount) || + !this->checkZeros(h_data.AntY, ehcoprfcount) || + !this->checkZeros(h_data.AntZ, ehcoprfcount) + ) { + printf("the ant position is not zeros!!!"); + } + + h_data.Np = ehcoprfcount; h_data.phdata= (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * ehcoprfcount * echofreqpoint); @@ -401,24 +455,49 @@ ErrorCode TBPImageAlgCls::ProcessGPU() shared_complexPtrToHostCuComplex(echoArr.get(), h_data.phdata, ehcoprfcount * echofreqpoint); - // BP - GPUDATA d_data; - initGPUData(h_data, d_data); - bpBasic0CUDA(d_data, 0); - printf("BP finished!!!\n"); - DeviceToHost(h_data.im_final, d_data.im_final, sizeof(cuComplex) * h_data.nx * h_data.ny); - DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * ehcoprfcount* echofreqpoint); - freeGPUData(d_data); - - std::shared_ptr < std::complex> echofftdata(new std::complex[ehcoprfcount * echofreqpoint], delArrPtr); - qDebug() << QString(" image block PRF:[%1] \t").arg((img_rid + imrowcount) * 100.0 / imHeight) << QString(" echo [%1] prfidx:\t ").arg((prfid + ehcoprfcount)*100.0/ PRFCount)<< prfid << "\t-\t" << prfid + ehcoprfcount; - HostCuComplexToshared_complexPtr(h_data.phdata, echofftdata.get(), ehcoprfcount* echofreqpoint); - outTimeEchoImg.saveImage(echofftdata, prfid, 0, ehcoprfcount, echofreqpoint,1); + + { + // BP + GPUDATA d_data; + initGPUData(h_data, d_data); - HostCuComplexToshared_complexPtr(h_data.im_final, imgArr.get(), imrowcount* imcolcount); + qDebug() << "------------------------------------------------------"; + double min_r_vec = 0, max_r_cev = 0; + minValueInArr(h_data.r_vec, h_data.Nfft, min_r_vec); + maxValueInArr(h_data.r_vec, h_data.Nfft, max_r_cev); + qDebug() << "r_vec:\t" << min_r_vec << " ~ " << max_r_cev; + + minValueInArr(h_data.Freq, h_data.Nfft, min_r_vec); + maxValueInArr(h_data.Freq, h_data.Nfft, max_r_cev); + qDebug() << "Freq:\t" << min_r_vec << " ~ " << max_r_cev; + // 打印参数 + qDebug() << "R0\t" << h_data.R0; + qDebug() << "deltaF\t" << h_data.deltaF; + qDebug() << "Nfft\t" << h_data.Nfft; + qDebug() << "R0\t" << h_data.R0; + qDebug() << "K\t" << h_data.K; + qDebug() << "Np\t" << h_data.Np; + qDebug() << "nx\t" << h_data.nx; + qDebug() << "ny\t" << h_data.ny; + qDebug() << "------------------------------------------------------"; + + bpBasic0CUDA(d_data, 0); + + printf("BP finished!!!\n"); + DeviceToHost(h_data.im_final, d_data.im_final, sizeof(cuComplex) * h_data.nx * h_data.ny); + DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * h_data.Np * h_data.Nfft); + gdalImageComplex outTimeEchoImg = CreategdalImageComplexNoProj(this->TimeEchoDataPath, h_data.Np, h_data.Nfft, 1); + std::shared_ptr> data_time(new std::complex[h_data.Np * h_data.Nfft], delArrPtr); + for (long i = 0; i < h_data.Np * h_data.Nfft; i++) { + data_time.get()[i] = std::complex(h_data.phdata[i].x, h_data.phdata[i].y); + } + outTimeEchoImg.saveImage(data_time, 0, 0, h_data.Np, h_data.Nfft, 1); + + HostCuComplexToshared_complexPtr(h_data.im_final, imgArr.get(), imrowcount * imcolcount); + freeGPUData(d_data); + } this->L1ds->saveImageRaster(imgArr, img_rid, imrowcount); - } freeHostData(h_data); @@ -525,3 +604,16 @@ void TBPImageAlgCls::EchoFreqToTime() } +bool TBPImageAlgCls::checkZeros(double* data, long long len) +{ + bool flag = true; +#pragma omp parallel for + for (long long i = 0; i < len; i++) { + if (abs(data[i]) < 1e-7) { + flag= false; + break; + } + } + return flag; +} + diff --git a/Toolbox/SimulationSARTool/SimulationSAR/TBPImageAlgCls.h b/Toolbox/SimulationSARTool/SimulationSAR/TBPImageAlgCls.h index 0b7ab9a..3222ec3 100644 --- a/Toolbox/SimulationSARTool/SimulationSAR/TBPImageAlgCls.h +++ b/Toolbox/SimulationSARTool/SimulationSAR/TBPImageAlgCls.h @@ -66,10 +66,13 @@ private:// void EchoFreqToTime(); +private: + bool checkZeros(double* data, long long len); }; + void CreatePixelXYZ(std::shared_ptr echoL0ds,QString outPixelXYZPath);