function data = runSimulation(caseName)

if isstruct(caseName)
    opt = obtainData('baseDefinition');
    opt = mergeStruct(caseName, opt);
else
    opt = obtainData(caseName);
end

%% preparation
% initialize output
data.rrVector = opt.rrVector;
data.mwMatrix = zeros(length(opt.rrVector), size(opt.dirMatrix, 2));
data.lamMatrix = zeros(length(opt.rrVector), size(opt.dirMatrix, 2));

% scale input hrep
if opt.aggregate == 1
    baseCenterMatrix = [0    0    -0.2;...
                        0.05 0.15  0.4; ...
                        0    0.4  -0.25];
    hrepBase = opt.hrepBase;
    hrepBase.h = hrepBase.h / computeMeanWidthHrep(hrepBase);
    hrepBase(2) = moveHrep(hrepBase(1), baseCenterMatrix(2,:));
    hrepBase(3) = moveHrep(hrepBase(1), baseCenterMatrix(3,:));
    hrepBase(1) = moveHrep(hrepBase(1), baseCenterMatrix(1,:));
%     close all
%     figure,hold on
%     viewPolytope(hrepBase(1))
%     viewPolytope(hrepBase(2))
%     viewPolytope(hrepBase(3))
%     return
elseif opt.aggregate == 2
        baseCenterMatrix = [0    0    -0.2;...
                            0.05 0.15  0.2];
    hrepBase = opt.hrepBase;
    hrepBase.h = hrepBase.h / computeMeanWidthHrep(hrepBase);
    hrepBase(2) = moveHrep(hrepBase(1), baseCenterMatrix(2,:));
    hrepBase(1) = moveHrep(hrepBase(1), baseCenterMatrix(1,:));
%     close all
%     figure,hold on
%     viewPolytope(hrepBase(1))
%     viewPolytope(hrepBase(2))
%     return
elseif opt.aggregate == 3
        baseCenterMatrix = [0    0    0;...
                            0    0    0.05];
    hrepBase = opt.hrepBase;
    hrepBase.h = hrepBase.h / computeMeanWidthHrep(hrepBase);
    hrepBase(2) = moveHrep(...
        stretchHrep(rotateHrep3d(hrepBase(1), 'vector', [1 1 4]), ...
                    0.9), baseCenterMatrix(2,:));
    hrepBase(1) = moveHrep(hrepBase(1), baseCenterMatrix(1,:));
%     close all
%     figure,hold on
%     viewPolytope(hrepBase(1))
%     viewPolytope(hrepBase(2))
%     return
else
    baseCenterMatrix = [0 0 0];
    hrepBase = opt.hrepBase;
    hrepBase.h = hrepBase.h / computeMeanWidthHrep(hrepBase);
end

% generate image matrix
image = false(opt.particleSize*2);
[xGrid yGrid] = meshgrid(linspace(...
    -opt.particleSize + 0.5, opt.particleSize - 0.5, ...
    2*opt.particleSize));
positionMatrix = [xGrid(:)'; yGrid(:)'] ./ opt.particleSize;

% generate structuring elements
se = cell(1, length(opt.rrVector));
seLow = cell(1, length(opt.rrVector));
seHigh = cell(1, length(opt.rrVector));
rrDiff = opt.rrDiff;
for iRR = 1:length(opt.rrVector)
    se{iRR} = strel('disk', ...
        round(0.5*opt.rrVector(iRR)*opt.particleSize), 0);
    if rrDiff ~= 0
        seLow{iRR} = strel('disk', max(0, ...
            round(0.5*(opt.rrVector(iRR)-rrDiff)*opt.particleSize)), 0);
        seHigh{iRR} = strel('disk', min(0.5*opt.particleSize, ...
            round(0.5*(opt.rrVector(iRR)+rrDiff)*opt.particleSize)), 0);
    end
end

% generate blurr filter
if opt.blurrStd ~= 0
    blurrFilter = fspecial('gaussian', opt.blurrSize, ...
        opt.blurrStd);
end

% load background image source
bgImage = imread('backgroundImage.png');

%% main loop
for iDir = 1:size(opt.dirMatrix, 2)
    % project polytope
    for iPoly = 1:length(hrepBase)
        hrep(iPoly) = rotateHrep3d(hrepBase(iPoly), 'vector', opt.dirMatrix(:, iDir));
        vrep(iPoly) = convertHrepToVrep(hrep(iPoly));
        vrep(iPoly).V = vrep(iPoly).V(:, 1:2);
        hrep(iPoly) = convertVrepToHrep(vrep(iPoly));     
    end
    % rotate center points
    vrepCenter.V = baseCenterMatrix;
    vrepCenter = rotateVrep3d(vrepCenter, 'vector', opt.dirMatrix(:, iDir));
    centerMatrix = vrepCenter.V(:,1:2);
    
for iRR = 1:length(opt.rrVector)
    disp(['iDir = ' num2str(iDir) ', iRR = ' num2str(iRR)]);
    % clear image
    image(:) = false;
    
    % write kernel polytope(s)
    for iPoly = 1:length(hrep)
    hrepIn(iPoly) = hrep(iPoly);
    hrepIn(iPoly) = stretchHrep(hrep(iPoly), 1-opt.rrVector(iRR), ...
        centerMatrix(iPoly, :)');
%     hrepIn(iPoly).h = hrepIn(iPoly).h * (1-opt.rrVector(iRR));
%     hrepIn(iPoly) = reduceHrep(hrepIn(iPoly));

    filterIn = all(bsxfun(@le, hrepIn(iPoly).A * positionMatrix, hrepIn(iPoly).h));
    if any(filterIn)
        image(filterIn) = true;
    else
        % switch a pixel on to allow imdilate - if the polytope is too
        % small to contain a single pixel.. ..then we could take the
        % closest pixel to one of the vertices
        distanceVector = sum(bsxfun(@minus, ...
            centerMatrix(iPoly,:)', positionMatrix).^2, 1);
        [dump thisIndex] = min(distanceVector);
        image(thisIndex(1)) = true;
%         image(round(end/2), round(end/2)) = true;
    end
    end
    
    % dilate with se
    image = imdilate(image, se{iRR});

    %% add differing additive sphericity
    if rrDiff ~= 0
        % lower mode (copied from above)
        imageLow = image & false;
        hrepIn = hrep;
        hrepIn.h = hrepIn.h * (1-max(opt.rrVector(iRR)-rrDiff, 0));
        hrepIn = reduceHrep(hrepIn);
        filterIn = all(bsxfun(@le, hrepIn.A * positionMatrix, hrepIn.h));
        if any(filterIn)
            imageLow(filterIn) = true;
        else
            imageLow(round(end/2), round(end/2)) = true;
        end
        imageLow = imdilate(imageLow, seLow{iRR});
           
        % higher image
        imageHigh = image & false;
        hrepIn = hrep;
        hrepIn.h = hrepIn.h * (1-min(opt.rrVector(iRR)+rrDiff, 1));
        hrepIn = reduceHrep(hrepIn);
        filterIn = all(bsxfun(@le, hrepIn.A * positionMatrix, hrepIn.h));
        if any(filterIn)
            imageHigh(filterIn) = true;
        else
            imageHigh(round(end/2), round(end/2)) = true;
        end
        imageHigh = imdilate(imageHigh, seHigh{iRR});
        
%         % transfer to image
%         thisSub = 1:optionStruct.particleSize;
%         image(thisSub, thisSub) = imageLow(thisSub, thisSub);
%         thisSub = 1:optionStruct.particleSize + optionStruct.particleSize;
%         image(thisSub, thisSub) = imageHigh(thisSub, thisSub);
    end
   
    %transfer to grayscale
    image = uint8(image);
    image(image == 1) = 30;
    % transfer to image
    if rrDiff ~= 0
        thisSub = 1:opt.particleSize;
        image(thisSub, thisSub) = imageLow(thisSub, thisSub)*35;
        thisSub = (1:opt.particleSize) + opt.particleSize;
        image(thisSub, thisSub) = imageHigh(thisSub, thisSub)*40;
    end
    
    % background
%     if all(bgInfo{1} == 0)
        image(image == 0) = 170;
%     else
%         xStart = round((size(bgImage, 1) - 5 - size(image, 1)) * ...
%                        bgInfo{1}(iRR, iDir));
%         yStart = round((size(bgImage, 2) - 5 - size(image, 2)) * ...
%                        bgInfo{2}(iRR, iDir));
%         thisBgImage = bgImage(...
%             xStart + (1:2*optionStruct.particleSize), ...
%             yStart + (1:2*optionStruct.particleSize));
%         
%         image(image == 0) = thisBgImage(image == 0);
%     end
    
    % write image if requested
    if opt.writeBaseImage
        imwrite(image, ['imageOriginal_' num2str(iDir) '_' num2str(iRR) ...
            '.png'], 'PNG');
    end
    
    % blurr image
    if opt.blurrStd ~= 0
        image = imfilter(image, blurrFilter, 'symmetric', 'conv');
    end
    
    % add noise
    if opt.noise ~= 0
        image = imnoise(image, 'gaussian', 0, opt.noise);
    end
    
    % background
    if all(opt.bgInfo{1} == 0)
%         image(image == 0) = 170;
    else
        xStart = round((size(bgImage, 1) - 5 - size(image, 1)) * ...
                       opt.bgInfo{1}(iRR, iDir));
        yStart = round((size(bgImage, 2) - 5 - size(image, 2)) * ...
                       opt.bgInfo{2}(iRR, iDir));
        thisBgImage = bgImage(...
            xStart + (1:2*opt.particleSize), ...
            yStart + (1:2*opt.particleSize));
        
%         image(image > 110) = thisBgImage(image > 110);
        image = min(image, thisBgImage);
    end
        
    % write image if requested
    if opt.writeBlurrImage
        imwrite(image, ['imageBlurr_' num2str(iDir) '_' num2str(iRR) ...
            '.png'], 'PNG');
    end
    
    % find contour lines
    image = edge(image > 110, 'sobel');
    [xBlur yBlur] = find(image == 1);
    
    % write image if requested
    if opt.writeEdgeImage
        imwrite(image, ['imageEdge_' num2str(iDir) '_' num2str(iRR) ...
            '.png'], 'PNG');
    end

    % identify
    bnd = obtainBnd([xBlur'; yBlur']);
    fitShape = fitShapeBnd(bnd, 'lambdaLocal', opt.lambdaLocal, ...
        'lambdaError', opt.lambdaError);
    data.mwMatrix(iRR, iDir) = fitShape.meanWidth;
    data.lamMatrix(iRR, iDir) = fitShape.lambda;
end
end

% refined results
data.rrMatrix = 2*data.lamMatrix ./ data.mwMatrix;

data.rrMeanVector = mean(data.rrMatrix, 2);
data.rrStdVector = std(data.rrMatrix, 0, 2);
data.rrMaxVector = max(data.rrMatrix, [], 2);
data.rrMinVector = min(data.rrMatrix, [], 2);

data.mwMeanVector = mean(data.mwMatrix, 2);
data.mwStdVector = std(data.mwMatrix, 0, 2);
data.mwMaxVector = max(data.mwMatrix, [], 2);
data.mwMinVector = min(data.mwMatrix, [], 2);

disp(['Save as: ' 'data' filesep opt.caseName '.mat'])
save(['data' filesep opt.caseName '.mat'], 'data', 'opt');
plotSimulationSummary(data, opt.caseName)
obtainData('init', 'database', 'data', 'type', 'mat')