Cycles: add Fisheye Lens Polynomial camera model

This allows real world cameras to be modeled by specifying the coordinates of a
4th degree polynomial that relates a pixels distance (in mm) from the optical
center on the sensor to the angle (in radians) of the world ray that is
projected onto that pixel.

This is available as part of the panoramic lens type, however it can also be
used to model lens distortions in projective cameras for example.

Differential Revision: https://developer.blender.org/D12691

9 files changed, 214 insertions, 1 deletions

diff --git a/intern/cycles/blender/CMakeLists.txt b/intern/cycles/blender/CMakeLists.txt
index b4a4d487355..1f05b1aa234 100644
--- a/intern/cycles/blender/CMakeLists.txt
+++ b/intern/cycles/blender/CMakeLists.txt
@@ -87,6 +87,7 @@ endif()
 
 set(ADDON_FILES
   addon/__init__.py
+  addon/camera.py
   addon/engine.py
   addon/operators.py
   addon/osl.py
diff --git a/intern/cycles/blender/addon/camera.py b/intern/cycles/blender/addon/camera.py
new file mode 100644
index 00000000000..d4133796875
--- /dev/null
+++ b/intern/cycles/blender/addon/camera.py
@@ -0,0 +1,84 @@
+#
+# Copyright 2011-2021 Blender Foundation
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+# <pep8 compliant>
+
+# Fit to match default projective camera with focal_length 50 and sensor_width 36.
+default_fisheye_polynomial = [-1.1735143712967577e-05,
+                              -0.019988736953434998,
+                              -3.3525322965709175e-06,
+                              3.099275275886036e-06,
+                              -2.6064646454854524e-08]
+
+# Utilities to generate lens polynomials to match built-in camera types, only here
+# for reference at the moment, not used by the code.
+def create_grid(sensor_height, sensor_width):
+    import numpy as np
+    if sensor_height is None:
+        sensor_height = sensor_width / (16 / 9)  # Default aspect ration 16:9
+    uu, vv = np.meshgrid(np.linspace(0, 1, 100), np.linspace(0, 1, 100))
+    uu = (uu - 0.5) * sensor_width
+    vv = (vv - 0.5) * sensor_height
+    rr = np.sqrt(uu ** 2 + vv ** 2)
+    return rr
+
+
+def fisheye_lens_polynomial_from_projective(focal_length=50, sensor_width=36, sensor_height=None):
+    import numpy as np
+    rr = create_grid(sensor_height, sensor_width)
+    polynomial = np.polyfit(rr.flat, (-np.arctan(rr / focal_length)).flat, 4)
+    return list(reversed(polynomial))
+
+
+def fisheye_lens_polynomial_from_projective_fov(fov, sensor_width=36, sensor_height=None):
+    import numpy as np
+    f = sensor_width / 2 / np.tan(fov / 2)
+    return fisheye_lens_polynomial_from_projective(f, sensor_width, sensor_height)
+
+
+def fisheye_lens_polynomial_from_equisolid(lens=10.5, sensor_width=36, sensor_height=None):
+    import numpy as np
+    rr = create_grid(sensor_height, sensor_width)
+    x = rr.reshape(-1)
+    x = np.stack([x**i for i in [1, 2, 3, 4]])
+    y = (-2 * np.arcsin(rr / (2 * lens))).reshape(-1)
+    polynomial = np.linalg.lstsq(x.T, y.T, rcond=None)[0]
+    return [0] + list(polynomial)
+
+
+def fisheye_lens_polynomial_from_equidistant(fov=180, sensor_width=36, sensor_height=None):
+    import numpy as np
+    return [0, -np.radians(fov) / sensor_width, 0, 0, 0]
+
+
+def fisheye_lens_polynomial_from_distorted_projective_polynomial(k1, k2, k3, focal_length=50, sensor_width=36, sensor_height=None):
+    import numpy as np
+    rr = create_grid(sensor_height, sensor_width)
+    r2 = (rr / focal_length) ** 2
+    r4 = r2 * r2
+    r6 = r4 * r2
+    r_coeff = 1 + k1 * r2 + k2 * r4 + k3 * r6
+    polynomial = np.polyfit(rr.flat, (-np.arctan(rr / focal_length * r_coeff)).flat, 4)
+    return list(reversed(polynomial))
+
+def fisheye_lens_polynomial_from_distorted_projective_divisions(k1, k2, focal_length=50, sensor_width=36, sensor_height=None):
+    import numpy as np
+    rr = create_grid(sensor_height, sensor_width)
+    r2 = (rr / focal_length) ** 2
+    r4 = r2 * r2
+    r_coeff = 1 + k1 * r2 + k2 * r4
+    polynomial = np.polyfit(rr.flat, (-np.arctan(rr / focal_length / r_coeff)).flat, 4)
+    return list(reversed(polynomial))
diff --git a/intern/cycles/blender/addon/properties.py b/intern/cycles/blender/addon/properties.py
index 3dd6a02946b..42f1e8f6f1a 100644
--- a/intern/cycles/blender/addon/properties.py
+++ b/intern/cycles/blender/addon/properties.py
@@ -33,6 +33,7 @@ from math import pi
 # enums
 
 from . import engine
+from . import camera
 
 enum_devices = (
     ('CPU', "CPU", "Use CPU for rendering"),
@@ -72,6 +73,8 @@ enum_panorama_types = (
     ('FISHEYE_EQUISOLID', "Fisheye Equisolid",
                           "Similar to most fisheye modern lens, takes sensor dimensions into consideration"),
     ('MIRRORBALL', "Mirror Ball", "Uses the mirror ball mapping"),
+    ('FISHEYE_LENS_POLYNOMIAL', "Fisheye Lens Polynomial",
+     "Defines the lens projection as polynomial to allow real world camera lenses to be mimicked."),
 )
 
 enum_curve_shape = (
@@ -891,6 +894,32 @@ class CyclesCameraSettings(bpy.types.PropertyGroup):
         default=pi,
     )
 
+    fisheye_polynomial_k0: FloatProperty(
+        name="Fisheye Polynomial K0",
+        description="Coefficient K0 of the lens polinomial",
+        default=camera.default_fisheye_polynomial[0], precision=6, step=0.1, subtype='ANGLE',
+    )
+    fisheye_polynomial_k1: FloatProperty(
+        name="Fisheye Polynomial K1",
+        description="Coefficient K1 of the lens polinomial",
+        default=camera.default_fisheye_polynomial[1], precision=6, step=0.1, subtype='ANGLE',
+    )
+    fisheye_polynomial_k2: FloatProperty(
+        name="Fisheye Polynomial K2",
+        description="Coefficient K2 of the lens polinomial",
+        default=camera.default_fisheye_polynomial[2], precision=6, step=0.1, subtype='ANGLE',
+    )
+    fisheye_polynomial_k3: FloatProperty(
+        name="Fisheye Polynomial K3",
+        description="Coefficient K3 of the lens polinomial",
+        default=camera.default_fisheye_polynomial[3], precision=6, step=0.1, subtype='ANGLE',
+    )
+    fisheye_polynomial_k4: FloatProperty(
+        name="Fisheye Polynomial K4",
+        description="Coefficient K4 of the lens polinomial",
+        default=camera.default_fisheye_polynomial[4], precision=6, step=0.1, subtype='ANGLE',
+    )
+
     @classmethod
     def register(cls):
         bpy.types.Camera.cycles = PointerProperty(
diff --git a/intern/cycles/blender/camera.cpp b/intern/cycles/blender/camera.cpp
index b5afcfa7fda..dffcceed42a 100644
--- a/intern/cycles/blender/camera.cpp
+++ b/intern/cycles/blender/camera.cpp
@@ -69,6 +69,12 @@ struct BlenderCamera {
   float pole_merge_angle_from;
   float pole_merge_angle_to;
 
+  float fisheye_polynomial_k0;
+  float fisheye_polynomial_k1;
+  float fisheye_polynomial_k2;
+  float fisheye_polynomial_k3;
+  float fisheye_polynomial_k4;
+
   enum { AUTO, HORIZONTAL, VERTICAL } sensor_fit;
   float sensor_width;
   float sensor_height;
@@ -200,6 +206,12 @@ static void blender_camera_from_object(BlenderCamera *bcam,
     bcam->longitude_min = RNA_float_get(&ccamera, "longitude_min");
     bcam->longitude_max = RNA_float_get(&ccamera, "longitude_max");
 
+    bcam->fisheye_polynomial_k0 = RNA_float_get(&ccamera, "fisheye_polynomial_k0");
+    bcam->fisheye_polynomial_k1 = RNA_float_get(&ccamera, "fisheye_polynomial_k1");
+    bcam->fisheye_polynomial_k2 = RNA_float_get(&ccamera, "fisheye_polynomial_k2");
+    bcam->fisheye_polynomial_k3 = RNA_float_get(&ccamera, "fisheye_polynomial_k3");
+    bcam->fisheye_polynomial_k4 = RNA_float_get(&ccamera, "fisheye_polynomial_k4");
+
     bcam->interocular_distance = b_camera.stereo().interocular_distance();
     if (b_camera.stereo().convergence_mode() == BL::CameraStereoData::convergence_mode_PARALLEL) {
       bcam->convergence_distance = FLT_MAX;
@@ -422,7 +434,8 @@ static void blender_camera_sync(Camera *cam,
   cam->set_full_height(height);
 
   /* panorama sensor */
-  if (bcam->type == CAMERA_PANORAMA && bcam->panorama_type == PANORAMA_FISHEYE_EQUISOLID) {
+  if (bcam->type == CAMERA_PANORAMA && (bcam->panorama_type == PANORAMA_FISHEYE_EQUISOLID ||
+                                        bcam->panorama_type == PANORAMA_FISHEYE_LENS_POLYNOMIAL)) {
     float fit_xratio = (float)bcam->render_width * bcam->pixelaspect.x;
     float fit_yratio = (float)bcam->render_height * bcam->pixelaspect.y;
     bool horizontal_fit;
@@ -465,6 +478,12 @@ static void blender_camera_sync(Camera *cam,
   cam->set_latitude_min(bcam->latitude_min);
   cam->set_latitude_max(bcam->latitude_max);
 
+  cam->set_fisheye_polynomial_k0(bcam->fisheye_polynomial_k0);
+  cam->set_fisheye_polynomial_k1(bcam->fisheye_polynomial_k1);
+  cam->set_fisheye_polynomial_k2(bcam->fisheye_polynomial_k2);
+  cam->set_fisheye_polynomial_k3(bcam->fisheye_polynomial_k3);
+  cam->set_fisheye_polynomial_k4(bcam->fisheye_polynomial_k4);
+
   cam->set_longitude_min(bcam->longitude_min);
   cam->set_longitude_max(bcam->longitude_max);
 
diff --git a/intern/cycles/kernel/camera/projection.h b/intern/cycles/kernel/camera/projection.h
index 0aea82fa812..ddbe8ad5748 100644
--- a/intern/cycles/kernel/camera/projection.h
+++ b/intern/cycles/kernel/camera/projection.h
@@ -146,6 +146,51 @@ fisheye_equisolid_to_direction(float u, float v, float lens, float fov, float wi
   return make_float3(cosf(theta), -cosf(phi) * sinf(theta), sinf(phi) * sinf(theta));
 }
 
+ccl_device_inline float3
+fisheye_lens_polynomial_to_direction(float u, float v, float coeff0, float4 coeffs,
+                                     float fov, float width, float height)
+{
+  u = (u - 0.5f) * width;
+  v = (v - 0.5f) * height;
+
+  float r = sqrtf(u * u + v * v);
+  float r2 = r*r;
+  float4 rr = make_float4(r, r2, r2*r, r2*r2);
+  float theta = -(coeff0 + dot(coeffs, rr));
+
+  if (fabsf(theta) > 0.5f * fov)
+    return zero_float3();
+
+  float phi = safe_acosf((r != 0.0f) ? u / r : 0.0f);
+
+  if (v < 0.0f)
+    phi = -phi;
+
+  return make_float3(cosf(theta), -cosf(phi) * sinf(theta), sinf(phi) * sinf(theta));
+}
+
+ccl_device float2 direction_to_fisheye_lens_polynomial(float3 dir, float coeff0, float4 coeffs,
+                                                       float width, float height)
+{
+  float theta = -safe_acosf(dir.x);
+
+  float r = (theta - coeff0) / coeffs.x;
+
+  for (int i=0; i<20; i++) {
+    float r2 = r*r;
+    float4 rr = make_float4(r, r2, r2*r, r2*r2);
+    r = (theta - (coeff0 + dot(coeffs, rr))) / coeffs.x;
+  }
+
+  float phi = atan2f(dir.z, dir.y);
+
+  float u = r * cosf(phi) / width + 0.5f;
+  float v = r * sinf(phi) / height + 0.5f;
+
+  return make_float2(u, v);
+}
+
+
 /* Mirror Ball <-> Cartesion direction */
 
 ccl_device float3 mirrorball_to_direction(float u, float v)
@@ -191,6 +236,10 @@ ccl_device_inline float3 panorama_to_direction(ccl_constant KernelCamera *cam, f
       return mirrorball_to_direction(u, v);
     case PANORAMA_FISHEYE_EQUIDISTANT:
       return fisheye_to_direction(u, v, cam->fisheye_fov);
+    case PANORAMA_FISHEYE_LENS_POLYNOMIAL:
+      return fisheye_lens_polynomial_to_direction(
+        u, v, cam->fisheye_lens_polynomial_bias, cam->fisheye_lens_polynomial_coefficients, cam->fisheye_fov,
+        cam->sensorwidth, cam->sensorheight);
     case PANORAMA_FISHEYE_EQUISOLID:
     default:
       return fisheye_equisolid_to_direction(
@@ -207,6 +256,10 @@ ccl_device_inline float2 direction_to_panorama(ccl_constant KernelCamera *cam, f
       return direction_to_mirrorball(dir);
     case PANORAMA_FISHEYE_EQUIDISTANT:
       return direction_to_fisheye(dir, cam->fisheye_fov);
+    case PANORAMA_FISHEYE_LENS_POLYNOMIAL:
+      return direction_to_fisheye_lens_polynomial(
+          dir, cam->fisheye_lens_polynomial_bias, cam->fisheye_lens_polynomial_coefficients,
+          cam->sensorwidth, cam->sensorheight);
     case PANORAMA_FISHEYE_EQUISOLID:
     default:
       return direction_to_fisheye_equisolid(
diff --git a/intern/cycles/kernel/types.h b/intern/cycles/kernel/types.h
index b44386e0da0..c39289224ad 100644
--- a/intern/cycles/kernel/types.h
+++ b/intern/cycles/kernel/types.h
@@ -478,6 +478,7 @@ enum PanoramaType {
   PANORAMA_FISHEYE_EQUIDISTANT = 1,
   PANORAMA_FISHEYE_EQUISOLID = 2,
   PANORAMA_MIRRORBALL = 3,
+  PANORAMA_FISHEYE_LENS_POLYNOMIAL = 4,
 
   PANORAMA_NUM_TYPES,
 };
@@ -922,6 +923,8 @@ typedef struct KernelCamera {
   float fisheye_fov;
   float fisheye_lens;
   float4 equirectangular_range;
+  float fisheye_lens_polynomial_bias;
+  float4 fisheye_lens_polynomial_coefficients;
 
   /* stereo */
   float interocular_offset;
diff --git a/intern/cycles/scene/camera.cpp b/intern/cycles/scene/camera.cpp
index 5bafe736fb5..12e106779e3 100644
--- a/intern/cycles/scene/camera.cpp
+++ b/intern/cycles/scene/camera.cpp
@@ -100,6 +100,7 @@ NODE_DEFINE(Camera)
   panorama_type_enum.insert("mirrorball", PANORAMA_MIRRORBALL);
   panorama_type_enum.insert("fisheye_equidistant", PANORAMA_FISHEYE_EQUIDISTANT);
   panorama_type_enum.insert("fisheye_equisolid", PANORAMA_FISHEYE_EQUISOLID);
+  panorama_type_enum.insert("fisheye_lens_polynomial", PANORAMA_FISHEYE_LENS_POLYNOMIAL);
   SOCKET_ENUM(panorama_type, "Panorama Type", panorama_type_enum, PANORAMA_EQUIRECTANGULAR);
 
   SOCKET_FLOAT(fisheye_fov, "Fisheye FOV", M_PI_F);
@@ -112,6 +113,12 @@ NODE_DEFINE(Camera)
   SOCKET_FLOAT(fov_pre, "FOV Pre", M_PI_4_F);
   SOCKET_FLOAT(fov_post, "FOV Post", M_PI_4_F);
 
+  SOCKET_FLOAT(fisheye_polynomial_k0, "Fisheye Polynomial K0", 0.0f);
+  SOCKET_FLOAT(fisheye_polynomial_k1, "Fisheye Polynomial K1", 0.0f);
+  SOCKET_FLOAT(fisheye_polynomial_k2, "Fisheye Polynomial K2", 0.0f);
+  SOCKET_FLOAT(fisheye_polynomial_k3, "Fisheye Polynomial K3", 0.0f);
+  SOCKET_FLOAT(fisheye_polynomial_k4, "Fisheye Polynomial K4", 0.0f);
+
   static NodeEnum stereo_eye_enum;
   stereo_eye_enum.insert("none", STEREO_NONE);
   stereo_eye_enum.insert("left", STEREO_LEFT);
@@ -418,6 +425,9 @@ void Camera::update(Scene *scene)
                                             -longitude_min,
                                             latitude_min - latitude_max,
                                             -latitude_min + M_PI_2_F);
+  kcam->fisheye_lens_polynomial_bias = fisheye_polynomial_k0;
+  kcam->fisheye_lens_polynomial_coefficients = make_float4(fisheye_polynomial_k1, fisheye_polynomial_k2,
+                                                           fisheye_polynomial_k3, fisheye_polynomial_k4);
 
   switch (stereo_eye) {
     case STEREO_LEFT:
diff --git a/intern/cycles/scene/camera.h b/intern/cycles/scene/camera.h
index 58e39599267..50586287bea 100644
--- a/intern/cycles/scene/camera.h
+++ b/intern/cycles/scene/camera.h
@@ -105,6 +105,12 @@ class Camera : public Node {
   NODE_SOCKET_API(float, longitude_min)
   NODE_SOCKET_API(float, longitude_max)
 
+  NODE_SOCKET_API(float, fisheye_polynomial_k0)
+  NODE_SOCKET_API(float, fisheye_polynomial_k1)
+  NODE_SOCKET_API(float, fisheye_polynomial_k2)
+  NODE_SOCKET_API(float, fisheye_polynomial_k3)
+  NODE_SOCKET_API(float, fisheye_polynomial_k4)
+
   /* panorama stereo */
   NODE_SOCKET_API(StereoEye, stereo_eye)
   NODE_SOCKET_API(bool, use_spherical_stereo)
diff --git a/release/scripts/startup/bl_ui/properties_data_camera.py b/release/scripts/startup/bl_ui/properties_data_camera.py
index b56182bb637..edd0623d8fe 100644
--- a/release/scripts/startup/bl_ui/properties_data_camera.py
+++ b/release/scripts/startup/bl_ui/properties_data_camera.py
@@ -110,6 +110,14 @@ class DATA_PT_lens(CameraButtonsPanel, Panel):
                     sub = col.column(align=True)
                     sub.prop(ccam, "longitude_min", text="Longitude Min")
                     sub.prop(ccam, "longitude_max", text="Max")
+                elif ccam.panorama_type == 'FISHEYE_LENS_POLYNOMIAL':
+                    col.prop(ccam, "fisheye_fov")
+                    col.prop(ccam, "fisheye_polynomial_k0", text="K0")
+                    col.prop(ccam, "fisheye_polynomial_k1", text="K1")
+                    col.prop(ccam, "fisheye_polynomial_k2", text="K2")
+                    col.prop(ccam, "fisheye_polynomial_k3", text="K3")
+                    col.prop(ccam, "fisheye_polynomial_k4", text="K4")
+
             elif engine in {'BLENDER_RENDER', 'BLENDER_EEVEE', 'BLENDER_WORKBENCH'}:
                 if cam.lens_unit == 'MILLIMETERS':
                     col.prop(cam, "lens")