def METHOD_NAME(point,cell,normal,gamma=1):
    n = normal[cell[:,0],:] + normal[cell[:,1],:] + normal[cell[:,2],:]
    v01 = point[cell[:,1],:] - point[cell[:,0],:]
    v02 = point[cell[:,2],:] - point[cell[:,0],:]
    nv = np.cross(v01,v02,1)
    length = np.sqrt(np.square(nv).sum(axis=1))
    area = length/2.0
    nv = nv/length.reshape((-1, 1))
    cellRho = np.abs(9-np.sum(n**2,1))/area
    cellRho = cellRho + np.sqrt(np.spacing(1))
    for i in range(12):
        sub = cell.flatten(0)
        cellRho_area = cellRho*area
        val = np.matlib.repmat(cellRho_area, 3, 1).flatten(0)
        pointRho = np.bincount(sub, weights=val, minlength=point.shape[0])
        sarea = np.matlib.repmat(area, 3, 1).flatten(0)
        b = np.bincount(sub, weights=sarea, minlength=point.shape[0])
        pointRho = pointRho/b
        cellRho = (pointRho[cell[:,0]] + pointRho[cell[:,1]] + pointRho[cell[:,2]])/3
    cur = pointRho
    pointRho = (cellRho/np.max(cellRho))**(gamma)
    return pointRho, cellRho, cur, area