function W2 = transition_weights(Twm,Trans,Cwm,W)

numM = size(Twm,2);

[i,j] = find(W);
ind = sub2ind(size(W),i,j);

Tij = zeros(3,numM,numM);
Tij(:,ind) = diff_config(Twm(:,i),Twm(:,j));
%Tij(3,ind) = mod(Tij(3,ind)+pi,2*pi)-pi;
%Trans(3,ind) = mod(Trans(3,ind)+pi,2*pi)-pi;

A = mul_config(Twm(:,i),Trans(:,ind));
B = mul_config(Twm(:,i),Tij(:,ind));
Res = A-B;
Res(3,:) = mod(Res(3,:)+pi,2*pi)-pi;

Sigma = Cwm(:,:,i)-Cwm(:,:,j);
Sigma = Sigma * (10)^2;
neg = find(Sigma(1,1,:)<0);
Sigma(:,:,neg) = -Sigma(:,:,neg);
[a,b,c] = ndgrid(1:3,1:3,0:length(ind)-1);
SigmaS = sparse(a(:)+3*c(:),b(:)+3*c(:),Sigma(:));
iS = inv(SigmaS);

iSRes = reshape(full(iS*Res(:)),size(Res));
Nres2 = abs(sum(Res.*iSRes,1));

W2 = zeros(numM);
W2(ind) = exp(-Nres2/2)+0.0;

%make weights symmetric
W2 = .5*(W2+W2');

if 0
%invert all 3x3 matrices
det = ...
   -Sigma(1,3,:).*Sigma(2,2,:).*Sigma(3,1) ...
   +Sigma(1,2,:).*Sigma(2,3,:).*Sigma(3,1) ...
   +Sigma(1,3,:).*Sigma(2,1,:).*Sigma(3,2) ...
   -Sigma(1,1,:).*Sigma(2,3,:).*Sigma(3,2) ...
   -Sigma(1,2,:).*Sigma(2,1,:).*Sigma(3,3) ...
   +Sigma(1,1,:).*Sigma(2,2,:).*Sigma(3,3);

iS = zeros(size(Sigma));
iS(1,1,:) = -Sigma(2,3,:).*Sigma(3,2,:) + Sigma(2,2,:).*Sigma(3,3,:);
iS(1,2,:) = -Sigma(1,3,:).*Sigma(3,2,:) - Sigma(1,2,:).*Sigma(3,3,:);
iS(1,3,:) = -Sigma(1,3,:).*Sigma(2,2,:) + Sigma(1,2,:).*Sigma(2,3,:);
iS(2,1,:) = -Sigma(2,3,:).*Sigma(3,1,:) - Sigma(2,1,:).*Sigma(3,3,:);
iS(2,2,:) = -Sigma(1,3,:).*Sigma(3,1,:) + Sigma(1,1,:).*Sigma(3,3,:);
iS(2,3,:) = -Sigma(1,3,:).*Sigma(2,1,:) - Sigma(1,1,:).*Sigma(2,3,:);
iS(3,1,:) = -Sigma(2,2,:).*Sigma(3,1,:) + Sigma(2,1,:).*Sigma(3,2,:);
iS(3,2,:) = -Sigma(1,2,:).*Sigma(3,1,:) - Sigma(1,1,:).*Sigma(3,2,:);
iS(3,3,:) = -Sigma(1,2,:).*Sigma(2,1,:) + Sigma(1,1,:).*Sigma(2,2,:);
iS(:,:,:) = iS(:,:,:)./det;

end
   
