tierneytim / OPM

Tools for OPM processing using SPM
GNU General Public License v2.0
17 stars 7 forks source link

OPM array simulation and source reconstruction #19

Open qishengjie9 opened 1 year ago

qishengjie9 commented 1 year ago

hello, I find the OPM sensor number is different from yours in your paper titled Pragmatic spatial sampling for wearable MEG arrays, such as 940 for 10mm space vs 784 in your paper. The codes I use are as follows. In addition, I got a problem when I try to use the sensor number as the spatial number in the module of source inversion iterative with the error information of not many this lead field.

for i = 1:Narray S =[]; S.space = space(i); S.meshres = 3; S.offset = 6.5; S.wholehead = 0; S.fname = sprintf('simopm%dmm',S.space); D = spm_opm_sim(S); newfile=fullfile(out_path, sprintf('simopm%dmm.mat',S.space)); save(newfile,'D'); end

tierneytim commented 1 year ago

Hi, the algorithm for sensor placement is not deterministic and has changed since the publication so the number of sensors will change when you run it. As for the iterative source inversion, I did not write that code and am not familiar with the issue you have described. It might be best to contact the author directly

qishengjie9 commented 1 year ago

Thanks for your kind help and advice. I will continue to try to solve the problem!

georgeoneill commented 1 year ago

Hello @qishengjie9 if you could copy in the code you use to perform the source inversion, and copy the error messages you get, we should be able to identify what the problem is.

qishengjie9 commented 1 year ago

Hello @georgeoneill I have copied the code and error in my email. I guess the problem is that the rank of the lead field matrix decomposed by svd is less than the sensor number when the number is 940 in my situation. You can check the code to give me inspiration and suggestions. Thanks for your kind help.

qishengjie9 commented 1 year ago

Hello,these are the snippet of my code where the error occurs and the error message. I find the error occured after performing the svd of the lead field matrix. image

9a5bae1d7f1349963d2bc63da6313e5
georgeoneill commented 1 year ago

I think this hsould be easy to fix.

You want to avoid the spatial decomposotion step entirely,

To do this, if n_channels is the number of channels of the dataset, set the code do be;

matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.nsmodes = n_channels;
matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.umodes= eye(n_channels);

the umodes option may need to be specified as a cell e.g. umodes= {eye(n_channels)};

qishengjie9 commented 1 year ago

The umodes I have uesd is the same as yours. And I found I can not avoid the spatial decomposition. The error message is the same as before. The spatial number available is 596 which is less than the sensor number of 940 in my simulation.

georgeoneill commented 1 year ago

Okay, can I think we need to check that everything communicated from the batch script is getting to the function.

1) can you show me the entire section of the matlabbatch you are using to call the function? 2) can you put a breakpoint at the line of the error. Run the code and it should pause at just before it throws an error. When it does, can you type inverse and report back what the command line returns?

qishengjie9 commented 1 year ago

Hello,I have done as you suggested. The matlabbatch code are as follows.

matlabbatch = []; matlabbatch{1}.spm.meeg.source.invertiter.D = {simfile}; matlabbatch{1}.spm.meeg.source.invertiter.val = 1; matlabbatch{1}.spm.meeg.source.invertiter.whatconditions.all = 1; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.invfunc = 'Classic'; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.invtype = algorithm; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.woi = [0 1000]; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.foi = [0 40]; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.hanning = 1; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.isfixedpatch.randpatch.npatches = 512; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.isfixedpatch.randpatch.niter = 1; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.patchfwhm = 0.6; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.mselect = 0; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.nsmodes = Nmodes; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.umodes = spatialmodesname; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.ntmodes = []; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.priors.priorsmask = {''}; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.priors.space = 1; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.restrict.locs = zeros(0, 3); matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.restrict.radius = 32; matlabbatch{1}.spm.meeg.source.invertiter.isstandard.custom.outinv = ''; matlabbatch{1}.spm.meeg.source.invertiter.modality = {'All'}; matlabbatch{1}.spm.meeg.source.invertiter.crossval = [0 1]; [a2,~]=spm_jobman('run', matlabbatch);

And the report is shown below after adding a breakpoint and type inverse. `Contains the following fields struct:

   priors: []
     type: 'EBB'
      woi: [0 999]
      Han: 1
      lpf: 0
      hpf: 40
mergeflag: 0
       Np: 512
       Nm: 940
       Nt: 4
   smooth: 0.6000
 modality: 'MEG'
     allF: 0
  PostMax: [20484×1 double]
        A: {}
       Ip: []

`

georgeoneill commented 1 year ago

It looks like the spatial modes are not being transferred to the function (it would be in inverse.A), the question I have is what is spatialmodesname? Is it a string? I beleive it should be either an array or an array ecapsulated in a cell. (e.g. umodes= {eye(n_channels)};) can you check please?

qishengjie9 commented 1 year ago

Thanks for your suggestions. And the spatialmodesnameis a string pointing to the spatial modes which contains the {eye(n_channels)}. The inverse.A is also empty after specifying umodes= {eye(n_channels)}. I think the inverse.A is dispensable if I choose all the sensor data(ie the number of spatial modes is equal to the number of sensors). Additionally, the inverse process will run successfully when the number of sensors is not that many in my simulation(eg about 200 sensors). I don't know if I understand correctly. Could you please give me some advice.