function [heMatrix, groupVector] = resampleData(pbeDef, result, ...
    absoluteSizeSteps, nSample)


nDim = size(result(1).hcMatrix, 1);


%% get heMean matrix
hcMeanMatrix = nan(nDim, length(result));
for iTime = 1:length(result)
    hcMeanMatrix(:, iTime) = mean(result(iTime).hcMatrix, 2);
end

%% construct time filter
% first, I need the mean width of the mean shapes over time
widthVector = computeMeasureDec(pbeDef.dec, 'hC', hcMeanMatrix, ...
    'meanWidth', [], [], 1);
minSize = min(widthVector);
maxSize = max(widthVector);
sizeVector = minSize:absoluteSizeSteps:maxSize;
nTime = length(sizeVector);
timeIndexVector = zeros(1, nTime);
for iTime = 1:nTime
    timeIndexVector(iTime) = find(widthVector >= sizeVector(iTime), 1, 'first');
end
if ~any(timeIndexVector == length(result))
    timeIndexVector(end+1) = length(result);
end

heMatrix = nan(nDim, nTime*nSample);
groupVector = repmat(1:nTime, nSample, 1);
groupVector = groupVector(:)';
for iTime = 1:nTime
    heMatrix(:, (iTime-1)*nSample + (1:nSample)) = ...
        result(timeIndexVector(iTime)).hcMatrix(:, 1:nSample);
end

% nSamplePerTime = 1+2*nDim;
% heMatrix = nan(nDim, nTime*nSamplePerTime);
% groupVector = repmat(1:nTime, nSamplePerTime, 1);
% groupVector = groupVector(:)';
% for iTime = 1:nTime
%     heMatrix(:, (iTime-1)*nSamplePerTime + 1) = hcMeanMatrix(...
%         :, timeIndexVector(iTime));
%     for iDim = 1:nDim
%         coordinateVector = result(timeIndexVector(iTime)).hcMatrix(iDim, :);
%         % minimum index
%         thisIndex = find(min(coordinateVector) == coordinateVector, 1);
%         heMatrix(:, (iTime-1)*nSamplePerTime + (iDim-1)*2 + 2) = ...
%             result(timeIndexVector(iTime)).hcMatrix(:, thisIndex);
%         % maximum index
%         thisIndex = find(max(coordinateVector) == coordinateVector, 1);
%         heMatrix(:, (iTime-1)*nSamplePerTime + (iDim-1)*2 + 3) = ...
%             result(timeIndexVector(iTime)).hcMatrix(:, thisIndex);
%     end
% end