From 4eb12578d275c4d15a73ee7733a156f77dee9b7f Mon Sep 17 00:00:00 2001 From: Leandro Farias Estrozi Date: Thu, 26 Sep 2024 18:02:35 +0200 Subject: [PATCH] Add average of power spectra similar to that of cryosparc job-average-power-spectra --- src/displayer.cpp | 126 ++++++++++++++++++++++++++++++++++++++++------ src/displayer.h | 9 +++- src/fftw.cpp | 17 ++++--- 3 files changed, 128 insertions(+), 24 deletions(-) diff --git a/src/displayer.cpp b/src/displayer.cpp index c0acfe939..f44cfe643 100644 --- a/src/displayer.cpp +++ b/src/displayer.cpp @@ -19,6 +19,7 @@ ***************************************************************************/ #include "src/displayer.h" +#include "src/ctf.h" //#define DEBUG // #ifdef HAVE_PNG @@ -301,7 +302,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 +333,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 +582,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 +613,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 +731,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 +775,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(); @@ -1024,22 +1044,67 @@ void multiViewerCanvas::showAverage(bool selected, bool show_stddev) MultidimArray sum2(ysize, xsize); int nn = 0; + + //I don't know how to properly retrive the pixel size here. Trying something that works in a simple case. + RFLOAT angpix = -1.0; + long int i,j; + int halfsizex,halfsizey; + RFLOAT x,y; + RFLOAT xs,ys; + CTF ctf; + MultidimArray sumN(ysize, xsize); + if(this->show_fourier_amplitudes || this->show_fourier_phase_angles) + { + int optics_group = 0; + if (obsModel->opticsMdt.containsLabel(EMDL_IMAGE_PIXEL_SIZE)) + obsModel->opticsMdt.getValue(EMDL_IMAGE_PIXEL_SIZE, angpix, optics_group); + halfsizex = xsize/2; + halfsizey = ysize/2; + xs = angpix*xsize; + ys = angpix*ysize; + //Initialized with ones to avoid division by zero later. + FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(sum) + { + DIRECT_MULTIDIM_ELEM(sumN, n) = 1.0; + } + } for (long int ipos = 0; ipos < boxes.size(); ipos++) { if (boxes[ipos]->selected == selected) { + if(this->show_fourier_amplitudes || this->show_fourier_phase_angles) ctf.readByGroup(boxes[ipos]->MDimg,obsModel); FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(sum) { int ival = boxes[ipos]->img_data[n]; if (ival < 0) ival += 256; - DIRECT_MULTIDIM_ELEM(sum, n) += ival; - DIRECT_MULTIDIM_ELEM(sum2, n) += ival * ival; + if(this->show_fourier_amplitudes || this->show_fourier_phase_angles) + { + i = n % xsize - halfsizex; + j = n / xsize - halfsizey; + x = (RFLOAT)i / xs; + y = (RFLOAT)j / ys; + if( fabs(ctf.getCTF(x, y, false, false, false, false)) > 0.3333) + { + DIRECT_MULTIDIM_ELEM(sum, n) += ival; + DIRECT_MULTIDIM_ELEM(sum2, n) += ival * ival; + DIRECT_MULTIDIM_ELEM(sumN, n)++; + } + } else { + DIRECT_MULTIDIM_ELEM(sum, n) += ival; + DIRECT_MULTIDIM_ELEM(sum2, n) += ival * ival; + } } nn++; } } - sum /= nn; - sum2 /= nn; + if(this->show_fourier_amplitudes || this->show_fourier_phase_angles) + { + sum /= sumN; + sum2 /= sumN; + } else { + sum /= nn; + sum2 /= nn; + } FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(sum) { @@ -1432,6 +1497,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/"; @@ -2771,14 +2868,13 @@ void Displayer::initialise() { if (do_pick || do_pick_startend) REPORT_ERROR("Displayer::initialise ERROR: cannot pick particles from Fourier maps!"); - if (fn_in.isStarFile()) - REPORT_ERROR("Displayer::initialise ERROR: use single 2D image files as input!"); - Image img; - img.read(fn_in, false); // dont read data yet: only header to get size - if ( (ZSIZE(img()) > 1) || (NSIZE(img()) > 1) ) - REPORT_ERROR("Displayer::initialise ERROR: cannot display Fourier maps for 3D images or stacks!"); +// if (fn_in.isStarFile()) +// REPORT_ERROR("Displayer::initialise ERROR: use single 2D image files as input!"); +// Image img; +// img.read(fn_in, false); // dont read data yet: only header to get size +// 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 +3090,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..9383dacda 100644 --- a/src/fftw.cpp +++ b/src/fftw.cpp @@ -1306,7 +1306,7 @@ void LoGFilterMap(MultidimArray &img, RFLOAT sigma, RFLOAT angpix) } else { - REPORT_ERROR("lowPassFilterMap: filtering of non-cube maps is not implemented..."); + REPORT_ERROR("LoGFilterMap: filtering of non-cube maps is not implemented..."); } } transformer.FourierTransform(img, FT, false); @@ -1322,7 +1322,7 @@ void LoGFilterMap(MultidimArray &img, RFLOAT sigma, RFLOAT angpix) } else { - REPORT_ERROR("lowPassFilterMap: filtering of non-cube maps is not implemented..."); + REPORT_ERROR("LoGFilterMap: filtering of non-cube maps is not implemented..."); } } @@ -1547,7 +1547,7 @@ void directionalFilterMap(MultidimArray &img, RFLOAT low_pass, RFLOAT a } else { - REPORT_ERROR("lowPassFilterMap: filtering of non-cube maps is not implemented..."); + REPORT_ERROR("directionalFilterMap: filtering of non-cube maps is not implemented..."); } } transformer.FourierTransform(img, FT, false); @@ -1563,7 +1563,7 @@ void directionalFilterMap(MultidimArray &img, RFLOAT low_pass, RFLOAT a } else { - REPORT_ERROR("lowPassFilterMap: filtering of non-cube maps is not implemented..."); + REPORT_ERROR("directionalFilterMap: filtering of non-cube maps is not implemented..."); } } @@ -1702,21 +1702,22 @@ void amplitudeOrPhaseMap(const MultidimArray &v, MultidimArray maxr2 = (XYdim - 1) * (XYdim - 1) / 4; FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM2D(Faux) { + long int r2 = ip * ip + jp * jp; if ( (ip > STARTINGY(out)) && (ip < FINISHINGY(out)) && (jp > STARTINGX(out)) && (jp < FINISHINGX(out)) - && ((ip * ip + jp * jp) < maxr2) ) + && (r2 < maxr2) ) { if (output_map_type == AMPLITUDE_MAP) - val = FFTW2D_ELEM(Faux, ip, jp).abs(); + // Down-weight the 10 pixels around origin for better visualization + val = (1.0-exp(-r2/100.0))*FFTW2D_ELEM(Faux, ip, jp).abs(); else if (output_map_type == PHASE_MAP) val = (180.) * (FFTW2D_ELEM(Faux, ip, jp).arg()) / PI; else REPORT_ERROR("fftw.cpp::amplitudeOrPhaseMap(): ERROR Unknown type of output map."); - A2D_ELEM(out, -ip, -jp) = A2D_ELEM(out, ip, jp) = val; } } - A2D_ELEM(out, 0, 0) = 0.; + //A2D_ELEM(out, 0, 0) = 0.; amp.clear(); amp = out; }