function x = ifdct_wrapping(C, is_real, M, N)
% ifdct_wrapping.m - Inverse Fast Discrete Curvelet Transform via wedge wrapping - Version 1.0
% This is in fact the adjoint, also the pseudo-inverse
%
% Inputs
% C Cell array containing curvelet coefficients (see
% description in fdct_wrapping.m)
% is_real As used in fdct_wrapping.m
% M, N Size of the image to be recovered (not necessary if finest
% = 2)
%
% Outputs
% x M-by-N matrix
%
% See also fdct_wrapping.m
%
% By Laurent Demanet, 2004
% Initialization
nbscales = length(C);
nbangles_coarse = length(C{2});
nbangles = [1, nbangles_coarse .* 2.^(ceil((nbscales-(nbscales:-1:2))/2))];
if length(C{end}) == 1, finest = 2; else finest = 1; end;
if finest == 2, nbangles(nbscales) = 1; end;
if nargin < 2, is_real = 0; end;
if nargin < 4,
if finest == 1, error('Syntax: IFCT_wrapping(C,M,N) where the matrix to be recovered is M-by-N'); end;
[N1,N2] = size(C{end}{1});
else
N1 = M;
N2 = N;
end;
M1 = N1/3;
M2 = N2/3;
if finest == 1;
bigN1 = 2*floor(2*M1)+1;
bigN2 = 2*floor(2*M2)+1;
X = zeros(bigN1,bigN2);
% Initialization: preparing the lowpass filter at finest scale
window_length_1 = floor(2*M1) - floor(M1) - 1 - (mod(N1,3)==0);
window_length_2 = floor(2*M2) - floor(M2) - 1 - (mod(N2,3)==0);
coord_1 = 0:(1/window_length_1):1;
coord_2 = 0:(1/window_length_2):1;
[wl_1,wr_1] = fdct_wrapping_window(coord_1);
[wl_2,wr_2] = fdct_wrapping_window(coord_2);
lowpass_1 = [wl_1, ones(1,2*floor(M1)+1), wr_1];
if mod(N1,3)==0, lowpass_1 = [0, lowpass_1, 0]; end;
lowpass_2 = [wl_2, ones(1,2*floor(M2)+1), wr_2];
if mod(N2,3)==0, lowpass_2 = [0, lowpass_2, 0]; end;
lowpass = lowpass_1'*lowpass_2;
scales = nbscales:-1:2;
else
M1 = M1/2;
M2 = M2/2;
bigN1 = 2*floor(2*M1)+1;
bigN2 = 2*floor(2*M2)+1;
X = zeros(bigN1,bigN2);
window_length_1 = floor(2*M1) - floor(M1) - 1;
window_length_2 = floor(2*M2) - floor(M2) - 1;
coord_1 = 0:(1/window_length_1):1;
coord_2 = 0:(1/window_length_2):1;
[wl_1,wr_1] = fdct_wrapping_window(coord_1);
[wl_2,wr_2] = fdct_wrapping_window(coord_2);
lowpass_1 = [wl_1, ones(1,2*floor(M1)+1), wr_1];
lowpass_2 = [wl_2, ones(1,2*floor(M2)+1), wr_2];
lowpass = lowpass_1'*lowpass_2;
hipass_finest = sqrt(1 - lowpass.^2);
scales = (nbscales-1):-1:2;
end;
% Loop: pyramidal reconstruction
Xj_topleft_1 = 1;
Xj_topleft_2 = 1;
for j = scales,
M1 = M1/2;
M2 = M2/2;
window_length_1 = floor(2*M1) - floor(M1) - 1;
window_length_2 = floor(2*M2) - floor(M2) - 1;
coord_1 = 0:(1/window_length_1):1;
coord_2 = 0:(1/window_length_2):1;
[wl_1,wr_1] = fdct_wrapping_window(coord_1);
[wl_2,wr_2] = fdct_wrapping_window(coord_2);
lowpass_1 = [wl_1, ones(1,2*floor(M1)+1), wr_1];
lowpass_2 = [wl_2, ones(1,2*floor(M2)+1), wr_2];
lowpass_next = lowpass_1'*lowpass_2;
hipass = sqrt(1 - lowpass_next.^2);
Xj = zeros(2*floor(4*M1)+1,2*floor(4*M2)+1);
% Loop: angles
l = 0;
nbquadrants = 2 + 2*(~is_real);
nbangles_perquad = nbangles(j)/4;
for quadrant = 1:nbquadrants
M_horiz = M2 * (mod(quadrant,2)==1) + M1 * (mod(quadrant,2)==0);
M_vert = M1 * (mod(quadrant,2)==1) + M2 * (mod(quadrant,2)==0);
if mod(nbangles_perquad,2),
wedge_ticks_left = round((0:(1/(2*nbangles_perquad)):.5)*2*floor(4*M_horiz) + 1);
wedge_ticks_right = 2*floor(4*M_horiz) + 2 - wedge_ticks_left;
wedge_ticks = [wedge_ticks_left, wedge_ticks_right(end:-1:1)];
else
wedge_ticks_left = round((0:(1/(2*nbangles_perquad)):.5)*2*floor(4*M_horiz) + 1);
wedge_ticks_right = 2*floor(4*M_horiz) + 2 - wedge_ticks_left;
wedge_ticks = [wedge_ticks_left, wedge_ticks_right((end-1):-1:1)];
end;
wedge_endpoints = wedge_ticks(2:2:(end-1)); % integers
wedge_midpoints = (wedge_endpoints(1:(end-1)) + wedge_endpoints(2:end))/2;
% Left corner wedge
l = l+1;
first_wedge_endpoint_vert = round(2*floor(4*M_vert)/(2*nbangles_perquad) + 1);
length_corner_wedge = floor(4*M_vert) - floor(M_vert) + ceil(first_wedge_endpoint_vert/4);
Y_corner = 1:length_corner_wedge;
[XX,YY] = meshgrid(1:(2*floor(4*M_horiz)+1),Y_corner);
width_wedge = wedge_endpoints(2) + wedge_endpoints(1) - 1;
slope_wedge = (floor(4*M_horiz) + 1 - wedge_endpoints(1))/floor(4*M_vert);
left_line = round(2 - wedge_endpoints(1) + slope_wedge*(Y_corner - 1));
[wrapped_XX, wrapped_YY] = deal(zeros(length_corner_wedge,width_wedge));
first_row = floor(4*M_vert)+2-ceil((length_corner_wedge+1)/2)+...
mod(length_corner_wedge+1,2)*(quadrant-2 == mod(quadrant-2,2));
first_col = floor(4*M_horiz)+2-ceil((width_wedge+1)/2)+...
mod(width_wedge+1,2)*(quadrant-3 == mod(quadrant-3,2));
for row = Y_corner
cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
new_row = 1 + mod(row - first_row, length_corner_wedge);
admissible_cols = round(1/2*(cols+1+abs(cols-1)));
wrapped_XX(new_row,:) = XX(row,admissible_cols);
wrapped_YY(new_row,:) = YY(row,admissible_cols);
end;
slope_wedge_right = (floor(4*M_horiz)+1 - wedge_midpoints(1))/floor(4*M_vert);
mid_line_right = wedge_midpoints(1) + slope_wedge_right*(wrapped_YY - 1);
% not integers
% in general
coord_right = 1/2 + floor(4*M_vert)/(wedge_endpoints(2) - wedge_endpoints(1)) * ...
(wrapped_XX - mid_line_right)./(floor(4*M_vert)+1 - wrapped_YY);
C2 = 1/(1/(2*(floor(4*M_horiz))/(wedge_endpoints(1) - 1) - 1) + 1/(2*(floor(4*M_vert))/(first_wedge_endpoint_vert - 1) - 1));
C1 = C2 / (2*(floor(4*M_vert))/(first_wedge_endpoint_vert - 1) - 1);
wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) + (wrapped_YY-1)/floor(4*M_vert) == 2) = ...
wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) + (wrapped_YY-1)/floor(4*M_vert) == 2) + 1;
coord_corner = C1 + C2 * ((wrapped_XX - 1)/(floor(4*M_horiz)) - (wrapped_YY - 1)/(floor(4*M_vert))) ./ ...
(2-((wrapped_XX - 1)/(floor(4*M_horiz)) + (wrapped_YY - 1)/(floor(4*M_vert))));
wl_left = fdct_wrapping_window(coord_corner);
[wl_right,wr_right] = fdct_wrapping_window(coord_right);
switch is_real
case 0
wrapped_data = fftshift(fft2(ifftshift(C{j}{l})))/sqrt(prod(size(C{j}{l})));
wrapped_data = rot90(wrapped_data,(quadrant-1));
case 1
x = C{j}{l} + sqrt(-1)*C{j}{l+nbangles(j)/2};
wrapped_data = fftshift(fft2(ifftshift(x)))/sqrt(prod(size(x)))/sqrt(2);
wrapped_data = rot90(wrapped_data,(quadrant-1));
end;
wrapped_data = wrapped_data .* (wl_left .* wr_right);
% Unwrapping data
for row = Y_corner
cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
admissible_cols = round(1/2*(cols+1+abs(cols-1)));
new_row = 1 + mod(row - first_row, length_corner_wedge);
Xj(row,admissible_cols) = Xj(row,admissible_cols) + wrapped_data(new_row,:);
% We use the following property: in an assignment
% A(B) = C where B and C are vectors, if
% some value x repeats in B, then the
% last occurrence of x is the one
% corresponding to the eventual assignment.
end;
% Regular wedges
length_wedge = floor(4*M_vert) - floor(M_vert);
没有合适的资源?快使用搜索试试~ 我知道了~
曲波去噪matlab代码
共14个文件
m:14个
需积分: 50 75 下载量 83 浏览量
2017-11-22
11:20:19
上传
评论 12
收藏 23KB ZIP 举报
温馨提示
曲波去噪,包括软阈值、硬阈值、软硬折中去噪三种参数选择
资源推荐
资源详情
资源评论
收起资源包目录
curvelet.zip (14个子文件)
curvelet
fdct_wrapping_param.m 8KB
fdct_wrapping_pos2idx.m 785B
fdct_wrapping_demo_wave.m 2KB
fdct_wrapping_demo_denoise.m 4KB
ifdct_wrapping.m 16KB
fdct_wrapping_demo_basic.m 1KB
fdct_wrapping_dispcoef.m 2KB
fdct_wrapping_window.m 751B
demo_fdct_wrapping_denoise.m 8KB
fdct_wrapping_demo_recon.m 1KB
fdct_wrapping_demo_disp.m 1KB
SoftThreshCurv.m 10KB
fdct_wrapping_demo_denoise_enhanced.m 8KB
fdct_wrapping.m 15KB
共 14 条
- 1
资源评论
hanwxa
- 粉丝: 0
- 资源: 1
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功