You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: dba_get_avg_vertarea.m
    path of file: ./repos/brainstorm3/external/dba
    
    the code of the file until where you have to start completion: function nbVert = dba_get_avg_vertarea(mesh_ref, mesh_2)
% 
% Compute nbVert for subsampling mesh#2 according to the average vertices area of mesh ref
% 
% Yohan Attal - HM-TC project 2013

disp = 0;

[FaceArea_ref, VertexArea_ref] = tess_area(mesh_ref.Vertices, mesh_ref.Faces);
avg_ref = mean(VertexArea_ref);
std_ref = std(VertexArea_ref);
if disp
    figure('color','w'),
    subplot(131), title('mesh ref')
    hold on
    plot(VertexArea_ref), ylabel('Vertex area')
    plot(ones(1,length(VertexArea_ref))*avg_ref,'r')
    plot(ones(1,length(VertexArea_ref))*(avg_ref+std_ref),'g--')
    plot(ones(1,length(VertexArea_ref))*(avg_ref-std_ref),'g--')
end
 
[FaceArea_2, VertexArea_2] = tess_area(mesh_2.Vertices, mesh_2.Faces);
avg_2 = mean(VertexArea_2);
std_2 = std(VertexArea_2);
if disp
    subplot(132), title('mesh 2')
    hold on
    plot(VertexArea_2), ylabel('Vertex area')
    plot(ones(1,length(VertexArea_2))*avg_2,'r')
    plot(ones(1,length(VertexArea_2))*(avg_2+std_2),'g--')
    plot(ones(1,length(VertexArea_2))*(avg_2-std_2),'g--')
end

if avg_2 > avg_ref
    disp('DBA> Error : The order of mesh are not consistent with for subsampling. Try to invert the meshes entry positions.');
    nbVert = -1;
else
    r = 1; 
    cst = 0.025; % subsampling by steps of 2.5%
    while(avg_2 <= avg_ref)
        r = r-cst; 
        [mesh_2_tmp.Faces, mesh_2_tmp.Vertices] = reducepatch(mesh_2.Faces, mesh_2.Vertices, r);
        % [fv_2_tmp.faces, fv_2_tmp.vertices] = reducepatch_subdiv
    
function [faces, vertices] = reducepatch_subdiv(TessMat, dsFactor)

newNbVertices       = length(TessMat.Vertices);

% Reduce number of vertices
[NewTessMat.Faces, NewTessMat.Vertices] = reducepatch(TessMat.Faces, TessMat.Vertices, dsFactor * 0.94);
% Find the vertices that were kept by reducepatch
[tmp, I, J] = intersect(TessMat.Vertices, NewTessMat.Vertices, 'rows');
% Save initial number for vertices
nVertLowInit = size(NewTessMat.Vertices,1);

% Progress bar
bst_progress('start', 'Resample surface', 'Analyzing surface...');
% Calulate face areas and perimeter
FaceArea  = tess_area(NewTessMat.Vertices, NewTessMat.Faces);
% Vertex connectivity, normals, Curvature
VertConn    = tess_vertconn(NewTessMat.Vertices, NewTessMat.Faces);
[VertNormals, FaceNormals] = tess_normals(NewTessMat.Vertices, NewTessMat.Faces, VertConn);
Curvature   = tess_curvature(NewTessMat.Vertices, VertConn, VertNormals);
% Get center of each face
FaceCenter = (NewTessMat.Vertices(NewTessMat.Faces(:,1),:) + NewTessMat.Vertices(NewTessMat.Faces(:,2),:) + NewTessMat.Vertices(NewTessMat.Faces(:,3),:)) ./ 3;
% Get center of mass of the vertices
SurfCenter = mean(NewTessMat.Vertices, 1);
% Get large faces to subdivide (perimeter or area)
iBigFaces = find((sum(FaceNormals .* bst_bsxfun(@minus, FaceCenter, [0 0 SurfCenter(3)]), 2) > 0.04) & ...  % Faces pointing outwards (normal in the same direction as position vector)
    (sum(Curvature(NewTessMat.Faces) > 0, 2) >= 2) & ...       % Curvature has to be > 0
    (FaceArea > mean(FaceArea) + 1*std(FaceArea)));            % Face area threshold
% If there are not enough points to add (white matter): perform search on all the surface
if (length(iBigFaces) < .75 * (newNbVertices - size(NewTessMat.Vertices,1)))
    iBigFaces = find(FaceArea > mean(FaceArea) + 2.5 * std(FaceArea));
end