vcycle.m

采用matlab编写的数字图像恢复程序

  function [uh] = vcycle(uh,bh,dcoef_xh,dcoef_yh,level,p0,nu)

%  [uh] = vcycle(uh,bh,dcoef_xh,dcoef_yh,level,p0,nu)
%
%  Apply one Multigrid V-cycle to the system
%    Lh*uh = bh.
%  Lh is the CCFD discretization of the diffusion operator 
%    L*v = -div (dcoef grad u).
 
  [nx,ny] = size(uh);
  N = nx*ny;
  Lh = get_L(dcoef_xh,dcoef_yh);
  L = tril(Lh,-1);
  D = spdiags(diag(Lh),0,N,N); 
  U = triu(Lh,1);
  bh = Qproject(bh,p0);
  
%  Perform nu steps of Gauss-Seidel presmoothing.
 
  for gsiter = 1:nu
    uh = uh(:) + (D + L)\(bh(:) - Lh*uh(:));
    uh = reshape(uh,nx,ny);
%    uh = Qproject(reshape(uh,nx,ny),p0);  %  Project onto orthocomp. of VH.
  end
 
%  Recursive call to vcycle.

  if level > p0 + 1
    level = level-1;
    rh = reshape( bh(:) - Lh*uh(:) , nx,ny);
    [r2h] = restrict(rh);
    [dcoef_xH,dcoef_yH] = restrict_dcoef(dcoef_xh,dcoef_yh);
    u2h = zeros(size(r2h));
    [u2h] = vcycle(u2h,r2h,dcoef_xH,dcoef_yH,level,p0,nu);
    uh = uh + prolong(u2h);
  else
%%%    uh = pinv(full(Lh)) * bh;
  end

%  Perform nu steps of Gauss-Seidel postsmoothing.  Note that pre- and 
%  post-smoothing operators are adjoints.
 
  for gsiter = 1:nu
    uh(:) = uh(:) + (D + U)\(bh(:) - Lh*uh(:));
%    uh = reshape( Qproject(uh,p0) , nx,ny);  %  Project onto orthocomp. of V0.
  end

uh = reshape( Qproject(uh,p0) , nx,ny);  %  Project onto orthocomp. of V0.