Translation Invariant Denoising with Cycle Spinning

Cycle spinning compensates for the lack of shift invariance in the critically-sampled wavelet transform by averaging over denoised cyclically-shifted versions of the signal or image. The appropriate inverse circulant shift operator is applied to the denoised signal/image and the results are averaged together to obtain the final denoised signal/image.

There are N unique cyclically-shifted versions of a signal of length, N. For an M-by-N image, there are MN versions. This makes using all possible shifted versions computationally prohibitive. However, in practice, good results can be obtained by using a small subset of the possible circular shifts.

1-D Cycle Spinning

This example shows how to denoise a 1-D signal using cycle spinning and the shift-variant orthogonal nonredundant wavelet transform. The example compares the results of the two denoising methods.

Create a noisy 1-D bumps signal with a signal-to-noise ratio of 6. The signal-to-noise ratio is defined as

N||X||22σ2

where N is the length of the signal, ||X||22 is the squared ℓ2 norm, and σ2 is the variance of the noise.

rng default;
[X,XN] = wnoise('bumps',10,sqrt(6));
subplot(211)
plot(X); title('Original Signal');
subplot(212)
plot(XN); title('Noisy Signal');

Denoise the signal using cycle spinning with 15 shifts, 7 to the left and 7 to the right, including the zero-shifted signal. Use Daubechies' least-asymmetric wavelet with 4 vanishing moments (sym4) and denoise the signal down to level 4 using soft thresholding and the universal threshold estimated from the level-1 detail coefficients.

ydenoise = zeros(length(XN),15);
for nn = -7:7
    yshift = circshift(XN,[0 nn]);
    [yd,cyd] = wden(yshift,'sqtwolog' ,'s','sln',4,'sym4');
    ydenoise(:,nn+8) = circshift(yd,[0, -nn]);
end
ydenoise = mean(ydenoise,2);

Denoise the signal using the orthogonal nonredundant discrete wavelet transform (DWT) with the same parameters. Compare the orthogonal DWT with cycle spinning.

xd = wden(XN,'sqtwolog','s','sln',4,'sym4');
subplot(211)
plot(ydenoise,'b','linewidth',2);
hold on;
plot(X,'r')
axis([1 1024 -10 10]);
legend('Denoised Signal','Original Signal','Location','SouthEast');
ylabel('Amplitude');
title('Cycle Spinning Denoising');
subplot(212)
plot(xd,'b','linewidth',2);
hold on;
plot(X,'r');
axis([1 1024 -10 10]);
legend('Denoised Signal','Original Signal','Location','SouthEast');
xlabel('Sample'); ylabel('Amplitude');
title('Standard Orthogonal Denoising');
absDiffDWT = norm(X-xd,2)
absDiffCycleSpin = norm(X-ydenoise',2)
absDiffDWT =

   18.0428

absDiffCycleSpin =

   15.4778

Cycle spinning with only 15 shifts has reduced the approximation error.

2-D Cycle Spinning

This example shows how to denoise an image using cycle spinning with 82=64 shifts.

Load the sine image and add zero-mean white Gaussian noise with a variance of 5.

load sinsin;
rng default;
Xnoisy = X+sqrt(5)*randn(size(X));
subplot(211)
imagesc(X); colormap(jet);
title('Original Image');
subplot(212)
imagesc(Xnoisy); title('Noisy Image');

Determine the universal threshold from the level-1 detail coefficients. Use the B-spline biorthogonal wavelet with 3 vanishing moments in the reconstruction wavelet and 5 vanishing moments in the decomposition wavelet.

wname = 'bior3.5';
[C,S] = wavedec2(Xnoisy,1,wname);
Cdet = C(4097:end);
THR = thselect(Cdet,'sqtwolog');

Create a grid of 8 shifts in both the X and Y directions. This results in a total of 64 shifts.

N = 8;
[deltaX, deltaY] = ndgrid(0:N-1,0:N-1);

Allocate a matrix of zeros the size of the image for the cycle spinning result. Specify soft thresholding and set the level to 3.

Xspin = zeros(size(X));
sorh = 's';
level = 3;

Use cycle spinning denoising and display the result.

for nn =1:N^2

    Xshift = circshift(Xnoisy, [deltaX(nn) deltaY(nn)]);
    [coefs,sizes] = wavedec2(Xshift,level,wname);
    [XDEN,cfsDEN,dimCFS] = wdencmp('gbl',coefs,sizes, ...
    wname,level,THR,sorh,1);
    XDEN = circshift(XDEN, -[deltaX(nn) deltaY(nn)]);
    Xspin = Xspin*(nn-1)/nn+XDEN/nn;

end

subplot(211)
imagesc(X); colormap(jet);
title('Original Image');
subplot(212)
imagesc(Xspin); title('Cycle Spinning');

Denoise the image using the identical parameters with the nonredundant DWT. Compare the peak signal-to-noise (PSNR), mean square error, and energy ratios obtained with cycle spinning and the nonredundant DWT.

[coefs,sizes] = wavedec2(Xnoisy,level,wname);
[XDEN,cfsDEN,dimCFS] = wdencmp('gbl',coefs,sizes,wname,3,THR,'s',1);
[PSNRcs,MSEcs,~,L2RATcs] = measerr(X,Xspin)
[PSNR,MSE,~,L2RAT] = measerr(X,XDEN)

The error measures show that cycle spinning has improved the image approximation. The PSNR, mean square error, and energy ratio are all better in the image denoised with cycle spinning.

Was this topic helpful?