function [grid, steps] = setup_grid(image_width, unwarped_image, unwarped_points, placement_method, n_points)
% SETUP_GRID: Create a grid for points (shapes) to be embedded
%             in.
%
% Code written by Katherine Smith, 2003
%
%    GENERAL
%
%      [grid, steps] = setup_grid(image_width, unwarped_image,
%                  unwarped_points, placement_method, n_points)
%
%    INPUT/S
%
%      -image_width:
%           The width of the input image.
%
%      -unwarped_image:
%           The input image.
%
%      -unwarped_points:
%           The control points of this image.
%
%      -placement_method:
%           One of: 'overlap grid', 'grid', 'random',
%           'random_and_scale', 'edge', 'edge_and_scale'.
%
%      -n_points:
%           Number of points in grid?
%           
%    OUTPUT/S
%
%      -grid:
%           The grid in 1-D
%
%      -steps:
%           The steps on the grid in 1-D
%
%    PENDING WORK
%
%      -
%
%    KNOWN BUG/S
%
%      -
%
%    COMMENT/S
%
%      -steps must be specified for single-point warps but not for
%      multi-point warps
%
%    RELATED FUNCTION/S
%
%      
%
%    ABOUT
%
%      -Created:     November 23rd, 2003
%      -Last update: Novermber 25th, 2003
%      -Revision:    0.0.2
%      -Author:      R. S. Schestowitz, University of Manchester
% ==============================================================

steps = [];

switch placement_method,
    
  case 'overlap_grid'
	ctr = 1;
	for i=1:5
        step = 2 / (i + 1); % 2/2=1, 2/3, 2/4=1/2, 2/5...
        for j=1:i % will execute 1,2,3... times
          grid(ctr)= -1 + j * step; % 0 + stepsize*number of steps
          steps(ctr) = step; % CTR STANDS FOR COUNTER, NOT CENTRE - save corresponding step size
          ctr = ctr + 1; % increment step index
        end
	end
    
  case 'grid'
    grid = -0.9:0.2:0.9; % divide into 10 segments (10 lines or knot-points)
  
  case 'alternating_grid'
%     random_number = rand;
%     if (rand < 0.3)
%        grid = -0.9:0.2:0.9; % divide into 10 segments
%     elseif (rand > 0.7)
%        grid = -0.9:0.3:0.9; % divide into 7 segments
%     else
%        grid = -0.8:0.4:0.8; % divide into 5 segments
%     end
   grid = -0.45:0.1:0.45;
      % ...still to be modified...
   
  case 'random'
    for ctr = 1:n_points
      grid(ctr) = rand * 2 - 1; % choose one to begin with... do-while loop would fit better here
      while(min(abs(grid(1:ctr-1) - grid(ctr))) < 0.01) % find the minimum of the differences and ensure it is very small
        grid(ctr) = rand * 2 - 1; % generate the next random value for a line margin in the grid until it's a good choice
      end
    end

  case 'random_and_scale'
      % width is only used here, but it is not really image width!!
    for ctr = 1:n_points
      steps(ctr) = image_width + 1;
      while (steps(ctr) > image_width / 2) % iterate until the value is good enough, i.e. small enough
        steps(ctr) = abs(randn(1)) * (image_width / 2) * 0.7; % some mysterious formulation
      end
      grid(ctr) = -1 + steps(ctr) + rand * (image_width - 2 * steps(ctr)); % set grid according to step above and apply some randomisation to step size too
%     grid(ctr) = steps(ctr) + 1 + rand*(image_width - 2*steps(ctr));
    end

  case 'edge'
    % Smith: put knotpoints on n_knots strongest edges
    % RSS: steps however is not set, hence no 'scale'
    n_knots = n_points;
      % Smith set this to 2
    edges = diff(unwarped_image);
    [dummy,indices] = sort(abs(edges));
    indices = flipud(indices);
    which = indices(1:n_knots);
    grid = unwarped_points(which);
  
  case 'edge_and_scale'
      % put knotpoints on n_knots strongest edges
    n_knots = n_points; 
      % Smith set this to 2 for bump data
    edges = diff(unwarped_image);
    [dummy,indices] = sort(abs(edges));
      % absolute value as we do not account for "mirroring" effects
    indices = flipud(indices);
      % sort them by magnitude
    which = indices(1:n_knots); 
      % choose these which are the strongest n edges
    grid = unwarped_points(which);
    steps = min([1 - abs(grid), 1 + abs(grid)], [], 2);
  
  otherwise
    error(['Placement method ',placement_method,' is unknown']);
end