BBTools

bbempoisson

This function is not yet fully documented. This is a transcript of the text-formatted help.

bbempoisson  Solve a non-negative system with Poisson noise
   X=bbempoisson(A,B) use the iterative Expectation Maximization (EM)
   algorithm of K. Lange and R. Carson to find a Maximum Likelihood (ML)
   solution of the linear system A*X=L, where A is a non-negative matrix,
   X is an unknown non-negative vector, and L is the average of the
   quantity being measured.

   Instead of L we are provided with a vector, B, containing independent
   samples from a Poisson-distribution. The mean of each element (or,
   equivalently, the parameter lambda) is given by L. This model is
   suitable for emission tomography and images measured in low-light
   conditions (e.g. astronomy).

   The result is a local ML estimate, and may not correspond to the global
   maximum. Furthermore, the solution may amplify noise and result in
   solutions with a large variance.

   NOISE. Several strategies are available to reduce noise. An easy and
   popular solution is post filtering, which may be applied after
   bbempoisson (e.g. using bbconvn).

   Another possibility is early stopping, which is supported via the
   options 'maxi' and 'all'. The first limit the number of iterations,
   the second returns all iterates allowing a selection strategy.

   SIEVES. Sieves have several purposes, and may be specified as an
   operator or a function handle. An operator must be N-by-N, N=size(A,2),
   and represent a non-negative transform. Similarly, a function must
   transform a non-negative column-vector of length N to another, but may
   be non-linear.

   Sieves are commonly used to control noise, using a blurring operation
   (e.g. bbconvn). However, they can also be used to enforce conditions
   such as constant values, known values, and limits. (The latter was
   even suggested in the original paper by Lange and Carson.)

   OPTIONS. Options are specified as [X,OPT2]=bbempoisson(A,B,'opt1',VAL1).
   Flags are not followed by a value. OPT2 is a struct with informational
   about the computation.

   The following options are available:
     'init':   Followed by a starting-vector, X0. This is the first
               approximation in the procedure (default: vector of ones).
     'maxi':   Limit on the number of iterations. This can be used to
               implement early stopping (default: 50).
     'tol':    Tolerance criteria. The iteration is stopped when 2-norm
               of two iterates has a relative difference smaller than
               the tolerance (default: 1e-3).
     'all':    If this flag is specified, X will be returned as a matrix
               where each column is an iterate.
     'sieve':  Operator, or function handle, applied on each iteration.
               Extra arguments may be passed to a function by putting
               the function handle and arguments in a cell-array.

   The struct OPT2 contains 'maxi' and 'tol' and the following fields:
     'iter':   The number of iterations actually used
     'norm':   A vector containing the 2-norm of each iterate


   Example:
     dim1=[150,150]; I0=ones(dim1); I0(25:75,40:100)=5;
     figure(1); imagesc(I0); axis image; colormap(gray);

     [bb,dim2,opt]=bbtomoemit(dim1,0:2:358);
     S=reshape(bb*I0(:),dim2);
     S=bbrand('poisson',S);
     figure(2); imagesc(S); colormap(gray);

     [xx,yy]=meshgrid(-8:8,-8:8); K=exp(-(xx.*xx+yy.*yy)/.5);
     Sieve=bbconvn(dim1,K,'same');
     X=bbempoisson(bb,S(:),'sieve',Sieve);
     figure(3); imagesc(reshape(X,dim1)); axis image; colormap(gray);

   See also bbradon, bbtomoemit, bbfbp, bbrand.

bbdocbbeye