Transform mesh to line up with empties, based on other child empties
In the end I decided to just learn how Blender's scripting works and wrote a script. I know there would be simpler ways to calculate the transformations but I'm no maths expert and I wanted to spend my time debugging my logic errors rather my maths errors.

Basically what this script does is take 3 input points (in the form of empties) relative to the photogrammetry mesh and 3 input points in the same position relative to the tracking data, then attempts to line them up with each other, bringing the photogrammetry mesh along with it. This actually works pretty well in my testing. It's obviously not the fastest workflow but it does work and that's the important part.

import bpy
import bmesh
from mathutils import Matrix

#####################################
# TO USE:
#
# 1) create empties at three easily-identified positions in tracking data, called ea1, ea2 and ea3
# 2) create empties at the same three positions but in photogrammetry data, called eb1, eb2 and eb3
#    NOTE: eb1, eb2 and eb3 must be positioned IN THE SAME ORDER as ea1, ea2 and ea3!
# 3) name the photogrammetry mesh "photogrammetry"
# 4) run the script and be amazed
#
#####################################

# from https://blender.stackexchange.com/a/118783/41192
def createLookRotation(forward, up):
    rot = Matrix.Identity(3)
    rot[0] = up                 # x
    rot[1] = forward.cross(up)  # y
    rot[2] = forward            # z
    rot = rot.transposed()
    return rot.to_quaternion()

def doMeshThing(obj, mesh, empty1, empty2, empty3):
    bm = bmesh.new()
    bm.verts.new(empty1.location)
    bm.verts.new(empty2.location)
    bm.verts.new(empty3.location)
    bm.verts.ensure_lookup_table()
    bm.faces.new([bm.verts[0], bm.verts[1], bm.verts[2]])
    bm.normal_update()
    bm.to_mesh(mesh)
    bpy.ops.object.select_all(action='DESELECT')
    obj.select = True
    bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY')
    bm.faces.ensure_lookup_table()
    normal = bm.faces[0].normal
    guide = empty2.location - empty1.location
    quatRot = createLookRotation(normal, guide)
    obj.rotation_mode = 'QUATERNION'
    obj.rotation_quaternion = quatRot.inverted()
    bpy.ops.object.transform_apply(location=False, rotation=True, scale=False)
    obj.rotation_quaternion = quatRot
    bm.free()

ea1 = bpy.data.objects['ea1']
ea2 = bpy.data.objects['ea2']
ea3 = bpy.data.objects['ea3']
eb1 = bpy.data.objects['eb1']
eb2 = bpy.data.objects['eb2']
eb3 = bpy.data.objects['eb3']

mesha = bpy.data.meshes.new('mesha')
meshb = bpy.data.meshes.new('meshb')

obja = bpy.data.objects.new('obja', mesha)
objb = bpy.data.objects.new('objb', meshb)

scene = bpy.context.scene
scene.objects.link(obja)
scene.objects.link(objb)

doMeshThing(obja, mesha, ea1, ea2, ea3)
doMeshThing(objb, meshb, eb1, eb2, eb3)

bpy.context.scene.update()

# position and rotate objb (photogrammetry) to match position and rotation of obja (tracking)
photogrammetryMesh = bpy.data.objects['photogrammetry']
photogrammetryMesh.parent = objb
photogrammetryMesh.matrix_parent_inverse = objb.matrix_world.inverted()

objb.location = obja.location
objb.rotation_quaternion = obja.rotation_quaternion

# calculate average dimension difference here and scale
scaleDiffX = obja.dimensions.x / objb.dimensions.x
scaleDiffY = obja.dimensions.y / objb.dimensions.y
scaleDiffZ = obja.dimensions.z / objb.dimensions.z
#averageScaleDiff = (scaleDiffX + scaleDiffY + scaleDiffZ) / 3
averageScaleDiff = (scaleDiffX + scaleDiffY) / 2 # seems to be in local space, so z is useless
objb.scale.x = averageScaleDiff
objb.scale.y = averageScaleDiff
objb.scale.z = averageScaleDiff

bpy.context.scene.update()

parented_wm = photogrammetryMesh.matrix_world.copy()
photogrammetryMesh.parent = None
photogrammetryMesh.matrix_world = parented_wm

bpy.ops.object.select_all(action='DESELECT')
obja.select = True
bpy.ops.object.delete()
bpy.ops.object.select_all(action='DESELECT')
objb.select = True
bpy.ops.object.delete()