function [x0, r, w, maxW] = fitCircleRobust(pts, initx0, initr, normals, ...
                                            sigmaGM, beta)
%
% function [x0, r, w, maxW] = fitCircleRobust(pts, initx0, initr,
%                                  normals, sigmaGM, beta)
%
%  minimize sum_i  rho([n_i^T (x_i-x_0) - r)]^2 +  beta [1+n_i' nrml_i])
%  w.r.t. x_0 and r.
% Input:
%  pts: an Nx2 matrix containing 2D x_i, i=1...N
%  initx0: 2-vector initial guess for centre
%  initr: initial guess for radius 
%  normals: an Nx2 matrix of corresponding 2D edge normal directions nrml_i  
%  beta is a weight on normal energy component of the objective function
% Output
%  x0 - 2x1 fitted circle's center position
%  r  - fitted circle's radius
%  w  - N x 1 robust weights for edgels.
%  maxW  - maximum possible robust weight (i.e. for zero error)
%          This may be used by the calling code to scale the weights 
%          to be between 0 and 1.
  
x0 = initx0(:)';  % Make it a row vector.

% FINISH THIS CODE