function computeApicobasal(o)

if ~o.available.rotational
    o.computeRotational;
end
o.printStatus('Computing apicobasal coordinate...');
t = toc;

% parameters
numTmVal = 10; % number of contour lines in transmural direction
numRtVal = 48; % number of contour lines in rotational direction, ideally a multiple of 6
numClPoints = 100; % number of points per contour line after resampling using cubic splines

% compute contour surfaces of tm (cs)
tmp = vtkDeleteDataArrays(o.m1.vol);
tmp.pointData.tm = o.m1.tm;
tmVal = linspace(0.5/numTmVal, 1-0.5/numTmVal, numTmVal);
cs = vtkContourFilter(tmp, 'points', 'tm', tmVal);

% assign regions to cs
cs = vtkConnectivityFilter(cs);
cs.pointData.csRegion = cs.pointData.RegionId;
cs = vtkDeleteDataArrays(cs, {'pointData.csRegion'});

% interpolate abLaplace from o.m2.vol to cs
Mcs = baryInterpMat(o.m2.vol.points, o.m2.vol.cells, cs.points);
cs.pointData.abLaplace = Mcs * o.m2.abLaplace;

if o.cfg.exportLevel > 2
    vtkWrite(cs, sprintf('%sabContourSurfaces.vtp', o.cfg.outPrefix));
end

% find apex point for each csRegion
csRegions = unique(cs.pointData.csRegion);
apexPoints = NaN(numel(csRegions),3);
for i = 1:numel(csRegions)
    pointIds = find(cs.pointData.csRegion==csRegions(i));
    [~,minId] = min(cs.pointData.abLaplace(pointIds));
    apexPoints(i,:) = cs.points(pointIds(minId),:);
end

% compute contour lines of rtSin (rtSinCl)
tmp = cs;
tmp.pointData.rtSin = Mcs * o.m2.rtSin;
rtSinVal = sin(pi*linspace(-0.25, 0.25, round(numRtVal/2)+1));
rtSinCl = vtkContourFilter(tmp, 'points', 'rtSin', rtSinVal);
rtSinCl.pointData = rmfield(rtSinCl.pointData, 'rtSin');

% compute contour lines of rtCos (rtCosCl)
tmp = cs;
tmp.pointData.rtCos = Mcs * o.m2.rtCos;
rtCosVal = rtSinVal(2:end-1);
rtCosCl = vtkContourFilter(tmp, 'points', 'rtCos', rtCosVal);
rtCosCl.pointData = rmfield(rtCosCl.pointData, 'rtCos');

% combine rtSin and rtCos contour lines (cl)
cl = vtkAppendPolyData({rtSinCl, rtCosCl});

% exclude points close to the apex to split the contour lines into two parts
% apexRadius needs to be large enough for each contour line to be splitted
apexRadius = 3*o.m1.meanEdgLen;
cl.pointData.clRegion = zeros(size(cl.points,1),1);
for i = 1:numel(csRegions)
    pointIds = find(cl.pointData.csRegion==csRegions(i));
    ids = rangesearch(cl.points(pointIds,:), apexPoints(i,:), apexRadius);
    cl.pointData.clRegion(pointIds(ids{1})) = -1;
end
cl = vtkThreshold(cl, 'points', 'clRegion', [0 inf]);

% identify separate contour lines --> clRegion
cl = vtkConnectivityFilter(cl);
cl.pointData.clRegion = cl.pointData.RegionId;
cl.pointData = rmfield(cl.pointData, 'RegionId');
cl = rmfield(cl, 'cellData');

% exclude too short contour lines based on values of abLaplace
clRegions = unique(cl.pointData.clRegion);
abLaplaceMin = NaN(numel(clRegions),1);
abLaplaceMax = NaN(numel(clRegions),1);
pointIds = cell(numel(clRegions),1);
for i = 1:numel(clRegions)
    pointIds{i} = find(cl.pointData.clRegion==clRegions(i));
    abLaplaceVal = cl.pointData.abLaplace(pointIds{i});
    abLaplaceMin(i) = min(abLaplaceVal);
    abLaplaceMax(i) = max(abLaplaceVal);
end
abLaplaceMinThresh = min(abLaplaceMin) + 2*(median(abLaplaceMin)-min(abLaplaceMin));
abLaplaceMaxThresh = 0.99;
for i = 1:numel(clRegions)
    if abLaplaceMin(i) > abLaplaceMinThresh || abLaplaceMax(i) < abLaplaceMaxThresh
        cl.pointData.clRegion(pointIds{i}) = -1;
    end
end
cl = vtkThreshold(cl, 'points', 'clRegion', [0 inf]);

% reassign clRegion
cl = vtkConnectivityFilter(cl);
cl.pointData.clRegion = cl.pointData.RegionId;
cl.pointData = rmfield(cl.pointData, 'RegionId');
cl = rmfield(cl, 'cellData');

if o.cfg.exportLevel > 2
    vtkWrite(cl, sprintf('%sabContourLines.vtp', o.cfg.outPrefix));
end

% use cubic smoothing spline to smooth and resample the contour lines
% and to compute the normalized distance along the contour lines (abSmoothCl)
clRegions = unique(cl.pointData.clRegion);
splineMisfit = 5e-4 * norm(mean(double(o.m1.sur.points(o.m1.sur.pointData.class==1,:)),1)-mean(apexPoints,1));
smoothClPoints = NaN(numel(clRegions)*numClPoints, 3);
abSmoothCl = zeros(numel(clRegions)*numClPoints, 1);
abLength = abSmoothCl;
for i = 1:numel(clRegions)
    pointIds = find(cl.pointData.clRegion==clRegions(i));
    [~,sortInd] = sort(cl.pointData.abLaplace(pointIds));
    apexPoint = apexPoints(csRegions==cl.pointData.csRegion(pointIds(sortInd(1))),:);
    P = [apexPoint; double(cl.points(pointIds(sortInd),:))];
    d = sqrt(sum(diff(P,1,1).^2,2));
    ind = d < 1e-4*o.m1.meanEdgLen;
    P(ind,:) = [];
    d(ind) = [];
    d = [0; cumsum(d)];
    w = [100; ones(size(P,1)-1,1)];
    tol = numel(d)*splineMisfit^2;
    sp = spaps(d, P', tol, w);
    P = fnval(sp, linspace(min(d), max(d), numClPoints))';
    d = [0; cumsum(sqrt(sum(diff(P,1,1).^2,2)))];
    ind = (i-1)*numClPoints+1:i*numClPoints;
    smoothClPoints(ind,:) = P;
    abSmoothCl(ind) = d/max(d);
    abLength(ind) = max(d);
end

if o.cfg.exportLevel > 2
    clSmooth.points = single(smoothClPoints);
    clSmooth.pointData.ab = single(abSmoothCl);
    clSmooth.cells = int32(repmat(reshape(repmat(0:numClPoints:numel(clRegions)*numClPoints-1,numClPoints-1,1),[],1),1,2) + repmat([(1:numClPoints-1)' (2:numClPoints)'],numel(clRegions),1));
    clSmooth.cellTypes = repmat(uint8(3), size(clSmooth.cells,1), 1);
    vtkWrite(clSmooth, sprintf('%sabContourLinesSmooth.vtp', o.cfg.outPrefix));
end

% BEGIN: Laplacian extrapolation of abSmoothCl to o.m1.vol

% ||M ab - abSmoothCl|| forces extrapolated values to fit to contour values
M = baryInterpMat(o.m1.vol.points, o.m1.vol.cells, smoothClPoints);

% ||E ab - 1|| forces extrapolated values to 1 at the base
baseIds = double(o.m1.surToVol(o.m1.sur.pointData.class==1));
baseVal = ones(numel(baseIds),1);
E = sparse(1:numel(baseIds), baseIds, ones(size(baseIds)), numel(baseIds), size(o.m1.vol.points,1));
eta = (numel(abSmoothCl)/(numel(baseIds)))^2;

% ||L ab|| forces extrapolated values to be smooth
L = o.m1.massMat \ o.m1.L;

% Initial guess for ab using nearest-neighbor interpolation
ab = abSmoothCl(knnsearch(smoothClPoints, o.m1.vol.points, 'NSMethod','kdtree'));

% Initial guess for lambda based on the relation between lambda and the 
% half-width at half-maximum of the point spread function of the operator
% inv(speye(size(L))+lambda*(L'*L))
hwhm = mean(abLength)/100;
lambda = 1.58*hwhm^3.57;
lambda = numel(abSmoothCl)/numel(ab)*lambda;
fprintf('\n0\t%.3e\n', lambda);

b = M'*abSmoothCl + eta*E'*baseVal;
MM = M'*M + eta*(E'*E);
LL = L'*L;
extrapMisfit = o.cfg.abExtrapSmooth/100;

try
    [~,flag,iter] = secant(@objFun, 1e-1*lambda, lambda, 1e-2*extrapMisfit);
    if flag
        warning('Secant stopped at iteration %i without converging.', iter);
    end
catch err
    disp(getReport(err, 'extended', 'hyperlinks', 'off'));
    warning('Optimization of lambda failed. Using default value instead.');
    objFun(lambda);
end

function objVal = objFun(lambda)
    A = MM + lambda*LL;
    icMat = ichol_autocomp(A, struct('michol','on'));
    [ab, flag, relres, iter] = pcg(A, b, o.cfg.tol, o.cfg.maxit, icMat, icMat', ab);
    if flag
        error('pcg failed at iteration %i with flag %i and relative residual %.1e.', iter, flag, relres);
    end
    objVal = rms(M*ab-abSmoothCl)-extrapMisfit;
end

% END: Laplacian extrapolation

o.m1.ab = min(max(ab,0),1);
o.m0.ab = min(max(o.m1.M*o.m1.ab,0),1);
o.result.pointData.ab = single(o.m0.ab);

if o.cfg.exportLevel > 1
    o.m1.debug.pointData.ab = single(o.m1.ab);
    vtkWrite(o.m1.debug, sprintf('%sdebug1.vtu', o.cfg.outPrefix));
    
    o.m0.debug.pointData.ab = single(o.m0.ab);
    vtkWrite(o.m0.debug, sprintf('%sdebug0.vtu', o.cfg.outPrefix));
end

o.printStatus(sprintf('%.1f seconds\n', toc-t), true);
o.available.apicobasal = true;

end