[MRG; ENH] Add "add_estimated_fiducials" to DigMontage (#9118)

* Adding updated unit test.

* Fixing unit test.

* Fixing flake.

* Update mne/channels/montage.py

Co-authored-by: Robert Luke <748691+rob-luke@users.noreply.github.com>

* Update mne/channels/montage.py

Co-authored-by: Robert Luke <748691+rob-luke@users.noreply.github.com>

* Update mne/channels/montage.py

Co-authored-by: Robert Luke <748691+rob-luke@users.noreply.github.com>

* Update mne/channels/montage.py

Co-authored-by: Eric Larson <larson.eric.d@gmail.com>

* Update mne/channels/tests/test_montage.py

Co-authored-by: Eric Larson <larson.eric.d@gmail.com>

* Address mne comments.

* Adding updated docstring.

* Adding docstring for test montage.

* ADding crossref.

* Adding changelog.

* Fix order of docstring.

* Reference tutorial.'

* Does this work.

* Update mne/channels/montage.py

Co-authored-by: Eric Larson <larson.eric.d@gmail.com>

* BUG: Link

* Fixing a docstring.g

Co-authored-by: Robert Luke <748691+rob-luke@users.noreply.github.com>
Co-authored-by: Eric Larson <larson.eric.d@gmail.com>

doc/changes/latest.inc

@@ -127,7 +127,9 @@ Enhancements
 
 - `mne.viz.plot_drop_log` and :meth:`mne.Epochs.plot_drop_log` now omit displaying the subject name in the title if ``subject=None`` is passed (:gh:`9015` by `Richard Höchenberger`_)
 
-- :func:`mne.io.anonymize_info` now anonymizes also sex and hand fields when ``keep_his`` is ``False`` (:gh:`9103` by |Rotem Falach|_)`
+- Add :func:`mne.channels.DigMontage.add_estimated_fiducials` which will add LPA, RPA and Nasion fiducial points to the ``DigMontage`` object in ``mri`` coordinate frame (:gh:`9118` by `Adam Li`_)
+
+  - :func:`mne.io.anonymize_info` now anonymizes also sex and hand fields when ``keep_his`` is ``False`` (:gh:`9103` by |Rotem Falach|_)`
 
 Bugs
 ~~~~

mne/channels/montage.py

@@ -322,6 +322,7 @@ def get_positions(self):
                 {
                     'ch_pos': {'EEG061': [0, 0, 0]},
                     'nasion': [0, 0, 1],
+                    'coord_frame': 'mni_tal',
                     'lpa': [0, 1, 0],
                     'rpa': [1, 0, 0],
                     'hsp': None,
@@ -331,7 +332,6 @@ def get_positions(self):
         # get channel positions as dict
         ch_pos = self._get_ch_pos()
 
-        # _get_fid_coords(self.dig)
         # get coordframe and fiducial coordinates
         montage_bunch = _get_data_as_dict_from_dig(self.dig)
         coord_frame = _frame_to_str.get(montage_bunch.coord_frame)
@@ -348,6 +348,60 @@ def get_positions(self):
         )
         return positions
 
+    @verbose
+    def add_estimated_fiducials(self, subject, subjects_dir=None,
+                                verbose=None):
+        """Estimate fiducials based on FreeSurfer ``fsaverage`` subject.
+
+        This takes a montage with the ``mri`` coordinate frame,
+        corresponding to the FreeSurfer RAS (xyz in the volume) T1w
+        image of the specific subject. It will call
+        :func:`mne.coreg.get_mni_fiducials` to estimate LPA, RPA and
+        Nasion fiducial points.
+
+        Parameters
+        ----------
+        %(subject)s
+        %(subjects_dir)s
+        %(verbose)s
+
+        Returns
+        -------
+        inst : instance of DigMontage
+            The instance, modified in-place.
+
+        See Also
+        --------
+        :ref:`plot_source_alignment`
+
+        Notes
+        -----
+        Since MNE uses the FIF data structure, it relies on the ``head``
+        coordinate frame. Any coordinate frame can be transformed
+        to ``head`` if the fiducials (i.e. LPA, RPA and Nasion) are
+        defined. One can use this function to estimate those fiducials
+        and then use ``montage.get_native_head_t()`` to get the
+        head <-> MRI transform.
+        """
+        from ..coreg import get_mni_fiducials
+
+        # get coordframe and fiducial coordinates
+        montage_bunch = _get_data_as_dict_from_dig(self.dig)
+
+        # get the coordinate frame as a string and check that it's MRI
+        if montage_bunch.coord_frame != FIFF.FIFFV_COORD_MRI:
+            raise RuntimeError(
+                f'Montage should be in mri coordinate frame to call '
+                f'`add_estimated_fiducials`. The current coordinate '
+                f'frame is {montage_bunch.coord_frame}')
+
+        # estimate LPA, nasion, RPA from FreeSurfer fsaverage
+        fids_mri = list(get_mni_fiducials(subject, subjects_dir))
+
+        # add those digpoints to front of montage
+        self.dig = fids_mri + self.dig
+        return self
+
 
 VALID_SCALES = dict(mm=1e-3, cm=1e-2, m=1)
 

mne/channels/tests/test_montage.py

@@ -1458,4 +1458,36 @@ def test_plot_montage():
     plt.close('all')
 
 
+@testing.requires_testing_data
+def test_montage_add_estimated_fiducials():
+    """Test montage can add estimated fiducials for rpa, lpa, nas."""
+    # get the fiducials from test file
+    subjects_dir = op.join(data_path, 'subjects')
+    subject = 'sample'
+    fid_fname = op.join(subjects_dir, subject, 'bem',
+                        'sample-fiducials.fif')
+    test_fids, test_coord_frame = read_fiducials(fid_fname)
+    test_fids = np.array([f['r'] for f in test_fids])
+
+    # create test montage and add estimated fiducials
+    test_ch_pos = {'A1': [0, 0, 0]}
+    montage = make_dig_montage(ch_pos=test_ch_pos, coord_frame='mri')
+    montage.add_estimated_fiducials(subject=subject, subjects_dir=subjects_dir)
+
+    # check that these fiducials are close to the estimated fiducials
+    ch_pos = montage.get_positions()
+    fids_est = [ch_pos['lpa'], ch_pos['nasion'], ch_pos['rpa']]
+
+    dists = np.linalg.norm(test_fids - fids_est, axis=-1) * 1000.  # -> mm
+    assert (dists < 8).all(), dists
+
+    # an error should be raised if the montage is not in `mri` coord_frame
+    # which is the FreeSurfer RAS
+    montage = make_dig_montage(ch_pos=test_ch_pos, coord_frame='mni_tal')
+    with pytest.raises(RuntimeError, match='Montage should be in mri '
+                                           'coordinate frame'):
+        montage.add_estimated_fiducials(subject=subject,
+                                        subjects_dir=subjects_dir)
+
+
 run_tests_if_main()