for each v in vertices {
    n_accum = vec3(0,0,0);
    for each f in faces {
        count = 0;
        if(v is in f){
            n_accum += f.normal;
            count++;
        }
    }
    n = n_accum / count;
    v.normal = n;
}

//with theshold:

thredshold = degtorad(60.0)
for each f1 in faces {
    for each v in f1 {
        n_accum = vec3(0,0,0);
        for each f2 in faces {
            count = 0;
            if(v is in f2 && dot(f2.normal, f1.normal) < thredshold){
                n_accum += f2.normal;
                count++;
            }
        }
        n = n_accum / count;
        v.normal = n;
    }
}