Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

SampledDensityMap - Which constructor is chosen affects results. #1036

Open
ichem001 opened this issue Nov 9, 2020 · 3 comments
Open

SampledDensityMap - Which constructor is chosen affects results. #1036

ichem001 opened this issue Nov 9, 2020 · 3 comments
Assignees
Labels
bug IMP.em Related to modules/em

Comments

@ichem001
Copy link
Contributor

ichem001 commented Nov 9, 2020

SampledDensityMap has 3 constructors for its class.

(i)   SampledDensityMap(KernelType kt = GAUSSIAN);

(ii)  SampledDensityMap(const DensityHeader &header, KernelType kt = GAUSSIAN);

(iii) SampledDensityMap(const ParticlesTemp &ps, emreal resolution,
                        emreal voxel_size,
                        IMP::FloatKey mass_key = IMP::atom::Mass::get_mass_key(),
                        int sig_cutoff = 3, KernelType kt = GAUSSIAN);

When creating a map from a PDB using (iii), maps are usually aligned with the PDB as expected.
@saltzberg found a few examples where when using (ii) in combination with

void set_particles(const ParticlesTemp &ps,
                   IMP::FloatKey mass_key = IMP::atom::Mass::get_mass_key());

the generated density has a different bouding box, and the density is translated to the (0,0,0) corners. This behavior is unexpected. See Figure at the bottom.

The test test_sample_particles.py does not display this difference in behavior when writing both xxx.mrc and yyy.mrc.

Code Snippet to reproduce error:

### Using (ii) in red in the figure
sampled_input_density = IMP.em.SampledDensityMap(mrc.get_header())
sampled_input_density.set_particles(ps)
sampled_input_density.resample()
sampled_input_density.calcRMS()
 
### Using (iii) in Green in the figure
sampled_input_density2 = IMP.em.SampledDensityMap(ps,
                                                  mrc.get_header().get_resolution(),
                                                  mrc.get_header().get_spacing())
    
IMP.em.write_map(sampled_input_density, "./5610.sid1.mrc", IMP.em.MRCReaderWriter())
IMP.em.write_map(sampled_input_density2, "./5610.sid2.mrc", IMP.em.MRCReaderWriter())

Screen Shot 2020-11-09 at 1 22 42 PM

@ichem001 ichem001 added bug IMP.em Related to modules/em labels Nov 9, 2020
@ichem001 ichem001 self-assigned this Nov 9, 2020
@ichem001
Copy link
Contributor Author

ichem001 commented Nov 10, 2020

Update 1:

  • Difference in box is because constructor (ii) based its box size on the header while constructor (iii) based its box size on the particles passed.
  • Changing the origin before resampling does not fix the issue.

@ichem001
Copy link
Contributor Author

ichem001 commented Nov 16, 2020

Update 2:

  • To fix constructor (ii), the origin needs to be updated before adding particles, e.g.:
sampled_input_density = IMP.em.SampledDensityMap(mrc.get_header())
    sampled_input_density.set_origin(p2d_mrc.get_origin())
    sampled_input_density.set_particles(ps)
    sampled_input_density.resample()
    sampled_input_density.calcRMS()

Just need to compute the origin when using constructor (ii) a priori following em::SampleddensityMap::determine_grid_size(*) if experimental map has origin (0,0,0).

@benmwebb is that enough to close this? I do not think it is a bug but rather a choice, using by default the origin in the experimental EM map header -

@ichem001
Copy link
Contributor Author

ichem001 commented Dec 13, 2020

Update 3:

  • Work around that seems to be working but need further testing:

Aligning the centroid of the model with the centroid of the bounding box seems to have fixed the issue for me. By default, it seems that the (0,0,0) corner of the box align itself with the centroid of the PDB reference used. There is some behavior arising due the bounding box not being cubic and there is a spacing /2 shift (this shift was a coding choice I assume to place center of voxel at density point).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug IMP.em Related to modules/em
Projects
None yet
Development

No branches or pull requests

1 participant