How to compute the centroid of a mesh with triangular faces without using Blender's built-in functions?
Bmesh approach

The bmesh module has for faces the methods to calculate centre and area.

There are three methods to calculate the centre of the face.


  calc_center_bounds()

Return bounds center of the face.
Returns:  a 3D vector.
Return type:  mathutils.Vector

  
  calc_center_median()

Return median center of the face.
Returns:  a 3D vector.
Return type:  mathutils.Vector

  
  calc_center_median_weighted()

Return median center of the face weighted by edge lengths.
Returns:  a 3D vector.
Return type:  mathutils.Vector



I've used calc_center_median() for $x_i$ and BMFace.calc_area() for the weight $w_i$ for each face of index $i$ the weighted mean is

$$\tilde x = \frac {w_0x_0 + w_1x_1 + ... + w_nx_n}{w_0 + w_1 + ... + w_n}$$

import bpy
import bmesh
from mathutils import Vector
context = bpy.context

ob = context.object
me = ob.data
bm = bmesh.new()
bm.from_mesh(me)

weights = [f.calc_area() for f in bm.faces]
weighted_centres = [f.calc_area() * f.calc_center_median() 
        for f in bm.faces]

local_weighted_mean = sum(weighted_centres, Vector()) / sum(weights)
print(local_weighted_mean)

# set as new origin
for v in bm.verts:
    v.co -= local_weighted_mean
mw = ob.matrix_world
mw.translation = mw * local_weighted_mean
bm.to_mesh(me)
me.update()



Result of running on newly added Suzanne. Scene cursor shows old origin.

Note: if modifiers are used on the object to triangulate, or subdivide, can (for example) use

bm.from_object(ob, context.scene, render=True)


to calculate weighted origin of modified mesh.