diff --git a/src/xmipp/libraries/reconstruction_adapt_cuda11/forward_art_zernike3d_gpu.cpp b/src/xmipp/libraries/reconstruction_adapt_cuda11/forward_art_zernike3d_gpu.cpp index ae96f3b6d..b88722e8c 100644 --- a/src/xmipp/libraries/reconstruction_adapt_cuda11/forward_art_zernike3d_gpu.cpp +++ b/src/xmipp/libraries/reconstruction_adapt_cuda11/forward_art_zernike3d_gpu.cpp @@ -33,6 +33,7 @@ #include #include #include +#include #include #include "data/cpu.h" @@ -76,6 +77,8 @@ void ProgForwardArtZernike3DGPU::readParams() debug_iter = checkParam("--debug_iter"); save_iter = getIntParam("--save_iter"); sort_last_N = getIntParam("--sort_last"); + sort_random = checkParam("--sort_random"); + fnSym = getParam("--sym"); FileName outPath = getParam("-o"); outPath = outPath.afterLastOf("/"); fnVolO = fnOutDir + "/" + outPath; @@ -107,6 +110,7 @@ void ProgForwardArtZernike3DGPU::show() const << "Zernike Degree: " << L1 << std::endl << "SH Degree: " << L2 << std::endl << "Step: " << loop_step << std::endl + << "Symmetry group: " << fnSym << std::endl << "Correct CTF: " << useCTF << std::endl << "Correct heretogeneity: " << useZernike << std::endl << "Remove negative values: " << removeNegValues << std::endl @@ -142,6 +146,7 @@ void ProgForwardArtZernike3DGPU::defineParams() addParamsLine(" [--ltk ] : Tikhonov regualrization"); addParamsLine(" [--ll1 ] : L1 regualrization"); addParamsLine(" [--lst ] : Soft threshold regualrization"); + addParamsLine(" [--sym ] : Symmetry to be considered during the reconstruction"); addParamsLine(" [--useZernike] : Correct heterogeneity with Zernike3D coefficients"); addParamsLine(" [--useCTF] : Correct CTF during ART reconstruction"); addParamsLine(" [--phaseFlipped] : Input images have been phase flipped"); @@ -152,6 +157,7 @@ void ProgForwardArtZernike3DGPU::defineParams() addParamsLine(" [--save_iter ] : Save intermidiate volume after #save_iter iterations"); addParamsLine( " [--sort_last ] : The algorithm sorts projections in the most orthogonally possible way. "); + addParamsLine(" [--sort_random] : Random sort of projections"); addParamsLine( " : The most orthogonal way is defined as choosing the projection which " "maximizes the "); @@ -277,8 +283,8 @@ void ProgForwardArtZernike3DGPU::preProcess() mask2D.setXmippOrigin(); vecSize = 0; - numCoefficients(L1, L2, vecSize); - fillVectorTerms(L1, L2, vL1, vN, vL2, vM); + numCoefficients(L1, L2); + fillVectorTerms(L1, L2); initX = STARTINGX(Vrefined()); endX = FINISHINGX(Vrefined()); @@ -336,6 +342,20 @@ void ProgForwardArtZernike3DGPU::preProcess() } + // Symmetry list + Matrix2D L(3, 3), R(3, 3); + L.initIdentity(3); + R.initIdentity(3); + SymList SL_aux; + SL_aux.readSymmetryFile(fnSym); + LV.push_back(L); + RV.push_back(R); + for (int isym = 0; isym < SL_aux.symsNo(); isym++) { + SL_aux.getMatrices(isym, L, R, false); + LV.push_back(L); + RV.push_back(R); + } + // Create GPU interface const cuda_forward_art_zernike3D::Program::ConstantParameters parameters = { .VRecMaskF = VRecMaskF, @@ -415,14 +435,29 @@ void ProgForwardArtZernike3DGPU::processImage(const FileName &fnImg, I.read(fnImg); I().setXmippOrigin(); - // Forward Model - artModel(); + auto rot_dbl = static_cast(rot); + auto tilt_dbl = static_cast(tilt); + auto psi_dbl = static_cast(psi); + double new_rot, new_tilt, new_psi; + + for (int isym = 0; isym < RV.size(); isym++) + { + + Euler_apply_transf(LV[isym], RV[isym], rot_dbl, tilt_dbl, psi_dbl, new_rot, new_tilt, new_psi); + + rot = static_cast(new_rot); + tilt = static_cast(new_tilt); + psi = static_cast(new_psi); - // ART update - artModel(); + // Forward Model + artModel(); + + // ART update + artModel(); + } } -void ProgForwardArtZernike3DGPU::numCoefficients(int l1, int l2, int &vecSize) +void ProgForwardArtZernike3DGPU::numCoefficients(int l1, int l2) { for (int h = 0; h <= l2; h++) { int numSPH = 2 * h + 1; @@ -436,12 +471,7 @@ void ProgForwardArtZernike3DGPU::numCoefficients(int l1, int l2, int &vecSize) } } -void ProgForwardArtZernike3DGPU::fillVectorTerms(int l1, - int l2, - Matrix1D &vL1, - Matrix1D &vN, - Matrix1D &vL2, - Matrix1D &vM) +void ProgForwardArtZernike3DGPU::fillVectorTerms(int l1, int l2) { int idx = 0; vL1.initZeros(vecSize); @@ -497,7 +527,18 @@ void ProgForwardArtZernike3DGPU::run() startProcessing(); - sortOrthogonal(); + if (sort_random) + { + std::vector imgs_id(getInputMd()->size()); + std::iota(imgs_id.begin(), imgs_id.end(), 0); + auto rng = std::default_random_engine {}; + std::shuffle(std::begin(imgs_id), std::end(imgs_id), rng); + ordered_list.initZeros(getInputMd()->size()); + for (size_t i=0; i < XSIZE(ordered_list); i++) + A1D_ELEM(ordered_list, i) = imgs_id[i]; + } + else + sortOrthogonal(); if (!oroot.empty()) { if (oext.empty()) @@ -593,12 +634,12 @@ void ProgForwardArtZernike3DGPU::sortOrthogonal() double min_prod = MAXFLOAT; ; int min_prod_proj = 0; - std::vector rot; - std::vector tilt; + std::vector rotV; + std::vector tiltV; Matrix2D v(numIMG, 3); Matrix2D euler; - getInputMd()->getColumnValues(MDL_ANGLE_ROT, rot); - getInputMd()->getColumnValues(MDL_ANGLE_TILT, tilt); + getInputMd()->getColumnValues(MDL_ANGLE_ROT, rotV); + getInputMd()->getColumnValues(MDL_ANGLE_TILT, tiltV); // Initialization ordered_list.resize(numIMG); @@ -609,7 +650,7 @@ void ProgForwardArtZernike3DGPU::sortOrthogonal() // Compute the Euler matrix for each image and keep only // the third row of each one - Euler_angles2matrix(rot[i], tilt[i], 0., euler); + Euler_angles2matrix(rotV[i], tiltV[i], 0., euler); euler.getRow(2, z); v.setRow(i, z); } diff --git a/src/xmipp/libraries/reconstruction_adapt_cuda11/forward_art_zernike3d_gpu.h b/src/xmipp/libraries/reconstruction_adapt_cuda11/forward_art_zernike3d_gpu.h index b0d309e76..60ea1c572 100644 --- a/src/xmipp/libraries/reconstruction_adapt_cuda11/forward_art_zernike3d_gpu.h +++ b/src/xmipp/libraries/reconstruction_adapt_cuda11/forward_art_zernike3d_gpu.h @@ -28,6 +28,9 @@ #define _PROG_FORWARD_ART_ZERNIKE3D_GPU #include +#include +#include +#include #include #include #include @@ -35,6 +38,7 @@ #include #include #include +#include #include #include @@ -85,10 +89,14 @@ class ProgForwardArtZernike3DGPU : public XmippMetadataProgram { double ltv, ltk, ll1, lst, lmr; // Remove negative values bool removeNegValues; + // Symmetry + FileName fnSym; + // Symmetry list + std::vector> LV, RV; public: /** Resume computations */ - bool resume; + bool resume = false; // Number of ART iterations int niter; // Sort last N projections @@ -121,6 +129,8 @@ class ProgForwardArtZernike3DGPU : public XmippMetadataProgram { bool flip; // CTF Check bool hasCTF; + // Random image sorting + bool sort_random; // Original defocus double defocusU, defocusV, defocusAngle; // CTF @@ -134,7 +144,7 @@ class ProgForwardArtZernike3DGPU : public XmippMetadataProgram { // Vector containing the degree of the spherical harmonics std::vector clnm; // Show optimization - bool showOptimization; + bool showOptimization = false; // Row ids ordered in a orthogonal fashion MultidimArray ordered_list; // Save iter counter @@ -182,23 +192,16 @@ class ProgForwardArtZernike3DGPU : public XmippMetadataProgram { An exception is thrown if any of the files is not found*/ void preProcess(); - /** Create the processing working files. - * The working files are: - * nmaTodo.xmd for images to process (nmaTodo = mdIn - nmaDone) - * nmaDone.xmd image already processed (could exists from a previous run) - */ - // virtual void createWorkFiles(); - /** Predict angles and shift. At the input the pose parameters must have an initial guess of the parameters. At the output they have the estimated pose.*/ void processImage(const FileName &fnImg, const FileName &fnImgOut, const MDRow &rowIn, MDRow &rowOut); /// Length of coefficients vector - void numCoefficients(int l1, int l2, int &vecSize); + void numCoefficients(int l1, int l2); /// Zernike and SPH coefficients allocation - void fillVectorTerms(int l1, int l2, Matrix1D &vL1, Matrix1D &vN, Matrix1D &vL2, Matrix1D &vM); + void fillVectorTerms(int l1, int l2); ///Deform a volumen using Zernike-Spherical harmonic basis void deformVol(MultidimArray &mP,