make particle ppp generation more robust

synpivimage/core.py

@@ -48,6 +48,13 @@ def take_image(laser: Laser,
         The peak count of the particles
     kwargs : dict
         Additional parameters
+
+    Returns
+    -------
+    img: np.ndarray
+        The created image
+    particles:
+        The updated Particles object
     """
     logger = kwargs.get('logger', LOGGER)
     # compute the particle intensity factor in order to reach particle_peak_count

synpivimage/io.py

@@ -366,12 +366,16 @@ def __exit__(self, exc_type, exc_val, exc_tb):
         img_idx_ds.make_scale()
 
         ds.dims[0].attach_scale(img_idx_ds)
-        nx_ds = ds.parent.create_dataset('nx', data=np.arange(self.camera.nx))
+        nx_ds = ds.parent.create_dataset('ix', data=np.arange(self.camera.nx))
+        ds.attrs['long_name'] = 'x pixel coordinate'
+        ds.attrs['standard_name'] = 'x_pixel_coordinate'  # https://matthiasprobst.github.io/pivmeta/#x_pixel_coordinate
         nx_ds.make_scale()
         ds.dims[2].attach_scale(nx_ds)
 
-        ny_ds = ds.parent.create_dataset('ny', data=np.arange(self.camera.ny))
+        ny_ds = ds.parent.create_dataset('iy', data=np.arange(self.camera.ny))
         ny_ds.make_scale()
+        ds.attrs['long_name'] = 'y pixel coordinate'
+        ds.attrs['standard_name'] = 'y_pixel_coordinate'  # https://matthiasprobst.github.io/pivmeta/#x_pixel_coordinate
         ds.dims[1].attach_scale(ny_ds)
 
         for cmp_name, cmp, onto_model in zip(('laser', 'camera'),

synpivimage/laser.py

@@ -77,7 +77,7 @@ def illuminate(self,
         Returns
         -------
         Particles
-            The illuminated particles (new object!)
+            The illuminated particles (same object!)
         """
         logger = kwargs.get('logger', LOGGER)
 
@@ -92,6 +92,8 @@ def illuminate(self,
             laser_intensity = tophat(self.width)
         else:
             laser_intensity = real(dz0, s)
+
+        particles.reset()
         particles.source_intensity = laser_intensity(particles.z)
 
         inside_laser = particles.source_intensity > np.exp(-2)
@@ -113,7 +115,9 @@ def illuminate(self,
             plt.ylabel('Normalized particle intensity in beam / -')
             plt.grid()
             plt.show()
-        return Particles(**particles.dict())
+
+        return particles
+        # return Particles(**particles.dict())
 
     def model_dump_jsonld(self) -> str:
         """Return JSON-LD string representation"""

synpivimage/particles.py

@@ -116,6 +116,9 @@ def _parse_array(_arr, _n: int, dtype=None):
         self._image_electrons = _parse_array(image_electrons, N)
         self._image_quantized_electrons = _parse_array(image_quantized_electrons, N)
         self._flag = _parse_array(flag, N, dtype=int)
+        self._xlim = None
+        self._ylim = None
+        self._zlim = None
 
     @property
     def x(self):
@@ -212,32 +215,44 @@ def reset(self):
         self._flag = np.zeros_like(self.x, dtype=int)
 
     @classmethod
-    def generate(cls, ppp: float, dx_max, dy_max, dz_max, camera, laser) -> "Particles":
+    def generate(cls, ppp: float, dx_max, dy_max, dz_max, size, camera, laser) -> "Particles":
         """Generate particles based on a certain ppp (particles per pixel). With
         dx, dy, dz the maximum displacement of the particles can be set. The camera and laser
         are used to determine the sensor size and the laser width."""
         from .utils import generate_particles
-        return generate_particles(ppp=ppp, dx_max=dx_max, dy_max=dy_max, dz_max=dz_max, camera=camera, laser=laser)
+        return generate_particles(ppp=ppp, dx_max=dx_max, dy_max=dy_max, dz_max=dz_max, size=size, camera=camera, laser=laser)
 
     def get_ppp(self, camera_size: int) -> float:
         """Return the particles per pixel"""
         return self.active.sum() / camera_size
 
     def regenerate(self) -> "Particles":
-        """Regenerate the particles. Does return a new object and does not touch the current object!
+        """Regenerate the particles of this object.
+        The locations (x, y, z) of the particles will change and the intensities
+        will be reset
 
-        The number of particles and hence the ppp will remain constant (within a statistical deviation)
-        but the position and size of the particles will change."""
-        new_part = self.copy()
+        .. note::
+
+            Does NOT create a new object!
+            The returned object is the same
+        """
+        self.reset()
         N = len(self.x)
-        x = np.random.uniform(min(self.x), max(self.x), N)
-        y = np.random.uniform(min(self.y), max(self.y), N)
-        z = np.random.uniform(min(self.z), max(self.z), N)
-        new_part._x = x
-        new_part._y = y
-        new_part._z = z
-        new_part.reset()
-        return new_part
+        if self._xlim is None:
+            self._x = np.random.uniform(min(self.x), max(self.x), N)
+        else:
+            self._x = np.random.uniform(*self._xlim, N)
+
+        if self._ylim is None:
+            self._y = np.random.uniform(min(self.y), max(self.y), N)
+        else:
+            self._y = np.random.uniform(*self._ylim, N)
+
+        if self._ylim is None:
+            self._z = np.random.uniform(min(self.z), max(self.z), N)
+        else:
+            self._z = np.random.uniform(*self._zlim, N)
+        return self
 
     def __len__(self):
         return self.x.size
@@ -508,7 +523,7 @@ def compute_intensity_distribution(
     xp : float
         x-coordinate of the particles on the sensor in pixels
     yp : float
-        y-coordinate of the particles on thesensor in pixels
+        y-coordinate of the particles on the sensor in pixels
     dp : float
         particle image diameter (in pixels)
     sigmax : float

synpivimage/utils.py

@@ -1,4 +1,5 @@
 import logging
+import warnings
 from typing import Tuple
 
 import numpy as np
@@ -16,9 +17,10 @@ def generate_particles(ppp: float,
                        dx_max: Tuple[float, float],
                        dy_max: Tuple[float, float],
                        dz_max: Tuple[float, float],
+                       size: float,
                        camera: Camera,
                        laser: Laser,
-                       iter_max: int = 40):
+                       **kwargs) -> Particles:
     """Generates a particle class based on the current setup and given max displacements
 
     Parameters
@@ -31,7 +33,27 @@ def generate_particles(ppp: float,
         Maximal displacement in y-direction [lower, upper]
     dz_max: Tuple[float, float]
         Maximal displacement in z-direction [lower, upper]
+    camera: Camera
+        camera
+    laser: Laser
+        laser
+    size: float
+        Particle size
+    kwargs:
+        iter_max: int=40
+            Max. iteration for particle number determination algorithm
+        N_max: int=10**7
+            Max. number of particles
+
+    Returns
+    -------
+    particles: Particles
+        The generated particles
     """
+
+    iter_max = kwargs.get('iter_max', 40)
+    N_max = kwargs.get('N_ma', 10 ** 7)
+
     logger.debug(f'Generating particles with a ppp of {ppp}')
     assert 0 < ppp < 1, f"Expected ppp to be between 0 and 1, got {ppp}"
 
@@ -43,49 +65,62 @@ def generate_particles(ppp: float,
         zmax = laser.width
     area = (camera.nx + (dx_max[1] - dx_max[0])) * (camera.ny + (dy_max[1] - dy_max[0]))
     N = int(area * ppp)
+    if N < 2:
+        raise ValueError(
+            'Number of particles is too low. For generation of only one particle,'
+            ' please use discrete generation by providing the coordinates'
+        )
+    logger.debug(f'Initial particle count: {N}')
     # Ntarget = ppp * camera.size
 
     curr_ppp = 0
     # rel_dev = abs((curr_ppp - ppp) / ppp)
 
+    _xlim = min(-dx_max[1], 0), max(camera.nx, camera.nx - dx_max[0])
+    _ylim = min(-dy_max[1], 0), max(camera.ny, camera.ny - dy_max[0])
+    _zlim = min(zmin, zmin - dz_max[1]), max(zmax, zmax + dz_max[0])
+
+    x_extent = _xlim[1] - _xlim[0]
+    y_extent = _ylim[1] - _ylim[0]
+
+    n_too_much_noise = 0
     i = 0
     while abs((curr_ppp - ppp) / ppp) > 0.01:
         i += 1
         if i == 1:
             # generate a initial random distribution. From this. particles are either added or removed
             logger.debug(f'Generate initial random distribution of {N} particles')
-            xe = np.random.uniform(min(-dx_max[1], 0), max(camera.nx, camera.nx - dx_max[0]), N)
-            ye = np.random.uniform(min(-dy_max[1], 0), max(camera.ny, camera.ny - dy_max[0]), N)
-            ze = np.random.uniform(min(zmin, zmin - dz_max[1]), max(zmax, zmax + dz_max[0]), N)
+            xe = np.random.uniform(*_xlim, N)
+            ye = np.random.uniform(*_ylim, N)
+            ze = np.random.uniform(*_zlim, N)
 
         particles = Particles(
             x=xe,
             y=ye,
             z=ze,
-            size=np.ones_like(xe) * 2
+            size=size * np.ones_like(xe)
         )
         _img, _part = take_image(particles=particles,
                                  camera=camera,
                                  laser=laser,
                                  particle_peak_count=1000)
-        curr_ppp = _part.active.sum() / camera.size
-        # print(f'curr ppp: {curr_ppp:.5f}')
+        curr_ppp = _part.get_ppp(camera.size)  # _part.active.sum() / camera.size
+
         diff_ppp = ppp - curr_ppp
-        # print(f'diff ppp: {diff_ppp:.5f}')
-        Nadd = int(diff_ppp * camera.size)  # dampen the addition by 10%
+
+        Nadd = int(diff_ppp * x_extent * y_extent)
 
         if Nadd == 0:
             logger.debug('Stopping early because no new particles to be added.')
             break
 
         logger.debug(f'Generate particles. Iteration {i}/{iter_max}: curr ppp: {curr_ppp:.5f}. '
                      f'diff ppp: {diff_ppp:.5f}. Adding {Nadd} particles')
-
         if Nadd > 0:
             # generate new particles within the given range
-            xe_new = np.random.uniform(min(-dx_max[1], 0), max(camera.nx, camera.nx - dx_max[0]), Nadd)
-            ye_new = np.random.uniform(min(-dy_max[1], 0), max(camera.ny, camera.ny - dy_max[0]), Nadd)
-            ze_new = np.random.uniform(min(zmin, zmin - dz_max[1]), max(zmax, zmax + dz_max[0]), Nadd)
+            xe_new = np.random.uniform(*_xlim, Nadd)
+            ye_new = np.random.uniform(*_ylim, Nadd)
+            ze_new = np.random.uniform(*_zlim, Nadd)
             # add them to the existing particles
             _n_new = len(xe) + Nadd
             xe = np.concatenate([xe, xe_new])
@@ -101,6 +136,18 @@ def generate_particles(ppp: float,
             assert len(xe) == N + Nadd
 
         N += Nadd
+
+        if N > N_max:
+            raise ValueError(f'Number of particles exceeded maximum of {N_max}. '
+                             f'Consider increasing the number of iterations or the particle size')
+
+        if curr_ppp == 0 and _part.in_fov.sum() > 0:
+            n_too_much_noise += 1
+            warnings.warn('No particles in the field of view are illuminated beyond the laser noise. Consider using a'
+                          'smaller noise level or a larger particle size.')
+            if n_too_much_noise > 10:
+                raise ValueError('Too many iterations without any particles in the field of view. Stopping.')
+
         err = abs((curr_ppp - ppp) / ppp)
 
         if err < 0.01:
@@ -110,4 +157,7 @@ def generate_particles(ppp: float,
             logger.debug(f'Reached max iteration of {iter_max}')
             break
 
+    _part._xlim = _xlim
+    _part._ylim = _ylim
+    _part._zlim = _zlim
     return _part