Frontier

I have no experience with line-fitting in 2-D and 3D spectra. I have never studied it and don't know what people have already done and which is the current approach. This is a blog and I want to use it like a blog, to express my ideas. I hope you'll also want to treat it like a blog, not as a textbook. If one day I am going to put my hands into this matter, my ideas will certainly evolve, as it always has been when passing from planning to real problem-solving.
nD spectra are more problematic than 1D spectra:

1. There are more parameters to estimate (the double in 2D and the triple in 3D).
2. There are less experimental points.
3. The shape is not Lorentzian, because the FID is multiplied by weighting functions.

A direct solution is to simulate the FID with exponentially damped sinusoids, apply the same weighting that has been applied to the spectrum, and compare the synthetic spectrum with the real counterpart. The fitting algorithms, alas, require the partial derivatives. I can't know their analytical formula (the weighting functions are not know at compile time, and even if they were...), but I can simply use the ratio Δf(x)/Δp, where f(x) is the spectrum and p the parameter. This direct solution requires more experimental points that is common to find. Requires quite a lot of programmatic work. It's purely computational, while I like interactive adjustments, where the user can manually change the parameters at any stage through a GUI.
Alternatively, I would start approximating. I'd assume, first, that the signal has not decayed too much or that the peaks have similar widths, therefore the shapes are common to all peaks (they would mostly depends on the weighting functions that have been applied, but it doesn't matter a lot). Second assumption: the common shape can be adequately approximated with a mixed Lorentzian/Gaussian function, namely:
f(ν,δ,I) = α L(ν,δ,I) + (1-α) G(ν,δ,I)
where 0 ≤ α ≤ 1, ν is the central frequency, δ the line-width at half height, I the area. Many peak shapes can be approximated with this formula. With this approach, The number of parameters decreases, because all the peaks would have the same δ and α along f₁ and another couple of shared values along f₂.
That could work, but: who needs it?
Here's a practical case: a ROESY with weak cross-peaks. They are dominated or polluted by DQF-like cross peaks that not only hinder the normal integration process, but are typically sharper than the genuine peaks, therefore my model can't work (not all the peaks have the same shape).