Added computation of new frequency-domain and time-domain Jacobians

[HirviP]

MCX originally computes the Jacobians for CW and TD intensity with respect to voxel-wise absorption and scattering coefficients, and for radio-frequency (RF/FD) intensity-modulated light with respect to absorption. We added the computation of the RF scattering Jacobians, and the Jacobians for photon mean time-of-flight (TOF) with respect to absorption and scattering. Specifically, we added three new output types: 1) cosine- and sine weighted average scattering counts in each voxel (similarly to current rf replay for absorption), 2) weighted average of "total TOF x total path length" in each voxel, and 3) weighted average of "total TOF x scattering count" in each voxel. These can be combined with the other replay modes to obtain the Jacobians for FD amplitude and phase shift with respect to scattering, and TD mean TOF with respect to absorption and scattering. See the following document for the required formulas: Jacobian_Formulas. We also removed the normalization with cfg.unitinmm from replay outputs that do not use path lengths in voxel grid units.

Check List

Before you submit your pull-request, please verify and check all below items

• You have only modified the minimum number of lines that are necessary for the update; you should never add/remove white spaces to other lines that are irrelevant to the desired change.
• You have run make pretty (requires astyle in the command line) under the src/ folder and formatted your C/C++/CUDA source codes before every commit; similarly, you should run python3 -m black *.py (pip install black first) to reformat all modified Python codes, or run mh_style --fix . (pip install miss-hit first) at the top-folder to format all MATLAB scripts.
• Add sufficient in-code comments following the doxygen C format
• In addition to source code changes, you should also update the documentation (README.md, mcx_utils.c and/or mcxlab.m) if any command line flag was added or changed.

If your commits included in this PR contain changes that did not follow the above guidelines, you are strongly recommended to create a clean patch using git rebase and git cherry-pick to prevent in-compliant history from appearing in the upstream code.

Moreover, you are highly recommended to

• Add a test in the mcx/test/testmcx.sh script, following existing examples, to test the newly added feature; or add a MATLAB script under mcxlab/examples to gives examples of the desired outputs
• MCX's simulation speed is currently limited by the number of GPU registers. In your change, please consider minimizing the use of registers by reusing existing ones. CUDA compilers may not be able to optimize register counts, thus require manual optimization.
• Please create a github Issue first with detailed descriptions of the problem, testing script and proposed fixes, and link the issue with this pull-request

[fangq]

@HirviP, so sorry for the long delay in getting back to you on this - I have been overwhelmed by a progress report for another parallel project.

I appreciate very much for the rigorous implementation and careful formatting the code.

before accepting this PR, one of my questions is related to the Test_mcxlab_rfmusreplay_PH.sh script included in this PR

do you need it to run the example? if not, I suggest to remove it.

[rantahar]

Build seems to fail because of the minimum required CMake version. I will try updating it.

[rantahar]

I think Test_mcxlab_rfmusreplay_PH.sh is specific to the HPC system at Aalto University. It's probably better to remove it.

[HirviP]

Dear @fangq ,

No worries at all, great if you have some time for this project now! I agree that the .sh script was extra, so thanks @rantahar for removing it. I also just added demo script 'mcxlab/examples/demo_replay_all_jacobian.m' for computing any of the original or new Jacobian types (CW/TD/FD and mua/mus). This could also be formatted as a function with input parameter defining the wanted Jacobian type.

I can make a separate GitHub issue on this matter, but note that the original 'wp' or 'nscat' replay can return NaN, very large values and negative values close to the source. One can observe this from matrix 'meanP' in the new demo script (or 'dlogIdmus' and 'dTOFdmus' Jacobians that utilize 'meanP'). I'm looking into this, and we are also adding validation results to the Overleaf document.