function R = rtsBaselineCorrect(A, b, dim)

%This function baselines an n-dimensional array given a series of points b
%B = baseline-corrected data
%A = raw data
%b = baseline points in the format [b0 b1 ... bn]
%dim = dimension along which the baseline correction will occur


%allocate space for the sub-component indices
sz = size(A);
inds = repmat({1},1,ndims(A));

%allocate space for the result
R = zeros(size(A));

%for each baseline point
for bp = 1:length(b) - 1
    
    %extract the data at the first baseline point
    n = b(bp);
    inds = getIndices(A, n, dim);
    BP1 = A(inds{:});

    %extract the data at the second baseline point
    n = b(bp+1);
    inds = getIndices(A, n, dim);
    BP2 = A(inds{:});

    %for each point in between the baseline points
    for p = b(bp):b(bp+1)

        %compute the weighting of the linear function between BP1 and BP2
        a2 = (p - b(bp))/(b(bp+1) - b(bp));
        a1 = 1.0 - a2;

        n = p;
        inds = getIndices(A, n, dim);

        baseFunc = a1 * BP1 + a2 * BP2;
        R(inds{:}) = A(inds{:}) - baseFunc;
    end
end


function I = getIndices(A, n, dim)

%this function returns indices for all elements in S at
%point n along dimension dim

sz = size(A);
I = repmat({1},1,ndims(A));

for d = 1:ndims(A)
    if d ~= dim
        I{d} = 1:sz(d);
    else
        I{d} = n;
    end
end