diff --git a/BaseConstVariable.h b/BaseConstVariable.h index 1aaa76b..2339c16 100644 --- a/BaseConstVariable.h +++ b/BaseConstVariable.h @@ -3,7 +3,7 @@ #ifndef BASECONSTVARIABLE_H #define BASECONSTVARIABLE_H #define EIGEN_USE_MKL_ALL -#define EIGEN_NO_DEBUG +//#define EIGEN_NO_DEBUG #define EIGEN_USE_BLAS @@ -14,6 +14,7 @@ #define __CUDANVCC___ // 定义CUDA函数 +//#define __PRFDEBUG__ @@ -29,7 +30,7 @@ #define r2d 180/3.141592653589793238462643383279 #define d2r 3.141592653589793238462643383279/180 #define LIGHTSPEED 299792458 -#define PRECISIONTOLERANCE 1e-9 +#define PRECISIONTOLERANCE 1e-6 #define Radians2Degrees(Radians) Radians*PI_180 #define Degrees2Radians(Degrees) Degrees*T180_PI #define EARTHWE 0.000072292115 diff --git a/BaseLibraryCPP.rar b/BaseLibraryCPP.rar deleted file mode 100644 index 8bf6eaa..0000000 Binary files a/BaseLibraryCPP.rar and /dev/null differ diff --git a/BaseTool.cpp b/BaseTool.cpp index 9e9fca4..57c343b 100644 --- a/BaseTool.cpp +++ b/BaseTool.cpp @@ -33,6 +33,7 @@ #include /* 提供了向量结构*/ #include #include +#include @@ -500,7 +501,7 @@ long FindValueInStdVectorLast(std::vector& list, double& insertValue) long high = list.size() - 1; long mid = -1; // 处理边界 - if (list[high] <= insertValue) + if (list[high]+ PRECISIONTOLERANCE < insertValue) { return -1; } @@ -584,3 +585,45 @@ QVector SatelliteAntPos2SatellitePos(QVector pose } return antposes; } + + +QString getDebugDataPath(QString filename) +{ + QString folderName = "debugdata"; + QString appDir = QCoreApplication::applicationDirPath(); + QString folderpath = JoinPath(appDir, folderName); + if (!QDir(folderpath).exists()) { + QDir(appDir).mkdir(folderName); + } + QString datapath = JoinPath(folderpath, filename); + QFile datafile(datapath); + if (datafile.exists()) { + datafile.remove(); + } + return datapath; +} + + +std::vector split(const std::string& str, char delimiter) { + std::vector tokens; + std::string token; + std::istringstream tokenStream(str); + + while (std::getline(tokenStream, token, delimiter)) { + tokens.push_back(token); + } + + return tokens; +} + + +Eigen::VectorXd linspace(double start, double stop, int num) { + Eigen::VectorXd result(num); + + double step = (stop - start) / (num - 1); // 计算步长 + for (int i = 0; i < num; ++i) { + result[i] = start + i * step; // 生成等间距数值 + } + + return result; +} \ No newline at end of file diff --git a/BaseTool.h b/BaseTool.h index ccb0c1b..a08e3bf 100644 --- a/BaseTool.h +++ b/BaseTool.h @@ -111,7 +111,7 @@ QVector SatellitePos2SatelliteAntPos(QVector pose QVector SatelliteAntPos2SatellitePos(QVector poses); - - - +QString getDebugDataPath(QString filename); +std::vector split(const std::string& str, char delimiter); +Eigen::VectorXd linspace(double start, double stop, int num); #endif \ No newline at end of file diff --git a/EchoDataFormat.cpp b/EchoDataFormat.cpp index b6ceecf..9437f61 100644 --- a/EchoDataFormat.cpp +++ b/EchoDataFormat.cpp @@ -249,6 +249,24 @@ QString EchoL0Dataset::getEchoDataFilename() return GPSPointFilePath; } +void EchoL0Dataset::initEchoArr(std::complex init0) +{ + long blockline = Memory1MB * 2000 / 8 / 2 / this->PlusePoints; + + long start = 0; + for (start = 0; start < this->PluseCount; start = start + blockline) { + long templine = start + blockline < this->PluseCount ? blockline : this->PluseCount - start; + std::shared_ptr> echotemp = this->getEchoArr(start, templine); + for (long i = 0; i < templine; i++) { + for (long j = 0; j < this->PlusePoints; j++) { + echotemp.get()[i * this->PlusePoints + j] = init0; + } + } + this->saveEchoArr(echotemp, start, templine); + qDebug() << "echo init col : " << start << "\t-\t" << start + templine << "\t/\t" << this->PluseCount; + } +} + // Getter Setter ʵ long EchoL0Dataset::getPluseCount() { return this->PluseCount; } void EchoL0Dataset::setPluseCount(long pulseCount) { this->PluseCount = pulseCount; } diff --git a/EchoDataFormat.h b/EchoDataFormat.h index 887baa4..f276bdd 100644 --- a/EchoDataFormat.h +++ b/EchoDataFormat.h @@ -108,6 +108,9 @@ public: QString getGPSPointFilename(); QString getEchoDataFilename(); + void initEchoArr(std::complex init0); + + private: // Ʒ QString folder; QString filename; diff --git a/GPUTool.cu b/GPUTool.cu index 18fce48..8eaa011 100644 --- a/GPUTool.cu +++ b/GPUTool.cu @@ -29,7 +29,8 @@ __device__ cuComplex cuCexpf(cuComplex x) // __device__ float GPU_getSigma0dB(CUDASigmaParam param, float theta) { - return param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6); + float sigma= param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6); + return sigma; } __device__ CUDAVector GPU_VectorAB(CUDAVector A, CUDAVector B) { @@ -96,6 +97,12 @@ __device__ float GPU_BillerInterpAntPattern(float* antpattern, float searththeta, float searchphi) { float stheta = searththeta; float sphi = searchphi; + if (stheta > 90) { + return 0; + } + else {} + + float pthetaid = (stheta - starttheta) / dtheta;// float pphiid = (sphi - startphi) / dphi; @@ -103,6 +110,7 @@ __device__ float GPU_BillerInterpAntPattern(float* antpattern, long nextTheta = lasttheta + 1; long lastphi = floorf(pphiid); long nextPhi = lastphi + 1; + if (lasttheta < 0 || nextTheta < 0 || lastphi < 0 || nextPhi < 0 || lasttheta >= thetapoints || nextTheta >= thetapoints || lastphi >= phipoints || nextPhi >= phipoints) @@ -124,10 +132,10 @@ __device__ float GPU_BillerInterpAntPattern(float* antpattern, float z22 = antpattern[nextTheta * phipoints + nextPhi]; - z11 = powf(10, z11 / 10); - z12 = powf(10, z12 / 10); - z21 = powf(10, z21 / 10); - z22 = powf(10, z22 / 10); + //z11 = powf(10, z11 / 10); // dB-> + //z12 = powf(10, z12 / 10); + //z21 = powf(10, z21 / 10); + //z22 = powf(10, z22 / 10); float GainValue = (z11 * (x2 - x) * (y2 - y) + z21 * (x - x1) * (y2 - y) @@ -170,7 +178,7 @@ __global__ void CUDA_B_DistanceA(float* Ax, float* Ay, float* Az, float Bx, floa __global__ void CUDA_make_VectorA_B(float sX, float sY, float sZ, float* tX, float* tY, float* tZ, float* RstX, float* RstY, float* RstZ, long len) { long idx = blockIdx.x * blockDim.x + threadIdx.x; if (idx < len) { - RstX[idx] = sX - tX[idx]; + RstX[idx] = sX - tX[idx]; // -> RstY[idx] = sY - tY[idx]; RstZ[idx] = sZ - tZ[idx]; } @@ -209,9 +217,9 @@ __global__ void CUDA_SatelliteAntDirectNormal(float* RstX, float* RstY, float* R , long len) { long idx = blockIdx.x * blockDim.x + threadIdx.x; if (idx < len) { - float Xst = -1 * RstX[idx]; // --> - float Yst = -1 * RstY[idx]; - float Zst = -1 * RstZ[idx]; + float Xst = -1 * RstX[idx]; // --> + float Yst = -1 * RstY[idx]; + float Zst = -1 * RstZ[idx]; float AntXaxisX = antXaxisX; float AntXaxisY = antXaxisY; float AntXaxisZ = antXaxisZ; @@ -221,17 +229,80 @@ __global__ void CUDA_SatelliteAntDirectNormal(float* RstX, float* RstY, float* R float AntZaxisX = antZaxisX; float AntZaxisY = antZaxisY; float AntZaxisZ = antZaxisZ; - // ָϵµֵ - float Xant = (Xst * (AntYaxisY * AntZaxisZ - AntYaxisZ * AntZaxisY) + Xst * (AntXaxisZ * AntZaxisY - AntXaxisY * AntZaxisZ) + Xst * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY)) / (AntXaxisX * (AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * (AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY)); - float Yant = (Yst * (AntYaxisZ * AntZaxisX - AntYaxisX * AntZaxisZ) + Yst * (AntXaxisX * AntZaxisZ - AntXaxisZ * AntZaxisX) + Yst * (AntYaxisX * AntXaxisZ - AntXaxisX * AntYaxisZ)) / (AntXaxisX * (AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * (AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY)); - float Zant = (Zst * (AntYaxisX * AntZaxisY - AntYaxisY * AntZaxisX) + Zst * (AntXaxisY * AntZaxisX - AntXaxisX * AntZaxisY) + Zst * (AntXaxisX * AntYaxisY - AntYaxisX * AntXaxisY)) / (AntXaxisX * (AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * (AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY)); + + // һ + float RstNorm = sqrtf(Xst * Xst + Yst * Yst + Zst * Zst); + float AntXaxisNorm = sqrtf(AntXaxisX * AntXaxisX + AntXaxisY * AntXaxisY + AntXaxisZ * AntXaxisZ); + float AntYaxisNorm = sqrtf(AntYaxisX * AntYaxisX + AntYaxisY * AntYaxisY + AntYaxisZ * AntYaxisZ); + float AntZaxisNorm = sqrtf(AntZaxisX * AntZaxisX + AntZaxisY * AntZaxisY + AntZaxisZ * AntZaxisZ); + + + float Rx = Xst / RstNorm; + float Ry = Yst / RstNorm; + float Rz = Zst / RstNorm; + float Xx = AntXaxisX / AntXaxisNorm; + float Xy = AntXaxisY / AntXaxisNorm; + float Xz = AntXaxisZ / AntXaxisNorm; + float Yx = AntYaxisX / AntYaxisNorm; + float Yy = AntYaxisY / AntYaxisNorm; + float Yz = AntYaxisZ / AntYaxisNorm; + float Zx = AntZaxisX / AntZaxisNorm; + float Zy = AntZaxisY / AntZaxisNorm; + float Zz = AntZaxisZ / AntZaxisNorm; + + float Xant = (Rx * Yy * Zz - Rx * Yz * Zy - Ry * Yx * Zz + Ry * Yz * Zx + Rz * Yx * Zy - Rz * Yy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx); + float Yant = -(Rx * Xy * Zz - Rx * Xz * Zy - Ry * Xx * Zz + Ry * Xz * Zx + Rz * Xx * Zy - Rz * Xy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx); + float Zant = (Rx * Xy * Yz - Rx * Xz * Yy - Ry * Xx * Yz + Ry * Xz * Yx + Rz * Xx * Yy - Rz * Xy * Yx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx); + + // theta phi float Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // pho float ThetaAnt = acosf(Zant / Norm); // theta Zļн - float YsinTheta = Yant / sinf(ThetaAnt); - float PhiAnt = (YsinTheta / abs(YsinTheta)) * acosf(Xant / (Norm * sinf(ThetaAnt))); - thetaAnt[idx] = ThetaAnt; - phiAnt[idx] = PhiAnt; + float PhiAnt = atanf(Yant / Xant); // -pi/2 ~pi/2 + + + if (abs(Yant) < PRECISIONTOLERANCE) { // X + PhiAnt = 0; + } + else if (abs(Xant) < PRECISIONTOLERANCE) { // Yϣԭ + if (Yant > 0) { + PhiAnt = PI / 2; + } + else { + PhiAnt = -PI / 2; + } + } + else if (Xant < 0) { + if (Yant > 0) { + PhiAnt = PI + PhiAnt; + } + else { + PhiAnt = -PI+PhiAnt ; + } + } + else { // Xant>0 X + + } + + if (isnan(PhiAnt)) { + printf("V=[%f,%f,%f];norm=%f;thetaAnt=%f;phiAnt=%f;\n", Xant, Yant, Zant,Norm, ThetaAnt, PhiAnt); + } + + //if (abs(ThetaAnt - 0) < PRECISIONTOLERANCE) { + // PhiAnt = 0; + //} + //else {} + + + thetaAnt[idx] = ThetaAnt*r2d; + phiAnt[idx] = PhiAnt*r2d; + //printf("Rst=[%f,%f,%f];AntXaxis = [%f, %f, %f];AntYaxis=[%f,%f,%f];AntZaxis=[%f,%f,%f];phiAnt=%f;thetaAnt=%f;\n", Xst, Yst, Zst + // , AntXaxisX, AntXaxisY, AntXaxisZ + // , AntYaxisX, AntYaxisY, AntYaxisZ + // , AntZaxisX, AntZaxisY, AntZaxisZ + // , phiAnt[idx] + // , thetaAnt[idx] + //); } } @@ -331,12 +402,12 @@ __global__ void CUDA_RTPC( //printf("\ntheta: %f\t,%f ,%f ,%f ,%f ,%f ,%f \n", localangle * r2d, sigma0Paramslist[clsid].p1, sigma0Paramslist[clsid].p2, sigma0Paramslist[clsid].p3, // sigma0Paramslist[clsid].p4, sigma0Paramslist[clsid].p5, sigma0Paramslist[clsid].p6); // ䷽ͼ - float transPattern = GPU_BillerInterpAntPattern(Tantpattern, + float transPattern = GPU_BillerInterpAntPattern(Tantpattern, Tstarttheta, Tstartphi, Tdtheta, Tdphi, Tthetapoints, Tphipoints, Rtanttheta.theta, Rtanttheta.phi) * r2d; // շͼ - float receivePattern = GPU_BillerInterpAntPattern(Rantpattern, + float receivePattern = GPU_BillerInterpAntPattern(Rantpattern, Rstarttheta, Rstartphi, Rdtheta, Rdphi, Rthetapoints, Rphipoints, Rtanttheta.theta, Rtanttheta.phi) * r2d; // λ @@ -352,7 +423,6 @@ __global__ void CUDA_RTPC( if (timeID < 0 || timeID >= Freqnumbers) { timeID = 0; amp = 0; - } else {} @@ -367,6 +437,127 @@ __global__ void CUDA_RTPC( } +__global__ void CUDA_TBPImage( + float* antPx, float* antPy, float* antPz, + float* imgx, float* imgy, float* imgz, + cuComplex* echoArr, cuComplex* imgArr, + float freq, float fs, float Rnear, float Rfar, + long rowcount, long colcount, + long prfid, long freqcount +) { + int idx = blockIdx.x * blockDim.x + threadIdx.x; + //printf("\nidx:\t %d %d %d\n", idx, linecount, plusepoint); + if (idx < rowcount * colcount) { + float R = sqrtf(powf(antPx[prfid] - imgx[idx], 2) + powf(antPy[prfid] - imgy[idx], 2) + powf(antPz[prfid] - imgz[idx], 2)); + float Ridf = ((R - Rnear) * 2 / LIGHTSPEED) * fs; + long Rid = floorf(Ridf); + if(Rid <0|| Rid >= freqcount){} + else { + float factorj = freq * 4 * PI / LIGHTSPEED; + cuComplex Rphi =cuCexpf(make_cuComplex(0, factorj * R));// У + imgArr[idx] = cuCaddf(imgArr[idx], cuCmulf(echoArr[Rid] , Rphi));// + } + } +} + + +__global__ void CUDA_calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt, + float* localangle, float* R, float* slopeangle, + float nearRange, float Fs, float Pt, float lamda, long FreqIDmax, + cuComplex* echoArr, long* FreqID, + long len) { + long idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < len) { + float r = R[idx]; + float amp = Pt * TransAnt[idx] * ReciveAnt[idx]; + amp = amp * sigma0[idx]; + amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(r, 4)); // ǿ + + // λ + float phi = (-4 * LAMP_CUDA_PI / lamda) * r; + cuComplex echophi = make_cuComplex(0, phi); + cuComplex echophiexp = cuCexpf(echophi); + + float timeR = 2 * (r - nearRange) / LIGHTSPEED * Fs; + long timeID = floorf(timeR); + if (timeID < 0 || timeID >= FreqIDmax) { + timeID = 0; + amp = 0; + } + + cuComplex echo; + echo.x = echophiexp.x * amp; + echo.y = echophiexp.y * amp; + echoArr[idx] = echo; + FreqID[idx] = timeID; + } +} + + +__global__ void CUDA_AntPatternInterpGain(float* anttheta, float* antphi, float* gain, + float* antpattern, float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints, long len) { + int idx = blockIdx.x * blockDim.x + threadIdx.x; + + if (idx < len) { + + float temptheta = anttheta[idx]; + float tempphi = antphi[idx]; + + + float antPatternGain = GPU_BillerInterpAntPattern(antpattern, + starttheta, startphi, dtheta, dphi, thetapoints, phipoints, + temptheta, tempphi) ; + gain[idx] = antPatternGain; + } +} + +//__global__ void Sigma0InterpPixel(long* demcls, float* demslopeangle, CUDASigmaParam* sigma0Paramslist, float* localangle, float* sigma0list, long sigmaparamslistlen, long len) +//{ +// long idx = blockIdx.x * blockDim.x + threadIdx.x; +// if (idx < len) { +// long clsid = demcls[idx]; +// if(clsid<=) +// sigma0list[idx] = 0; +// } +//} + + +__global__ void CUDA_InterpSigma( + long* demcls, float* sigmaAmp, float* localanglearr, long len, + CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen) { + long idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < len) { + long clsid = demcls[idx]; + float localangle = localanglearr[idx] * r2d; + CUDASigmaParam tempsigma = sigma0Paramslist[clsid]; + //printf("cls:%d;localangle=%f;\n",clsid, localangle); + + if (localangle < 0 || localangle >= 90) { + sigmaAmp[idx] = 0; + } + else {} + + if (abs(tempsigma.p1)< PRECISIONTOLERANCE&& + abs(tempsigma.p2) < PRECISIONTOLERANCE && + abs(tempsigma.p3) < PRECISIONTOLERANCE && + abs(tempsigma.p4) < PRECISIONTOLERANCE&& + abs(tempsigma.p5) < PRECISIONTOLERANCE&& + abs(tempsigma.p6) < PRECISIONTOLERANCE + ) { + sigmaAmp[idx] = 0; + } + else { + float sigma = GPU_getSigma0dB(tempsigma, localangle); + sigma = powf(10.0, sigma / 10.0);// ɢϵ + //printf("cls:%d;localangle=%f;sigma0=%f;\n", clsid, localangle, sigma); + sigmaAmp[idx] = sigma; + } + } +} + + + + //ʾ void checkCudaError(cudaError_t err, const char* msg) { if (err != cudaSuccess) { @@ -380,6 +571,7 @@ void checkCudaError(cudaError_t err, const char* msg) { extern "C" void* mallocCUDAHost(long memsize) { void* ptr; cudaMallocHost(&ptr, memsize); + #ifdef __CUDADEBUG__ cudaError_t err = cudaGetLastError(); if (err != cudaSuccess) { @@ -387,6 +579,7 @@ extern "C" void* mallocCUDAHost(long memsize) { // Possibly: exit(-1) if program cannot continue.... } #endif // __CUDADEBUG__ + cudaDeviceSynchronize(); return ptr; } @@ -400,7 +593,7 @@ extern "C" void FreeCUDAHost(void* ptr) { // Possibly: exit(-1) if program cannot continue.... } #endif // __CUDADEBUG__ - + cudaDeviceSynchronize(); } // GPUڴ @@ -414,7 +607,7 @@ extern "C" void* mallocCUDADevice(long memsize) { // Possibly: exit(-1) if program cannot continue.... } #endif // __CUDADEBUG__ - + cudaDeviceSynchronize(); return ptr; } @@ -428,7 +621,7 @@ extern "C" void FreeCUDADevice(void* ptr) { // Possibly: exit(-1) if program cannot continue.... } #endif // __CUDADEBUG__ - + cudaDeviceSynchronize(); } // GPU ڴת @@ -443,7 +636,7 @@ extern "C" void HostToDevice(void* hostptr, void* deviceptr, long memsize) { } #endif // __CUDADEBUG__ - + cudaDeviceSynchronize(); } extern "C" void DeviceToHost(void* hostptr, void* deviceptr, long memsize) { @@ -455,7 +648,7 @@ extern "C" void DeviceToHost(void* hostptr, void* deviceptr, long memsize) { // Possibly: exit(-1) if program cannot continue.... } #endif // __CUDADEBUG__ - + cudaDeviceSynchronize(); } extern "C" void CUDATestHelloWorld(float a,long len) { @@ -463,7 +656,7 @@ extern "C" void CUDATestHelloWorld(float a,long len) { int blockSize = 256; // ÿ߳ int numBlocks = (len + blockSize - 1) / blockSize; // pixelcount С // CUDA ˺ - CUDA_Test_HelloWorld << > > (a, len); + CUDA_Test_HelloWorld << > > (a, len); #ifdef __CUDADEBUG__ cudaError_t err = cudaGetLastError(); if (err != cudaSuccess) { @@ -474,12 +667,53 @@ extern "C" void CUDATestHelloWorld(float a,long len) { cudaDeviceSynchronize(); } +void CUDATBPImage(float* antPx, float* antPy, float* antPz, + float* imgx, float* imgy, float* imgz, + cuComplex* echoArr, cuComplex* imgArr, + float freq, float fs, float Rnear, float Rfar, + long rowcount, long colcount, + long prfid, long freqcount) +{ + int blockSize = 256; // ÿ߳ + int numBlocks = (rowcount * colcount + blockSize - 1) / blockSize; // pixelcount С + //printf("\nCUDA_RTPC_SiglePRF blockSize:%d ,numBlock:%d\n",blockSize,numBlocks); + // CUDA ˺ CUDA_RTPC_Kernel + + CUDA_TBPImage << > > ( + antPx, antPy, antPz, + imgx, imgy, imgz, + echoArr, imgArr, + freq, fs, Rnear, Rfar, + rowcount, colcount, + prfid, freqcount + ); + + +#ifdef __CUDADEBUG__ + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { + printf("CUDATBPImage CUDA Error: %s\n", cudaGetErrorString(err)); + // Possibly: exit(-1) if program cannot continue.... + } +#endif // __CUDADEBUG__ + cudaDeviceSynchronize(); +} + extern "C" void distanceAB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* R, long len) { // CUDA ˺Ϳijߴ int blockSize = 256; // ÿ߳ int numBlocks = (len + blockSize - 1) / blockSize; // pixelcount С // CUDA ˺ - CUDA_DistanceAB << > > (Ax, Ay, Az, Bx, By, Bz, R, len); + CUDA_DistanceAB << > > (Ax, Ay, Az, Bx, By, Bz, R, len); + +#ifdef __CUDADEBUG__ + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { + printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err)); + // Possibly: exit(-1) if program cannot continue.... + } +#endif // __CUDADEBUG__ + cudaDeviceSynchronize(); } extern "C" void BdistanceAs(float* Ax, float* Ay, float* Az, float Bx, float By, float Bz, float* R, long len) { @@ -487,7 +721,15 @@ extern "C" void BdistanceAs(float* Ax, float* Ay, float* Az, float Bx, float By, int blockSize = 256; // ÿ߳ int numBlocks = (len + blockSize - 1) / blockSize; // pixelcount С // CUDA ˺ - CUDA_B_DistanceA << > > (Ax, Ay, Az, Bx, By, Bz, R, len); + CUDA_B_DistanceA << > > (Ax, Ay, Az, Bx, By, Bz, R, len); + +#ifdef __CUDADEBUG__ + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { + printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err)); + // Possibly: exit(-1) if program cannot continue.... + } +#endif // __CUDADEBUG__ cudaDeviceSynchronize(); } @@ -496,7 +738,15 @@ extern "C" void make_VectorA_B(float sX, float sY, float sZ, float* tX, float* t int blockSize = 256; // ÿ߳ int numBlocks = (len + blockSize - 1) / blockSize; // pixelcount С // CUDA ˺ - CUDA_make_VectorA_B << > > (sX, sY, sZ, tX, tY, tZ, RstX, RstY, RstZ, len); + CUDA_make_VectorA_B << > > (sX, sY, sZ, tX, tY, tZ, RstX, RstY, RstZ, len); + +#ifdef __CUDADEBUG__ + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { + printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err)); + // Possibly: exit(-1) if program cannot continue.... + } +#endif // __CUDADEBUG__ cudaDeviceSynchronize(); } @@ -505,7 +755,15 @@ extern "C" void Norm_Vector(float* Vx, float* Vy, float* Vz, float* R, long len) int blockSize = 256; // ÿ߳ int numBlocks = (len + blockSize - 1) / blockSize; // pixelcount С // CUDA ˺ - CUDA_Norm_Vector << > > (Vx, Vy, Vz, R, len); + CUDA_Norm_Vector << > > (Vx, Vy, Vz, R, len); + +#ifdef __CUDADEBUG__ + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { + printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err)); + // Possibly: exit(-1) if program cannot continue.... + } +#endif // __CUDADEBUG__ cudaDeviceSynchronize(); } @@ -513,7 +771,15 @@ extern "C" void cosAngle_VA_AB(float* Ax, float* Ay, float* Az, float* Bx, float int blockSize = 256; // ÿ߳ int numBlocks = (len + blockSize - 1) / blockSize; // pixelcount С // CUDA ˺ - CUDA_cosAngle_VA_AB << > > (Ax, Ay, Az, Bx, By, Bz, anglecos, len); + CUDA_cosAngle_VA_AB << > > (Ax, Ay, Az, Bx, By, Bz, anglecos, len); + +#ifdef __CUDADEBUG__ + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { + printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err)); + // Possibly: exit(-1) if program cannot continue.... + } +#endif // __CUDADEBUG__ cudaDeviceSynchronize(); } @@ -528,25 +794,76 @@ extern "C" void SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ, int blockSize = 256; // ÿ߳ int numBlocks = (len + blockSize - 1) / blockSize; // pixelcount С // CUDA ˺ - CUDA_SatelliteAntDirectNormal << > > (RstX, RstY, RstZ, + CUDA_SatelliteAntDirectNormal << > > (RstX, RstY, RstZ, antXaxisX, antXaxisY, antXaxisZ, antYaxisX, antYaxisY, antYaxisZ, antZaxisX, antZaxisY, antZaxisZ, antDirectX, antDirectY, antDirectZ, thetaAnt, phiAnt , len); + +#ifdef __CUDADEBUG__ + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { + printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err)); + // Possibly: exit(-1) if program cannot continue.... + } +#endif // __CUDADEBUG__ cudaDeviceSynchronize(); } +extern "C" void AntPatternInterpGain(float* anttheta, float* antphi, float* gain, + float* antpattern, float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints, long len) { + int blockSize = 256; // ÿ߳ + int numBlocks = (len + blockSize - 1) / blockSize; // pixelcount С + //printf("\nCUDA_RTPC_SiglePRF blockSize:%d ,numBlock:%d\n", blockSize, numBlocks); + + CUDA_AntPatternInterpGain << > > ( anttheta,antphi, gain, + antpattern, + starttheta, startphi, dtheta, dphi, thetapoints, phipoints, + len); +#ifdef __CUDADEBUG__ + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { + printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err)); + // Possibly: exit(-1) if program cannot continue.... + } +#endif // __CUDADEBUG__ + cudaDeviceSynchronize(); +} + extern "C" void CUDARTPCPRF(float antPx, long len) { int blockSize = 256; // ÿ߳ int numBlocks = (len + blockSize - 1) / blockSize; // pixelcount С printf("\nCUDA_RTPC_SiglePRF blockSize:%d ,numBlock:%d\n", blockSize, numBlocks); - CUDA_Test_HelloWorld << > > (antPx, len); + CUDA_Test_HelloWorld << > > (antPx, len); +#ifdef __CUDADEBUG__ + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { + printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err)); + // Possibly: exit(-1) if program cannot continue.... + } +#endif // __CUDADEBUG__ + cudaDeviceSynchronize(); +} +extern "C" void calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt, + float* localangle, float* R, float* slopeangle, + float nearRange, float Fs, float pt, float lamda, long FreqIDmax, + cuComplex* echoAmp, long* FreqID, + long len) +{ + int blockSize = 256; // ÿ߳ + int numBlocks = (len + blockSize - 1) / blockSize; // pixelcount С + // CUDA ˺ + CUDA_calculationEcho << > > (sigma0, TransAnt, ReciveAnt, + localangle, R, slopeangle, + nearRange, Fs, pt, lamda, FreqIDmax, + echoAmp, FreqID, + len); #ifdef __CUDADEBUG__ cudaError_t err = cudaGetLastError(); if (err != cudaSuccess) { @@ -610,7 +927,25 @@ CUDA_RTPC << > > ( } - +extern "C" void CUDAInterpSigma( + long* demcls,float* sigmaAmp, float* localanglearr,long len, + CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen) {// ر-sigmaֵӦ-ulaby + int blockSize = 256; // ÿ߳ + int numBlocks = (len + blockSize - 1) / blockSize; // pixelcount С + // CUDA ˺ + CUDA_InterpSigma << > > ( + demcls, sigmaAmp, localanglearr, len, + sigma0Paramslist, sigmaparamslistlen + ); +#ifdef __CUDADEBUG__ + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { + printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err)); + // Possibly: exit(-1) if program cannot continue.... + } +#endif // __CUDADEBUG__ + cudaDeviceSynchronize(); +} #endif diff --git a/GPUTool.cuh b/GPUTool.cuh index dd19e48..5d0ad17 100644 --- a/GPUTool.cuh +++ b/GPUTool.cuh @@ -58,6 +58,9 @@ extern "C" void make_VectorA_B(float sX, float sY, float sZ, float* tX, float* t extern "C" void Norm_Vector(float* Vx, float* Vy, float* Vz, float* R, long member); extern "C" void cosAngle_VA_AB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* anglecos, long len); extern "C" void SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ, float antXaxisX, float antXaxisY, float antXaxisZ, float antYaxisX, float antYaxisY, float antYaxisZ, float antZaxisX, float antZaxisY, float antZaxisZ, float antDirectX, float antDirectY, float antDirectZ, float* thetaAnt, float* phiAnt, long len); + +extern "C" void AntPatternInterpGain(float* anttheta, float* antphi, float* gain, float* antpattern, float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints,long len); + extern "C" void CUDA_RTPC_SiglePRF( float antPx, float antPy, float antPZ, float antXaxisX, float antXaxisY, float antXaxisZ, @@ -75,22 +78,41 @@ extern "C" void CUDA_RTPC_SiglePRF( ); extern "C" void CUDARTPCPRF(float antPx, long len); - - - extern "C" void CUDATestHelloWorld(float a, long len); + +extern "C" void CUDATBPImage( + float* antPx, + float* antPy, + float* antPz, + float* imgx, + float* imgy, + float* imgz, + cuComplex* echoArr, + cuComplex* imgArr, + float freq, float fs, float Rnear, float Rfar, + long rowcount, long colcount, + long prfid, long freqcount +); + + +extern "C" void calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt, + float* localangle, float* R, float* slopeangle, + float nearRange, float Fs, float pt, float lamda, long FreqIDmax, + cuComplex* echoAmp, long* FreqID, + long len); + + +extern "C" void CUDAInterpSigma( + long* demcls, float* sigmaAmp, float* localanglearr, long len, + CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen); + #endif #endif - - - - - /** * diff --git a/ImageOperatorBase.cpp b/ImageOperatorBase.cpp index 45156ee..fac1fdd 100644 --- a/ImageOperatorBase.cpp +++ b/ImageOperatorBase.cpp @@ -819,8 +819,6 @@ void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col } - - int datarows = data.rows(); int datacols = data.cols(); @@ -1174,7 +1172,7 @@ RasterExtend gdalImage::getExtend() } gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, - Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite) + Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite, bool isEnvi) { if(exists_test(img_path.toUtf8().constData())) { if(overwrite) { @@ -1187,11 +1185,21 @@ gdalImage CreategdalImage(const QString& img_path, int height, int width, int ba } GDALAllRegister(); CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7 - GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); + GDALDriver* poDriver = nullptr; + if (isEnvi) { + poDriver = GetGDALDriverManager()->GetDriverByName("ENVI"); + } + else { + poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); + } + + GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_Float32, NULL); // 锟斤拷锟斤拷锟斤拷 if(need_gt) { - poDstDS->SetProjection(projection.toUtf8().constData()); + if (!projection.isEmpty()) { + poDstDS->SetProjection(projection.toUtf8().constData()); + } double gt_ptr[6] = { 0 }; for(int i = 0; i < gt.rows(); i++) { for(int j = 0; j < gt.cols(); j++) { @@ -1932,6 +1940,18 @@ ErrorCode MergeRasterInGeoCoding(QVector imgdslist, gdalImage resulti return ErrorCode::SUCCESS; } +bool saveEigenMatrixXd2Bin(Eigen::MatrixXd data, QString dataStrPath) +{ + + Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3); + + gdalImage img=CreategdalImage(dataStrPath, data.rows(), data.cols(), 1, gt, "", false, false,true); + + img.saveImage(data, 0,0,1); + + return true; +} + gdalImageComplex::gdalImageComplex(const QString& raster_path) { omp_lock_t lock; @@ -1942,7 +1962,7 @@ gdalImageComplex::gdalImageComplex(const QString& raster_path) GDALAllRegister(); CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen( - raster_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 + raster_path.toUtf8().constData(), GA_ReadOnly)); this->width = rasterDataset->GetRasterXSize(); this->height = rasterDataset->GetRasterYSize(); this->band_num = rasterDataset->GetRasterCount(); diff --git a/ImageOperatorBase.h b/ImageOperatorBase.h index 3e2d33a..c26ad1d 100644 --- a/ImageOperatorBase.h +++ b/ImageOperatorBase.h @@ -226,7 +226,7 @@ public: }; // 创建影像 -gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false); +gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false, bool isEnvi = false); gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, long espgcode, GDALDataType eType=GDT_Float32, bool need_gt = true, bool overwrite = false,bool isENVI=false); gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false); @@ -268,12 +268,16 @@ ErrorCode MergeRasterProcess(QVector filepath, QString outfileptah, QSt ErrorCode MergeRasterInGeoCoding(QVector inimgs, gdalImage resultimg,gdalImage maskimg, ShowProessAbstract* dia = nullptr); +// 保存矩阵转换为envi文件，默认数据格式为double +bool saveEigenMatrixXd2Bin(Eigen::MatrixXd data, QString dataStrPath); + + + template std::shared_ptr readDataArr(gdalImage& imgds, int start_row, int start_col, int rows_count, int cols_count, int band_ids, GDALREADARRCOPYMETHOD method) { std::shared_ptr result = nullptr; - omp_lock_t lock; omp_init_lock(&lock); omp_set_lock(&lock); @@ -303,7 +307,6 @@ std::shared_ptr readDataArr(gdalImage& imgds, int start_row, int start_col, i } } - delete[] temp; } else if (gdal_datatype == GDT_UInt16) { @@ -416,40 +419,50 @@ std::shared_ptr readDataArr(gdalImage& imgds, int start_row, int start_col, i } delete[] temp; } - else if (gdal_datatype == GDT_CFloat32) { - float* temp = new float[rows_count * cols_count*2]; - demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); + //else if (gdal_datatype == GDT_CFloat32) { + // if (std::is_same>::value || std::is_same>::value) { + // float* temp = new float[rows_count * cols_count * 2]; + // demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); - result = std::shared_ptr(new T[rows_count * cols_count], delArrPtr); - if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { - std::memcpy(result.get(), temp, rows_count * cols_count); - } - else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { - long count = rows_count * cols_count; - for (long i = 0; i < count; i++) { - result.get()[i] = std::complex(temp[i * 2 ], - temp[i * 2 + 1]); - } - } - delete[] temp; - } - else if (gdal_datatype == GDT_CFloat64) { - double* temp = new double[rows_count * cols_count*2]; - demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); + // result = std::shared_ptr(new T[rows_count * cols_count], delArrPtr); + // if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { + // std::memcpy(result.get(), temp, rows_count * cols_count); + // } + // else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { + // long count = rows_count * cols_count; + // for (long i = 0; i < count; i++) { + // result.get()[i] = T(temp[i * 2], + // temp[i * 2 + 1]); + // } + // } + // delete[] temp; + // } + // else { + // result = nullptr; + // } + //} + //else if (gdal_datatype == GDT_CFloat64 ) { + // if (std::is_same>::value || std::is_same>::value) { + // double* temp = new double[rows_count * cols_count * 2]; + // demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); - result = std::shared_ptr(new T[rows_count * cols_count], delArrPtr); - if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { - std::memcpy(result.get(), temp, rows_count * cols_count); - } - else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { - long count = rows_count * cols_count; - for (long i = 0; i < count; i++) { - result.get()[i] = std::complex(temp[i * 2], - temp[i * 2 + 1]); - } - } - delete[] temp; - } + // result = std::shared_ptr(new T[rows_count * cols_count], delArrPtr); + // if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { + // std::memcpy(result.get(), temp, rows_count * cols_count); + // } + // else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { + // long count = rows_count * cols_count; + // for (long i = 0; i < count; i++) { + // result.get()[i] = T(temp[i * 2], + // temp[i * 2 + 1]); + // } + // } + // delete[] temp; + // } + // else { + // result = nullptr; + // } + //} else { } GDALClose((GDALDatasetH)rasterDataset); @@ -459,6 +472,78 @@ std::shared_ptr readDataArr(gdalImage& imgds, int start_row, int start_col, i return result; } +template +std::shared_ptr readDataArrComplex(gdalImageComplex& imgds, int start_row, int start_col, int rows_count, int cols_count, int band_ids, GDALREADARRCOPYMETHOD method) +{ + std::shared_ptr result = nullptr; + + omp_lock_t lock; + omp_init_lock(&lock); + omp_set_lock(&lock); + GDALAllRegister(); + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); + GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(imgds.img_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 + + GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); + GDALRasterBand* demBand = rasterDataset->GetRasterBand(band_ids); + + rows_count = start_row + rows_count <= imgds.height ? rows_count : imgds.height - start_row; + cols_count = start_col + cols_count <= imgds.width ? cols_count : imgds.width - start_col; + + //Eigen::MatrixXd datamatrix(rows_count, cols_count); + + if (gdal_datatype == GDT_CFloat32) { + if (std::is_same>::value || std::is_same>::value) { + float* temp = new float[rows_count * cols_count * 2]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); + + result = std::shared_ptr(new T[rows_count * cols_count], delArrPtr); + if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { + std::memcpy(result.get(), temp, rows_count * cols_count); + } + else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { + long count = rows_count * cols_count; + for (long i = 0; i < count; i++) { + result.get()[i] = T(temp[i * 2], + temp[i * 2 + 1]); + } + } + delete[] temp; + } + else { + result = nullptr; + } + } + else if (gdal_datatype == GDT_CFloat64) { + if (std::is_same>::value || std::is_same>::value) { + double* temp = new double[rows_count * cols_count * 2]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); + + result = std::shared_ptr(new T[rows_count * cols_count], delArrPtr); + if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { + std::memcpy(result.get(), temp, rows_count * cols_count); + } + else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { + long count = rows_count * cols_count; + for (long i = 0; i < count; i++) { + result.get()[i] = T(temp[i * 2], + temp[i * 2 + 1]); + } + } + delete[] temp; + } + else { + result = nullptr; + } + } + else { + } + GDALClose((GDALDatasetH)rasterDataset); + omp_unset_lock(&lock); // 锟酵放伙拷斤拷 + omp_destroy_lock(&lock); // 劫伙拷斤拷 + // GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH + return result; +} //--------------------- 图像分块 ------------------------------