How to correct surface normal for translated object
There are two problems in the script, the location and the rotation of the cone objects aren't correct.

The location is incorrect because the script currently sets it to the position of the randomly chosen point, not the location of the closest point on the surface. These coordinates can be retrieved through the return values of closest_point_on_mesh():

(result, location, normal, index) = emitter.closest_point_on_mesh(point)


The global coordinates are therefore:

emitter.matrix_world * location




If the goal is to align an arbitrarily chosen axis of an object with the normal, then the problem that needs to be solved is to find a rotation matrix. This matrix has to rotate the vector of the local axis to match the normal vector exactly, assuming that the vector is normalized.

One way of determining the rotation matrix is through the Rodrigues' rotation formula. It requires an angle and an axis of rotation. The angle $\theta$ between two vectors can be determined through the dot product, which is defined as:

$$a \cdot b =  \|a\| \|b\| \cos(\theta)$$

It follows

$$\theta = \arccos\left(\frac{a \cdot b}{\|a\| \|b\|}\right) $$

Any rotation angle in range $[0, \pi)$ radians or $[0°, 180°)$ can be used unambiguously. If $a$ and $b$ are unit vectors, then their magnitude is one and the division isn't necessary.

In order to rotate the vector, an axis for the rotation is required. Two vectors that a linearly independent define a plane and the cross product $a \times b$ is the normal to that plane, which is perpendicular to the two vectors.  

Blender already implements the Rodrigues' formula in mathutils.Matrix.Rotation(angle, size, axis). Assuming we have a normalized vector $a$ that we want to align to vector $b$, the rotation matrix can be determined by the following code:

cross_prod = a.cross(b) # Axis used for rotation
dot_prod = np.arccos(a.dot(b)) # Angle in radians (ambiguous if &gt;= pi)
rotation_matrix = mathutils.Matrix.Rotation(dot_prod, 4, cross_prod)


This could actually be optimized and not use arccos altogether as explained by Kevin Moran.

You can chose an arbitrary vector for $a$ as long as it isn't linearly dependent to $b$. For your specific application it makes sense that $a$ is either one of the global axis or one of the local axis of the cone object in the global coordinate system. For the latter you will have to decompose the world matrix of the cone object in order to determine it's rotation matrix and multiply with the desired axis.

# Determine the vector that will be aligned with the surface normal
if axis == 'X':
    axis_vector = mathutils.Vector((1.0, 0.0, 0.0))
elif axis == 'Y':
    axis_vector = mathutils.Vector((0.0, 1.0, 0.0))
else:
    axis_vector = mathutils.Vector((0.0, 0.0, 1.0))

# Apply the existing rotation of the object when the local axis is used
if use_local_axis:
    (old_location, old_rotation, old_scale) = new_obj.matrix_world.decompose()
    axis_vector = old_rotation * axis_vector


The adjustment of both rotation and location can be combined by creating a new world matrix for the target object.

new_obj.matrix_world = (mathutils.Matrix.Translation(emitter.matrix_world * location) * mathutils.Matrix.Rotation(angle, 4, rotation_axis) * mathutils.Matrix.Scale(1.0, 4))




Below you can find an implementation of the described solution as an add-on for Blender 2.79.



bl_info = {
    "name": "Create and Align Objects with Surface Normal",
    "author": "Robert Guetzkow",
    "version": (1, 0),
    "blender": (2, 79, 0),
    "location": "3D View &gt; Tools &gt; Create and Align",
    "description": "Create objects that are aligned with the surface normal of another object",
    "warning": "",
    "wiki_url": "",
    "category": "3D View"}

import bpy
import bpy_extras
import mathutils
import numpy as np
from bpy_extras import mesh_utils


class CreateAlignSettings(bpy.types.PropertyGroup):
    emitter = bpy.props.PointerProperty(type=bpy.types.Object,
                                        name="Emitter",
                                        description="The emitters surface is used to create"
                                                    " and align objects on it, according to its surface normal")

    dupli_object = bpy.props.PointerProperty(type=bpy.types.Object,
                                             name="Dupli Object",
                                             description="The object duplicated along the emitters surface")

    axis = bpy.props.EnumProperty(name="Axis",
                                  description="The axis of the duplicated object that will be aligned with the surface "
                                              "normal",
                                  items=[('X', 'X', 'X-axis', 0),
                                         ('Y', 'Y', 'Y-axis', 1),
                                         ('Z', 'Z', 'Z-axis', 2)],
                                  default='Z')

    use_local_axis = bpy.props.BoolProperty(name="Local Axis",
                                            description="Use the possibly rotated local coordinate system of the "
                                                        "duplicated object, to determine the axis orientation in "
                                                        "global space. If disable the global axis is used and the "
                                                        "object's rotation is ignored.")

    number_duplicates = bpy.props.IntProperty(name="Number of Duplicates",
                                              default=0)


class CreateAlign:
    @staticmethod
    def create(emitter, dupli_object, axis, use_local_axis, number_duplicates):
        # Remove previously created duplicate objects
        CreateAlign.remove_dupli_objects(emitter, dupli_object)

        emitter.data.calc_tessface()
        surface_points = bpy_extras.mesh_utils.face_random_points(1, emitter.data.tessfaces)

        for idx, point in enumerate(surface_points):
            if idx == number_duplicates:
                break

            CreateAlign.duplicate_and_align(emitter, dupli_object, axis, use_local_axis, point)

    @staticmethod
    def duplicate_and_align(emitter, dupli_object, axis, use_local_axis, point):
        new_obj = dupli_object.copy()
        new_obj.data = dupli_object.data.copy()
        bpy.context.scene.objects.link(new_obj)
        new_obj.animation_data_clear()
        new_obj.parent = emitter

        (result, location, normal, index) = emitter.closest_point_on_mesh(point)

        # Determine the vector that will be aligned with the surface normal
        if axis == 'X':
            axis_vector = mathutils.Vector((1.0, 0.0, 0.0))
        elif axis == 'Y':
            axis_vector = mathutils.Vector((0.0, 1.0, 0.0))
        else:
            axis_vector = mathutils.Vector((0.0, 0.0, 1.0))

        # Apply the existing rotation of the object when the local axis is used
        if use_local_axis:
            (old_location, old_rotation, old_scale) = new_obj.matrix_world.decompose()
            axis_vector = old_rotation * axis_vector

        normalized_axis_vector = axis_vector.normalized()
        rotation_axis = normalized_axis_vector.cross(normal)
        angle = np.arccos(normalized_axis_vector.dot(normal))

        # Create a new world matrix with the position of the closests surface point
        # and the rotation to align with the surface normal
        new_obj.matrix_world = (mathutils.Matrix.Translation(emitter.matrix_world * location) *
                                mathutils.Matrix.Rotation(angle, 4, rotation_axis) *
                                mathutils.Matrix.Scale(1.0, 4))

    @staticmethod
    def remove_dupli_objects(emitter, dupli_object):
        children = emitter.children

        for child in children:
            if child.name.startswith(dupli_object.name):
                bpy.data.objects.remove(child, True)


class CREATE_ALIGN_OT_create_and_align(bpy.types.Operator):
    bl_idname = "create_align.create_and_align"
    bl_label = "Create"
    bl_description = "Create duplicated objects that are aligned with the surface normal of the emitter"
    bl_options = {'REGISTER', 'UNDO'}

    def execute(self, context):
        settings = context.scene.create_align_settings
        CreateAlign.create(settings.emitter,
                           settings.dupli_object,
                           settings.axis,
                           settings.use_local_axis,
                           settings.number_duplicates)
        return {'FINISHED'}


class CREATE_ALIGN_OT_remove(bpy.types.Operator):
    bl_idname = "create_align.remove"
    bl_label = "Remove"
    bl_description = "Remove duplicated objects"
    bl_options = {'REGISTER', 'UNDO'}

    def execute(self, context):
        settings = context.scene.create_align_settings
        CreateAlign.remove_dupli_objects(settings.emitter, settings.dupli_object)
        return {'FINISHED'}


class CREATE_ALIGN_PT_surface_normal(bpy.types.Panel):
    bl_label = "Create and Align with Surface Normal"
    bl_category = "Create and Align"
    bl_space_type = "VIEW_3D"
    bl_region_type = "TOOLS"

    def draw(self, context):
        layout = self.layout
        layout.prop(context.scene.create_align_settings, "emitter")
        layout.prop(context.scene.create_align_settings, "dupli_object")
        layout.prop(context.scene.create_align_settings, "axis")
        layout.prop(context.scene.create_align_settings, "use_local_axis")
        layout.prop(context.scene.create_align_settings, "number_duplicates")
        layout.separator()
        layout.operator(CREATE_ALIGN_OT_create_and_align.bl_idname)
        layout.operator(CREATE_ALIGN_OT_remove.bl_idname)


classes = (CreateAlignSettings,
           CREATE_ALIGN_OT_create_and_align,
           CREATE_ALIGN_OT_remove,
           CREATE_ALIGN_PT_surface_normal)


def register():
    for cls in classes:
        bpy.utils.register_class(cls)
    bpy.types.Scene.create_align_settings = bpy.props.PointerProperty(type=CreateAlignSettings)


def unregister():
    for cls in classes:
        bpy.utils.unregister_class(cls)
    del bpy.types.Scene.create_align_settings


if __name__ == "__main__":
    register()


Porting to Blender 2.80:
Note that I haven't ported this to 2.80, because calc_tessface and related functions have been removed and I didn't know whether you had special requirements for the point selection. You would need to use a different approach to determine random face points.

Other necessary changes include:


Matrix multiplication is @ instead of *.
Annotation would be used instead of assignment for the properties.
The tool sidebar doesn't exist anymore, therefore the panel should have bl_region_type = "UI"




Instead of sending me the offered money, please consider donating it to either one of these organisations:


Médecins Sans Frontières / Doctors Without Borders
Blender Foundation