|
You only need the one shear, so that you have 6 parameters in total.
You can test this by looking for affine transforms that can not be
represented by those six parameters:
for i=1:10000
M0=[randn(2,3)*10; 0 0 1];
P1=spm_imatrix2D(M0);
M1=spm_matrix2D(P1);
t=sum(sum((M1-M0).^2));
if t>1e-9, disp('Problem'); break; end
end
Best regards,
-John
On 6 May 2011 21:48, Siddharth Srivastava <[log in to unmask]> wrote:
> Hi John, please ignore the first part of my previous mail: made some errors
> in
> cut and paste. It runs fine now. The question about the shears, however,
> still lingers: do I have
> to somehow use both Sx and Sy, or as your implementation suggests, only one
> would suffice?
>
> thanks again, and sorry for being such a bother...
> sid.
>
> On Fri, May 6, 2011 at 12:23 PM, Siddharth Srivastava <[log in to unmask]>
> wrote:
>>
>> Thanks John, for the answers. The code, however exits with "Problem". M0
>> and M1 differ.
>> Also, why is there just one shear component in the matrix? I have seen
>> general shears modeled as [1,h;g,1]: is this superfluous, and just one "h"
>> will be sufficient for 2d case?
>>
>> Thanks,
>> Sid.
>>
>>
>>
>> On Fri, May 6, 2011 at 11:01 AM, John Ashburner <[log in to unmask]>
>> wrote:
>>>
>>> Your model uses 7 parameters, but there are only 6 elements in
>>> A(1:2,1:3). Also, you need to be able to deal with a broader range of
>>> angles, so atan2(s,c) is used to compute the angle. I've tested it
>>> with the following code, so hopefully it works for all cases.
>>>
>>> for i=1:10000
>>> P=randn(1,6)*10;
>>> M0=spm_matrix2D(P);
>>> P1=spm_imatrix2D(M0);
>>> M1=spm_matrix2D(P1);
>>> t=sum(sum((M1-M0).^2));
>>> if t>1e-9, disp('Problem'); break; end
>>> end
>>>
>>> Note that spm_imatrix2D(spm_matrix2D(P)) does not necessarily have to
>>> equal P, but the matrices generated from the two parameter sets should
>>> be the same.
>>>
>>> Best regards,
>>> -John
>>>
>>> On 6 May 2011 17:11, Siddharth Srivastava <[log in to unmask]> wrote:
>>> > Hi everyone,
>>> > Can someone help me check the 2D versions of the spm_matrix and
>>> > spm_imatrix? My
>>> > implementation currently is as follows:
>>> >
>>> > ---2D spm_matrix ---
>>> > function [A] = spm_matrix2D(P)
>>> > % P = [Tx, Ty, R, Zx, Zy, Sx,Sy]
>>> > p0 = [0,0,0,1,1,0,0];
>>> > P = [P, p0((length(P)+1):7)]
>>> >
>>> > A = eye(3);
>>> > % T -> R -> Zoom -> Shear
>>> > A = A * [1,0,P(1); ...
>>> > 0,1,P(2); ...
>>> > 0,0, 1];
>>> > A = A * [cos(P(3)), sin(P(3)), 0; ...
>>> > -sin(P(3)), cos(P(3)), 0; ...
>>> > 0, 0, 1];
>>> > A = A * [P(4), 0, 0; ...
>>> > 0, P(5), 0; ...
>>> > 0, 0, 1];
>>> > A = A * [1, P(6), 0; ...
>>> > P(7) , 1, 0; ...
>>> > 0, 0, 1];
>>> > --------------------------------------------------
>>> > -----2D spm_imatrix ------------------------
>>> > function P = spm_imatrix2D(M);
>>> >
>>> >
>>> >
>>> > R = M(1:2,1:2);
>>> > C = chol(R' * R);
>>> > P = [M(1:2,3)',0,diag(C)',0,0];
>>> > if(det(R) < 0) P(4) = -P(4); end
>>> > C = diag(diag(C))\C;
>>> > P(6:7) = C([2,3]);
>>> >
>>> > R0 = spm_matrix2D([0,0,0,P(4:end)]);
>>> > R0 = R0(1:2,1:2);
>>> > R1 = R/R0;
>>> >
>>> > th = asin(rang(R1(1,2)));
>>> > P(3) = th;
>>> >
>>> > % There may be slight rounding errors making b>1 or b<-1.
>>> > function a = rang(b)
>>> > a = min(max(b, -1), 1);
>>> > return;
>>> >
>>> > ---------------------------------------------------------------------------
>>> >
>>> > I know there is something not correct, since when, for a parameter set
>>> > 'p',
>>> > i do spm_imatrix2D(spm_matrix2D(p)), I do not recover P for all
>>> > entries,
>>> > specially
>>> > rotation and shears. C(2) above always evaluates to zero. The result is
>>> > that
>>> > these errors
>>> > accumulate (or are not accounted for) on repeated application of this
>>> > sequence of
>>> > operations in the registration routine.
>>> >
>>> > Thanks in anticipation,
>>> > Sid.
>>> >
>>> >
>>> >
>>> > On Fri, Apr 22, 2011 at 9:13 AM, Siddharth Srivastava
>>> > <[log in to unmask]>
>>> > wrote:
>>> >>
>>> >> Hi John,
>>> >> Thanks a lot for your answer, and it has been very helpful for
>>> >> my
>>> >> understanding
>>> >> of the details of the implementation, and my own development efforts.
>>> >> I
>>> >> have
>>> >> 2 very quick questions, though:
>>> >>
>>> >> 1) the parameter set x(:) that is returned from spm_mireg, and the
>>> >> "params" that
>>> >> spm_affsub3 returns (in spm99, which is also the version that I am
>>> >> looking
>>> >> at) ,
>>> >> are they equivalent? In other words, If I am collating the parameters
>>> >> from
>>> >> the mireg
>>> >> routine, would the distribution x(7:12) be used to estimate the
>>> >> icovar0
>>> >> entries in affsub3?
>>> >>
>>> >> 2) also would VG.mat\spm_matrix(x)*VF.mat transform G to F similar to
>>> >> VG.mat\spm_matrix(params)*VF.mat , or do I have to do
>>> >> M = VG.mat\spm_matrix(x(:)')*VF.mat
>>> >> pp = spm_imatrix(M)
>>> >> and use pp(7:12) ?
>>> >>
>>> >>
>>> >>
>>> >> Thanks again.
>>> >> Sid.
>>> >>
>>> >> On Thu, Apr 14, 2011 at 2:13 PM, John Ashburner <[log in to unmask]>
>>> >> wrote:
>>> >>>
>>> >>> > I had a look at the code in spm_maff, and
>>> >>> > it looks like it also uses the values of isig and mu in the call to
>>> >>> > affreg(). My guess is
>>> >>> > that during the initial registration on a large data set for
>>> >>> > estimating
>>> >>> > the
>>> >>> > prior
>>> >>> > distribution of parameters, this will be called with typ = 'none' .
>>> >>> > Is
>>> >>> > this
>>> >>> > correct?
>>> >>>
>>> >>> Once the registration has locked in on a solution, the priors have
>>> >>> less influence. Their main influence is in the early stages of the
>>> >>> registration, or when there are relatively few slices in the data. I
>>> >>> don't actually remember if the registration was regularised when I
>>> >>> computed their values.
>>> >>>
>>> >>> >
>>> >>> > Further, do we have to do a non-rigid registration for the
>>> >>> > estimation
>>> >>> > (since, in the comments the selected
>>> >>> > files are filtered as *seg_inv_sn.mat"? Since only p.Affine is
>>> >>> > being
>>> >>> > used,
>>> >>> > can this work with only
>>> >>> > an affine / rigid transformation?
>>> >>>
>>> >>> Only the affine transform in the files is used, so basically yes.
>>> >>> You
>>> >>> could just do affine registration.
>>> >>>
>>> >>> >
>>> >>> > I would also be happy if you could also answer my other question
>>> >>> > about
>>> >>> > examining the
>>> >>> > effects of the values in isig and mu (what are the units of mu?) on
>>> >>> > the
>>> >>> > transformed image.
>>> >>>
>>> >>> They are unit-less, and are just a measure of the zooms and shears.
>>> >>>
>>> >>> > For example,
>>> >>> > if i set a mu of [mu1, mu2, ...] and a isig of [sig1, 0, 0...; 0,
>>> >>> > sig2,
>>> >>> > 0,
>>> >>> > 0..; ...] (with only diagonal components),
>>> >>> > based on the Hencky tensor penalty, it seems that large deviations
>>> >>> > of
>>> >>> > parameter 1 estimate (T - mu) from
>>> >>> > mu1 will be penalized. How is the value of isig1 modifying the
>>> >>> > penalty?
>>> >>>
>>> >>> Its assuming that the elements of the matrix log of the part that
>>> >>> does
>>> >>> shears and zooms has a multivariate normal distribution, so the
>>> >>> penalty is a kind of negative log likelhood. The isig is essentially
>>> >>> a precision matrix, which is the inverse of a covariance matrix.
>>> >>>
>>> >>> p(x|mu,S) = 1/sqrt(det(2*pi*S)) * exp(-0.5*(x-mu)'*inv(S)*(x-mu))
>>> >>>
>>> >>> Taking logs and negating, you get 0.5*(x-mu)'*inv(S)*(x-mu) + const
>>> >>>
>>> >>> > How
>>> >>> > do the off diagonal terms affect
>>> >>> > the penalty?
>>> >>>
>>> >>> The inverse of isig would just encode the covariance between the
>>> >>> various shears and zooms.
>>> >>>
>>> >>> >
>>> >>> > Could you also point me to a reference where I can learn more about
>>> >>> > Hencky
>>> >>> > strain tensor (specifically
>>> >>> > what information it carries in the context of image processing
>>> >>> > etc.)
>>> >>>
>>> >>> I don't know if this is of any use:
>>> >>>
>>> >>> http://en.wikipedia.org/wiki/Deformation_(mechanics)
>>> >>>
>>> >>> Best regards,
>>> >>> -John
>>> >>>
>>> >>>
>>> >>> > On Thu, Apr 14, 2011 at 10:17 AM, John Ashburner
>>> >>> > <[log in to unmask]>
>>> >>> > wrote:
>>> >>> >>
>>> >>> >> They were computed by affine registering loads of (227) images to
>>> >>> >> MNI
>>> >>> >> space. This gave the affine transforms, which were then used to
>>> >>> >> build
>>> >>> >> up a mean and inverse covariance matrix. In the comments of
>>> >>> >> spm_maff.m, you should see the following...
>>> >>> >>
>>> >>> >> %% Values can be derived by...
>>> >>> >> %sn = spm_select(Inf,'.*seg_inv_sn.mat$');
>>> >>> >> %X = zeros(size(sn,1),12);
>>> >>> >> %for i=1:size(sn,1),
>>> >>> >> % p = load(deblank(sn(i,:)));
>>> >>> >> % M = p.VF(1).mat*p.Affine/p.VG(1).mat;
>>> >>> >> % J = M(1:3,1:3);
>>> >>> >> % V = sqrtm(J*J');
>>> >>> >> % R = V\J;
>>> >>> >> % lV = logm(V);
>>> >>> >> % lR = -logm(R);
>>> >>> >> % P = zeros(12,1);
>>> >>> >> % P(1:3) = M(1:3,4);
>>> >>> >> % P(4:6) = lR([2 3 6]);
>>> >>> >> % P(7:12) = lV([1 2 3 5 6 9]);
>>> >>> >> % X(i,:) = P';
>>> >>> >> %end;
>>> >>> >> %mu = mean(X(:,7:12));
>>> >>> >> %XR = X(:,7:12) - repmat(mu,[size(X,1),1]);
>>> >>> >> %isig = inv(XR'*XR/(size(X,1)-1))
>>> >>> >>
>>> >>> >> You may be able to re-code this to work with your 2D transforms.
>>> >>> >> Note
>>> >>> >> that I only use the components that describe shears and zooms.
>>> >>> >> Also,
>>> >>> >> if I was re-coding it all again, I would probably do things
>>> >>> >> slightly
>>> >>> >> differently, using the following decomposition:
>>> >>> >>
>>> >>> >> J=M(1:3,1:3);
>>> >>> >> L=real(logm(J));
>>> >>> >> S=(L+L')/2; % Symmetric part (for zooms and shears, which should
>>> >>> >> be
>>> >>> >> penalised)
>>> >>> >> A=(L-L')/2; % Anti-symmetric part (for rotations)
>>> >>> >>
>>> >>> >> I might even base the penalty on lambda(2)*trace(S)^2 +
>>> >>> >> lambda(1)*trace(S*S'), which is often used as a measure of
>>> >>> >> "elastic
>>> >>> >> energy" for penalising non-linear registration. Of course, there
>>> >>> >> would also be a bit of annoying extra stuff to deal with due to
>>> >>> >> MNI
>>> >>> >> space being a bit bigger than average brains are. Figuring out
>>> >>> >> the
>>> >>> >> mean would require (I think) an iterative process similar to the
>>> >>> >> one
>>> >>> >> described in "Characterizing volume and surface deformations in an
>>> >>> >> atlas framework: theory, applications, and implementation", by
>>> >>> >> Woods
>>> >>> >> in NeuroImage (2003).
>>> >>> >>
>>> >>> >> Best regards,
>>> >>> >> -John
>>> >>> >>
>>> >>> >> On 14 April 2011 17:15, Siddharth Srivastava <[log in to unmask]>
>>> >>> >> wrote:
>>> >>> >> > Hi all,
>>> >>> >> > My question pertains specifically to the values of isig and
>>> >>> >> > mu
>>> >>> >> > that
>>> >>> >> > are
>>> >>> >> > incorporated
>>> >>> >> > in the code, and I wish to disentangle the effect that these
>>> >>> >> > numbers
>>> >>> >> > have on
>>> >>> >> > individual parameters
>>> >>> >> > during the optimization stage. I am trying to map it to a 2D
>>> >>> >> > case,
>>> >>> >> > and
>>> >>> >> > these
>>> >>> >> > numbers
>>> >>> >> > will have to be calculated ab-initio for my case. Before I begin
>>> >>> >> > with
>>> >>> >> > estimating the prior distribution
>>> >>> >> > for these parameters, I need to check the effect that they have
>>> >>> >> > on
>>> >>> >> > the
>>> >>> >> > registered images. Regarding this
>>> >>> >> > 1) Can someone suggest a way I can examine the effect on each
>>> >>> >> > parameter
>>> >>> >> > separately, if it is
>>> >>> >> > at all possible.
>>> >>> >> > 2) Some advise on how to estimate these parameters, if it has
>>> >>> >> > been
>>> >>> >> > done
>>> >>> >> > for
>>> >>> >> > images other
>>> >>> >> > than brain images, will be really helpful for me.
>>> >>> >> >
>>> >>> >> > Thanks,
>>> >>> >> > Sid.
>>> >>> >> >
>>> >>> >> > On Mon, Mar 28, 2011 at 7:04 PM, Siddharth Srivastava
>>> >>> >> > <[log in to unmask]>
>>> >>> >> > wrote:
>>> >>> >> >>
>>> >>> >> >> Dear John,
>>> >>> >> >> Thanks for your answer. I was able to correctly run
>>> >>> >> >> my 2D
>>> >>> >> >> implementation
>>> >>> >> >> based on your advise. However, I am am not very confident of my
>>> >>> >> >> results
>>> >>> >> >> and wanted to
>>> >>> >> >> consult you on certain aspects of my mapping from 3D to 2D.
>>> >>> >> >> 1) In function reg(..), the comment says that the
>>> >>> >> >> derivatives
>>> >>> >> >> w.r.t rotations and
>>> >>> >> >> translations is zero. Consequently, I ran the indices from 4:7.
>>> >>> >> >> According
>>> >>> >> >> to your previous
>>> >>> >> >> reply, rotations will be lumped together with scaling and
>>> >>> >> >> affine in
>>> >>> >> >> 4
>>> >>> >> >> parameters, and
>>> >>> >> >> translations, independently, in the last two. I understand that
>>> >>> >> >> in
>>> >>> >> >> this
>>> >>> >> >> case, the indices
>>> >>> >> >> run over parameters that have been separated into individual
>>> >>> >> >> components
>>> >>> >> >> (similar to
>>> >>> >> >> the form accepted by spm_matrix). Is this correct for the
>>> >>> >> >> function
>>> >>> >> >> reg(..)? In fact, I went back to the
>>> >>> >> >> old spm code (spm99), and found parameters pdesc and free, and
>>> >>> >> >> got
>>> >>> >> >> more
>>> >>> >> >> confused as to when I should
>>> >>> >> >> use the 6 parameter, and when to use separate parameters (as if
>>> >>> >> >> I
>>> >>> >> >> was
>>> >>> >> >> passing them to spm_matrix).
>>> >>> >> >> So, when the loop runs over p1 and perturbs p1(i) by ds, does
>>> >>> >> >> it
>>> >>> >> >> run
>>> >>> >> >> over
>>> >>> >> >> translations/rotations etc
>>> >>> >> >> or does it run over the martix elements that define these
>>> >>> >> >> transformations
>>> >>> >> >> (2x2 + 1x2)?
>>> >>> >> >>
>>> >>> >> >> 2) In function penalty(..), I have used unique
>>> >>> >> >> elements
>>> >>> >> >> as
>>> >>> >> >> [1,3,4], and my code looks like
>>> >>> >> >> T = my_matrix(p)*M;
>>> >>> >> >> T = T(1:2,1:2);
>>> >>> >> >> T = 0.5*logm(T'*T);
>>> >>> >> >> T
>>> >>> >> >> T = T(els)' - mu;
>>> >>> >> >> h = T'*isig*T
>>> >>> >> >>
>>> >>> >> >> Assuming an application specific isig, is this correct?
>>> >>> >> >> 3) I have not yet incorporated the prior information
>>> >>> >> >> in my
>>> >>> >> >> code,
>>> >>> >> >> But for testing purposes,
>>> >>> >> >> can you suggest some neutral values for mu and isig that I can
>>> >>> >> >> try
>>> >>> >> >> to
>>> >>> >> >> concentrate on the
>>> >>> >> >> accuracy of the implementation minus the priors for now?
>>> >>> >> >> 4) My solution curently seems to oscillate around an
>>> >>> >> >> optimum.
>>> >>> >> >> I
>>> >>> >> >> am hoping that after
>>> >>> >> >> I have removed any errors after your replies, it will converge
>>> >>> >> >> gracefully.
>>> >>> >> >>
>>> >>> >> >> Thanks again for all the help.
>>> >>> >> >> Siddharth.
>>> >>> >> >>
>>> >>> >> >>
>>> >>> >> >> On Fri, Mar 18, 2011 at 2:05 AM, John Ashburner
>>> >>> >> >> <[log in to unmask]>
>>> >>> >> >> wrote:
>>> >>> >> >>>
>>> >>> >> >>> A 2D affine transform would use 6 parameters, rather than 7: a
>>> >>> >> >>> 2x2
>>> >>> >> >>> matrix to encode rotations, zooms and shears, and a 2x1 vector
>>> >>> >> >>> for
>>> >>> >> >>> translations.
>>> >>> >> >>>
>>> >>> >> >>> Best regards,
>>> >>> >> >>> -John
>>> >>> >> >>>
>>> >>> >> >>> On 18 March 2011 05:05, Siddharth Srivastava
>>> >>> >> >>> <[log in to unmask]>
>>> >>> >> >>> wrote:
>>> >>> >> >>> > Dear list users,
>>> >>> >> >>> > I am currently trying to understand the
>>> >>> >> >>> > implementation
>>> >>> >> >>> > in
>>> >>> >> >>> > spm_affreg, and
>>> >>> >> >>> > to make matters simple, I tried to work out a 2D version of
>>> >>> >> >>> > the
>>> >>> >> >>> > registration
>>> >>> >> >>> > procedure. I struck a roadblock, however...
>>> >>> >> >>> >
>>> >>> >> >>> > The function make_A creates part of the design matrix. For
>>> >>> >> >>> > the
>>> >>> >> >>> > 3D
>>> >>> >> >>> > case,
>>> >>> >> >>> > the
>>> >>> >> >>> > 13 parameters
>>> >>> >> >>> > are encoded in columns of A1, which then is used to generate
>>> >>> >> >>> > AA
>>> >>> >> >>> > +
>>> >>> >> >>> > A1'
>>> >>> >> >>> > A1. AA
>>> >>> >> >>> > is 13x13, and everything
>>> >>> >> >>> > works. For the 2D case, however, there will be 7+1
>>> >>> >> >>> > parameters( 2
>>> >>> >> >>> > translations, 1 rotation, 2 zoom, 2 shears and 1 intensity
>>> >>> >> >>> > scaling).
>>> >>> >> >>> > The A1 matrix for 2D, then, will be (in make_A)
>>> >>> >> >>> > A = [x1.*dG1 x1.*dG2 ...
>>> >>> >> >>> > x2.*dG1 x2.*dG2 ...
>>> >>> >> >>> > dG1 dG2 t];
>>> >>> >> >>> >
>>> >>> >> >>> > which is (number of sampled points) x 7. The AA matrix (in
>>> >>> >> >>> > function
>>> >>> >> >>> > costfun)
>>> >>> >> >>> > is 8x8
>>> >>> >> >>> > so, which is the parameter I am missing in modeling a 2D
>>> >>> >> >>> > example? I
>>> >>> >> >>> > hope I
>>> >>> >> >>> > have got the number of parameters
>>> >>> >> >>> > right for the 2D case? Is there an extra column that I have
>>> >>> >> >>> > to
>>> >>> >> >>> > add
>>> >>> >> >>> > to
>>> >>> >> >>> > A1 to
>>> >>> >> >>> > make the rest of the matrix computation
>>> >>> >> >>> > consistent?
>>> >>> >> >>> >
>>> >>> >> >>> > Thanks for the help
>>> >>> >> >>> >
>>> >>> >> >>
>>> >>> >> >
>>> >>> >> >
>>> >>> >
>>> >>> >
>>> >>
>>> >
>>> >
>>
>
>
|
