function iridge = extractRidges(tfm,penalty,nr,nbins)
% tfm - time frequency matrix
% penalty - penalty term (0 or greater)
% nr - number of ridges
% nbins - number of bins to remove from each ridge

gamma = 1e-8;
M = size(tfm,1);
N = size(tfm,2);
c = zeros(N,nr);
% Obtain energy transform
E = abs(tfm);
E = E/sum(E(:));
E = -log(E+gamma);
% If the number of ridges is 1.
% both tfridge and wsstridge check that nr is nonempty and positive
% nr = 1 is the default in both functions
if nr == 1
  % Call into compiled code for ridge extraction
  c = signalwavelet.extractcurve(E,penalty);
  % Make sure result is column vector
  c = c(:);
  
elseif nr > 1

  
  for ni = 1:nr
    
    ctmp = signalwavelet.extractcurve(E,penalty);
    ctmp = ctmp(:);
    c(:,ni) = ctmp;

    % Remove this curve from the synchrosqueezed transform prior
    % to computing determing the next curve.
    for ii=0:nbins
      % Find lower bound for index. Do not go below 1.
      lb = max(1, c(:,ni)-ii);
      % Find upper bound. Do not exceed M
      ub = min(M, c(:,ni)+ii);
      % Replace extracted ridge with -log(realmin) to ensure that the 
      % extracted ridge has the minimum energy
      E((0:N-1)'*M+lb) = -log(realmin);
      if ii>0, E((0:N-1)'*M+ub) = -log(realmin); end
    end
  end

    
end

iridge = c;

end
