Edit:
Now that you edited your question, I should add that the procedure of Tonry & Davis that I describe below is pretty general. In fact, appart from reducing the squared difference between the convoluted, shifted and scaled template and the observed spectrum (which can be thought of as a Maximum Likelihood estimator of the parameters with equal errors, which can be accounted for [see my comments]) the only part where they make strong assumptions is on the shape of the convolution of the template and on the shape of the largest peak in the CCF. You described that your objects may have different shapes; have you tried actually plotting the shape of the CCF? Maybe an extended gaussian or a Moffat function are more suitable choices in your case in order to model the 2D CCF peak. Analytically tracking the shape of the CCF, I think, is the easier path; this would allow you to analytically track the errors on the CCFs.
My original answer:
At least in astronomy we "have" a standard way of doing this, and it's the method of Tonry & Davis (1979). The part that may interest you in that paper is Section III, but let me explain to you a little of the background of their method in order to have an easier reading :-).
The idea in astronomy is that you have a source that emits certain flux (energy per meter$^2$ per second) at different wavelengths, say, $f_{\text{source}}(\lambda)$: this is called the spectrum of the source. However, because the objects that Tonry & Davis measure (galaxies, by the way) are receding from our point of view in our galaxy, all the known features of this object (say, bumps at certain wavelengths because of atomic absorption or emission) are generally shifted towards the red (i.e., we see them at longer wavelengths. For example, if we expected to see a bump at $\lambda=4500$, we may actually see it at $\lambda_{\text{source}}=4510$). Their work, then, is focused on calculating this wavelength shift and, of course, measuring the error on this shift.
Note that in their work, Tonry & Davis make a conversion between wavelength and bins (because we bin this flux as a function in wavelength in pixels; both to have higher signal to noise ratio and because a CCD camera is the best instrument to date for measuring flux) and from there they measure this wavelength shift. Maybe in your case you don't actually need this conversion between the bins Tonry & Davis talk about and wavelength, so you may want to change equation (1) in their paper to suit your needs. Another important feature of this paper is that they weight any deviations from zero quadratically, because these "bumps" I was talking about in the first paragraph of this answer are clearly more important. However, you can modify their $\chi^2$ reduction scheme to suit your needs ;-).
Finally, note that "all they do" is to approximate the largest peak in the Cross-Correlation Function (CCF) by a gaussian, and from there do the error analysis. In practice, I've seen this approximation to work pretty well, but always check the shape of your CCF, just in case.
PS: I didn't post this answer on that post you cited because JBWhitmore appeared to work on astronomy (or be an astronomer), and almost every astronomer I know knows the paper of Tonry & Davis, so I thought he was searching for something else.