At a quick glance:

• One of the first hurdles is the fact that you're basing the output modelled image on the input header, which could include b=0 volumes (or even other b-value shells technically), whereas the internal CSD structure will only include those volumes actually being utilised. So you'd need to either construct a new template header for that output image based on the number of volumes in the DW shell being utilised (and make sure dw_scheme is updated in that header also), or you'd need to zero-fill in the output image the volumes not utilised in CSD, and additionally know the mapping from rows in the CSD arithmetic to image volumes.

• Your forward convolution is using "HR_dirs", which I think are used for the non-negativity regularisation, and therefore don't reflect the input DWI data sensitisation directions.

Hi Tom! Hope all's well.

Just thought I'd mention before you invest too much time into this that the functionality you're after probably already exists with a combination of shconv & sh2amp, following the latest changes in #1617. This is available as of version 3.0.0. Probably not as convenient as an option in dwifod csd, but not all that tricky either...

@jdtournier: Confirm you agree with https://github.com/tclose/mrtrix3/blob/4716aacf1f53529ed89ad1089853ad7097b5aa04/cpp/cmd/dwi2fod.cpp#L162?