Skip to content

Commit

Permalink
Add average of power spectra similar to that of cryosparc job-average…
Browse files Browse the repository at this point in the history 
…-power-spectra
  • Loading branch information
@estrozi
estrozi committed Sep 26, 2024
1 parent 12cf15d commit 4eb1257
Show file tree
Hide file tree
Showing 3 changed files with 128 additions and 24 deletions.
126 changes: 111 additions & 15 deletions src/displayer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
***************************************************************************/

#include "src/displayer.h"
#include "src/ctf.h"
//#define DEBUG
//
#ifdef HAVE_PNG
Expand Down Expand Up @@ -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());
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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))
{
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -1024,22 +1044,67 @@ void multiViewerCanvas::showAverage(bool selected, bool show_stddev)
MultidimArray<RFLOAT> 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<RFLOAT> 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)
{
Expand Down Expand Up @@ -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/";
Expand Down Expand Up @@ -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<RFLOAT> 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<RFLOAT> 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()
Expand Down Expand Up @@ -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
{
Expand Down
9 changes: 8 additions & 1 deletion src/displayer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<RFLOAT> image, RFLOAT _minval, RFLOAT _maxval, RFLOAT _sigma_contrast, RFLOAT _scale);
int fillPickerViewerCanvas(MultidimArray<RFLOAT> image, RFLOAT _minval, RFLOAT _maxval, RFLOAT _sigma_contrast, RFLOAT _scale, RFLOAT _coord_scale,
int _particle_radius, bool do_startend = false, FileName _fn_coords = "",
Expand All @@ -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;
Expand Down Expand Up @@ -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();
Expand Down
17 changes: 9 additions & 8 deletions src/fftw.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1306,7 +1306,7 @@ void LoGFilterMap(MultidimArray<RFLOAT > &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);
Expand All @@ -1322,7 +1322,7 @@ void LoGFilterMap(MultidimArray<RFLOAT > &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...");
}
}

Expand Down Expand Up @@ -1547,7 +1547,7 @@ void directionalFilterMap(MultidimArray<RFLOAT > &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);
Expand All @@ -1563,7 +1563,7 @@ void directionalFilterMap(MultidimArray<RFLOAT > &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...");
}
}

Expand Down Expand Up @@ -1702,21 +1702,22 @@ void amplitudeOrPhaseMap(const MultidimArray<RFLOAT > &v, MultidimArray<RFLOAT >
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;
}
Expand Down

0 comments on commit 4eb1257

Please sign in to comment.