function pbeDefinition = growthDensitySampleMetroVolume

%% call common definition
pbeDefinition = growthDensityCommon;


%% Initial distribution
rng(1);
% discrete sampling
%! ToDo: I need to understand cov to lay out the normal distribution
%        properly.
% Design code:
% stdX=1e-1;stdY=5e-1;corr=1/sqrt(2);
% disD = obtainDistributionDiscrete('normal', 'mean', [3 3], 'cov', [stdX^2 corr*stdX*stdY; corr*stdX*stdY stdY^2], 'nSample', 1e3)
% plot(disD.positionMatrix(1,:), disD.positionMatrix(2,:), 'r.')
%%
% xlim([0 2]),ylim([0 2])
std = max([pbeDefinition.seeding.std1(:); ...
           pbeDefinition.seeding.std2(:)]);
% densityHandle = @(h)(pbeDefinition.seeding.fhandle(h));
densityHandle = @(h)(pbeDefinition.seeding.fhandle(h) .* ...
    computeMeasureDec(pbeDefinition.dec, 'hC', h, 'volume', [], [], 1));
disD = obtainDistributionDiscrete('metro', ...
    'nSample', pbeDefinition.seeding.nSample, ...
    'metroHandle', densityHandle, ...
    'metroInit', pbeDefinition.seeding.mean1, ...  %m
    'metroCov', [std^2 0; 0 std^2], ...
    'metroRelTol', 0.01);
%%
% application to pbeDefinition
pbeDefinition.seeding.positionMatrix = disD.positionMatrix;       % m
pbeDefinition.seeding.numberVector = ...
    pbeDefinition.seeding.fhandle(disD.positionMatrix) ./ ...
    disD.probabilityVector; % #
pbeDefinition.seeding.probabilityVector = disD.probabilityVector;    

% close all
% plot3(disD.positionMatrix(1,:), disD.positionMatrix(2,:), ...
%     pbeDefinition.seeding.fhandle(disD.positionMatrix), 'r.')