Problem with casting rays from the camera into the scene for each pixel

Question

To retrieve information about the objects that lie behind each pixel, I have to cast rays into the scene for each pixel of a rendered frame.

My problem is that the ray directions which I calculate are not correct.

The eye rays that I generate do not match my image frame vertically (they either start/end outside of it, or too far inside).

For example here is the first direction(top left of the image space) that i calculate

And here is the last direction (bottom right of the image space)

I think that maybe my calculation of the vertical field of view might be off. But so far I have not been able to fix it.

#calculates the direction the ray has to take for a u, v coordinate
def get_ray_direction(u, v, position, direction, up, right, width, height):
    #mapping the u and v coordinates into normalised image coordinates
    x = ((2*u - width)/width) 
    y = ((2*v - height)/height) 

    direction = x * right + y * up + direction + position
    return direction

s = bpy.data.scenes["Scene"]    
cam = s.camera
camData = cam.data

view_direction = cam.matrix_world.to_quaternion() * Vector((0.0, 0.0, -1.0))
up = cam.matrix_world.to_quaternion() * Vector((0.0, 1.0, 0.0))
aspect_ratio = s.render.resolution_x / s.render.resolution_y
right = mathutils.Vector.cross(view_direction, up)

up = up.normalized()
right = right.normalized()

#For the FOV I have to use the focal length and sensor size as base values
#since I will be given these values and have to set the camera accordingly
#I took the formula from: http://paulbourke.net/miscellaneous/lens/
fov_v = 2 * math.atan((camData.sensor_height * 0.5) / camData.lens)
fov_h = 2 * math.atan((camData.sensor_width * 0.5) / camData.lens)

up = up * math.tan(fov_v/2)     
right = right * math.tan(fov_h/2)

height = s.render.resolution_y/1.0
width = s.render.resolution_x/1.0

#Starting top left, going row by row 
for i in reversed(range(s.render.resolution_y)):
    for j in range(s.render.resolution_x):       
        dir = get_ray_direction(j, i, pos, view_direction, up, right, width, height)
        ray_dir = dir * 10000
        ray_result = s.ray_cast(cam.location, ray_dir);       

This is the code that I currently have in total, which contains some debugging stuff like empties that are placed at the ray direction

Here is the primary source I used (among others, but I found this one to be most comprehensible) for the generation of eye rays

Answer 1

A solution when you render in 3D View. Maybe something in here will be handy for you. This comes from transforming mouse coords to 3d so maybe its a bit more general than you need.

Getting camera frame corners in screen space (2d):

import bpy
from bpy_extras.view3d_utils import location_3d_to_region_2d

cam = bpy.context.scene.camera
frame = cam.data.view_frame(bpy.context.scene)

# Its in local space so transform to global
frame = [cam.matrix_world * corner for corner in frame]

# Transform into screenspace
region = bpy.context.region
rv3d = bpy.context.region_data
frame_px = [location_3d_to_region_2d(region, rv3d, corner) for corner in frame]

Now you know where exactly the camera frame is and how its panned or zoomed etc. Now you can start firing rays from every pixel in that frame like so:

from bpy_extras.view3d_utils import region_2d_to_origin_3d
from bpy_extras.view3d_utils import region_2d_to_vector_3d

def 2d_to_ray(context, point_px):
    region = context.region
    rv3d = context.region_data

    ray_origin = region_2d_to_origin_3d(region, rv3d, point_px)
    ray_vector = region_2d_to_vector_3d(region, rv3d, point_px)

    return ray_origin, ray_vector

Now you have a ray vector with origin all in global space for every pixel of camera frame. Here is also how to intersect such ray with a plane (I am using a parallel plane to camera view located at 3d cursor):

def ray_intersect_plane(plane_normal, plane_point, ray_vector, ray_origin):
    '''Does not work in Ortho view'''
    t = plane_normal.dot(plane_point - ray_origin) / normal.dot(ray_vector)
    return line_point + t*ray_vector

# This camera location will work also when there is no camera
# It will work for view-port rendering too
rv3d = context.region_data
cam_loc = rv3d.view_matrix.inverted().translation
cam_look_at = rv3d.view_location
plane_normal = cam_loc - cam_look_at

location = ray_intersect_plane(plane_normal,
                               bpy.context.scene.cursor_location,
                               ray_vector,
                               ray_origin)

If you want to know the pixel dimensions of current area (where the mouse is) for the view-port rendering you do it like this:

regions = bpy.context.area.regions

width = next(region.width for region in regions if region.type=='WINDOW')
height = next(region.height for region in regions if region.type=='WINDOW')

---

A solution when you render with F12

Don't transform the camera frame points to 2d space. Use the 3d global frame corners and render resolution to interpolate all the render pixels in 3d space.

From camera location and these 3d pixels get your rays. Else is similar.

Answer 2

In case anyone needs a full working example, here is a minimal code snippet which casts rays from the camera into the scene.

import bpy
from mathutils import Vector, Quaternion
import numpy as np
import bmesh
objects to consider
obj = bpy.data.objects['suzanne']
background = bpy.data.objects['plane']
targets = [background, obj]
camera object which defines ray source
cam = bpy.data.objects['camera']
save current view mode
mode = bpy.context.area.type
set view mode to 3D to have all needed variables available
bpy.context.area.type = "VIEW_3D"
get vectors which define view frustum of camera
frame = cam.data.view_frame(scene=bpy.context.scene)
topRight = frame[0]
bottomRight = frame[1]
bottomLeft = frame[2]
topLeft = frame[3]
number of pixels in X/Y direction
resolutionX = int(bpy.context.scene.render.resolution_x * (bpy.context.scene.render.resolution_percentage / 100))
resolutionY = int(bpy.context.scene.render.resolution_y * (bpy.context.scene.render.resolution_percentage / 100))
setup vectors to match pixels
xRange = np.linspace(topLeft[0], topRight[0], resolutionX)
yRange = np.linspace(topLeft[1], bottomLeft[1], resolutionY)
array to store hit information
values = np.empty((xRange.size, yRange.size), dtype=object)
indices for array mapping
indexX = 0
indexY = 0
filling array with None
for x in xRange:
    for y in yRange:
        values[indexX,indexY] = None
        indexY += 1
    indexX += 1
    indexY = 0
iterate over all targets
for target in targets:
    # calculate origin
    matrixWorld = target.matrix_world
    matrixWorldInverted = matrixWorld.inverted()
    origin = matrixWorldInverted @ cam.matrix_world.translation
# reset indices
indexX = 0
indexY = 0

# iterate over all X/Y coordinates
for x in xRange:
    for y in yRange:
        # get current pixel vector from camera center to pixel
        pixelVector = Vector((x, y, topLeft[2]))

        # rotate that vector according to camera rotation
        pixelVector.rotate(cam.matrix_world.to_quaternion())

        # calculate direction vector
        destination = matrixWorldInverted @ (pixelVector + cam.matrix_world.translation) 
        direction = (destination - origin).normalized()

        # perform the actual ray casting
        hit, location, norm, face =  target.ray_cast(origin, direction)

        if hit:
            values[indexX,indexY] = (matrixWorld @ location)

        # update indices
        indexY += 1

    indexX += 1
    indexY = 0


create new mesh
source: https://devtalk.blender.org/t/alternative-in-2-80-to-create-meshes-from-python-using-the-tessfaces-api/7445/3
mesh = bpy.data.meshes.new(name='created mesh')
bm = bmesh.new()
iterate over all possible hits
for index, location in np.ndenumerate(values):
    # no hit at this position
    if location is None:
        continue
# add new vertex
bm.verts.new((location[0], location[1], location[2]))  


make the bmesh the object's mesh
bm.to_mesh(mesh)

bm.free()  # always do this when finished
We're done setting up the mesh values, update mesh object and
let Blender do some checks on it
mesh.update()
mesh.validate()
Create Object whose Object Data is our new mesh
obj = bpy.data.objects.new('created object', mesh)
Add Object to the scene, not the mesh
scene = bpy.context.scene
scene.collection.objects.link(obj)
Select the new object and make it active
bpy.ops.object.select_all(action='DESELECT')
obj.select_set(True)
bpy.context.view_layer.objects.active = obj
reset view mode
bpy.context.area.type = mode
print("Done.")

Answer 3

You can render a depth buffer of your scene into a bpy.types.Image / OpenGL texture and sample it later on for depth data. I think this approach will be more effective than raycasting to each (!!!) pixel.

This can be done this way:

Example 1: Rendering to Image. Only good for visualizing depth as this way we loose actual depth value due to linearization which is important to perform since depth values are not linear.

import bpy
import bgl

def linearize(depth, znear, zfar):
    return (*([(-zfar * znear / (depth * (zfar - znear) - zfar)) / zfar] * 3), 1.0)

def render_depth_to_image(pos_x: int, pos_y: int, width: int, height: int, space_view_3d: bpy.types.SpaceView3D)
    """
    Render depth to bpy.types.Image
    width, height could be region height/width or camera in your case.
    pos_x, pos_y is a starting position to render the rectangle from (0,0) for whole viewport render.
    space_view_3d is context.space_data or whatever viewport you are rendering from
    """

    depth_buf = Buffer(GL_FLOAT, width * height)
    glReadPixels(pos_x, pos_y, width, height, GL_DEPTH_COMPONENT, GL_FLOAT, depth_buf)

    image_name = "scene_depth"
    if not image_name in bpy.data.images:
        bpy.data.images.new(image_name, width, height)

    znear = space_view_3d.clip_start
    zfar = space_view_3d.clip_end

    image = bpy.data.images[image_name]
    image.scale(width, height)

    image.pixels = [y for x in [linearize(v, znear, zfar) for v in self.depth_buf] for y in x]

Example 2: Using the rendered Buffer directly. No depth data loss.

import bgl
from typing import List, Callable

class DepthRenderAndProcess:

    """
    Render depth to bpy.types.Image
    width, height could be region height/width or camera in your case.
    pos_x, pos_y is a starting position to render the rectangle from (0,0) for whole viewport render.
    processing_callbacks is a list of functions/other callables that expect one argument which is a flat list of depth values per pixel
    """

    def __init__(self, pos_x, pos_y, width, height, processing_callbacks: List[Callable]):
        self.pos_x = pos_x
        self.pos_y = pos_y
        self.width = width
        self.height = height
        self.depth_buf = self.render_depth_to_opengl_texture()

        for callback in processing_callbacks:
            callback(depth_buf.to_list())


    def render_depth_to_opengl_texture(self):
        depth_buf = Buffer(GL_FLOAT, width * height)
        glReadPixels(self.pos_x, self.pos_y, self.width, self.height, GL_DEPTH_COMPONENT, GL_FLOAT, depth_buf)

        return depth_buf