Costantini相位解包裹MATLAB代码，最小费用最大流

function PhiUnwrap = cunwrap（Psi options）
% PhiUnwrap = cunwrap（Psi options）
%
% Purpose: Costantini‘s 2D unwrapping based on Minimum Network Flow
%
% INPUTS:
%   - Psi: 2D-array wrapped phase （radian）
%   - options: optional structure used to modify the behavior of the 
%              unwrapping algorithm. All fields are optional.
%              The fieldname does not need to exactly case-matched.
%       Weight: array of same dimension as the Psi with is used as
%               the positive weight of the L1 norm of the residual （to be
%               minimized）. For example user can provide WEIGHT as function
%               of the interference amplitude or function of phase
%               measurement quality （e.g. fringe contrast）.
%               By default: all internal pixels have the same weight of 1
%               0.5 for edge-pixels and 0.25 for corner-pixels. 
%       CutSize: even integer scalar default [4]: the width （in pixel） of
%                the Gaussian kernel use to lower the weight around pixels
%                that do not fulfill rotational relation.
%                If CutSize is 0 no rotational-based weighting will be
%                carried out.
%       RoundK: Boolean [false] by default. When ROUNDK is true CUNWRAP
%               forces the partial derivative residuals K1/K2 to be
%               integers.
%       MaxBlockSize: scalar default [125] target linear blocksize.
%                     Set to Inf for single-block global unwrapping.
%                     Note: network flow is costly in CPU for large size.
%       Overlap: scalar default [0.25] fractional overlapping between two
%                neighboring blocks.
%       Verbose: logical [true] control printout information.
%       LPoption: LINPROG option structure as return by OPTIMSET（...）.
% OUTPUT:
%   - Phi: 2D-array unwrapped phase （radian）
%
% Minimum network flow computation engine is Matlab LINPROG （optimization
%   toolbox is required）
%
% References:
%   - http://earth.esa.int/workshops/ers97/papers/costantini/
%   - Costantini M. （1998） A novel phase unwrapping method based on network
%   programming. IEEE Tran. on Geoscience and Remote Sensing 36 813-821.
%
% Author: Bruno Luong 
% History:
%   Orginal: 27-Aug-2009
%       28-Aug-2009: modification in block splitting default size changes
%                    to 125 x 125
%       28-Aug-2009: correct a serious indexing bug

if nargin<2
    options = struct（）;
end

if ndims（Psi）>2
    error（‘CUNWRAP: input Psi must be 2D array.‘）;
end

[ny nx] = size（Psi）;

if nx<2 || ny<2
    error（‘CUNWRAP: size of Psi must be larger than 2‘）;
end

roundK = getoption（options ‘roundk‘ false）;

verbose = getoption（options ‘verbose‘ true）;

% Default weight
w1 = ones（ny1）; w1（[1 end]） = 0.5;
w2 = ones（1nx）; w2（[1 end]） = 0.5;
weight = w1*w2; % tensorial product
weight = getoption（options ‘weight‘