Merge branch 'code-update/design-constraints-20260919' into dev

Author: Mihou

deeplens/geolens.py

@@ -10,7 +10,6 @@
     Xinge Yang, Qiang Fu, and Wolfgang Heidrich, "Curriculum learning for ab initio deep learned refractive optics," Nature Communications 2024.
 """
 
-import logging
 import math
 
 import numpy as np
@@ -972,8 +971,10 @@ def calc_fov(self):
             1. **Perspective projection** — from focal length and sensor size
                (effective FoV, ignoring distortion).
             2. **Forward ray tracing** — sweeps FOV angles from object side,
-               traces to sensor, and finds the angle whose centroid image height
-               matches the sensor half-diagonal. This avoids the failure of the
+               traces to sensor, and interpolates the first angle whose centroid
+               reaches the sensor half-diagonal. If a field loses its rays first,
+               uses the preceding surviving field (zero if the on-axis field
+               is unavailable). This avoids the failure of the
                old backward-tracing approach on wide-angle lenses where pupil
                aberration at full field leaves zero valid rays.
 
@@ -1002,33 +1003,46 @@ def calc_fov(self):
         # Sweep FOV angles from object side, trace to sensor, and find which
         # angle produces an image height matching r_sensor.
         num_fov = 64
-        fov_lo = float(np.rad2deg(self.rfov_eff)) * 0.5
         fov_hi = min(float(np.rad2deg(self.rfov_eff)) * 1.8, 89.0)
         fov_samples = torch.linspace(
-            fov_lo, fov_hi, num_fov, device=self.device, dtype=self.dtype
+            0.0, fov_hi, num_fov, device=self.device, dtype=self.dtype
         )
 
         ray = self.sample_from_fov(fov_x=0.0, fov_y=fov_samples.tolist(), num_rays=256)
         ray = self.trace2sensor(ray)
 
         # Centroid image height per FOV angle, shape [num_fov]
         valid = ray.is_valid > 0  # [num_fov, num_rays]
-        masked_y = ray.o[..., 1] * valid
+        masked_y = torch.where(valid, ray.o[..., 1], 0.0)
         n_valid = valid.sum(dim=-1).clamp(min=1)
         imgh = (masked_y.sum(dim=-1) / n_valid).abs()
 
-        # Find the FOV angle whose image height is closest to r_sensor
+        # Stop at the first crossing or lost field; a folded outer image must
+        # not be mistaken for a larger usable field of view.
         has_valid = valid.sum(dim=-1) > 10
         if has_valid.any():
-            imgh[~has_valid] = float("inf")
-            diff = (imgh - self.r_sensor).abs()
-            best_idx = diff.argmin().item()
-            rfov = fov_samples[best_idx].item() * math.pi / 180.0
+            beyond = (~has_valid) | (imgh >= self.r_sensor)
+            if not beyond.any():
+                best_deg = fov_samples[-1]
+            else:
+                idx = int(beyond.int().argmax())
+                if idx == 0:
+                    best_deg = fov_samples[0]
+                elif has_valid[idx]:
+                    fraction = (self.r_sensor - imgh[idx - 1]) / (
+                        imgh[idx] - imgh[idx - 1]
+                    ).clamp_min(EPSILON)
+                    best_deg = fov_samples[idx - 1] + fraction * (
+                        fov_samples[idx] - fov_samples[idx - 1]
+                    )
+                else:
+                    best_deg = fov_samples[idx - 1]
+            rfov = math.radians(float(best_deg))
             self.rfov = rfov
             self.real_dfov = 2 * rfov
         else:
-            self.rfov = self.rfov_eff
-            self.real_dfov = self.dfov
+            self.rfov = 0.0
+            self.real_dfov = 0.0
 
         # 3. Compute 35mm equivalent focal length. 35mm sensor: 36mm * 24mm
         self.eqfl = 21.63 / math.tan(self.rfov_eff)
@@ -1472,32 +1486,48 @@ def refocus(self, foc_dist=float("inf")):
 
     @torch.no_grad()
     def set_fnum(self, fnum):
-        """Set F-number and aperture radius using binary search.
+        """Set F-number by bracketing the pupil radius, then bisecting.
 
         Args:
             fnum (float): target F-number.
+
+        Raises:
+            ValueError: F-number or target pupil radius is not finite and positive.
+            RuntimeError: The pupil target cannot be reached. The original
+                aperture radius is restored on failure.
         """
+        if not math.isfinite(fnum) or fnum <= 0:
+            raise ValueError("F-number must be finite and positive.")
         target_pupil_r = self.foclen / fnum / 2
-        aper_r = self.surfaces[self.aper_idx].r
-        lo, hi = 0.1 * aper_r, 5.0 * aper_r
-
-        pupilr = None
-        for _ in range(40):
-            mid = 0.5 * (lo + hi)
-            self.surfaces[self.aper_idx].update_r(float(mid))
-            _, pupilr = self.calc_entrance_pupil_rayaiming()
-            if abs(pupilr - target_pupil_r) / target_pupil_r < 1e-3:
-                break
-            if pupilr > target_pupil_r:
-                hi = mid
-            else:
-                lo = mid
-        else:
-            logging.warning(
-                f"set_fnum: did not converge, pupil_r={pupilr:.4f}, target={target_pupil_r:.4f}"
+        if not math.isfinite(target_pupil_r) or target_pupil_r <= 0:
+            raise ValueError("Target pupil radius must be finite and positive.")
+        aperture = self.surfaces[self.aper_idx]
+        original_r = aperture.r
+        lo, hi = 0.0, None
+        radius = original_r
+        try:
+            for _ in range(80):
+                aperture.update_r(float(radius))
+                radius = aperture.r  # update_r enforces the surface's physical limit.
+                _, pupilr = self.calc_entrance_pupil_rayaiming()
+                if not math.isfinite(pupilr) or pupilr <= 0:
+                    break
+                if abs(pupilr - target_pupil_r) / target_pupil_r < 1e-3:
+                    self.calc_pupil()
+                    return
+                if pupilr > target_pupil_r:
+                    hi = radius
+                else:
+                    if radius <= lo:
+                        break
+                    lo = radius
+                radius = 2 * radius if hi is None else 0.5 * (lo + hi)
+            raise RuntimeError(
+                f"Cannot reach F/{fnum}: target pupil radius {target_pupil_r:.6g} mm."
             )
-
-        self.calc_pupil()
+        except Exception:
+            aperture.update_r(original_r)
+            raise
 
     @torch.no_grad()
     def set_target_fov_fnum(self, rfov, fnum):
@@ -1521,12 +1551,8 @@ def set_target_fov_fnum(self, rfov, fnum):
         self.real_dfov = 2 * self.rfov
         self.foclen = self.r_sensor / math.tan(self.rfov_eff)
         self.eqfl = 21.63 / math.tan(self.rfov_eff)
-        self.fnum = fnum
-        aper_r = self.foclen / fnum / 2
-        self.surfaces[self.aper_idx].update_r(float(aper_r))
-
-        # Update pupil after setting aperture radius
-        self.calc_pupil()
+        # F-number specifies the entrance pupil, not an internal stop's radius.
+        self.set_fnum(fnum)
 
     @torch.no_grad()
     def set_fov(self, rfov):

deeplens/geolens_pkg/eval.py

@@ -72,8 +72,8 @@
             field angles (infinite conjugates), or a traced bundle.
 
     Comprehensive Analysis:
-        analysis_spot: RMS and geometric spot radii averaged over RGB at
-            multiple field positions.
+        analysis_spot: RMS and geometric spot radii and ray throughput over RGB
+            at multiple field positions.
         analysis_rendering: Render a test image through the lens and report
             PSNR / SSIM.
         analysis: One-call entry point that chains layout drawing, spot
@@ -1549,7 +1549,7 @@ def analysis_rendering(
 
     @torch.no_grad()
     def analysis_spot(self, num_field=3, depth=float("inf")):
-        """Compute RMS and geometric spot radii at multiple field positions for RGB.
+        """Compute RGB spot radii and ray throughput at multiple field positions.
 
         Traces rays at ``num_field`` evenly-spaced field positions along the
         meridional direction for three wavelengths (R, G, B), and computes
@@ -1581,6 +1581,8 @@ def analysis_spot(self, num_field=3, depth=float("inf")):
                 ``'fov0.5'``, ``'fov1.0'``). Each value is a dict with:
                     - ``'rms'``: Polychromatic RMS spot radius in μm.
                     - ``'radius'``: Polychromatic geometric spot radius in μm.
+                    - ``'throughput'``: Fraction of sampled rays reaching the sensor;
+                      zero marks a lost field even when its spot radius is zero.
         """
         # Trace each wavelength and pool rays across wavelengths per field
         xy_list = []
@@ -1624,6 +1626,7 @@ def analysis_spot(self, num_field=3, depth=float("inf")):
         # Convert mm → μm
         avg_rms_radius_um = spot_rms * 1000.0
         avg_geo_radius_um = spot_radius * 1000.0
+        throughput = (valid_all > 0).float().mean(-1)
 
         # Print results
         print(f"Ray spot analysis results for depth {depth}:")
@@ -1633,6 +1636,9 @@ def analysis_spot(self, num_field=3, depth=float("inf")):
         print(
             f"Geo radius: FoV (0.0) {avg_geo_radius_um[0]:.3f} um, FoV (0.5) {avg_geo_radius_um[num_field // 2]:.3f} um, FoV (1.0) {avg_geo_radius_um[-1]:.3f} um"
         )
+        print(
+            f"Throughput: FoV (0.0) {throughput[0]:.3f}, FoV (0.5) {throughput[num_field // 2]:.3f}, FoV (1.0) {throughput[-1]:.3f}"
+        )
 
         # Save to dict
         rms_results = {}
@@ -1642,6 +1648,7 @@ def analysis_spot(self, num_field=3, depth=float("inf")):
             rms_results[f"fov{fov}"] = {
                 "rms": round(avg_rms_radius_um[i].item(), 4),
                 "radius": round(avg_geo_radius_um[i].item(), 4),
+                "throughput": round(throughput[i].item(), 4),
             }
 
         return rms_results

deeplens/geolens_pkg/optim.py

@@ -331,9 +331,10 @@ def loss_profile(self):
     def loss_bound(self):
         """Penalize geometry-bound violations in a single surface-sampling pass.
 
-        Each surface pair is sampled once and its distances feed both the
-        clearance (min) and envelope (max) relu penalties for air gaps,
-        glass thickness, BFL, and TTL.
+        Each surface pair is sampled over its shared physical radius and at
+        each surface's own rim. These distances feed both the clearance (min)
+        and envelope (max) relu penalties for air gaps, glass thickness, BFL,
+        and TTL, without extrapolating beyond either surface.
 
         Returns:
             loss_clearance (torch.Tensor): Scalar clearance penalty for parts
@@ -380,7 +381,9 @@ def loss_bound(self):
                 r_center, 0.0, valid_check=False, d=self.surf_d(i + 1)
             )
 
-            r_edge = torch.linspace(0.5, 1.0, 16, device=self.device) * current_surf.r
+            r_edge = torch.linspace(0.0, 1.0, 16, device=self.device) * min(
+                current_surf.r, next_surf.r
+            )
             z_prev_edge = current_surf.surface_with_offset(
                 r_edge, 0.0, valid_check=False, d=self.surf_d(i)
             )
@@ -390,6 +393,19 @@ def loss_bound(self):
 
             dist_center = z_next_center - z_prev_center
             dist_edges = z_next_edge - z_prev_edge
+            # Keep the same physical-rim sag bracket used by prune_surf.
+            rim_gap = next_surf.surface_with_offset(
+                torch.as_tensor(next_surf.r, device=self.device),
+                0.0,
+                valid_check=False,
+                d=self.surf_d(i + 1),
+            ) - current_surf.surface_with_offset(
+                torch.as_tensor(current_surf.r, device=self.device),
+                0.0,
+                valid_check=False,
+                d=self.surf_d(i),
+            )
+            dist_edges = torch.cat((dist_edges, rim_gap.reshape(1)))
             dist_edge_lo = torch.min(dist_edges)
             dist_edge_hi = torch.max(dist_edges)