From dc3664bbc659b11f12c9ab17912edd5983f79077 Mon Sep 17 00:00:00 2001 From: "Leandro F. Estrozi" Date: Thu, 15 May 2025 14:19:25 +0200 Subject: [PATCH 1/7] parallel_applyHelicaSymmetry --- src/backprojector.cpp | 280 ++++++++++++++++++++++-------------------- src/backprojector.h | 2 +- src/displayer.cpp | 60 ++++++++- src/displayer.h | 9 +- src/fftw.cpp | 15 ++- src/ml_optimiser.cpp | 4 +- 6 files changed, 224 insertions(+), 146 deletions(-) diff --git a/src/backprojector.cpp b/src/backprojector.cpp index e76ff72ad..0863dc7de 100644 --- a/src/backprojector.cpp +++ b/src/backprojector.cpp @@ -2140,7 +2140,7 @@ void BackProjector::symmetrise(int nr_helical_asu, RFLOAT helical_twist, RFLOAT enforceHermitianSymmetry(); // Then apply helical and point group symmetry (order irrelevant?) - applyHelicalSymmetry(nr_helical_asu, helical_twist, helical_rise); + applyHelicalSymmetry(nr_helical_asu, helical_twist, helical_rise, threads); applyPointGroupSymmetry(threads); } @@ -2164,158 +2164,166 @@ void BackProjector::enforceHermitianSymmetry() } } -void BackProjector::applyHelicalSymmetry(int nr_helical_asu, RFLOAT helical_twist, RFLOAT helical_rise) +void BackProjector::applyHelicalSymmetry(int nr_helical_asu, RFLOAT helical_twist, RFLOAT helical_rise, int threads) { - if ( (nr_helical_asu < 2) || (ref_dim != 3) ) + if ( (nr_helical_asu < 2) || (ref_dim != 3) || data.getDim() != 3) return; int rmax2 = ROUND(r_max * padding_factor) * ROUND(r_max * padding_factor); - Matrix2D R(4, 4); // A matrix from the list - MultidimArray sum_weight; - MultidimArray sum_data; - RFLOAT x, y, z, fx, fy, fz, xp, yp, zp, r2; - bool is_neg_x; - int x0, x1, y0, y1, z0, z1; - Complex d000, d001, d010, d011, d100, d101, d110, d111; - Complex dx00, dx01, dx10, dx11, dxy0, dxy1, ddd; - RFLOAT dd000, dd001, dd010, dd011, dd100, dd101, dd110, dd111; - RFLOAT ddx00, ddx01, ddx10, ddx11, ddxy0, ddxy1; + MultidimArray sum_weight(weight); + MultidimArray sum_data(data); // First symmetry operator (not stored in SL) is the identity matrix - sum_weight = weight; - sum_data = data; int h_min = -nr_helical_asu/2; int h_max = -h_min + nr_helical_asu%2; - for (int hh = h_min; hh < h_max; hh++) + #pragma omp parallel num_threads(threads) { - if (hh != 0) // h==0 is done before the for loop (where sum_data = data) + MultidimArray sum_weight_private(sum_weight); + MultidimArray sum_data_private(sum_data); + sum_weight_private.initZeros(); + sum_data_private.initZeros(); + Matrix2D R(4, 4); // A matrix from the list + + #pragma omp for schedule(dynamic) + for (int hh = h_min; hh < h_max; hh++) { - RFLOAT rot_ang = hh * (-helical_twist); - rotation3DMatrix(rot_ang, 'Z', R); - R.setSmallValuesToZero(); // TODO: invert rotation matrix? - - // Loop over all points in the output (i.e. rotated, or summed) array - FOR_ALL_ELEMENTS_IN_ARRAY3D(sum_weight) + if (hh != 0) // h==0 is done before the for loop (where sum_data = data) { - x = (RFLOAT)j; // STARTINGX(sum_weight) is zero! - y = (RFLOAT)i; - z = (RFLOAT)k; - r2 = x*x + y*y + z*z; - if (r2 <= rmax2) - { - // coords_output(x,y) = A * coords_input (xp,yp) - xp = x * R(0, 0) + y * R(0, 1) + z * R(0, 2); - yp = x * R(1, 0) + y * R(1, 1) + z * R(1, 2); - zp = x * R(2, 0) + y * R(2, 1) + z * R(2, 2); + RFLOAT rot_ang = hh * (-helical_twist); + rotation3DMatrix(rot_ang, 'Z', R); + R.setSmallValuesToZero(); // TODO: invert rotation matrix? - // Only asymmetric half is stored - if (xp < 0) - { - // Get complex conjugated hermitian symmetry pair - xp = -xp; - yp = -yp; - zp = -zp; - is_neg_x = true; - } - else + // Loop over all points in the output (i.e. rotated, or summed) array + FOR_ALL_ELEMENTS_IN_ARRAY3D(sum_weight_private) + { + RFLOAT x = (RFLOAT)j; // STARTINGX(sum_weight) is zero! + RFLOAT y = (RFLOAT)i; + RFLOAT z = (RFLOAT)k; + RFLOAT r2 = x*x + y*y + z*z; + if (r2 <= rmax2) { - is_neg_x = false; - } - - // Trilinear interpolation (with physical coords) - // Subtract STARTINGY and STARTINGZ to accelerate access to data (STARTINGX=0) - // In that way use DIRECT_A3D_ELEM, rather than A3D_ELEM - x0 = FLOOR(xp); - fx = xp - x0; - x1 = x0 + 1; - - y0 = FLOOR(yp); - fy = yp - y0; - y0 -= STARTINGY(data); - y1 = y0 + 1; - - z0 = FLOOR(zp); - fz = zp - z0; - z0 -= STARTINGZ(data); - z1 = z0 + 1; + // coords_output(x,y) = A * coords_input (xp,yp) + RFLOAT xp = x * R(0, 0) + y * R(0, 1) + z * R(0, 2); + RFLOAT yp = x * R(1, 0) + y * R(1, 1) + z * R(1, 2); + RFLOAT zp = x * R(2, 0) + y * R(2, 1) + z * R(2, 2); + + // Only asymmetric half is stored + bool is_neg_x = false; + if (xp < 0) + { + // Get complex conjugated hermitian symmetry pair + xp = -xp; + yp = -yp; + zp = -zp; + is_neg_x = true; + } + + // Trilinear interpolation (with physical coords) + // Subtract STARTINGY and STARTINGZ to accelerate access to data (STARTINGX=0) + // In that way use DIRECT_A3D_ELEM, rather than A3D_ELEM + int x0 = FLOOR(xp); + RFLOAT fx = xp - x0; + int x1 = x0 + 1; + + int y0 = FLOOR(yp); + RFLOAT fy = yp - y0; + y0 -= STARTINGY(data); + int y1 = y0 + 1; + + int z0 = FLOOR(zp); + RFLOAT fz = zp - z0; + z0 -= STARTINGZ(data); + int z1 = z0 + 1; #ifdef CHECK_SIZE - if (x0 < 0 || y0 < 0 || z0 < 0 || - x1 < 0 || y1 < 0 || z1 < 0 || - x0 >= XSIZE(data) || y0 >= YSIZE(data) || z0 >= ZSIZE(data) || - x1 >= XSIZE(data) || y1 >= YSIZE(data) || z1 >= ZSIZE(data) ) - { - std::cerr << " x0= " << x0 << " y0= " << y0 << " z0= " << z0 << std::endl; - std::cerr << " x1= " << x1 << " y1= " << y1 << " z1= " << z1 << std::endl; - data.printShape(); - REPORT_ERROR("BackProjector::applyPointGroupSymmetry: checksize!!!"); - } + if (x0 < 0 || y0 < 0 || z0 < 0 || + x1 < 0 || y1 < 0 || z1 < 0 || + x0 >= XSIZE(data) || y0 >= YSIZE(data) || z0 >= ZSIZE(data) || + x1 >= XSIZE(data) || y1 >= YSIZE(data) || z1 >= ZSIZE(data) ) + { + std::cerr << " x0= " << x0 << " y0= " << y0 << " z0= " << z0 << std::endl; + std::cerr << " x1= " << x1 << " y1= " << y1 << " z1= " << z1 << std::endl; + data.printShape(); + REPORT_ERROR("BackProjector::applyPointGroupSymmetry: checksize!!!"); + } #endif - // First interpolate (complex) data - d000 = DIRECT_A3D_ELEM(data, z0, y0, x0); - d001 = DIRECT_A3D_ELEM(data, z0, y0, x1); - d010 = DIRECT_A3D_ELEM(data, z0, y1, x0); - d011 = DIRECT_A3D_ELEM(data, z0, y1, x1); - d100 = DIRECT_A3D_ELEM(data, z1, y0, x0); - d101 = DIRECT_A3D_ELEM(data, z1, y0, x1); - d110 = DIRECT_A3D_ELEM(data, z1, y1, x0); - d111 = DIRECT_A3D_ELEM(data, z1, y1, x1); - - dx00 = LIN_INTERP(fx, d000, d001); - dx01 = LIN_INTERP(fx, d100, d101); - dx10 = LIN_INTERP(fx, d010, d011); - dx11 = LIN_INTERP(fx, d110, d111); - dxy0 = LIN_INTERP(fy, dx00, dx10); - dxy1 = LIN_INTERP(fy, dx01, dx11); - - // Take complex conjugated for half with negative x - ddd = LIN_INTERP(fz, dxy0, dxy1); - - if (is_neg_x) - ddd = conj(ddd); - - // Also apply a phase shift for helical translation along Z - if (ABS(helical_rise) > 0.) - { - RFLOAT zshift = hh * helical_rise; - zshift /= - ori_size * (RFLOAT)padding_factor; - RFLOAT dotp = 2 * PI * (z * zshift); - RFLOAT a = cos(dotp); - RFLOAT b = sin(dotp); - RFLOAT c = ddd.real; - RFLOAT d = ddd.imag; - RFLOAT ac = a * c; - RFLOAT bd = b * d; - RFLOAT ab_cd = (a + b) * (c + d); - ddd = Complex(ac - bd, ab_cd - ac - bd); - } - // Accumulated sum of the data term - A3D_ELEM(sum_data, k, i, j) += ddd; - - // Then interpolate (real) weight - dd000 = DIRECT_A3D_ELEM(weight, z0, y0, x0); - dd001 = DIRECT_A3D_ELEM(weight, z0, y0, x1); - dd010 = DIRECT_A3D_ELEM(weight, z0, y1, x0); - dd011 = DIRECT_A3D_ELEM(weight, z0, y1, x1); - dd100 = DIRECT_A3D_ELEM(weight, z1, y0, x0); - dd101 = DIRECT_A3D_ELEM(weight, z1, y0, x1); - dd110 = DIRECT_A3D_ELEM(weight, z1, y1, x0); - dd111 = DIRECT_A3D_ELEM(weight, z1, y1, x1); - - ddx00 = LIN_INTERP(fx, dd000, dd001); - ddx01 = LIN_INTERP(fx, dd100, dd101); - ddx10 = LIN_INTERP(fx, dd010, dd011); - ddx11 = LIN_INTERP(fx, dd110, dd111); - ddxy0 = LIN_INTERP(fy, ddx00, ddx10); - ddxy1 = LIN_INTERP(fy, ddx01, ddx11); - - A3D_ELEM(sum_weight, k, i, j) += LIN_INTERP(fz, ddxy0, ddxy1); - - } // end if r2 <= rmax2 - } // end loop over all elements of sum_weight - } // end if hh!=0 - } // end loop over hh + // First interpolate (complex) data + Complex d000 = DIRECT_A3D_ELEM(data, z0, y0, x0); + Complex d001 = DIRECT_A3D_ELEM(data, z0, y0, x1); + Complex d010 = DIRECT_A3D_ELEM(data, z0, y1, x0); + Complex d011 = DIRECT_A3D_ELEM(data, z0, y1, x1); + Complex d100 = DIRECT_A3D_ELEM(data, z1, y0, x0); + Complex d101 = DIRECT_A3D_ELEM(data, z1, y0, x1); + Complex d110 = DIRECT_A3D_ELEM(data, z1, y1, x0); + Complex d111 = DIRECT_A3D_ELEM(data, z1, y1, x1); + + Complex dx00 = LIN_INTERP(fx, d000, d001); + Complex dx01 = LIN_INTERP(fx, d100, d101); + Complex dx10 = LIN_INTERP(fx, d010, d011); + Complex dx11 = LIN_INTERP(fx, d110, d111); + Complex dxy0 = LIN_INTERP(fy, dx00, dx10); + Complex dxy1 = LIN_INTERP(fy, dx01, dx11); + + // Take complex conjugated for half with negative x + Complex ddd = LIN_INTERP(fz, dxy0, dxy1); + + if (is_neg_x) + ddd = conj(ddd); + + // Also apply a phase shift for helical translation along Z + if (ABS(helical_rise) > 0.) + { + RFLOAT zshift = hh * helical_rise; + zshift /= - ori_size * (RFLOAT)padding_factor; + RFLOAT dotp = 2 * PI * (z * zshift); + RFLOAT a = cos(dotp); + RFLOAT b = sin(dotp); + RFLOAT c = ddd.real; + RFLOAT d = ddd.imag; + RFLOAT ac = a * c; + RFLOAT bd = b * d; + RFLOAT ab_cd = (a + b) * (c + d); + ddd = Complex(ac - bd, ab_cd - ac - bd); + } + // Accumulated sum of the data term + A3D_ELEM(sum_data_private, k, i, j) += ddd; + + // Then interpolate (real) weight + RFLOAT dd000 = DIRECT_A3D_ELEM(weight, z0, y0, x0); + RFLOAT dd001 = DIRECT_A3D_ELEM(weight, z0, y0, x1); + RFLOAT dd010 = DIRECT_A3D_ELEM(weight, z0, y1, x0); + RFLOAT dd011 = DIRECT_A3D_ELEM(weight, z0, y1, x1); + RFLOAT dd100 = DIRECT_A3D_ELEM(weight, z1, y0, x0); + RFLOAT dd101 = DIRECT_A3D_ELEM(weight, z1, y0, x1); + RFLOAT dd110 = DIRECT_A3D_ELEM(weight, z1, y1, x0); + RFLOAT dd111 = DIRECT_A3D_ELEM(weight, z1, y1, x1); + + RFLOAT ddx00 = LIN_INTERP(fx, dd000, dd001); + RFLOAT ddx01 = LIN_INTERP(fx, dd100, dd101); + RFLOAT ddx10 = LIN_INTERP(fx, dd010, dd011); + RFLOAT ddx11 = LIN_INTERP(fx, dd110, dd111); + RFLOAT ddxy0 = LIN_INTERP(fy, ddx00, ddx10); + RFLOAT ddxy1 = LIN_INTERP(fy, ddx01, ddx11); + + A3D_ELEM(sum_weight_private, k, i, j) += LIN_INTERP(fz, ddxy0, ddxy1); + + } // end if r2 <= rmax2 + } // end loop over all elements of sum_weight + } // end if hh!=0 + } // end loop over hh + #pragma omp critical + { + FOR_ALL_ELEMENTS_IN_ARRAY3D(sum_data) + { + A3D_ELEM(sum_data, k, i, j) += A3D_ELEM(sum_data_private, k, i, j); + } + FOR_ALL_ELEMENTS_IN_ARRAY3D(sum_weight) + { + A3D_ELEM(sum_weight, k, i, j) += A3D_ELEM(sum_weight_private, k, i, j); + } + } + } data = sum_data; weight = sum_weight; diff --git a/src/backprojector.h b/src/backprojector.h index 7651d5530..832fc3c9d 100644 --- a/src/backprojector.h +++ b/src/backprojector.h @@ -369,7 +369,7 @@ class BackProjector: public Projector /* Applies helical symmetry. Note that helical_rise is in PIXELS here, as BackProjector doesn't know angpix */ - void applyHelicalSymmetry(int nr_helical_asu = 1, RFLOAT helical_twist = 0., RFLOAT helical_rise = 0.); + void applyHelicalSymmetry(int nr_helical_asu = 1, RFLOAT helical_twist = 0., RFLOAT helical_rise = 0., int threads = 1); /* Applies the symmetry from the SymList object to the weight and the data array */ diff --git a/src/displayer.cpp b/src/displayer.cpp index c0acfe939..2d8e8e373 100644 --- a/src/displayer.cpp +++ b/src/displayer.cpp @@ -301,7 +301,7 @@ int DisplayBox::unSelect() int basisViewerWindow::fillCanvas(int viewer_type, MetaDataTable &MDin, ObservationModel *obsModel, EMDLabel display_label, EMDLabel text_label, bool _do_read_whole_stacks, bool _do_apply_orient, RFLOAT _minval, RFLOAT _maxval, RFLOAT _sigma_contrast, RFLOAT _scale, RFLOAT _ori_scale, int _ncol, long int max_nr_images, RFLOAT lowpass, RFLOAT highpass, bool _do_class, MetaDataTable *_MDdata, int _nr_regroup, bool _do_recenter, bool _is_data, MetaDataTable *_MDgroups, - bool do_allow_save, FileName fn_selected_imgs, FileName fn_selected_parts, int max_nr_parts_per_class) + bool do_allow_save, FileName fn_selected_imgs, FileName fn_selected_parts, int max_nr_parts_per_class, bool _show_fourier_amplitudes, bool _show_fourier_phase_angles) { // Scroll bars Fl_Scroll scroll(0, 0, w(), h()); @@ -332,6 +332,8 @@ int basisViewerWindow::fillCanvas(int viewer_type, MetaDataTable &MDin, Observat canvas.fn_selected_imgs= fn_selected_imgs; canvas.fn_selected_parts = fn_selected_parts; canvas.max_nr_parts_per_class = max_nr_parts_per_class; + canvas.show_fourier_amplitudes = _show_fourier_amplitudes; + canvas.show_fourier_phase_angles = _show_fourier_phase_angles; canvas.fill(MDin, obsModel, display_label, text_label, _do_apply_orient, _minval, _maxval, _sigma_contrast, _scale, _ncol, _do_recenter, max_nr_images, lowpass, highpass); canvas.nr_regroups = _nr_regroup; canvas.do_recenter = _do_recenter; @@ -579,6 +581,12 @@ void basisViewerCanvas::fill(MetaDataTable &MDin, ObservationModel *obsModel, EM if (highpass > 0. && have_optics_group) highPassFilterMap(img(), highpass, angpix); + if(this->show_fourier_amplitudes) + { + amplitudeOrPhaseMap(img(), img(), AMPLITUDE_MAP); + } else { + if(this->show_fourier_phase_angles) amplitudeOrPhaseMap(img(), img(), PHASE_MAP); + } // Dont change the user-provided _minval and _maxval in the getImageContrast routine! RFLOAT myminval = _minval; RFLOAT mymaxval = _maxval; @@ -604,6 +612,11 @@ void basisViewerCanvas::fill(MetaDataTable &MDin, ObservationModel *obsModel, EM int xcoor = icol * xsize_box; DisplayBox* my_box = new DisplayBox(xcoor, ycoor, xsize_box, ysize_box, ""); + if(this->show_fourier_phase_angles) + { + myminval = -180.0; + mymaxval = 180.0; + } my_box->setData(img(), MDin.getObject(my_ipos), my_ipos, myminval, mymaxval, _scale, false); if (MDin.containsLabel(text_label)) { @@ -717,15 +730,17 @@ int multiViewerCanvas::handle(int ev) { "Show Fourier phase angles (2x)" }, { "Show helical layer line profile" }, { "Show particles from selected classes" }, + { "Show Fourier amplitudes (2x) from selected classes" }, + { "Show Fourier phase angles (2x) from selected classes" }, { "Set selection type" }, - { "Save selected classes" }, // idx = 14; change below when re-ordered!! + { "Save selected classes" }, // idx = 16; change below when re-ordered!! { "Quit" }, { 0 } }; if (!do_allow_save) { - rclick_menu[14].deactivate(); + rclick_menu[16].deactivate(); } const Fl_Menu_Item *m = rclick_menu->popup(Fl::event_x(), Fl::event_y(), 0, 0, 0); @@ -759,6 +774,10 @@ int multiViewerCanvas::handle(int ev) setSelectionType(); else if ( strcmp(m->label(), "Show particles from selected classes") == 0 ) showSelectedParticles(current_selection_type); + else if ( strcmp(m->label(), "Show Fourier amplitudes (2x) from selected classes") == 0 ) + showSelectedFourierAmplitudes(current_selection_type); + else if ( strcmp(m->label(), "Show Fourier phase angles (2x) from selected classes") == 0 ) + showSelectedFourierPhaseAngles(current_selection_type); else if ( strcmp(m->label(), "Save selected classes") == 0 ) { saveBackupSelection(); @@ -1432,6 +1451,38 @@ void multiViewerCanvas::showSelectedParticles(int save_selected) std::cout <<" No classes selected. First select one or more classes..." << std::endl; } +void multiViewerCanvas::showSelectedFourierAmplitudes(int save_selected) +{ + MetaDataTable MDpart; + makeStarFileSelectedParticles(save_selected, MDpart); + int nparts = MDpart.numberOfObjects(); + if (nparts > 0) + { + basisViewerWindow win(MULTIVIEW_WINDOW_WIDTH, MULTIVIEW_WINDOW_HEIGHT, "Amplitudes in the selected classes"); + win.fillCanvas(MULTIVIEWER, MDpart, obsModel, EMDL_IMAGE_NAME, text_label, do_read_whole_stacks, do_apply_orient, 0., 0., 0., boxes[0]->scale, ori_scale, ncol, + multi_max_nr_images, 0, 0, false, MDdata, nr_regroups, do_recenter, false, MDgroups, + do_allow_save, fn_selected_imgs, fn_selected_parts, max_nr_parts_per_class, true); + } + else + std::cout <<" No classes selected. First select one or more classes..." << std::endl; +} + +void multiViewerCanvas::showSelectedFourierPhaseAngles(int save_selected) +{ + MetaDataTable MDpart; + makeStarFileSelectedParticles(save_selected, MDpart); + int nparts = MDpart.numberOfObjects(); + if (nparts > 0) + { + basisViewerWindow win(MULTIVIEW_WINDOW_WIDTH, MULTIVIEW_WINDOW_HEIGHT, "Phase angles in the selected classes"); + win.fillCanvas(MULTIVIEWER, MDpart, obsModel, EMDL_IMAGE_NAME, text_label, do_read_whole_stacks, do_apply_orient, 0., 0., 0., boxes[0]->scale, ori_scale, ncol, + multi_max_nr_images, 0, 0, false, MDdata, nr_regroups, do_recenter, false, MDgroups, + do_allow_save, fn_selected_imgs, fn_selected_parts, max_nr_parts_per_class, false, true); + } + else + std::cout <<" No classes selected. First select one or more classes..." << std::endl; +} + void multiViewerCanvas::saveTrainingSet() { FileName fn_rootdir = "/net/dstore1/teraraid3/scheres/trainingset/"; @@ -2778,7 +2829,6 @@ void Displayer::initialise() if ( (ZSIZE(img()) > 1) || (NSIZE(img()) > 1) ) REPORT_ERROR("Displayer::initialise ERROR: cannot display Fourier maps for 3D images or stacks!"); } - } int Displayer::runGui() @@ -2994,7 +3044,7 @@ void Displayer::run() basisViewerWindow win(MULTIVIEW_WINDOW_WIDTH, MULTIVIEW_WINDOW_HEIGHT, fn_in.c_str()); win.fillCanvas(MULTIVIEWER, MDin, &obsModel, display_label, text_label, do_read_whole_stacks, do_apply_orient, minval, maxval, sigma_contrast, scale, ori_scale, ncol, max_nr_images, lowpass, highpass, do_class, &MDdata, nr_regroups, do_recenter, fn_in.contains("_data.star"), &MDgroups, - do_allow_save, fn_selected_imgs, fn_selected_parts, max_nr_parts_per_class); + do_allow_save, fn_selected_imgs, fn_selected_parts, max_nr_parts_per_class, show_fourier_amplitudes, show_fourier_phase_angles); } else { diff --git a/src/displayer.h b/src/displayer.h index cf9ed5437..d82db2b6a 100644 --- a/src/displayer.h +++ b/src/displayer.h @@ -145,7 +145,7 @@ class basisViewerWindow : public Fl_Window RFLOAT _scale, RFLOAT _ori_scale, int _ncol, long int max_nr_images = -1, RFLOAT lowpass = -1.0 , RFLOAT highpass = -1.0, bool do_class = false, MetaDataTable *MDdata = NULL, int _nr_regroup = -1, bool do_recenter = false, bool _is_data = false, MetaDataTable *MDgroups = NULL, - bool do_allow_save = false, FileName fn_selected_imgs="", FileName fn_selected_parts="", int max_nr_parts_per_class = -1); + bool do_allow_save = false, FileName fn_selected_imgs="", FileName fn_selected_parts="", int max_nr_parts_per_class = -1, bool _show_fourier_amplitudes = false, bool _show_fourier_phase_angles = false); int fillSingleViewerCanvas(MultidimArray image, RFLOAT _minval, RFLOAT _maxval, RFLOAT _sigma_contrast, RFLOAT _scale); int fillPickerViewerCanvas(MultidimArray image, RFLOAT _minval, RFLOAT _maxval, RFLOAT _sigma_contrast, RFLOAT _scale, RFLOAT _coord_scale, int _particle_radius, bool do_startend = false, FileName _fn_coords = "", @@ -168,6 +168,11 @@ class basisViewerCanvas : public Fl_Widget int ysize_box; int xoff; int yoff; + // Show Fourier amplitudes? + bool show_fourier_amplitudes; + + // Show Fourier phase angles? + bool show_fourier_phase_angles; // To get positions in scrolled canvas... Fl_Scroll *scroll; @@ -277,6 +282,8 @@ class multiViewerCanvas : public basisViewerCanvas void makeStarFileSelectedParticles(int save_selected, MetaDataTable &MDpart); void saveSelectedParticles(int save_selected); void showSelectedParticles(int save_selected); + void showSelectedFourierAmplitudes(int save_selected); + void showSelectedFourierPhaseAngles(int save_selected); void saveTrainingSet(); void saveSelected(int save_selected); void saveBackupSelection(); diff --git a/src/fftw.cpp b/src/fftw.cpp index 1a7ce755d..dcbf3722a 100644 --- a/src/fftw.cpp +++ b/src/fftw.cpp @@ -1700,6 +1700,19 @@ void amplitudeOrPhaseMap(const MultidimArray &v, MultidimArray out.setXmippOrigin(); out.initZeros(XYdim, XYdim); maxr2 = (XYdim - 1) * (XYdim - 1) / 4; + //Get maximum value for PS log-normalization. The normalisation of the transform is undone. + RFLOAT val_max = -99999; + unsigned long int size = MULTIDIM_SIZE(out); + FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM2D(Faux) + { + if ( (ip > STARTINGY(out)) && (ip < FINISHINGY(out)) + && (jp > STARTINGX(out)) && (jp < FINISHINGX(out)) + && ((ip * ip + jp * jp) < maxr2) ) + { + if(FFTW2D_ELEM(Faux, ip, jp).abs()*size > val_max) val_max = FFTW2D_ELEM(Faux, ip, jp).abs()*size; + } + } + RFLOAT logOfMax = log10f(val_max+1.0f); FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM2D(Faux) { if ( (ip > STARTINGY(out)) && (ip < FINISHINGY(out)) @@ -1707,7 +1720,7 @@ void amplitudeOrPhaseMap(const MultidimArray &v, MultidimArray && ((ip * ip + jp * jp) < maxr2) ) { if (output_map_type == AMPLITUDE_MAP) - val = FFTW2D_ELEM(Faux, ip, jp).abs(); + val = log10f(FFTW2D_ELEM(Faux, ip, jp).abs()*size+1.0f)/logOfMax; else if (output_map_type == PHASE_MAP) val = (180.) * (FFTW2D_ELEM(Faux, ip, jp).arg()) / PI; else diff --git a/src/ml_optimiser.cpp b/src/ml_optimiser.cpp index c9deb509d..4cdcf920c 100644 --- a/src/ml_optimiser.cpp +++ b/src/ml_optimiser.cpp @@ -4743,7 +4743,7 @@ void MlOptimiser::symmetriseReconstructions() wsum_model.BPref[ith_recons].applyHelicalSymmetry( mymodel.helical_nr_asu, mymodel.helical_twist[ith_recons], - mymodel.helical_rise[ith_recons] / mymodel.pixel_size); + mymodel.helical_rise[ith_recons] / mymodel.pixel_size, nr_threads); if (fn_multi_sym.size() > ith_recons) // Always false if size=0 { @@ -4767,7 +4767,7 @@ void MlOptimiser::symmetriseReconstructions() wsum_model.BPref[iclass_half].applyHelicalSymmetry( mymodel.helical_nr_asu, mymodel.helical_twist[ith_recons], - mymodel.helical_rise[ith_recons] / mymodel.pixel_size); + mymodel.helical_rise[ith_recons] / mymodel.pixel_size, nr_threads); if (fn_multi_sym.size() > ith_recons) // Always false if size=0 { From 7651c761a701b8c22e3a1577b80ac4351e4899de Mon Sep 17 00:00:00 2001 From: "Leandro F. Estrozi" Date: Thu, 15 May 2025 16:34:45 +0200 Subject: [PATCH 2/7] parallel_applyHelicaSymmetry (witouht displayer changes) --- src/displayer.cpp | 60 ++++------------------------------------------- src/displayer.h | 9 +------ src/fftw.cpp | 15 +----------- 3 files changed, 7 insertions(+), 77 deletions(-) diff --git a/src/displayer.cpp b/src/displayer.cpp index 2d8e8e373..c0acfe939 100644 --- a/src/displayer.cpp +++ b/src/displayer.cpp @@ -301,7 +301,7 @@ int DisplayBox::unSelect() int basisViewerWindow::fillCanvas(int viewer_type, MetaDataTable &MDin, ObservationModel *obsModel, EMDLabel display_label, EMDLabel text_label, bool _do_read_whole_stacks, bool _do_apply_orient, RFLOAT _minval, RFLOAT _maxval, RFLOAT _sigma_contrast, RFLOAT _scale, RFLOAT _ori_scale, int _ncol, long int max_nr_images, RFLOAT lowpass, RFLOAT highpass, bool _do_class, MetaDataTable *_MDdata, int _nr_regroup, bool _do_recenter, bool _is_data, MetaDataTable *_MDgroups, - bool do_allow_save, FileName fn_selected_imgs, FileName fn_selected_parts, int max_nr_parts_per_class, bool _show_fourier_amplitudes, bool _show_fourier_phase_angles) + bool do_allow_save, FileName fn_selected_imgs, FileName fn_selected_parts, int max_nr_parts_per_class) { // Scroll bars Fl_Scroll scroll(0, 0, w(), h()); @@ -332,8 +332,6 @@ int basisViewerWindow::fillCanvas(int viewer_type, MetaDataTable &MDin, Observat canvas.fn_selected_imgs= fn_selected_imgs; canvas.fn_selected_parts = fn_selected_parts; canvas.max_nr_parts_per_class = max_nr_parts_per_class; - canvas.show_fourier_amplitudes = _show_fourier_amplitudes; - canvas.show_fourier_phase_angles = _show_fourier_phase_angles; canvas.fill(MDin, obsModel, display_label, text_label, _do_apply_orient, _minval, _maxval, _sigma_contrast, _scale, _ncol, _do_recenter, max_nr_images, lowpass, highpass); canvas.nr_regroups = _nr_regroup; canvas.do_recenter = _do_recenter; @@ -581,12 +579,6 @@ void basisViewerCanvas::fill(MetaDataTable &MDin, ObservationModel *obsModel, EM if (highpass > 0. && have_optics_group) highPassFilterMap(img(), highpass, angpix); - if(this->show_fourier_amplitudes) - { - amplitudeOrPhaseMap(img(), img(), AMPLITUDE_MAP); - } else { - if(this->show_fourier_phase_angles) amplitudeOrPhaseMap(img(), img(), PHASE_MAP); - } // Dont change the user-provided _minval and _maxval in the getImageContrast routine! RFLOAT myminval = _minval; RFLOAT mymaxval = _maxval; @@ -612,11 +604,6 @@ void basisViewerCanvas::fill(MetaDataTable &MDin, ObservationModel *obsModel, EM int xcoor = icol * xsize_box; DisplayBox* my_box = new DisplayBox(xcoor, ycoor, xsize_box, ysize_box, ""); - if(this->show_fourier_phase_angles) - { - myminval = -180.0; - mymaxval = 180.0; - } my_box->setData(img(), MDin.getObject(my_ipos), my_ipos, myminval, mymaxval, _scale, false); if (MDin.containsLabel(text_label)) { @@ -730,17 +717,15 @@ int multiViewerCanvas::handle(int ev) { "Show Fourier phase angles (2x)" }, { "Show helical layer line profile" }, { "Show particles from selected classes" }, - { "Show Fourier amplitudes (2x) from selected classes" }, - { "Show Fourier phase angles (2x) from selected classes" }, { "Set selection type" }, - { "Save selected classes" }, // idx = 16; change below when re-ordered!! + { "Save selected classes" }, // idx = 14; change below when re-ordered!! { "Quit" }, { 0 } }; if (!do_allow_save) { - rclick_menu[16].deactivate(); + rclick_menu[14].deactivate(); } const Fl_Menu_Item *m = rclick_menu->popup(Fl::event_x(), Fl::event_y(), 0, 0, 0); @@ -774,10 +759,6 @@ int multiViewerCanvas::handle(int ev) setSelectionType(); else if ( strcmp(m->label(), "Show particles from selected classes") == 0 ) showSelectedParticles(current_selection_type); - else if ( strcmp(m->label(), "Show Fourier amplitudes (2x) from selected classes") == 0 ) - showSelectedFourierAmplitudes(current_selection_type); - else if ( strcmp(m->label(), "Show Fourier phase angles (2x) from selected classes") == 0 ) - showSelectedFourierPhaseAngles(current_selection_type); else if ( strcmp(m->label(), "Save selected classes") == 0 ) { saveBackupSelection(); @@ -1451,38 +1432,6 @@ void multiViewerCanvas::showSelectedParticles(int save_selected) std::cout <<" No classes selected. First select one or more classes..." << std::endl; } -void multiViewerCanvas::showSelectedFourierAmplitudes(int save_selected) -{ - MetaDataTable MDpart; - makeStarFileSelectedParticles(save_selected, MDpart); - int nparts = MDpart.numberOfObjects(); - if (nparts > 0) - { - basisViewerWindow win(MULTIVIEW_WINDOW_WIDTH, MULTIVIEW_WINDOW_HEIGHT, "Amplitudes in the selected classes"); - win.fillCanvas(MULTIVIEWER, MDpart, obsModel, EMDL_IMAGE_NAME, text_label, do_read_whole_stacks, do_apply_orient, 0., 0., 0., boxes[0]->scale, ori_scale, ncol, - multi_max_nr_images, 0, 0, false, MDdata, nr_regroups, do_recenter, false, MDgroups, - do_allow_save, fn_selected_imgs, fn_selected_parts, max_nr_parts_per_class, true); - } - else - std::cout <<" No classes selected. First select one or more classes..." << std::endl; -} - -void multiViewerCanvas::showSelectedFourierPhaseAngles(int save_selected) -{ - MetaDataTable MDpart; - makeStarFileSelectedParticles(save_selected, MDpart); - int nparts = MDpart.numberOfObjects(); - if (nparts > 0) - { - basisViewerWindow win(MULTIVIEW_WINDOW_WIDTH, MULTIVIEW_WINDOW_HEIGHT, "Phase angles in the selected classes"); - win.fillCanvas(MULTIVIEWER, MDpart, obsModel, EMDL_IMAGE_NAME, text_label, do_read_whole_stacks, do_apply_orient, 0., 0., 0., boxes[0]->scale, ori_scale, ncol, - multi_max_nr_images, 0, 0, false, MDdata, nr_regroups, do_recenter, false, MDgroups, - do_allow_save, fn_selected_imgs, fn_selected_parts, max_nr_parts_per_class, false, true); - } - else - std::cout <<" No classes selected. First select one or more classes..." << std::endl; -} - void multiViewerCanvas::saveTrainingSet() { FileName fn_rootdir = "/net/dstore1/teraraid3/scheres/trainingset/"; @@ -2829,6 +2778,7 @@ void Displayer::initialise() if ( (ZSIZE(img()) > 1) || (NSIZE(img()) > 1) ) REPORT_ERROR("Displayer::initialise ERROR: cannot display Fourier maps for 3D images or stacks!"); } + } int Displayer::runGui() @@ -3044,7 +2994,7 @@ void Displayer::run() basisViewerWindow win(MULTIVIEW_WINDOW_WIDTH, MULTIVIEW_WINDOW_HEIGHT, fn_in.c_str()); win.fillCanvas(MULTIVIEWER, MDin, &obsModel, display_label, text_label, do_read_whole_stacks, do_apply_orient, minval, maxval, sigma_contrast, scale, ori_scale, ncol, max_nr_images, lowpass, highpass, do_class, &MDdata, nr_regroups, do_recenter, fn_in.contains("_data.star"), &MDgroups, - do_allow_save, fn_selected_imgs, fn_selected_parts, max_nr_parts_per_class, show_fourier_amplitudes, show_fourier_phase_angles); + do_allow_save, fn_selected_imgs, fn_selected_parts, max_nr_parts_per_class); } else { diff --git a/src/displayer.h b/src/displayer.h index d82db2b6a..cf9ed5437 100644 --- a/src/displayer.h +++ b/src/displayer.h @@ -145,7 +145,7 @@ class basisViewerWindow : public Fl_Window RFLOAT _scale, RFLOAT _ori_scale, int _ncol, long int max_nr_images = -1, RFLOAT lowpass = -1.0 , RFLOAT highpass = -1.0, bool do_class = false, MetaDataTable *MDdata = NULL, int _nr_regroup = -1, bool do_recenter = false, bool _is_data = false, MetaDataTable *MDgroups = NULL, - bool do_allow_save = false, FileName fn_selected_imgs="", FileName fn_selected_parts="", int max_nr_parts_per_class = -1, bool _show_fourier_amplitudes = false, bool _show_fourier_phase_angles = false); + bool do_allow_save = false, FileName fn_selected_imgs="", FileName fn_selected_parts="", int max_nr_parts_per_class = -1); int fillSingleViewerCanvas(MultidimArray image, RFLOAT _minval, RFLOAT _maxval, RFLOAT _sigma_contrast, RFLOAT _scale); int fillPickerViewerCanvas(MultidimArray image, RFLOAT _minval, RFLOAT _maxval, RFLOAT _sigma_contrast, RFLOAT _scale, RFLOAT _coord_scale, int _particle_radius, bool do_startend = false, FileName _fn_coords = "", @@ -168,11 +168,6 @@ class basisViewerCanvas : public Fl_Widget int ysize_box; int xoff; int yoff; - // Show Fourier amplitudes? - bool show_fourier_amplitudes; - - // Show Fourier phase angles? - bool show_fourier_phase_angles; // To get positions in scrolled canvas... Fl_Scroll *scroll; @@ -282,8 +277,6 @@ class multiViewerCanvas : public basisViewerCanvas void makeStarFileSelectedParticles(int save_selected, MetaDataTable &MDpart); void saveSelectedParticles(int save_selected); void showSelectedParticles(int save_selected); - void showSelectedFourierAmplitudes(int save_selected); - void showSelectedFourierPhaseAngles(int save_selected); void saveTrainingSet(); void saveSelected(int save_selected); void saveBackupSelection(); diff --git a/src/fftw.cpp b/src/fftw.cpp index dcbf3722a..1a7ce755d 100644 --- a/src/fftw.cpp +++ b/src/fftw.cpp @@ -1700,19 +1700,6 @@ void amplitudeOrPhaseMap(const MultidimArray &v, MultidimArray out.setXmippOrigin(); out.initZeros(XYdim, XYdim); maxr2 = (XYdim - 1) * (XYdim - 1) / 4; - //Get maximum value for PS log-normalization. The normalisation of the transform is undone. - RFLOAT val_max = -99999; - unsigned long int size = MULTIDIM_SIZE(out); - FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM2D(Faux) - { - if ( (ip > STARTINGY(out)) && (ip < FINISHINGY(out)) - && (jp > STARTINGX(out)) && (jp < FINISHINGX(out)) - && ((ip * ip + jp * jp) < maxr2) ) - { - if(FFTW2D_ELEM(Faux, ip, jp).abs()*size > val_max) val_max = FFTW2D_ELEM(Faux, ip, jp).abs()*size; - } - } - RFLOAT logOfMax = log10f(val_max+1.0f); FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM2D(Faux) { if ( (ip > STARTINGY(out)) && (ip < FINISHINGY(out)) @@ -1720,7 +1707,7 @@ void amplitudeOrPhaseMap(const MultidimArray &v, MultidimArray && ((ip * ip + jp * jp) < maxr2) ) { if (output_map_type == AMPLITUDE_MAP) - val = log10f(FFTW2D_ELEM(Faux, ip, jp).abs()*size+1.0f)/logOfMax; + val = FFTW2D_ELEM(Faux, ip, jp).abs(); else if (output_map_type == PHASE_MAP) val = (180.) * (FFTW2D_ELEM(Faux, ip, jp).arg()) / PI; else From 6b05507cb6ff85fb5af4c557511716e5a9a765fd Mon Sep 17 00:00:00 2001 From: "Leandro F. Estrozi" Date: Thu, 22 May 2025 18:16:10 +0200 Subject: [PATCH 3/7] parallel_applyHelicaSymmetry threads with 2D buffers --- src/backprojector.cpp | 250 +++++++++++++++++++----------------------- 1 file changed, 114 insertions(+), 136 deletions(-) diff --git a/src/backprojector.cpp b/src/backprojector.cpp index 0863dc7de..f6f49f783 100644 --- a/src/backprojector.cpp +++ b/src/backprojector.cpp @@ -2177,154 +2177,132 @@ void BackProjector::applyHelicalSymmetry(int nr_helical_asu, RFLOAT helical_twis // First symmetry operator (not stored in SL) is the identity matrix int h_min = -nr_helical_asu/2; int h_max = -h_min + nr_helical_asu%2; - #pragma omp parallel num_threads(threads) + + Matrix2D R(4, 4); // A matrix from the list + + for (int hh = h_min; hh < h_max; hh++) { - MultidimArray sum_weight_private(sum_weight); - MultidimArray sum_data_private(sum_data); - sum_weight_private.initZeros(); - sum_data_private.initZeros(); - Matrix2D R(4, 4); // A matrix from the list - - #pragma omp for schedule(dynamic) - for (int hh = h_min; hh < h_max; hh++) + if (hh == 0) continue; // h==0 is done before the for loop (where sum_data = data) + + RFLOAT rot_ang = hh * (-helical_twist); + rotation3DMatrix(rot_ang, 'Z', R); + R.setSmallValuesToZero(); // TODO: invert rotation matrix? + + #pragma omp parallel for schedule(dynamic) num_threads(threads) default(none) shared(data,weight,sum_data,sum_weight,rmax2,helical_twist,helical_rise,ori_size,padding_factor) firstprivate(hh,R) + for (long int k=STARTINGZ(sum_weight); k<=FINISHINGZ(sum_weight); k++) { - if (hh != 0) // h==0 is done before the for loop (where sum_data = data) + // Allocate minimal 2D slices per thread instead of full 3D arrays + MultidimArray slice_weight(YSIZE(data), XSIZE(data)); + MultidimArray slice_data(YSIZE(data), XSIZE(data)); + slice_weight.xdim = sum_weight.xdim; + slice_weight.ydim = sum_weight.ydim; + slice_weight.yxdim = sum_weight.yxdim; + slice_weight.xinit = sum_weight.xinit; + slice_weight.yinit = sum_weight.yinit; + slice_data.xdim = sum_data.xdim; + slice_data.ydim = sum_data.ydim; + slice_data.yxdim = sum_data.yxdim; + slice_data.xinit = sum_data.xinit; + slice_data.yinit = sum_data.yinit; + + slice_weight.initZeros(weight.ydim,weight.xdim); + slice_data.initZeros(data.ydim,data.xdim); + + for (long int i=STARTINGY(sum_weight); i<=FINISHINGY(sum_weight); i++) + for (long int j=STARTINGX(sum_weight); j<=FINISHINGX(sum_weight); j++) { - RFLOAT rot_ang = hh * (-helical_twist); - rotation3DMatrix(rot_ang, 'Z', R); - R.setSmallValuesToZero(); // TODO: invert rotation matrix? + RFLOAT x = (RFLOAT)j; // STARTINGX(sum_weight) is zero! + RFLOAT y = (RFLOAT)i; + RFLOAT z = (RFLOAT)k; + RFLOAT r2 = x*x + y*y + z*z; - // Loop over all points in the output (i.e. rotated, or summed) array - FOR_ALL_ELEMENTS_IN_ARRAY3D(sum_weight_private) + if (r2 <= rmax2) { - RFLOAT x = (RFLOAT)j; // STARTINGX(sum_weight) is zero! - RFLOAT y = (RFLOAT)i; - RFLOAT z = (RFLOAT)k; - RFLOAT r2 = x*x + y*y + z*z; - if (r2 <= rmax2) + // coords_output(x,y) = A * coords_input (xp,yp) + RFLOAT xp = x * R(0, 0) + y * R(0, 1); + RFLOAT yp = x * R(1, 0) + y * R(1, 1); + RFLOAT zp = z; // Z remains unchanged + + // Only asymmetric half is stored + bool is_neg_x = xp < 0; + if (is_neg_x) { - // coords_output(x,y) = A * coords_input (xp,yp) - RFLOAT xp = x * R(0, 0) + y * R(0, 1) + z * R(0, 2); - RFLOAT yp = x * R(1, 0) + y * R(1, 1) + z * R(1, 2); - RFLOAT zp = x * R(2, 0) + y * R(2, 1) + z * R(2, 2); - - // Only asymmetric half is stored - bool is_neg_x = false; - if (xp < 0) - { - // Get complex conjugated hermitian symmetry pair - xp = -xp; - yp = -yp; - zp = -zp; - is_neg_x = true; - } - - // Trilinear interpolation (with physical coords) - // Subtract STARTINGY and STARTINGZ to accelerate access to data (STARTINGX=0) - // In that way use DIRECT_A3D_ELEM, rather than A3D_ELEM - int x0 = FLOOR(xp); - RFLOAT fx = xp - x0; - int x1 = x0 + 1; - - int y0 = FLOOR(yp); - RFLOAT fy = yp - y0; - y0 -= STARTINGY(data); - int y1 = y0 + 1; - - int z0 = FLOOR(zp); - RFLOAT fz = zp - z0; - z0 -= STARTINGZ(data); - int z1 = z0 + 1; + xp = -xp; + yp = -yp; + zp = -zp; + } + // Trilinear interpolation (with physical coords) + // Subtract STARTINGY and STARTINGZ to accelerate access to data (STARTINGX=0) + // In that way use DIRECT_A3D_ELEM, rather than A3D_ELEM + int x0 = FLOOR(xp); + RFLOAT fx = xp - x0; + int x1 = x0 + 1; -#ifdef CHECK_SIZE - if (x0 < 0 || y0 < 0 || z0 < 0 || - x1 < 0 || y1 < 0 || z1 < 0 || - x0 >= XSIZE(data) || y0 >= YSIZE(data) || z0 >= ZSIZE(data) || - x1 >= XSIZE(data) || y1 >= YSIZE(data) || z1 >= ZSIZE(data) ) - { - std::cerr << " x0= " << x0 << " y0= " << y0 << " z0= " << z0 << std::endl; - std::cerr << " x1= " << x1 << " y1= " << y1 << " z1= " << z1 << std::endl; - data.printShape(); - REPORT_ERROR("BackProjector::applyPointGroupSymmetry: checksize!!!"); - } -#endif - // First interpolate (complex) data - Complex d000 = DIRECT_A3D_ELEM(data, z0, y0, x0); - Complex d001 = DIRECT_A3D_ELEM(data, z0, y0, x1); - Complex d010 = DIRECT_A3D_ELEM(data, z0, y1, x0); - Complex d011 = DIRECT_A3D_ELEM(data, z0, y1, x1); - Complex d100 = DIRECT_A3D_ELEM(data, z1, y0, x0); - Complex d101 = DIRECT_A3D_ELEM(data, z1, y0, x1); - Complex d110 = DIRECT_A3D_ELEM(data, z1, y1, x0); - Complex d111 = DIRECT_A3D_ELEM(data, z1, y1, x1); - - Complex dx00 = LIN_INTERP(fx, d000, d001); - Complex dx01 = LIN_INTERP(fx, d100, d101); - Complex dx10 = LIN_INTERP(fx, d010, d011); - Complex dx11 = LIN_INTERP(fx, d110, d111); - Complex dxy0 = LIN_INTERP(fy, dx00, dx10); - Complex dxy1 = LIN_INTERP(fy, dx01, dx11); - - // Take complex conjugated for half with negative x - Complex ddd = LIN_INTERP(fz, dxy0, dxy1); - - if (is_neg_x) - ddd = conj(ddd); - - // Also apply a phase shift for helical translation along Z - if (ABS(helical_rise) > 0.) - { - RFLOAT zshift = hh * helical_rise; - zshift /= - ori_size * (RFLOAT)padding_factor; - RFLOAT dotp = 2 * PI * (z * zshift); - RFLOAT a = cos(dotp); - RFLOAT b = sin(dotp); - RFLOAT c = ddd.real; - RFLOAT d = ddd.imag; - RFLOAT ac = a * c; - RFLOAT bd = b * d; - RFLOAT ab_cd = (a + b) * (c + d); - ddd = Complex(ac - bd, ab_cd - ac - bd); - } - // Accumulated sum of the data term - A3D_ELEM(sum_data_private, k, i, j) += ddd; - - // Then interpolate (real) weight - RFLOAT dd000 = DIRECT_A3D_ELEM(weight, z0, y0, x0); - RFLOAT dd001 = DIRECT_A3D_ELEM(weight, z0, y0, x1); - RFLOAT dd010 = DIRECT_A3D_ELEM(weight, z0, y1, x0); - RFLOAT dd011 = DIRECT_A3D_ELEM(weight, z0, y1, x1); - RFLOAT dd100 = DIRECT_A3D_ELEM(weight, z1, y0, x0); - RFLOAT dd101 = DIRECT_A3D_ELEM(weight, z1, y0, x1); - RFLOAT dd110 = DIRECT_A3D_ELEM(weight, z1, y1, x0); - RFLOAT dd111 = DIRECT_A3D_ELEM(weight, z1, y1, x1); - - RFLOAT ddx00 = LIN_INTERP(fx, dd000, dd001); - RFLOAT ddx01 = LIN_INTERP(fx, dd100, dd101); - RFLOAT ddx10 = LIN_INTERP(fx, dd010, dd011); - RFLOAT ddx11 = LIN_INTERP(fx, dd110, dd111); - RFLOAT ddxy0 = LIN_INTERP(fy, ddx00, ddx10); - RFLOAT ddxy1 = LIN_INTERP(fy, ddx01, ddx11); - - A3D_ELEM(sum_weight_private, k, i, j) += LIN_INTERP(fz, ddxy0, ddxy1); - - } // end if r2 <= rmax2 - } // end loop over all elements of sum_weight - } // end if hh!=0 - } // end loop over hh - #pragma omp critical - { - FOR_ALL_ELEMENTS_IN_ARRAY3D(sum_data) - { - A3D_ELEM(sum_data, k, i, j) += A3D_ELEM(sum_data_private, k, i, j); + int y0 = FLOOR(yp); + RFLOAT fy = yp - y0; + y0 -= STARTINGY(data); + int y1 = y0 + 1; + + int z0 = FLOOR(zp); + z0 -= STARTINGZ(data); + + // First interpolate (complex) data + Complex d00 = DIRECT_A3D_ELEM(data, z0, y0, x0); + Complex d01 = DIRECT_A3D_ELEM(data, z0, y0, x1); + Complex d10 = DIRECT_A3D_ELEM(data, z0, y1, x0); + Complex d11 = DIRECT_A3D_ELEM(data, z0, y1, x1); + + Complex dx00 = LIN_INTERP(fx ,d00, d01); + Complex dx10 = LIN_INTERP(fx ,d10, d11); + + // Take complex conjugated for half with negative x + Complex ddd = LIN_INTERP(fy, dx00, dx10); + + if (is_neg_x) + ddd = conj(ddd); + + if (ABS(helical_rise) > 0.) + { + RFLOAT zshift = hh * helical_rise; + zshift /= - ori_size * (RFLOAT)padding_factor; + RFLOAT dotp = 2 * PI * (z * zshift); + RFLOAT a = cos(dotp); + RFLOAT b = sin(dotp); + RFLOAT c = ddd.real; + RFLOAT d = ddd.imag; + RFLOAT ac = a * c; + RFLOAT bd = b * d; + RFLOAT ab_cd = (a + b) * (c + d); + ddd = Complex(ac - bd, ab_cd - ac - bd); + } + // Accumulated sum of the data term + A2D_ELEM(slice_data,i,j) += ddd; + + // Then interpolate (real) weight + RFLOAT dd00 = DIRECT_A3D_ELEM(weight, z0, y0, x0); + RFLOAT dd01 = DIRECT_A3D_ELEM(weight, z0, y0, x1); + RFLOAT dd10 = DIRECT_A3D_ELEM(weight, z0, y1, x0); + RFLOAT dd11 = DIRECT_A3D_ELEM(weight, z0, y1, x1); + + RFLOAT ddx00 = LIN_INTERP(fx, dd00, dd01); + RFLOAT ddx10 = LIN_INTERP(fx, dd10, dd11); + + A2D_ELEM(slice_weight,i,j) += LIN_INTERP(fy, ddx00, ddx10); + } } - FOR_ALL_ELEMENTS_IN_ARRAY3D(sum_weight) + #pragma omp critical { - A3D_ELEM(sum_weight, k, i, j) += A3D_ELEM(sum_weight_private, k, i, j); + for (long int i=STARTINGY(sum_weight); i<=FINISHINGY(sum_weight); i++) + for (long int j=STARTINGX(sum_weight); j<=FINISHINGX(sum_weight); j++) + { + A3D_ELEM(sum_data,k,i,j) += A2D_ELEM(slice_data,i,j); + A3D_ELEM(sum_weight,k,i,j) += A2D_ELEM(slice_weight,i,j); + } } } } - + + // Update original arrays data = sum_data; weight = sum_weight; } From 693b9548ac506b1162188bde3997dee7b2d63b13 Mon Sep 17 00:00:00 2001 From: "Leandro F. Estrozi" Date: Fri, 23 May 2025 14:54:18 +0200 Subject: [PATCH 4/7] bugfix in parallel_applyHelicaSymmetry threads with 2D buffers --- src/backprojector.cpp | 38 +++++++++++++++++++++++++++++++++++--- 1 file changed, 35 insertions(+), 3 deletions(-) diff --git a/src/backprojector.cpp b/src/backprojector.cpp index f6f49f783..5593e0e86 100644 --- a/src/backprojector.cpp +++ b/src/backprojector.cpp @@ -2188,8 +2188,10 @@ void BackProjector::applyHelicalSymmetry(int nr_helical_asu, RFLOAT helical_twis rotation3DMatrix(rot_ang, 'Z', R); R.setSmallValuesToZero(); // TODO: invert rotation matrix? + //Only the positive half of the z-range has to be calculated. + //The other half is calculated by the x<0 part of the x-loop. #pragma omp parallel for schedule(dynamic) num_threads(threads) default(none) shared(data,weight,sum_data,sum_weight,rmax2,helical_twist,helical_rise,ori_size,padding_factor) firstprivate(hh,R) - for (long int k=STARTINGZ(sum_weight); k<=FINISHINGZ(sum_weight); k++) + for (long int k=0; k<=FINISHINGZ(sum_weight); k++) { // Allocate minimal 2D slices per thread instead of full 3D arrays MultidimArray slice_weight(YSIZE(data), XSIZE(data)); @@ -2295,8 +2297,38 @@ void BackProjector::applyHelicalSymmetry(int nr_helical_asu, RFLOAT helical_twis for (long int i=STARTINGY(sum_weight); i<=FINISHINGY(sum_weight); i++) for (long int j=STARTINGX(sum_weight); j<=FINISHINGX(sum_weight); j++) { - A3D_ELEM(sum_data,k,i,j) += A2D_ELEM(slice_data,i,j); - A3D_ELEM(sum_weight,k,i,j) += A2D_ELEM(slice_weight,i,j); + RFLOAT x = (RFLOAT)j; // STARTINGX(sum_weight) is zero! + RFLOAT y = (RFLOAT)i; + RFLOAT z = (RFLOAT)k; + RFLOAT r2 = x*x + y*y + z*z; + if (r2 <= rmax2) + { + A3D_ELEM(sum_data,k,i,j) += A2D_ELEM(slice_data,i,j); + A3D_ELEM(sum_weight,k,i,j) += A2D_ELEM(slice_weight,i,j); + // coords_output(x,y) = A * coords_input (xp,yp) + RFLOAT xp = x * R(0, 0) + y * R(0, 1); + RFLOAT yp = x * R(1, 0) + y * R(1, 1); + RFLOAT zp = z; // Z remains unchanged + bool is_neg_x = xp < 0; + if (is_neg_x) + { + xp = -xp; + yp = -yp; + zp = -zp; + A3D_ELEM(sum_data,-k,i,j) += A2D_ELEM(slice_data,i,j); + A3D_ELEM(sum_weight,-k,i,j) += A2D_ELEM(slice_weight,i,j); + } else { +//In the x,y loop, for xp<0 the interpolations for -z were stored at +z and that's OK (perhaps a z-flip?). +//Specially considering we are only going through z=0..FINISHINGZ(sum_weight). +//But for the special case of xp == 0 and zp != 0, it is still necessary to have the values stored +//at both +z and -z slices. So the lines below. + if(xp == 0 && zp != 0) + { + A3D_ELEM(sum_data,-k,-i,j) += A2D_ELEM(slice_data,i,j); + A3D_ELEM(sum_weight,-k,-i,j) += A2D_ELEM(slice_weight,i,j); + } + } + } } } } From d6833754af5c0d4c889f422e5b97e9ce8ac225a2 Mon Sep 17 00:00:00 2001 From: "Leandro F. Estrozi" Date: Fri, 23 May 2025 15:01:26 +0200 Subject: [PATCH 5/7] little clean-up in parallel_applyHelicaSymmetry threads with 2D buffers --- src/backprojector.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/backprojector.cpp b/src/backprojector.cpp index 5593e0e86..fbc50ced5 100644 --- a/src/backprojector.cpp +++ b/src/backprojector.cpp @@ -2307,14 +2307,9 @@ void BackProjector::applyHelicalSymmetry(int nr_helical_asu, RFLOAT helical_twis A3D_ELEM(sum_weight,k,i,j) += A2D_ELEM(slice_weight,i,j); // coords_output(x,y) = A * coords_input (xp,yp) RFLOAT xp = x * R(0, 0) + y * R(0, 1); - RFLOAT yp = x * R(1, 0) + y * R(1, 1); - RFLOAT zp = z; // Z remains unchanged bool is_neg_x = xp < 0; if (is_neg_x) { - xp = -xp; - yp = -yp; - zp = -zp; A3D_ELEM(sum_data,-k,i,j) += A2D_ELEM(slice_data,i,j); A3D_ELEM(sum_weight,-k,i,j) += A2D_ELEM(slice_weight,i,j); } else { @@ -2322,7 +2317,7 @@ void BackProjector::applyHelicalSymmetry(int nr_helical_asu, RFLOAT helical_twis //Specially considering we are only going through z=0..FINISHINGZ(sum_weight). //But for the special case of xp == 0 and zp != 0, it is still necessary to have the values stored //at both +z and -z slices. So the lines below. - if(xp == 0 && zp != 0) + if(xp == 0 && z != 0) { A3D_ELEM(sum_data,-k,-i,j) += A2D_ELEM(slice_data,i,j); A3D_ELEM(sum_weight,-k,-i,j) += A2D_ELEM(slice_weight,i,j); From c2303b407f28b4db65cc6656e630ce532e0843da Mon Sep 17 00:00:00 2001 From: "Leandro F. Estrozi" Date: Fri, 23 May 2025 18:14:41 +0200 Subject: [PATCH 6/7] last fix in parallel_applyHelicaSymmetry threads with 2D buffers --- src/backprojector.cpp | 20 +------------------- 1 file changed, 1 insertion(+), 19 deletions(-) diff --git a/src/backprojector.cpp b/src/backprojector.cpp index fbc50ced5..220b472a8 100644 --- a/src/backprojector.cpp +++ b/src/backprojector.cpp @@ -2191,7 +2191,7 @@ void BackProjector::applyHelicalSymmetry(int nr_helical_asu, RFLOAT helical_twis //Only the positive half of the z-range has to be calculated. //The other half is calculated by the x<0 part of the x-loop. #pragma omp parallel for schedule(dynamic) num_threads(threads) default(none) shared(data,weight,sum_data,sum_weight,rmax2,helical_twist,helical_rise,ori_size,padding_factor) firstprivate(hh,R) - for (long int k=0; k<=FINISHINGZ(sum_weight); k++) + for (long int k=STARTINGZ(sum_weight); k<=FINISHINGZ(sum_weight); k++) { // Allocate minimal 2D slices per thread instead of full 3D arrays MultidimArray slice_weight(YSIZE(data), XSIZE(data)); @@ -2305,24 +2305,6 @@ void BackProjector::applyHelicalSymmetry(int nr_helical_asu, RFLOAT helical_twis { A3D_ELEM(sum_data,k,i,j) += A2D_ELEM(slice_data,i,j); A3D_ELEM(sum_weight,k,i,j) += A2D_ELEM(slice_weight,i,j); - // coords_output(x,y) = A * coords_input (xp,yp) - RFLOAT xp = x * R(0, 0) + y * R(0, 1); - bool is_neg_x = xp < 0; - if (is_neg_x) - { - A3D_ELEM(sum_data,-k,i,j) += A2D_ELEM(slice_data,i,j); - A3D_ELEM(sum_weight,-k,i,j) += A2D_ELEM(slice_weight,i,j); - } else { -//In the x,y loop, for xp<0 the interpolations for -z were stored at +z and that's OK (perhaps a z-flip?). -//Specially considering we are only going through z=0..FINISHINGZ(sum_weight). -//But for the special case of xp == 0 and zp != 0, it is still necessary to have the values stored -//at both +z and -z slices. So the lines below. - if(xp == 0 && z != 0) - { - A3D_ELEM(sum_data,-k,-i,j) += A2D_ELEM(slice_data,i,j); - A3D_ELEM(sum_weight,-k,-i,j) += A2D_ELEM(slice_weight,i,j); - } - } } } } From 9f1a240da56b6e7ce5bbd149a35fba1d5c2bf28a Mon Sep 17 00:00:00 2001 From: "Leandro F. Estrozi" Date: Wed, 28 May 2025 16:53:14 +0200 Subject: [PATCH 7/7] not returning when data.getDim()==0 --- src/backprojector.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/backprojector.cpp b/src/backprojector.cpp index 220b472a8..53ee68c13 100644 --- a/src/backprojector.cpp +++ b/src/backprojector.cpp @@ -2166,7 +2166,7 @@ void BackProjector::enforceHermitianSymmetry() void BackProjector::applyHelicalSymmetry(int nr_helical_asu, RFLOAT helical_twist, RFLOAT helical_rise, int threads) { - if ( (nr_helical_asu < 2) || (ref_dim != 3) || data.getDim() != 3) + if ( (nr_helical_asu < 2) || (ref_dim != 3) ) return; int rmax2 = ROUND(r_max * padding_factor) * ROUND(r_max * padding_factor);