SampledDensityMap - Which constructor is chosen affects results.

[ichem001]

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())

[ichem001]

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]

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]

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).