From 207c3163039ee42535799c88dc3f7e773b76efd9 Mon Sep 17 00:00:00 2001 From: "Leandro F. Estrozi" Date: Mon, 28 Jul 2025 15:59:16 +0200 Subject: [PATCH 1/4] applyHelicalSymmetry and PS averaging --- src/apps/image_handler.cpp | 96 +++++++++++++++++-- src/backprojector.cpp | 189 ++++++++++++++++++------------------- src/backprojector.h | 2 +- src/displayer.cpp | 53 ++++++++++- src/displayer.h | 1 + src/ml_optimiser.cpp | 4 +- 6 files changed, 236 insertions(+), 109 deletions(-) diff --git a/src/apps/image_handler.cpp b/src/apps/image_handler.cpp index 249d50e00..61655504d 100644 --- a/src/apps/image_handler.cpp +++ b/src/apps/image_handler.cpp @@ -37,7 +37,7 @@ class image_handler_parameters public: FileName fn_in, fn_out, fn_sel, fn_img, fn_sym, fn_sub, fn_mult, fn_div, fn_add, fn_subtract, fn_mask, fn_fsc, fn_adjust_power, fn_correct_ampl, fn_fourfilter, fn_cosDPhi; int bin_avg, avg_first, avg_last, edge_x0, edge_xF, edge_y0, edge_yF, filter_edge_width, new_box, minr_ampl_corr, my_new_box_size; - bool do_add_edge, do_invert_hand, do_flipXY, do_flipmXY, do_flipZ, do_flipX, do_flipY, do_shiftCOM, do_stats, do_calc_com, do_avg_ampl, do_avg_ampl2, do_avg_ampl2_ali, do_average, do_remove_nan, do_average_all_frames, do_power, do_guinier, do_ignore_optics, do_optimise_scale_subtract, write_float16; + bool do_add_edge, do_invert_hand, do_flipXY, do_flipmXY, do_flipZ, do_flipX, do_flipY, do_shiftCOM, do_stats, do_calc_com, do_avg_ampl, do_avg_ampl2, do_avg_ampl2_ali, do_avg_phase_ali, do_average, do_remove_nan, do_average_all_frames, do_power, do_guinier, do_ignore_optics, do_optimise_scale_subtract, write_float16; RFLOAT multiply_constant, divide_constant, add_constant, subtract_constant, threshold_above, threshold_below, angpix, requested_angpix, real_angpix, force_header_angpix, lowpass, highpass, logfilter, bfactor, shift_x, shift_y, shift_z, replace_nan, randomize_at, optimise_bfactor_subtract; // PNG options RFLOAT minval, maxval, sigma_contrast, guinier_fit_minres, guinier_fit_maxres; @@ -53,6 +53,7 @@ class image_handler_parameters Image Iop; Image Imask; MultidimArray avg_ampl; + MultidimArray avg_phase; MetaDataTable MD; FourierTransformer transformer; std::map n_images; @@ -127,6 +128,7 @@ class image_handler_parameters do_avg_ampl = parser.checkOption("--avg_ampl", "Calculate average amplitude spectrum for all images?"); do_avg_ampl2 = parser.checkOption("--avg_ampl2", "Calculate average amplitude spectrum for all images?"); do_avg_ampl2_ali = parser.checkOption("--avg_ampl2_ali", "Calculate average amplitude spectrum for all aligned images?"); + do_avg_phase_ali = parser.checkOption("--avg_phase_ali", "Calculate average phases for all aligned images?"); do_average = parser.checkOption("--average", "Calculate average of all images (without alignment)"); fn_correct_ampl = parser.getOption("--correct_avg_ampl", "Correct all images with this average amplitude spectrum", ""); minr_ampl_corr = textToInteger(parser.getOption("--minr_ampl_corr", "Minimum radius (in Fourier pixels) to apply average amplitudes", "0")); @@ -911,7 +913,11 @@ class image_handler_parameters if (XSIZE(Iop()) != xdim || YSIZE(Iop()) != ydim || ZSIZE(Iop()) != zdim) REPORT_ERROR("Error: operate-image is not of the correct size"); - if (do_avg_ampl || do_avg_ampl2 || do_avg_ampl2_ali) + if (do_avg_phase_ali) + { + avg_phase.initZeros(zdim, ydim, xdim/2+1); + } + else if (do_avg_ampl || do_avg_ampl2 || do_avg_ampl2_ali) { avg_ampl.initZeros(zdim, ydim, xdim/2+1); } @@ -941,11 +947,11 @@ class image_handler_parameters if (VEC_XSIZE(com) > 2) std::cout << " z " << ZZ(com); std::cout << std::endl; } - else if (do_avg_ampl || do_avg_ampl2 || do_avg_ampl2_ali) + else if (do_avg_ampl || do_avg_ampl2 || do_avg_ampl2_ali || do_avg_phase_ali) { Iin.read(fn_img); - if (do_avg_ampl2_ali) + if (do_avg_ampl2_ali || do_avg_phase_ali) { RFLOAT xoff = 0.; RFLOAT yoff = 0.; @@ -959,6 +965,20 @@ class image_handler_parameters MAT_ELEM(A,0, 2) = xoff; MAT_ELEM(A,1, 2) = yoff; selfApplyGeometry(Iin(), A, IS_NOT_INV, DONT_WRAP); + // +1 -1 chessboard pattern multiplication to shift the FFT origin to the middle of the array + bool xpair,ypair; + for (int i = 0; i < YSIZE(Iin()); i++) + { + ypair = (i % 2 == 0); + for (int j = 0; j < XSIZE(Iin()); j++) + { + xpair = (j % 2 == 0); + if ( ! ( (xpair && ypair) || (! xpair && ! ypair) ) ) + { + DIRECT_A2D_ELEM(Iin(), i, j) = -DIRECT_A2D_ELEM(Iin(), i, j); + } + } + } } MultidimArray FT; @@ -971,6 +991,13 @@ class image_handler_parameters DIRECT_MULTIDIM_ELEM(avg_ampl, n) += abs(DIRECT_MULTIDIM_ELEM(FT, n)); } } + else if (do_avg_phase_ali) + { + FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(FT) + { + DIRECT_MULTIDIM_ELEM(avg_phase, n) += arg(DIRECT_MULTIDIM_ELEM(FT, n)); + } + } else if (do_avg_ampl2 || do_avg_ampl2_ali) { FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(FT) @@ -1080,10 +1107,65 @@ class image_handler_parameters } - if (do_avg_ampl || do_avg_ampl2 || do_avg_ampl2_ali || do_average || do_average_all_frames) + if (do_avg_ampl || do_avg_ampl2 || do_avg_ampl2_ali || do_avg_phase_ali || do_average || do_average_all_frames) { - avg_ampl /= (RFLOAT)i_img; - Iout() = avg_ampl; + if(do_avg_phase_ali) + { + avg_phase /= 180. * (RFLOAT)i_img / PI; + } + else + { + avg_ampl /= (RFLOAT)i_img; + } + if(do_avg_ampl2_ali || do_avg_phase_ali) + { + //set the size of the output to NxN instead of the FFT Nx(N/2+1) and copy the missing half coefficients + Iout().resize(ydim,xdim); + } + if(do_avg_phase_ali) + { + avg_phase.resize(ydim,xdim); + for (int i = 0; i < YSIZE(avg_phase); i++) + { + for (int j = 0; j < XSIZE(avg_phase)/2; j++) + { + DIRECT_A2D_ELEM(avg_phase, YSIZE(avg_phase)-1-i, XSIZE(avg_phase)-1-j) = DIRECT_A2D_ELEM(avg_phase, i, j); + } + } + Iout() = avg_phase; + } + else if(do_avg_ampl2_ali) + { + avg_ampl.resize(ydim,xdim); + for (int i = 0; i < YSIZE(avg_ampl); i++) + { + for (int j = 0; j < XSIZE(avg_ampl)/2; j++) + { + DIRECT_A2D_ELEM(avg_ampl, YSIZE(avg_ampl)-1-i, XSIZE(avg_ampl)-1-j) = DIRECT_A2D_ELEM(avg_ampl, i, j); + } + } + //Get maximum value for PS log-normalization. The normalisation of the transform is undone. + RFLOAT val_max = -99999; + unsigned long int size = MULTIDIM_SIZE(avg_ampl); + for (int i = 0; i < YSIZE(avg_ampl); i++) + { + for (int j = 0; j < XSIZE(avg_ampl)/2; j++) + { + if(DIRECT_A2D_ELEM(avg_ampl, i, j)*size > val_max) val_max = DIRECT_A2D_ELEM(avg_ampl, i, j)*size; + } + } + RFLOAT logOfMax = log10f(val_max+1.0f); + for (int i = 0; i < YSIZE(avg_ampl); i++) + { + for (int j = 0; j < XSIZE(avg_ampl); j++) + { + DIRECT_A2D_ELEM(avg_ampl, i, j) = log10f(DIRECT_A2D_ELEM(avg_ampl, i, j)*size+1.0f)/logOfMax; + } + } + Iout() = avg_ampl; + } else { + Iout() = avg_ampl; + } Iout.write(fn_out, -1, false, WRITE_OVERWRITE, write_float16 ? Float16: Float); } diff --git a/src/backprojector.cpp b/src/backprojector.cpp index e76ff72ad..53ee68c13 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,118 +2164,105 @@ 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) ) 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; + + Matrix2D R(4, 4); // A matrix from the list + for (int hh = h_min; hh < h_max; hh++) { - if (hh != 0) // h==0 is done before the for loop (where sum_data = data) + 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? + + //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++) { - RFLOAT rot_ang = hh * (-helical_twist); - rotation3DMatrix(rot_ang, 'Z', R); - R.setSmallValuesToZero(); // TODO: invert rotation matrix? + // 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); - // Loop over all points in the output (i.e. rotated, or summed) array - FOR_ALL_ELEMENTS_IN_ARRAY3D(sum_weight) + for (long int i=STARTINGY(sum_weight); i<=FINISHINGY(sum_weight); i++) + for (long int j=STARTINGX(sum_weight); j<=FINISHINGX(sum_weight); j++) { - x = (RFLOAT)j; // STARTINGX(sum_weight) is zero! - y = (RFLOAT)i; - z = (RFLOAT)k; - r2 = x*x + y*y + z*z; + 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) - 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 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 - if (xp < 0) + bool is_neg_x = xp < 0; + if (is_neg_x) { - // Get complex conjugated hermitian symmetry pair xp = -xp; yp = -yp; zp = -zp; - is_neg_x = true; } - else - { - 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; + int x0 = FLOOR(xp); + RFLOAT fx = xp - x0; + int x1 = x0 + 1; - y0 = FLOOR(yp); - fy = yp - y0; - y0 -= STARTINGY(data); - y1 = y0 + 1; + int y0 = FLOOR(yp); + RFLOAT fy = yp - y0; + y0 -= STARTINGY(data); + int y1 = y0 + 1; - z0 = FLOOR(zp); - fz = zp - z0; + int z0 = FLOOR(zp); z0 -= STARTINGZ(data); - 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!!!"); - } -#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); + 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 - ddd = LIN_INTERP(fz, dxy0, dxy1); + Complex ddd = LIN_INTERP(fy, dx00, dx10); 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; @@ -2291,32 +2278,40 @@ void BackProjector::applyHelicalSymmetry(int nr_helical_asu, RFLOAT helical_twis ddd = Complex(ac - bd, ab_cd - ac - bd); } // Accumulated sum of the data term - A3D_ELEM(sum_data, k, i, j) += ddd; + A2D_ELEM(slice_data,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 + 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); + } + } + #pragma omp critical + { + 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 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); + } + } + } + } + } + + // Update original arrays 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..2469b40da 100644 --- a/src/displayer.cpp +++ b/src/displayer.cpp @@ -717,15 +717,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 +761,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 ) + showSelectedFourierAmplitudesOrPhases(current_selection_type,true); + else if ( strcmp(m->label(), "Show Fourier phase angles (2x) from selected classes") == 0 ) + showSelectedFourierAmplitudesOrPhases(current_selection_type,false); else if ( strcmp(m->label(), "Save selected classes") == 0 ) { saveBackupSelection(); @@ -1432,6 +1438,49 @@ void multiViewerCanvas::showSelectedParticles(int save_selected) std::cout <<" No classes selected. First select one or more classes..." << std::endl; } +void multiViewerCanvas::showSelectedFourierAmplitudesOrPhases(int save_selected, bool ampl) +{ + MetaDataTable MDpart; + makeStarFileSelectedParticles(save_selected, MDpart); + int nparts = MDpart.numberOfObjects(); + if (nparts > 0) + { + std::string command; + int res = 0; + if (exists("./PS_tmp.star")) + { + command = "rm -f ./PS_tmp.star"; + res = system(command.c_str()); + } + FileName fn = "./PS_tmp.star"; + MDpart.write(fn); + if(ampl) + command = "relion_image_handler --i ./PS_tmp.star --o PS_tmp.mrcs --avg_ampl2_ali"; + else + command = "relion_image_handler --i ./PS_tmp.star --o PS_tmp.mrcs --avg_phase_ali"; + res = system(command.c_str()); + command = "relion_display --i PS_tmp.mrcs --scale " + floatToString(ori_scale); + command += " --sigma_contrast " + floatToString(sigma_contrast); + command += " --black " + floatToString(minval); + command += " --white " + floatToString(maxval); + switch (colour_scheme) + { + case (BLACKGREYREDSCALE): { command += " --colour_fire"; break; } + case (BLUEGREYWHITESCALE): { command += " --colour_ice"; break; } + case (BLUEGREYREDSCALE): { command += " --colour_fire-n-ice"; break; } + case (RAINBOWSCALE): { command += " --colour_rainbow"; break; } + case (CYANBLACKYELLOWSCALE): { command += " --colour_difference"; break; } + } + // send job in the background + command += " &"; + res = system(command.c_str()); + command = "rm -f ./PS_tmp.star"; + res = system(command.c_str()); + } + 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/"; diff --git a/src/displayer.h b/src/displayer.h index cf9ed5437..f12892ebf 100644 --- a/src/displayer.h +++ b/src/displayer.h @@ -277,6 +277,7 @@ class multiViewerCanvas : public basisViewerCanvas void makeStarFileSelectedParticles(int save_selected, MetaDataTable &MDpart); void saveSelectedParticles(int save_selected); void showSelectedParticles(int save_selected); + void showSelectedFourierAmplitudesOrPhases(int save_selected, bool ampl=true); void saveTrainingSet(); void saveSelected(int save_selected); void saveBackupSelection(); 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 d9a904e04f86e222f4c9b4b79c195cde45380f8f Mon Sep 17 00:00:00 2001 From: "Leandro F. Estrozi" Date: Mon, 28 Jul 2025 16:45:04 +0200 Subject: [PATCH 2/4] added PS_angpix --- src/apps/image_handler.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/apps/image_handler.cpp b/src/apps/image_handler.cpp index 61655504d..78c947151 100644 --- a/src/apps/image_handler.cpp +++ b/src/apps/image_handler.cpp @@ -827,6 +827,8 @@ class image_handler_parameters if (verb > 0) init_progress_bar(MD.numberOfObjects()); + RFLOAT PS_angpix = -1.0; + bool do_md_out = false; FOR_ALL_OBJECTS_IN_METADATA_TABLE(MD) { @@ -852,6 +854,8 @@ class image_handler_parameters Image Ihead; Ihead.read(fn_img, false); Ihead.getDimensions(xdim, ydim, zdim, ndim); + if(do_avg_ampl2_ali || do_avg_phase_ali) + PS_angpix = Ihead.samplingRateX(); if (zdim > 1 && (do_add_edge || do_flipXY || do_flipmXY)) REPORT_ERROR("ERROR: you cannot perform 2D operations like --add_edge, --flipXY or --flipmXY on 3D maps. If you intended to operate on a movie, use .mrcs extensions for stacks!"); @@ -1166,6 +1170,7 @@ class image_handler_parameters } else { Iout() = avg_ampl; } + Iout.setSamplingRateInHeader(PS_angpix); Iout.write(fn_out, -1, false, WRITE_OVERWRITE, write_float16 ? Float16: Float); } From 22da7236132798a9c97de5f2997223b575f6722c Mon Sep 17 00:00:00 2001 From: "Leandro F. Estrozi" Date: Mon, 28 Jul 2025 19:00:01 +0200 Subject: [PATCH 3/4] added PS_padding_factor --- src/apps/image_handler.cpp | 57 +++++++++++++++++++------------------- 1 file changed, 28 insertions(+), 29 deletions(-) diff --git a/src/apps/image_handler.cpp b/src/apps/image_handler.cpp index 78c947151..2d811fb28 100644 --- a/src/apps/image_handler.cpp +++ b/src/apps/image_handler.cpp @@ -36,7 +36,7 @@ class image_handler_parameters { public: FileName fn_in, fn_out, fn_sel, fn_img, fn_sym, fn_sub, fn_mult, fn_div, fn_add, fn_subtract, fn_mask, fn_fsc, fn_adjust_power, fn_correct_ampl, fn_fourfilter, fn_cosDPhi; - int bin_avg, avg_first, avg_last, edge_x0, edge_xF, edge_y0, edge_yF, filter_edge_width, new_box, minr_ampl_corr, my_new_box_size; + int bin_avg, avg_first, avg_last, edge_x0, edge_xF, edge_y0, edge_yF, filter_edge_width, new_box, minr_ampl_corr, my_new_box_size, PS_padding_factor; bool do_add_edge, do_invert_hand, do_flipXY, do_flipmXY, do_flipZ, do_flipX, do_flipY, do_shiftCOM, do_stats, do_calc_com, do_avg_ampl, do_avg_ampl2, do_avg_ampl2_ali, do_avg_phase_ali, do_average, do_remove_nan, do_average_all_frames, do_power, do_guinier, do_ignore_optics, do_optimise_scale_subtract, write_float16; RFLOAT multiply_constant, divide_constant, add_constant, subtract_constant, threshold_above, threshold_below, angpix, requested_angpix, real_angpix, force_header_angpix, lowpass, highpass, logfilter, bfactor, shift_x, shift_y, shift_z, replace_nan, randomize_at, optimise_bfactor_subtract; // PNG options @@ -129,6 +129,7 @@ class image_handler_parameters do_avg_ampl2 = parser.checkOption("--avg_ampl2", "Calculate average amplitude spectrum for all images?"); do_avg_ampl2_ali = parser.checkOption("--avg_ampl2_ali", "Calculate average amplitude spectrum for all aligned images?"); do_avg_phase_ali = parser.checkOption("--avg_phase_ali", "Calculate average phases for all aligned images?"); + PS_padding_factor = textToInteger(parser.getOption("--PS_padding_factor", "Pad images into a this factor times bigger box before the FFT", "2")); do_average = parser.checkOption("--average", "Calculate average of all images (without alignment)"); fn_correct_ampl = parser.getOption("--correct_avg_ampl", "Correct all images with this average amplitude spectrum", ""); minr_ampl_corr = textToInteger(parser.getOption("--minr_ampl_corr", "Minimum radius (in Fourier pixels) to apply average amplitudes", "0")); @@ -919,11 +920,11 @@ class image_handler_parameters if (do_avg_phase_ali) { - avg_phase.initZeros(zdim, ydim, xdim/2+1); + avg_phase.initZeros(zdim, PS_padding_factor*ydim, (PS_padding_factor*xdim/2)+1); } else if (do_avg_ampl || do_avg_ampl2 || do_avg_ampl2_ali) { - avg_ampl.initZeros(zdim, ydim, xdim/2+1); + avg_ampl.initZeros(zdim, PS_padding_factor*ydim, (PS_padding_factor*xdim/2)+1); } else if (do_average || do_average_all_frames) { @@ -955,6 +956,7 @@ class image_handler_parameters { Iin.read(fn_img); + MultidimArray FT; if (do_avg_ampl2_ali || do_avg_phase_ali) { RFLOAT xoff = 0.; @@ -969,24 +971,20 @@ class image_handler_parameters MAT_ELEM(A,0, 2) = xoff; MAT_ELEM(A,1, 2) = yoff; selfApplyGeometry(Iin(), A, IS_NOT_INV, DONT_WRAP); - // +1 -1 chessboard pattern multiplication to shift the FFT origin to the middle of the array - bool xpair,ypair; - for (int i = 0; i < YSIZE(Iin()); i++) - { - ypair = (i % 2 == 0); - for (int j = 0; j < XSIZE(Iin()); j++) - { - xpair = (j % 2 == 0); - if ( ! ( (xpair && ypair) || (! xpair && ! ypair) ) ) - { - DIRECT_A2D_ELEM(Iin(), i, j) = -DIRECT_A2D_ELEM(Iin(), i, j); - } - } - } - } - MultidimArray FT; - transformer.FourierTransform(Iin(), FT); + MultidimArray out; + out.clear(); + CenterFFTbySign(Iin()); + padAndFloat2DMap(Iin(), out, PS_padding_factor); + long int XYdim = XSIZE(out); + transformer.FourierTransform(out, FT, false); + out.setXmippOrigin(); + out.initZeros(XYdim, XYdim); + } + else + { + transformer.FourierTransform(Iin(), FT); + } if (do_avg_ampl) { @@ -1123,12 +1121,12 @@ class image_handler_parameters } if(do_avg_ampl2_ali || do_avg_phase_ali) { - //set the size of the output to NxN instead of the FFT Nx(N/2+1) and copy the missing half coefficients - Iout().resize(ydim,xdim); + //set the output to squared dimensions N*padding_factor x N*padding_factor instead of the half-size FFTW N*padding_factor x ((N*padding_factor)/2+1) and copy/complete the missing half coefficients + Iout().resize(PS_padding_factor*ydim,PS_padding_factor*xdim); } if(do_avg_phase_ali) { - avg_phase.resize(ydim,xdim); + avg_phase.resize(PS_padding_factor*ydim,PS_padding_factor*xdim); for (int i = 0; i < YSIZE(avg_phase); i++) { for (int j = 0; j < XSIZE(avg_phase)/2; j++) @@ -1140,7 +1138,7 @@ class image_handler_parameters } else if(do_avg_ampl2_ali) { - avg_ampl.resize(ydim,xdim); + avg_ampl.resize(PS_padding_factor*ydim,PS_padding_factor*xdim); for (int i = 0; i < YSIZE(avg_ampl); i++) { for (int j = 0; j < XSIZE(avg_ampl)/2; j++) @@ -1148,22 +1146,23 @@ class image_handler_parameters DIRECT_A2D_ELEM(avg_ampl, YSIZE(avg_ampl)-1-i, XSIZE(avg_ampl)-1-j) = DIRECT_A2D_ELEM(avg_ampl, i, j); } } - //Get maximum value for PS log-normalization. The normalisation of the transform is undone. + //Get maximum value for PS log-normalization RFLOAT val_max = -99999; - unsigned long int size = MULTIDIM_SIZE(avg_ampl); for (int i = 0; i < YSIZE(avg_ampl); i++) { for (int j = 0; j < XSIZE(avg_ampl)/2; j++) { - if(DIRECT_A2D_ELEM(avg_ampl, i, j)*size > val_max) val_max = DIRECT_A2D_ELEM(avg_ampl, i, j)*size; + if(DIRECT_A2D_ELEM(avg_ampl, i, j) > val_max) val_max = DIRECT_A2D_ELEM(avg_ampl, i, j); } } - RFLOAT logOfMax = log10f(val_max+1.0f); + //The normalization of the FFT transform by fftw.cpp is undone here (multiplication by fftsize). + unsigned long int fftsize = ydim * PS_padding_factor * MULTIDIM_SIZE(avg_ampl); + RFLOAT logOfMax = log10f(val_max*fftsize+1.0f); for (int i = 0; i < YSIZE(avg_ampl); i++) { for (int j = 0; j < XSIZE(avg_ampl); j++) { - DIRECT_A2D_ELEM(avg_ampl, i, j) = log10f(DIRECT_A2D_ELEM(avg_ampl, i, j)*size+1.0f)/logOfMax; + DIRECT_A2D_ELEM(avg_ampl, i, j) = log10f(DIRECT_A2D_ELEM(avg_ampl, i, j)*fftsize+1.0f)/logOfMax; } } Iout() = avg_ampl; From 1edfcf917aed71af24d3dc3b223ff85e1fa1d9bf Mon Sep 17 00:00:00 2001 From: "Leandro F. Estrozi" Date: Tue, 29 Jul 2025 18:15:10 +0200 Subject: [PATCH 4/4] added PNG output filename indicating pixel size. Minor bugfixes --- src/apps/image_handler.cpp | 30 ++++++++++++++++++------------ src/displayer.cpp | 26 +++++++++++++++++++------- src/displayer.h | 5 ++++- 3 files changed, 41 insertions(+), 20 deletions(-) diff --git a/src/apps/image_handler.cpp b/src/apps/image_handler.cpp index 2d811fb28..73459ea46 100644 --- a/src/apps/image_handler.cpp +++ b/src/apps/image_handler.cpp @@ -957,32 +957,38 @@ class image_handler_parameters Iin.read(fn_img); MultidimArray FT; + MultidimArray out; if (do_avg_ampl2_ali || do_avg_phase_ali) { RFLOAT xoff = 0.; RFLOAT yoff = 0.; RFLOAT psi = 0.; - MD.getValue(EMDL_ORIENT_ORIGIN_X, xoff); - MD.getValue(EMDL_ORIENT_ORIGIN_Y, yoff); + MD.getValue(EMDL_ORIENT_ORIGIN_X_ANGSTROM, xoff); + MD.getValue(EMDL_ORIENT_ORIGIN_Y_ANGSTROM, yoff); MD.getValue(EMDL_ORIENT_PSI, psi); // Apply the actual transformation Matrix2D A; rotation2DMatrix(psi, A); - MAT_ELEM(A,0, 2) = xoff; - MAT_ELEM(A,1, 2) = yoff; + MAT_ELEM(A, 0, 2) = COSD(psi) * xoff - SIND(psi) * yoff; + MAT_ELEM(A, 1, 2) = COSD(psi) * yoff + SIND(psi) * xoff; selfApplyGeometry(Iin(), A, IS_NOT_INV, DONT_WRAP); - - MultidimArray out; + transformer.clear(); + FT.clear(); out.clear(); CenterFFTbySign(Iin()); padAndFloat2DMap(Iin(), out, PS_padding_factor); + if ( (XSIZE(out) != YSIZE(out)) || (ZSIZE(out) > 1) || (NSIZE(out) > 1) ) + REPORT_ERROR("fftw.cpp::amplitudeOrPhaseMap(): ERROR MultidimArray should be 2D square."); long int XYdim = XSIZE(out); transformer.FourierTransform(out, FT, false); + CenterFFTbySign(FT); out.setXmippOrigin(); out.initZeros(XYdim, XYdim); } else { + transformer.clear(); + FT.clear(); transformer.FourierTransform(Iin(), FT); } @@ -1113,7 +1119,7 @@ class image_handler_parameters { if(do_avg_phase_ali) { - avg_phase /= 180. * (RFLOAT)i_img / PI; + avg_phase /= (PI*(RFLOAT)i_img)/180.; } else { @@ -1122,28 +1128,28 @@ class image_handler_parameters if(do_avg_ampl2_ali || do_avg_phase_ali) { //set the output to squared dimensions N*padding_factor x N*padding_factor instead of the half-size FFTW N*padding_factor x ((N*padding_factor)/2+1) and copy/complete the missing half coefficients - Iout().resize(PS_padding_factor*ydim,PS_padding_factor*xdim); + Iout().resize(PS_padding_factor*ydim, PS_padding_factor*xdim); } if(do_avg_phase_ali) { - avg_phase.resize(PS_padding_factor*ydim,PS_padding_factor*xdim); + avg_phase.resize(zdim, PS_padding_factor*ydim, PS_padding_factor*xdim); for (int i = 0; i < YSIZE(avg_phase); i++) { for (int j = 0; j < XSIZE(avg_phase)/2; j++) { - DIRECT_A2D_ELEM(avg_phase, YSIZE(avg_phase)-1-i, XSIZE(avg_phase)-1-j) = DIRECT_A2D_ELEM(avg_phase, i, j); + DIRECT_A2D_ELEM(avg_phase, YSIZE(avg_phase)-i, XSIZE(avg_phase)-j) = DIRECT_A2D_ELEM(avg_phase, i, j); } } Iout() = avg_phase; } else if(do_avg_ampl2_ali) { - avg_ampl.resize(PS_padding_factor*ydim,PS_padding_factor*xdim); + avg_ampl.resize(zdim, PS_padding_factor*ydim, PS_padding_factor*xdim); for (int i = 0; i < YSIZE(avg_ampl); i++) { for (int j = 0; j < XSIZE(avg_ampl)/2; j++) { - DIRECT_A2D_ELEM(avg_ampl, YSIZE(avg_ampl)-1-i, XSIZE(avg_ampl)-1-j) = DIRECT_A2D_ELEM(avg_ampl, i, j); + DIRECT_A2D_ELEM(avg_ampl, YSIZE(avg_ampl)-i, XSIZE(avg_ampl)-j) = DIRECT_A2D_ELEM(avg_ampl, i, j); } } //Get maximum value for PS log-normalization diff --git a/src/displayer.cpp b/src/displayer.cpp index 2469b40da..73a8fd41b 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, RFLOAT angpix) { // Scroll bars Fl_Scroll scroll(0, 0, w(), h()); @@ -339,6 +339,7 @@ int basisViewerWindow::fillCanvas(int viewer_type, MetaDataTable &MDin, Observat canvas.obsModel = obsModel; canvas.text_label = text_label; canvas.metadata_table_name = MDin.getName(); + canvas.angpix = angpix; if (canvas.nr_regroups > 0) canvas.MDgroups = _MDgroups; if (_do_class) @@ -1455,11 +1456,13 @@ void multiViewerCanvas::showSelectedFourierAmplitudesOrPhases(int save_selected, FileName fn = "./PS_tmp.star"; MDpart.write(fn); if(ampl) - command = "relion_image_handler --i ./PS_tmp.star --o PS_tmp.mrcs --avg_ampl2_ali"; + command = "relion_image_handler --i ./PS_tmp.star --o ./PS_tmp_Apix_" + floatToString(angpix) + ".mrc --avg_ampl2_ali"; else - command = "relion_image_handler --i ./PS_tmp.star --o PS_tmp.mrcs --avg_phase_ali"; + command = "relion_image_handler --i ./PS_tmp.star --o ./PS_tmp_Apix_" + floatToString(angpix) + ".mrc --avg_phase_ali"; res = system(command.c_str()); - command = "relion_display --i PS_tmp.mrcs --scale " + floatToString(ori_scale); + command = "relion_image_handler --i ./PS_tmp_Apix_" + floatToString(angpix) + ".mrc --o ./PS_tmp_Apix_" + floatToString(angpix) + ".png"; + res = system(command.c_str()); + command = "relion_display --i PS_tmp_Apix_" + floatToString(angpix) + ".mrc --scale " + floatToString(ori_scale); command += " --sigma_contrast " + floatToString(sigma_contrast); command += " --black " + floatToString(minval); command += " --white " + floatToString(maxval); @@ -2782,6 +2785,17 @@ void Displayer::initialise() std::cout <<" Warning: cannot find model.star file for " << fn_in << " needed for regrouping..." << std::endl; } + if (fn_in.isStarFile() && angpix <= 0.) + { + if (MDdata.containsLabel(EMDL_IMAGE_OPTICS_GROUP)) + { + std::cout <<" Warning: setting angpix from 1st optics group..." << std::endl; + int optics_group; + MDdata.getValue(EMDL_IMAGE_OPTICS_GROUP, optics_group, 0); + optics_group--; + obsModel.opticsMdt.getValue(EMDL_IMAGE_PIXEL_SIZE, angpix, optics_group); + } + } // Check if input STAR file contains pixel-size information if (!do_class && (do_apply_orient || lowpass > 0 || highpass > 0)) @@ -2813,7 +2827,6 @@ void Displayer::initialise() } } - if (show_fourier_amplitudes && show_fourier_phase_angles) REPORT_ERROR("Displayer::initialise ERROR: cannot display Fourier amplitudes and phase angles at the same time!"); if (show_fourier_amplitudes || show_fourier_phase_angles) @@ -2827,7 +2840,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() @@ -3043,7 +3055,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, angpix); } else { diff --git a/src/displayer.h b/src/displayer.h index f12892ebf..b3403b207 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, RFLOAT angpix = -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 = "", @@ -238,6 +238,9 @@ class multiViewerCanvas : public basisViewerCanvas // Sjors 12mar18: read/write information-containing backup_selection for Liyi's project MetaDataTable MDbackup; + // Pixel size to calculate lowpass filter in Angstroms and translations in apply_orient + RFLOAT angpix; + // Scale for showing the original image RFLOAT ori_scale;