function [grain_r_vecs, grain_mask, grain_bb, grain_avg_r_vec] = gtDefEBSDSegmentGrainFromSeed(EBSD_r_map, EBSD_mask, seed, threshold)
    if (~exist('threshold', 'var') || isempty(threshold))
        threshold = 5;
    end
    if (isempty(EBSD_mask))
        EBSD_mask = any(EBSD_r_map > eps('single'), 3);
    end
    if (~EBSD_mask(seed(1), seed(2)))
        error('gtDefEBSDSegmentGrainFromSeed:wrong_seed', ...
            'The chosen seed (%d, %d) is outside the EBSD map valid region', seed);
    end

    ref_r_vec = reshape(EBSD_r_map(seed(1), seed(2), :), [], 1);
    EBSD_r_map_dmvol = permute(EBSD_r_map, [1 2 4 3]);
    [EBSD_r_map_gvdm, EBSD_r_map_dmvol_size] = gtDefDmvol2Gvdm(EBSD_r_map_dmvol);

    valid_indices = find(EBSD_mask);
    EBSD_r_map_gvdm = EBSD_r_map_gvdm(:, valid_indices);

    dists = EBSD_r_map_gvdm - ref_r_vec(:, ones(numel(valid_indices), 1));
    dists = sqrt(sum(dists .^ 2, 1));
    dists = 2 * atand(dists);
    full_grain_mask = inf(EBSD_r_map_dmvol_size(1:2));
    full_grain_mask(valid_indices) = dists;
    full_grain_mask = full_grain_mask < threshold;

    marker = false(EBSD_r_map_dmvol_size(1:2));
    marker(seed(1), seed(2)) = true;

    full_grain_mask = imreconstruct(marker, full_grain_mask);

    regprop = regionprops(full_grain_mask, 'BoundingBox', 'Area');
    if (length(regprop) > 1)
        % let's take the biggest! (usually happens without morpho)
        ind     = find([regprop.Area] == max([regprop.Area]), 1);
        regprop = regprop(ind);
    end

    grain_bb = ceil(regprop.BoundingBox([2 1 4 3]));
    crop_bb = [grain_bb(1:2), 1, grain_bb(3:4), 1];

    % ceil because Matlab regionprops yields "0.5" (instead of 1)
    % for bb start values
    grain_mask = gtCrop(full_grain_mask, crop_bb);

    grain_inds = find(grain_mask);
    grain_inds = cat(1, grain_inds, grain_inds + numel(grain_mask), ...
        grain_inds + 2 * numel(grain_mask));
    full_grain_inds = find(full_grain_mask);
    full_grain_inds = cat(1, full_grain_inds, ...
        full_grain_inds + numel(full_grain_mask), ...
        full_grain_inds + 2 * numel(full_grain_mask));

    grain_r_vecs = zeros([size(grain_mask), 3], 'single');
    grain_r_vecs(grain_inds) = EBSD_r_map(full_grain_inds);

    r_vecs = reshape(grain_r_vecs(grain_inds), [], 3);
    grain_avg_r_vec = sum(r_vecs, 1) / size(r_vecs, 1);
end
