function isosurface_ar = area_regularize(isosurface_vf, alpha, nbr, face_around_vertex)

n_vertex=size(nbr,1);

% calculate area of each triangle (area_all) and average area (avg_area_all)

area_all = zeros(size(isosurface_vf.faces,1),1);

for i=1:size(isosurface_vf.faces,1)
    triangle_vertices = isosurface_vf.vertices(isosurface_vf.faces(i,:),:);
    area_all(i,1) = triangle_area_3d(triangle_vertices');
end

avg_area_all = mean(area_all);

% area regularization

for i=1:n_vertex
    current_nbr=nbr(i,nbr(i,:)~=0);
    
    % calculate average area of each neighbor (avg_area)
    
    avg_area = zeros(size(current_nbr,2),1);
    
    for j=1:size(current_nbr,2)
        current_faces = face_around_vertex(current_nbr(1,j),:);
        area_faces = area_all(unique(current_faces(1,current_faces~=0)),1);
        avg_area(j,1) = mean(area_faces);
    end
    
    % regularize areas
    
    if (min(avg_area)<avg_area_all)
        [temp,max_area] = max(avg_area);
        isosurface_vf.vertices(i,:) = alpha*isosurface_vf.vertices(i,:)+(1-alpha)*isosurface_vf.vertices(current_nbr(1,max_area),:);
    end
end

isosurface_ar = isosurface_vf;