Add average of power spectra for helical samples.

[estrozi]

Dear Developers,

We have implemented a new feature that allows for the creation of an average power spectrum from the individual power spectra of images within selected 2D classes. This approach differs from calculating the power spectrum of a single 2D class average and is particularly useful for reducing artifacts when working with helical assemblies.

The procedure involves two steps:

First you need to check the "Apply orientations?" box in the "Relion display GUI" window before showing 2D classes then, right-click on a selected 2D class average and choose the new option "Show Fourier amplitudes (2x) from selected classes."

In the newly opened display window, select all the boxes containing the individual power spectra, right-click, and choose "Show average of selection" from the menu.

Best wishes,

Ambroise Desfosses and Leandro F. Estrozi.

[estrozi]

:-(

[biochem-fan]

It is not a good idea to change the behavior of an existing library function. Make this an optional behavior using a default argument r = 0.

[biochem-fan]

This has to be resolved.

[biochem-fan]

Thank you very much for your suggestion. I am not involved in filament codes, so I am not in a position to decide whether to merge this or not. Because others were not responding, I nonetheless reviewed and commented two obvious issues above. And most importantly:

is particularly useful for reducing artifacts when working with helical assemblies.

Could you please elaborate on this? Why does this approach reduce artifacts? Can you show examples (screenshots) of artifacts using the old approach and your new approach?

[estrozi]

Thank you very much for your suggestion. I am not involved in filament codes, so I am not in a position to decide whether to merge this or not. Because others were not responding, I nonetheless reviewed and commented two obvious issues above. And most importantly:

is particularly useful for reducing artifacts when working with helical assemblies.

Could you please elaborate on this? Why does this approach reduce artifacts? Can you show examples (screenshots) of artifacts using the old approach and your new approach?

Thanks for reviewing this.

In fact, mathematically, the FFT(PS) of the average image is different from the average of the images' FFT(PS). Both are rotation dependent but the latter is shift-invariant. Additionally, with helices, mixing up images with different rotations around the helix axis might kill some spots. I have one slide from Ambroise Desfosses illustrating that:

Could you help me please with some hints about your two other remarks? About the one about not changing the behavior of the library, I thought that my modifications didn't changed the behavior as the arguments bool _show_fourier_amplitudes and bool _show_fourier_phase_angles default to false. And about finding the pixel size correctly I trying to figure that out.

Thanks again !

Leandro.

[biochem-fan]

In fact, mathematically, the FFT(PS) of the average image is different from the average of the images' FFT(PS). Both are rotation dependent but the latter is shift-invariant.

OK, the situation is the same as "the PS of the sum of frames" vs "the sum of the PS of frames" in CTF estimation. If frame alignment is good, the former is better (coherent addition), while the latter is better when frame alignment is not great.

Additionally, with helices, mixing up images with different rotations around the helix axis might kill some spots. I have one slide from Ambroise Desfosses illustrating that:

So, the question is whether alignment by Class2D is good enough for coherent addition or not. Was your example produced by RELION? Do you see such differences only on difficult datasets or almost all datasets?

[biochem-fan]

About the one about not changing the behavior of the library

In amplitudeOrPhaseMap(const MultidimArray<RFLOAT > &v, MultidimArray<RFLOAT > &amp, int output_map_type), you added an arbitrary low pass filtering:

// Down-weight the 10 pixels around origin for better visualization
val = (1.0-exp(-r2/100.0))*FFTW2D_ELEM(Faux, ip, jp).abs();

This should be done elsewhere.

[estrozi]

In fact, mathematically, the FFT(PS) of the average image is different from the average of the images' FFT(PS). Both are rotation dependent but the latter is shift-invariant.

OK, the situation is the same as "the PS of the sum of frames" vs "the sum of the PS of frames" in CTF estimation. If frame alignment is good, the former is better (coherent addition), while the latter is better when frame alignment is not great.

Additionally, with helices, mixing up images with different rotations around the helix axis might kill some spots. I have one slide from Ambroise Desfosses illustrating that:

So, the question is whether alignment by Class2D is good enough for coherent addition or not. Was your example produced by RELION? Do you see such differences only on difficult datasets or almost all datasets?

Not really the same because it is common to observe the "mixing up images with different rotations around the helix axis" phenomenon. These rotations might be difficult to separate with 2D classification only (3D might be necessary but then knowing the correct helical parameters is important). So, Class2D is not always good enough for coherent addition. Only if you are "lucky" but we cannot count on that. These papers discuss this problem in some detail:

[ref1] Egelman,Helical reconstruction, again,Curr Opin Struct Biol, 2024 85:102788. doi: 10.1016/j.sbi.2024.102788
[ref2] Kreutzberger, Sonani, Egelman, Cryo-EM reconstruction of helical polymers: Beyond the simple cases, Q Rev Biophys.2024 57:e16. doi: 10.1017/S0033583524000155

We have observed this problem with RELION and cryosparc and it can happen with any kind of dataset as it has more to do with the underlying helical symmetry, helix radial density distribution of the structure and the presence of out-of-plane tilt. A common result of Class2D (RELION and other programs) is to nicely align helices on one side of the "tube" and make a blurred/fuzzy "other side" producing asymmetric PS with respect the meridian.

[estrozi]

About the one about not changing the behavior of the library

In amplitudeOrPhaseMap(const MultidimArray<RFLOAT > &v, MultidimArray<RFLOAT > &amp, int output_map_type), you added an arbitrary low pass filtering:

// Down-weight the 10 pixels around origin for better visualization
val = (1.0-exp(-r2/100.0))*FFTW2D_ELEM(Faux, ip, jp).abs();

This should be done elsewhere.

OK! I will move that somewhere else. Thanks!

[biochem-fan]

@estrozi Thank you very much for explanation. It makes sense.

@scheres Do you agree merging this feature once the two coding problems I pointed out above are resolved?

[biochem-fan]

I realized that this feature might be more useful in relion_image_handler than in relion_display. In this way, the sum of the power spectra can be calculated only once and saved in MRC. Users can display it later with different colors, conrast etc. If we implement the calculation in relion_display, all particles have to be loaded and summed again and again every time the user wants to adjust the contrast.

[estrozi]

I realized that this feature might be more useful in relion_image_handler than in relion_display. In this way, the sum of the power spectra can be calculated only once and saved in MRC. Users can display it later with different colors, conrast etc. If we implement the calculation in relion_display, all particles have to be loaded and summed again and again every time the user wants to adjust the contrast.

Hi, I have made the two corrections you suggested above. I'm going to update my pull request soon.
I totally agree with you about doing it also through relion_image_handler for practical reasons, however we also use relion for teaching helical indexing and doing the PS average through the GUI facilitates a lot when most students are still unfamiliar with command line tools. Can we have both? Thank a lot for the interest in the matter.

[estrozi]

PR updated :-)

[biochem-fan]

Optics groups: apparently your code does not support different pixel sizes (am I right?). This is fine but in this case you should add a check so that if some particles have a different pixel size, raise an error and stop.

Can we have both?

I noticed that you are calculating averages after images have been discretized to 8 bit (int ival = boxes[ipos]->img_data[n];). How about implementing the algorithm in floating points and call it from relion_image_handler and relion_display?

Please also avoid hard-coding ad-hoc thresholds (i.e. CTF 0.333 and 10 px radius). Make them optional arguments of the function or at least define them as constants. e.g.:

const RFLOAT CTF_CUTOFF = 0.333 // some comment to justify it

[estrozi]

Optics groups: apparently your code does not support different pixel sizes (am I right?). This is fine but in this case you should add a check so that if some particles have a different pixel size, raise an error and stop.

In principle I could make it support different pixel sizes but I'm not sure it's worth it. What the original (unmodified) showAverage function does in this case ? The first line seems to ignore such a possibility (int xsize = boxes[0]->xsize_data;). Am I right?

Can we have both?

I noticed that you are calculating averages after images have been discretized to 8 bit (int ival = boxes[ipos]->img_data[n];). How about implementing the algorithm in floating points and call it from relion_image_handler and relion_display?

It's a good idea but I would need some hints about how to do it.

Please also avoid hard-coding ad-hoc thresholds (i.e. CTF 0.333 and 10 px radius). Make them optional arguments of the function or at least define them as constants. e.g.:

const RFLOAT CTF_CUTOFF = 0.333 // some comment to justify it

OK.

[biochem-fan]

In principle I could make it support different pixel sizes but I'm not sure it's worth it. What the original (unmodified) showAverage function does in this case ? The first line seems to ignore such a possibility (int xsize = boxes[0]->xsize_data;). Am I right?

It is not worth implementing pixel size scaling; just adding a check is fine. (showAverage is an ancient function before the introduction of the optics group so check was not performed)

It's a good idea but I would need some hints about how to do it.

OK, then limit this to relion_display as a first attempt.

[estrozi]

I kept the changes to the absolute minimum by removing the ad-hoc thresholds etc. The result I get/need is pretty much the same. It should be fine now. Thanks for your help.

[biochem-fan]

@estrozi Thanks. It looks good to me.

@scheres Can we merge this?

[scheres]

Sorry, I am travelling and can't really look. If it looks good, please merge.On 3 Apr 2025 01:19, biochem_fan ***@***.***> wrote:CAUTION: This email originated from outside of the LMB: .***@***.*** Do not click links or open attachments unless you recognize the sender and know the content is safe.If you think this is a phishing email, please forward it to ***@***.*** -- @estrozi Thanks. It looks good to me. @scheres Can we merge this?—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***> biochem-fan left a comment (3dem/relion#1185) @estrozi Thanks. It looks good to me. @scheres Can we merge this? —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>

[biochem-fan]

I would like to test this myself. I will download extracted particles from https://www.ebi.ac.uk/empiar/EMPIAR-12121/ and run the patch.

[biochem-fan]

Hmm, the "bent particles" in EMPIAR-12121 did not give nice class averages, even after I fixed errors in the deposited STAR file.

Do you have clean Class2D results (or extracted particles) I can use for testing? I do not have time to process filaments from movies in EMPIAR.

[estrozi]

Hmm, the "bent particles" in EMPIAR-12121 did not give nice class averages, even after I fixed errors in the deposited STAR file.

Do you have clean Class2D results (or extracted particles) I can use for testing? I do not have time to process filaments from movies in EMPIAR.

Yes ! I will send you a link. Is there direct messages here in github ?

[biochem-fan]

Ideally the dataset should be small (up to a few GB) and public. If you want to share unpublished data, please write a mail to Takanori and Sjors (email address on the CCPEM archive).

[estrozi]

Ideally the dataset should be small (up to a few GB) and public. If you want to share unpublished data, please write a mail to Takanori and Sjors (email address on the CCPEM archive).

@biochem-fan : I need some time to find a nice example where difference between the two ways of averaging is clear. I will post it here.

[biochem-fan]

I tested your patch on your data.

I confirmed that your "the sum of power spectra" is much cleaner than "power spectrum of the sum". This is very good.

I continue to feel what I wrote above, that is, probably this functionality is better implemented as an option for relion_image_handler or relion_helix_toolbox. Doing this in relion_display is not ideal because one has to repeat Fourier transform and summation only to change the color scale. This is very wasteful, because the calculation takes ~ 100 seconds even for this tiny test data. I found getting the right scale was tricky (your suggested sigma contrast = 1 showed only the central part).

By doing this in a stand-alone program, we can implement this in a cleaner and better way. For example, "Apply orientations?" in relion_display is very crude; we should use BackProjector for proper interpolation. We can also parallelize it.

Another question is why you first show power spectra of individual particles and then ask a user to select "show average", instead of doing summation simultaneously. We have to keep all individual power spectra on memory and relion_display used > 10 GB of RAM in this test. Is this really necessary?

First you need to check the "Apply orientations?" box in the "Relion display GUI" window before showing 2D classes

Cannot you automatically enable this when a user select Show Fourier amplitudes (2x) from selected classes?

In summary, I like the idea of "the sum of power spectra" but am not sure if we should merge this in the current form. One view is that an imperfect implementation is still better than nothing. The other view is that incorporating a niche feature with many caveats is confusing to those who do not know how to use it. I leave the final decision to Sjors.

[estrozi]

I tested your patch on your data.

I confirmed that your "the sum of power spectra" is much cleaner than "power spectrum of the sum". This is very good.

I continue to feel what I wrote above, that is, probably this functionality is better implemented as an option for relion_image_handler or relion_helix_toolbox. Doing this in relion_display is not ideal because one has to repeat Fourier transform and summation only to change the color scale. This is very wasteful, because the calculation takes ~ 100 seconds even for this tiny test data. I found getting the right scale was tricky (your suggested sigma contrast = 1 showed only the central part).

By doing this in a stand-alone program, we can implement this in a cleaner and better way. For example, "Apply orientations?" in relion_display is very crude; we should use BackProjector for proper interpolation. We can also parallelize it.

Another question is why you first show power spectra of individual particles and then ask a user to select "show average", instead of doing summation simultaneously. We have to keep all individual power spectra on memory and relion_display used > 10 GB of RAM in this test. Is this really necessary?

First you need to check the "Apply orientations?" box in the "Relion display GUI" window before showing 2D classes

Cannot you automatically enable this when a user select Show Fourier amplitudes (2x) from selected classes?

In summary, I like the idea of "the sum of power spectra" but am not sure if we should merge this in the current form. One view is that an imperfect implementation is still better than nothing. The other view is that incorporating a niche feature with many caveats is confusing to those who do not know how to use it. I leave the final decision to Sjors.

@biochem-fan , I completely understand and agree with your perspective. The approach you suggested would definitely result in a better solution than what I’ve been able to achieve so far. However, I’m not yet familiar enough with the code to implement it that way. For now, I believe that an imperfect implementation is still better than no implementation at all. Hopefully, Sjors will also recognize the value and importance of having this feature, even in its current form.

I would also suggest the attached patch to make the PS in logarithm contrast. This makes much easier to see helical diffraction spots as shown in the picture. Thanks again for all your efforts!

fftw_log_contrast_patch.txt

[scheres]

How is this different from the option --avg_ampl2_ali in relion_image_handler? I agree with Takanori that doing this in relion_display is not an elegant solution.

[estrozi]

How is this different from the option --avg_ampl2_ali in relion_image_handler? I agree with Takanori that doing this in relion_display is not an elegant solution.

I didn't know that option! Thanks for the hint. I will probably be able to make something better based on it. I post here when/if I manage. ;-)

[estrozi]

Hi,

I've made two small modifications to relion_image_handler to enable generating a neatly squared average power spectrum with the (0,0) coefficient centered (rather than at the corner) using the --avg_ampl2_ali option (see details below).

I have two quick questions:

• Would this be an acceptable change to include?
• Do you think it would be considered 'elegant' to invoke relion_image_handler from displayer.cpp in order to display the average power spectrum of a selected 2D class?

Thanks in advance for your thoughts!

Best regards,

Leandro

[estrozi]

image_handler_patch.txt