Set transform of a child object to an arbitrary point using a script I've been working on a script to assign transforms to child objects, such that the transforms will be transferred to Unity ( relative to the origin set in Blender, along with their rotation ) when using the experimental .fbx options. After another user kindly shared a very helpful script with me, I refactored it to match my requirements. Below is the script.
import bpy
import math
from mathutils import Quaternion

context = bpy.context

bpy.ops.object.mode_set()
baseRotation = Quaternion()
for someObject in context.selected_objects:
    polygonsLength = len( someObject.data.polygons )
    if polygonsLength &lt; 1:
        continue

print(someObject)
objectData = someObject.data
objectWorldMatrix = someObject.matrix_world
objectWorldMatrixInverted = objectWorldMatrix.inverted()


polygonsAverage = [0.0, 0.0, 0.0]
for everyPolygon in someObject.data.polygons:
    for i in range(3):
        polygonsAverage[i] += everyPolygon.center[i]
for i in range(3):
    polygonsAverage[i] /= polygonsLength

smallestDistanceBetweenCenters = 1000000000
closestPolygonToCenter = someObject.data.polygons[0]
for everyPolygon in someObject.data.polygons:
    distanceToAveragePoint = math.sqrt( 
        pow( polygonsAverage[0] - everyPolygon.center[0], 2) + 
        pow( polygonsAverage[1] - everyPolygon.center[1], 2) +
        pow( polygonsAverage[2] - everyPolygon.center[2], 2) )
    if distanceToAveragePoint &lt; smallestDistanceBetweenCenters:
        smallestDistanceBetweenCenters = distanceToAveragePoint
        closestPolygonToCenter = everyPolygon
print(closestPolygonToCenter.center)

somePolygonInObject = closestPolygonToCenter


somePolygonCenter = objectWorldMatrix @ somePolygonInObject.center
somePolygonNormal = objectWorldMatrix @ somePolygonInObject.normal

rotationDifference = baseRotation.rotation_difference(somePolygonNormal.to_track_quat())
rotationDifferenceMatrix = rotationDifference.to_matrix().to_4x4()
rotationDifferenceMatrix.translation = somePolygonCenter
objectData.transform(objectWorldMatrixInverted @ rotationDifferenceMatrix.inverted())
someObject.matrix_world = rotationDifferenceMatrix

The script finds the mean of the object's constituent polygons' centers, and then picks the center of the closest polygon as the overall object's center. This worked reasonably well, however when checking the resultant transform center, it at times was not within the object or not close to the object's center, as the chosen polygon's center was not inside the object. Image below.

I would like to adjust this script to be able to assign purely the mean center of the polygons as the object center, provided that this mean center is within a  constituent polygon, relative to the origin. In thinking about this, it seems to me that the most reliable approach to checking for this would be some kind of raycasting from the origin to the mean center, and ensuring that one of the polygons collides with this ray, though I'm not sure if there may be a better / simpler approach to this.
How could I perform such a check, and how would I then assign the actual mean transform (which is currently being stored as a simple 3-element array) to the object's transform?
Any help would be greatly appreciated! :)
Use Vector maths.
One of the fundamentals of 3D programming is Vectors. If have two location vectors A and B,
A vector C describing the shortest distance from A to B
C = B - A

distance from A to B
C.length

there is almost never a need to calculate this using the Euclidean distance formula.
Making a small adjustment to answer to Calculating and exporting global rotation of faces created from splitting an isosphere
Instead of choosing the single face find the poly that matches the criteria of the question.

Average the face centers
Sort by distance from face center to calculated average
Use the closest

Since this is returning a face, there is no way (without weird ngons) this can be &quot;outside the faces&quot;
Note this will favour small tightly clumped faces.  Consider perhaps weighting by the face areas How to compute the centroid of a mesh with triangular faces without using Blender&#39;s built-in functions?
Also a point on another face may be closer to the median of face center medians.  Aside from raycast, Object.closest_point_on_mesh(...) will also return the face.
hit, loc, norm, index = ob.closest_point_on_mesh(o)

if hit:
    f = ob.data.polygons[index]

Test script: Run in object mode, sets origin to each selected object based on criteria above.
import bpy
from mathutils import Quaternion, Vector

context = bpy.context

z = Quaternion()

def closest_poly_to_centers_avg(polys):
    o = sum((p.center for p in polys), Vector()) / len(polys)
    return(sorted(polys, key=lambda p: (p.center - o).length).pop(0))

for ob in context.selected_objects:
    print(ob)
    me = ob.data
    if not me.polygons:
        continue
    mw = ob.matrix_world
    mwi = mw.inverted()
    f = closest_poly_to_centers_avg(me.polygons)
    c = mw @ f.center
    n = mw @ f.normal

    q = z.rotation_difference(n.to_track_quat())
    M =  q.to_matrix().to_4x4()
    M.translation = c
    me.transform(mwi @ M.inverted())
    ob.matrix_world = M

FWIW I  really loathe this
for everyObject in objects:

having a variable named every object that is one of each, may &quot;read well&quot;  but really? Equally with someObject, somePolygon opposed to simply ob and f