You have indeed found a bug. The effect is likely to be minimal though, as it
is usually called with values of epsilon and meth of 1,1. In the buggy code,
epsilon will therefore be 1. meth is a flag, with a value that is either
zero or non-zero. The bit of memory that it always accesses seems to always
have a value of zero. The difference between zero and one is that zero fixes
the boundaries, whereas a value of one leaves them to move freely.
If you fix this bug, then there are others (which I've just discovered) that
you should also fix in order to get freely moving boundaries. These are the
following:
if (meth & 1)
- tweek(0, x1, x2, y0, y1, y2,
dim_g, g, dim_f, vox_f, f, 1|msk2,
+ tweek(1, x1, x2, y0, y1, y2,
dim_g, g, dim_f, vox_f, f, 1|msk2,
*lambda, *epsilon,
*sigma2, *scale, &hf, &hp, &n, &sumf, &sumg, &cnt,sgn);
for(x0=2; x0<dim_g[0]; x0++)
tweek(x0, x1, x2, y0, y1, y2,
dim_g, g, dim_f, vox_f, f, msk2,
*lambda, *epsilon,
*sigma2, *scale, &hf, &hp, &n, &sumf, &sumg, &cnt,sgn);
if (meth & 1)
- tweek(dim_g[0]-1, x1, x2, y0,
y1, y2, dim_g, g, dim_f, vox_f, f, 2|msk2,
+ tweek(dim_g[0], x1, x2, y0,
y1, y2, dim_g, g, dim_f, vox_f, f, 2|msk2,
*lambda, *epsilon,
*sigma2, *scale, &hf, &hp, &n, &sumf, &sumg, &cnt,sgn);
if (meth & 1)
- tweek(dim_g[0]-1, x1, x2, y0,
y1, y2, dim_g, g, dim_f, vox_f, f, 2|msk2,
+ tweek(dim_g[0], x1, x2, y0,
y1, y2, dim_g, g, dim_f, vox_f, f, 2|msk2,
*lambda, *epsilon,
*sigma2, *scale, &hf, &hp, &n, &sumf, &sumg, &cnt,sgn);
for(x0=dim_g[0]-1; x0>1; x0--)
tweek(x0, x1, x2, y0, y1, y2,
dim_g, g, dim_f, vox_f, f, msk2,
*lambda, *epsilon,
*sigma2, *scale, &hf, &hp, &n, &sumf, &sumg, &cnt,sgn);
if (meth & 1)
- tweek(0, x1, x2, y0, y1, y2,
dim_g, g, dim_f, vox_f, f, 1|msk2,
+ tweek(1, x1, x2, y0, y1, y2,
dim_g, g, dim_f, vox_f, f, 1|msk2,
*lambda, *epsilon,
*sigma2, *scale, &hf, &hp, &n, &sumf, &sumg, &cnt,sgn);
A fixed version of the code will be in the next set of updates.
Best regards,
-John
On Thursday 03 August 2006 20:44, Matthew L. Senjem wrote:
> Hi John (and list),
>
> I am trying to understand the code in spm_warp.c, and I think I may have
> stumbled upon a bug:
>
> in file: $spm2/toolbox/Deformations/spm_warp.c
>
> At lines 798 - 801, it says:
> its = mxGetPr(prhs[6])[0];
> lambda = mxGetPr(prhs[6])[1];
> epsilon = mxGetPr(prhs[6])[3];
> meth = mxGetPr(prhs[6])[4];
>
> but I believe it should read:
> its = mxGetPr(prhs[6])[0];
> lambda = mxGetPr(prhs[6])[1];
> epsilon = mxGetPr(prhs[6])[2];
> meth = mxGetPr(prhs[6])[3];
>
> because in this case, prhs[6] has only 4 elements, and in C the indexing
> should go 0,1,2,3, correct?
>
>
> I believe the same bug is also in $spm5/toolbox/HDW/spm_warp.c
|