How to correct surface normal for translated object

Question

**Edit - I'm willing to pay a $15USD bounty via Paypal for the accepted answer to this question.

I have a script which generates random points on the surface of a sphere and then places cone objects at each point.

Using the emitter.closest_point_on_mesh call I can get the surface normal from the sphere at a spot very close to the generated point and then apply that normal to the rotation of the cone object, thereby pointing the cone out from the surface of the sphere at an angle that very closely approximates the normal from the sphere surface.

However this only works when the sphere is at the world origin point. When I move the sphere, it no longer works. The cones no longer point away from the sphere surface at the normal direction.

I've tried different techniques involved multiplying the sphere's matrix_world value by the normal, but I can't seem to find the right combo to make it work.

Any help would be GREATLY appreciated!

import bpy
from mathutils.bvhtree import BVHTree
import mathutils
import math
import bpy_extras
from bpy_extras import mesh_utils
import timeit
from random import random
from math import radians
from mathutils import Vector
import bmesh

scn = bpy.context.scene
objs = bpy.data.objects


def DeletePreviousParticles():

    for o in bpy.data.objects:
        if o.name.startswith("particle"):
            objs.remove(o, True)
    #select sphere to update window
    sphere = bpy.data.objects['Sphere']
    bpy.context.scene.objects.active = sphere



def PlaceParticle(emitter,src_obj, point):

    dupCount = 0
    maxObjectDimension = 2


    new_obj = src_obj.copy()
    new_obj.data = src_obj.data.copy()
    scn.objects.link(new_obj)
    new_obj.animation_data_clear()

    new_obj.location[0]=point[0]
    new_obj.location[1]=point[1]
    new_obj.location[2]=point[2]


    new_obj.name = "particle_" + str(dupCount)
    dupCount+=1


    #calculate normal from surface of emitter 
    partLoc = new_obj.location
    surfaceNormal = emitter.closest_point_on_mesh( partLoc )[1]
    print('surfaceNormal' + str(surfaceNormal))

    #the following line does not work when emitter object is not at world origin point
    new_obj.rotation_euler = surfaceNormal.to_track_quat('Z', 'Y').to_euler() 





def Main():

    emitter = bpy.data.objects['Sphere']
    src_obj = bpy.data.objects['sourceParticle']

    print("starting------------------------------------------------------------")

    DeletePreviousParticles()

    me = emitter.data
    me.calc_tessface() # recalculate tessfaces
    tessfaces_select = me.tessfaces #[f for f in me.tessfaces if f.select]
    pointList = bpy_extras.mesh_utils.face_random_points(1, tessfaces_select)

    totalPoints = 20
    pointCounter = 0

    for p in pointList:

            if(pointCounter<totalPoints):

                mat1 = emitter.matrix_world
                adjPoints = mat1 * p

                PlaceParticle(emitter, src_obj, adjPoints)

                pointCounter +=1
            else:
                break    

    print("done")



Main()

Here's the source file:

@rjg - Thanks again for your help! I've made a 15 Euro donation to the Blender Foundation as requested.

Answer 1

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 >= 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 > Tools > 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