Repository: magpylib/magpylib (language: python)

Issue:
Add sensor pixel aggregation
# Description
It would be a nice feature to be able to aggregate sensor pixel data directly from the `getB`/`getH` functions, and compute the `min`, `max` or `mean` values over the pixel axes.

Also, the aggregation could make it acceptable to have sensors with different pixel shapes to be evaluated, as it would be compatible with a `ndarray` after aggregation. 

# Expected behavior / Example

```python
iimport magpylib as magpy
import numpy as np
src1 = magpy.magnet.Cuboid((0,0,1000),(1,1,1)).move([[1,0,0]])
sens1 = magpy.Sensor(position=(0,0,1), pixel=np.zeros((4,5,3)))
sens2 = sens1.copy(position=(0,0,2))
sens3 = sens1.copy(position=(0,0,3))

sources = src1,
sensors = sens1, sens2, sens3
B = magpy.getB(sources, sensors, squeeze=False, pixel_agg=None)
print('expected: (1, 2, 3, 4, 5, 3)', 'out: ', B.shape)
B = magpy.getB(sources, sensors, squeeze=False, pixel_agg='mean')
print('expected: (1, 2, 3, 1, 3)', 'out: ', B.shape)
B = magpy.getB(sources, sensors, squeeze=True, pixel_agg=None)
print('expected: (2, 3, 4, 5, 3)', 'out: ', B.shape)
B = magpy.getB(sources, sensors, squeeze=True, pixel_agg='mean')
print('expected: (2, 3, 3)', 'out: ', B.shape)
```

# Possible Implementation
In  `getBH`/`getBH_level3` as  `pixel_agg` named keyword argument:
```python
"""
pixel_agg: str, default=`None`
        A compatible numpy aggregator string (e.g. `'min', 'max', 'mean'`) which applies on pixel
        output values.
"""
```

```python
 # sumup over sources
    if kwargs['sumup']:
        B = np.sum(B, axis=0, keepdims=True)

    # aggregate pixel values
    pixel_agg = kwargs['pixel_agg']
    if pixel_agg is not None:
        B = getattr(np, pixel_agg)(B, axis=tuple(range(3-B.ndim,-1)))
        if not kwargs['squeeze']:
            # add missing dimension since `pixel_agg` reduces pixel
            # dimensions to zero. Only needed if `squeeze is False``
            B = np.expand_dims(B, axis=-2)

    # reduce all size-1 levels
    if kwargs['squeeze']:
        B = np.squeeze(B)
```
Nested collections
# Related Issues
#488 

# Notes

This PR implements a new **nested collections** philosophy. 

### BREAKING CHANGES
- A child has a unique parent, and cannot be in multiple collections anymore.
- Adding up objects creates a new Collection -> new x1 + x2 +x3 behavior -> new c += x behavior
- No more Nesting lists/tuples as collection input (single flat list input still allowed)
- remove the `__sub__` method

### Base behavior
A Collection spans a local reference frame and holds object references to that frame. 
- [x] Children can be sensors, sources or other collections. Setters and getters apply only for the base tree level and not the nested one.
- [x] Moving/Rotating the Collection applies to the frame with all objects inside. Relative positions add orientations of nested object is kept, including paths. 
- [x] Duplicate objects are not allowed. **A child as a unique parent**.

### `__add__`  behavior
- [x] Simple rule: **addition always results in super collection for all type of objects**:
 - Example 1: `obj1 + obj2` -> `Collection(obj, obj2)` 
 - Example 2: `obj1 + obj2 + obj3` -> `Collection(Collection(obj, obj2), obj3)` 


### `getBH` behavior:
- [x] `getBH` of a collection gives the total field generated by this collection. It sums up all the sources found in the nested structure.
- [x] `getBH` of a sensor collection gives the field at every sensor found in the flatten structure

### Other features
- [x] `set_children_styles` walks down the nested tree (has recursive argument)
- [x] `remove(obj)` looks for the object in the nested structure (has recursive argument)
- [x] Add a `describe(...)` method for all objects and tree view for nested collections 

### Other tasks
- Adapt plotting backend for nested structure
  - [x] plotly
  - [x] matplotlib
- [x] Check consistency when building collections with `__add__`,  `__radd__` and `__sub__`
- Testing  @OrtnerMichael
  - [x] Fix current tests
  - [x] new behavior tests
- [x] Check if nested behavior is an acceptable overhead @OrtnerMichael 
- [x] Update docstrings
- [x] Go back to simple `__repr__` -> `Collection(id=..., label=...)`
- [x] Update changelog (specify breaking changes) @OrtnerMichael
- [x] Fix docs @OrtnerMichael

### Examples

```python
import magpylib as magpy
import numpy as np

magpy.defaults.display.style.magnet.magnetization.show = False

s1 = magpy.magnet.Cuboid((0, 0, 1000), (1, 1, 1), (0,0,0), style_label="cuboid1", style_color=None)
s2 = magpy.magnet.Cylinder((0, 0, 1000), (1, 1), (2,0,0), style_label="cylinder1")
s3 = magpy.magnet.Sphere((0, 0, 1000), 1, (4,0,0), style_label="sphere1")
sens1 = magpy.Sensor((1,0,2),style_label="sensor1")
sens2 = magpy.Sensor((3,0,2),style_label="sensor2")

src_col = magpy.Collection(s1, s2, style_label="src_col", style_color='orange')
sens_col = magpy.Collection(sens1, style_label="sens_col")
mixed_col = magpy.Collection(s3, sens2, style_label="mixed_col")
nested_col = magpy.Collection(src_col, sens_col, mixed_col, style_label="nested_col")
magpy.show(nested_col, backend='plotly')
nested_col.describe(properties=True)
```
Nested Collection tree:
```
Collection nested_col (id=140559898024592)
│   • position: [0. 0. 0.] mm
│   • orientation: [0. 0. 0.] degrees
├── Collection src_col (id=140559896365616)
│   │   • position: [0. 0. 0.] mm
│   │   • orientation: [0. 0. 0.] degrees
│   ├── Cuboid cuboid1 (id=140559896366672)
│   │       • position: [0. 0. 0.] mm
│   │       • orientation: [0. 0. 0.] degrees
│   │       • dimension: [1. 1. 1.] mm
│   │       • magnetization: [   0.    0. 1000.] mT
│   └── Cylinder cylinder1 (id=140559896365472)
│           • position: [2. 0. 0.] mm
│           • orientation: [0. 0. 0.] degrees
│           • dimension: [1. 1.] mm
│           • magnetization: [   0.    0. 1000.] mT
├── Collection sens_col (id=140559893170640)
│   │   • position: [0. 0. 0.] mm
│   │   • orientation: [0. 0. 0.] degrees
│   └── Sensor sensor1 (id=140559896364656)
│           • position: [1. 0. 2.] mm
│           • orientation: [0. 0. 0.] degrees
│           • pixel: [0. 0. 0.]
└── Collection mixed_col (id=140559898031056)
    │   • position: [0. 0. 0.] mm
    │   • orientation: [0. 0. 0.] degrees
    ├── Sphere sphere1 (id=140559896366192)
    │       • position: [4. 0. 0.] mm
    │       • orientation: [0. 0. 0.] degrees
    │       • diameter: 1.0 mm
    │       • magnetization: [   0.    0. 1000.] mT
    └── Sensor sensor2 (id=140559896367056)
            • position: [3. 0. 2.] mm
            • orientation: [0. 0. 0.] degrees
            • pixel: [0. 0. 0.]
```

<img width="313" alt="image" src="https://user-images.githubusercontent.com/29252289/160291550-32aef244-f7b6-4171-aa92-9d259ec90988.png">


Sensor rotation path bad overhead
```python
import magpylib as magpy
from time import perf_counter as pf
import numpy as np

NN = 100000

# case static sensor norot
src = magpy.magnet.Cuboid((1000,0,0),(1,1,1),(0,0,2))
sens = magpy.Sensor()
src.rotate_from_angax(np.linspace(0,360,NN),'z', start=0)

T0 = pf()
magpy.getB(src,sens)
T1 = pf()
print(T1-T0)


# case static sensor rot
src = magpy.magnet.Cuboid((1000,0,0),(1,1,1),(0,0,2))
sens = magpy.Sensor()
src.rotate_from_angax(np.linspace(0,360,NN),'z', start=0)
sens.rotate_from_angax(45,'z')

T0 = pf()
magpy.getB(src,sens)
T1 = pf()
print(T1-T0)


# case sensor rot path
src = magpy.magnet.Cuboid((1000,0,0),(1,1,1),(0,0,2))
sens = magpy.Sensor()
sens.rotate_from_angax(np.linspace(0,360,NN),'z', start=0)


T0 = pf()
magpy.getB(src,sens)
T1 = pf()
print(T1-T0)
```

The loop cycling through all path positions in a sensor-rotation path leads to an unexpected overhead !!!!

results: 
2.0304203999999997
2.0631880999999996
5.2061805


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 63b7cd6d4..ef3d52147 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -35,6 +35,7 @@ This is a major update that includes
   - Improved performance of `getB` and `getH`.
   - Fixed array dimension wrongly reduced when `sumup=True` and `squeeze=False` ind `getB` and `getH` functions ([#425](https://github.com/magpylib/magpylib/issues/425), [#426](https://github.com/magpylib/magpylib/pull/426))
   - Minimal non-squeeze output shape is (1,1,1,1,3), meaning that a single pixel is now also represented. ([#493](https://github.com/magpylib/magpylib/pull/493))
+- With the new kwarg `pixel_agg` it is now possible to apply a numpy function with reducing functionality (like `mean`, `min`, `average`) to the pixel output. In this case, it is allowed to provide `getB` and `getH` with different observer input shapes. ([#503](https://github.com/magpylib/magpylib/pull/503))
 
 ### Major graphic output overhaul:
 - Styles:
@@ -74,7 +75,7 @@ This is a major update that includes
 - Removed `increment` argument from `move` and `rotate` functions ([#438](https://github.com/magpylib/magpylib/discussions/438), [#444](https://github.com/magpylib/magpylib/issues/444))
 
 ### Modifications to the `Collection` class
-- Collections can now contain source, `Sensor` and other `Collection` objects and can function as source and observer inputs in `getB` and `getH`. ([#410](https://github.com/magpylib/magpylib/issues/410), [#415](https://github.com/magpylib/magpylib/pull/415), [#297](https://github.com/magpylib/magpylib/issues/297))
+- Collections can now contain source, `Sensor` and other `Collection` objects and can function as source and observer inputs in `getB` and `getH`. ([#502](https://github.com/magpylib/magpylib/pull/502), [#410](https://github.com/magpylib/magpylib/issues/410), [#415](https://github.com/magpylib/magpylib/pull/415), [#297](https://github.com/magpylib/magpylib/issues/297))
 - Instead of the property `Collection.sources` there are now the `Collection.children`, `Collection.sources`, `Collection.sensors` and `Collection.collections` properties. Setting these collection properties will automatically override parents. ([#446](https://github.com/magpylib/magpylib/issues/446), [#502](https://github.com/magpylib/magpylib/pull/502))
 - `Collection` has it's own `position`, `orientation` and `style`. ([#444](https://github.com/magpylib/magpylib/issues/444), [#461](https://github.com/magpylib/magpylib/issues/461))
 - All methods applied to a collection maintain relative child-positions in the local reference frame.
diff --git a/magpylib/_src/fields/field_wrap_BH_level2.py b/magpylib/_src/fields/field_wrap_BH_level2.py
index 7ad13256f..c24473a6c 100644
--- a/magpylib/_src/fields/field_wrap_BH_level2.py
+++ b/magpylib/_src/fields/field_wrap_BH_level2.py
@@ -12,68 +12,22 @@
     check_excitations,
     check_dimensions,
     check_format_input_observers,
+    check_format_pixel_agg,
 )
 
 
-def tile_mag(group: list, n_pp: int):
-    """ tile up magnetizations of shape (3,)
-    """
-    mags = np.array([src.magnetization for src in group])
-    magv = np.tile(mags, n_pp).reshape((-1, 3))
-    return magv
-
-
-def tile_dim_cuboid(group: list, n_pp: int):
-    """ tile up cuboid dimension
-    """
-    dims = np.array([src.dimension for src in group])
-    dimv = np.tile(dims, n_pp).reshape((-1, 3))
-    return dimv
-
-
-def tile_dim_cylinder(group: list, n_pp: int):
-    """ tile up cylinder dimensions.
-    """
-    dims = np.array([src.dimension for src in group])
-    dimv = np.tile(dims, n_pp).reshape((-1, 2))
-    return dimv
-
-
-def tile_dim_cylinder_segment(group: list, n_pp: int):
-    """ tile up cylinder segment dimensions.
-    """
-    dims = np.array([src.dimension for src in group])
-    dimv = np.tile(dims, n_pp).reshape((-1, 5))
-    return dimv
-
-
-def tile_dia(group: list, n_pp: int):
-    """ tile up diameter
-    """
-    dims = np.array([src.diameter for src in group])
-    dimv = np.tile(dims, n_pp).flatten()
-    return dimv
-
-
-def tile_moment(group: list, n_pp: int):
-    """ tile up moments of shape (3,)
-    """
-    moms = np.array([src.moment for src in group])
-    momv = np.tile(moms, n_pp).reshape((-1, 3))
-    return momv
-
-
-def tile_current(group: list, n_pp: int):
-    """ tile up current inputs
-    """
-    currs = np.array([src.current for src in group])
-    currv = np.tile(currs, n_pp).flatten()
-    return currv
+def tile_group_property(group: list, n_pp: int, prop_name: str):
+    """tile up group property"""
+    prop0 = getattr(group[0], prop_name)
+    out = np.array([getattr(src, prop_name) for src in group])
+    out = np.tile(out, n_pp)
+    if np.isscalar(prop0):
+        return out.flatten()
+    return out.reshape((-1, prop0.shape[0]))
 
 
 def get_src_dict(group: list, n_pix: int, n_pp: int, poso: np.ndarray) -> dict:
-    """ create dictionaries for level1 input
-    """
+    """create dictionaries for level1 input"""
     # pylint: disable=protected-access
     # pylint: disable=too-many-return-statements
 
@@ -88,56 +42,54 @@ def get_src_dict(group: list, n_pix: int, n_pp: int, poso: np.ndarray) -> dict:
     rotobj = R.from_quat(rotv)
 
     # pos_obs
-    posov = np.tile(poso, (len(group),1))
+    posov = np.tile(poso, (len(group), 1))
 
     # determine which group we are dealing with and tile up dim and excitation
     src_type = group[0]._object_type
 
-    kwargs = {'source_type': src_type, 'position': posv, 'observer': posov, 'orientation': rotobj}
-
-    if src_type == 'Sphere':
-        magv = tile_mag(group, n_pp)
-        diav = tile_dia(group, n_pp)
-        kwargs.update({'magnetization':magv, 'diameter':diav})
-
-    elif src_type == 'Cuboid':
-        magv = tile_mag(group, n_pp)
-        dimv = tile_dim_cuboid(group, n_pp)
-        kwargs.update({'magnetization':magv,  'dimension':dimv})
-
-    elif src_type == 'Cylinder':
-        magv = tile_mag(group, n_pp)
-        dimv = tile_dim_cylinder(group, n_pp)
-        kwargs.update({'magnetization':magv,  'dimension':dimv})
-
-    elif src_type == 'CylinderSegment':
-        magv = tile_mag(group, n_pp)
-        dimv = tile_dim_cylinder_segment(group, n_pp)
-        kwargs.update({'magnetization':magv,  'dimension':dimv})
+    kwargs = {
+        "source_type": src_type,
+        "position": posv,
+        "observer": posov,
+        "orientation": rotobj,
+    }
+
+    if src_type in ("Sphere", "Cuboid", "Cylinder", "CylinderSegment"):
+        magv = tile_group_property(group, n_pp, "magnetization")
+        kwargs.update(magnetization=magv)
+        if src_type == "Sphere":
+            diav = tile_group_property(group, n_pp, "diameter")
+            kwargs.update(diameter=diav)
+        else:
+            dimv = tile_group_property(group, n_pp, "dimension")
+            kwargs.update(dimension=dimv)
 
-    elif src_type == 'Dipole':
-        momv = tile_moment(group, n_pp)
-        kwargs.update({'moment':momv})
+    elif src_type == "Dipole":
+        momv = tile_group_property(group, n_pp, "moment")
+        kwargs.update({"moment": momv})
 
-    elif src_type == 'Loop':
-        currv = tile_current(group, n_pp)
-        diav = tile_dia(group, n_pp)
-        kwargs.update({'current':currv, 'diameter':diav})
+    elif src_type == "Loop":
+        currv = tile_group_property(group, n_pp, "current")
+        diav = tile_group_property(group, n_pp, "diameter")
+        kwargs.update({"current": currv, "diameter": diav})
 
-    elif src_type == 'Line':
+    elif src_type == "Line":
         # get_BH_line_from_vert function tiles internally !
-        #currv = tile_current(group, n_pp)
+        # currv = tile_current(group, n_pp)
         currv = np.array([src.current for src in group])
         vert_list = [src.vertices for src in group]
-        kwargs.update({'current':currv, 'vertices':vert_list})
-
-    elif src_type == 'CustomSource':
-        kwargs.update({
-            'field_B_lambda': group[0].field_B_lambda,
-            'field_H_lambda': group[0].field_H_lambda})
+        kwargs.update({"current": currv, "vertices": vert_list})
+
+    elif src_type == "CustomSource":
+        kwargs.update(
+            {
+                "field_B_lambda": group[0].field_B_lambda,
+                "field_H_lambda": group[0].field_H_lambda,
+            }
+        )
 
     else:
-        raise MagpylibInternalError('Bad source_type in get_src_dict')
+        raise MagpylibInternalError("Bad source_type in get_src_dict")
 
     return kwargs
 
@@ -156,6 +108,8 @@ def getBH_level2(sources, observers, **kwargs) -> np.ndarray:
     - 'sumup' (bool): False returns [B1,B2,...] for every source, True returns sum(Bi)
         for all sources.
     - 'squeeze' (bool): True output is squeezed (axes of length 1 are eliminated)
+    - 'pixel_agg' : str: A compatible numpy aggregator string (e.g. `'min', 'max', 'mean'`)
+       which applies on pixel output values.
     - 'field' (str): 'B' computes B field, 'H' computes H-field
     - getBH_dict inputs
 
@@ -179,16 +133,12 @@ def getBH_level2(sources, observers, **kwargs) -> np.ndarray:
 
     # CHECK AND FORMAT INPUT ---------------------------------------------------
     if isinstance(sources, str):
-        return getBH_dict_level2(
-            source_type=sources,
-            observer=observers,
-            **kwargs
-        )
+        return getBH_dict_level2(source_type=sources, observer=observers, **kwargs)
 
     # bad user inputs mixing getBH_dict kwargs with object oriented interface
     kwargs_check = kwargs.copy()
-    for popit in ['field', 'sumup', 'squeeze']:
-        kwargs_check.pop(popit)
+    for popit in ["field", "sumup", "squeeze", "pixel_agg"]:
+        kwargs_check.pop(popit, None)
     if kwargs_check:
         raise MagpylibBadUserInput(
             f"Keyword arguments {tuple(kwargs_check.keys())} are only allowed when the source "
@@ -205,22 +155,27 @@ def getBH_level2(sources, observers, **kwargs) -> np.ndarray:
     check_dimensions(sources)
     check_excitations(sources)
 
+
     # format observer inputs:
     #   allow only bare sensor, collection, pos_vec or list thereof
     #   transform input into an ordered list of sensors (pos_vec->pixel)
-    #   check if all pixel shapes are similar.
-    sensors = check_format_input_observers(observers)
-    pix_shape = sensors[0].pixel.shape
-    if pix_shape == (3,):
-        pix_shape = (1,3)
+    #   check if all pixel shapes are similar - or else if pixel_agg is given
+    pixel_agg = kwargs.get("pixel_agg", None)
+    pixel_agg_func = check_format_pixel_agg(pixel_agg)
+    sensors, pix_shapes = check_format_input_observers(observers, pixel_agg)
+    pix_nums = [
+        int(np.product(ps[:-1])) for ps in pix_shapes
+    ]  # number of pixel for each sensor
+    pix_inds = np.cumsum([0] + pix_nums)  # cummulative indices of pixel for each sensor
+    pix_all_same = len(set(pix_shapes)) == 1
 
     # check which sensors have unit roation
     #   so that they dont have to be rotated back later (performance issue)
     #   this check is made now when sensor paths are not yet tiled.
-    unitQ = np.array([0,0,0,1.])
-    unrotated_sensors = [all(all(r==unitQ)
-        for r in sens._orientation.as_quat())
-        for sens in sensors]
+    unitQ = np.array([0, 0, 0, 1.0])
+    unrotated_sensors = [
+        all(all(r == unitQ) for r in sens._orientation.as_quat()) for sens in sensors
+    ]
 
     # check which sensors have a static orientation
     #   either static sensor or translation path
@@ -228,111 +183,144 @@ def getBH_level2(sources, observers, **kwargs) -> np.ndarray:
     static_sensor_rot = check_static_sensor_orient(sensors)
 
     # some important quantities -------------------------------------------------
-    obj_list = set(src_list + sensors) # unique obj entries only !!!
-    l0 = len(sources)
-    l = len(src_list)
-    k = len(sensors)
+    obj_list = set(src_list + sensors)  # unique obj entries only !!!
+    num_of_sources = len(sources)
+    num_of_src_list = len(src_list)
+    num_of_sensors = len(sensors)
 
     # tile up paths -------------------------------------------------------------
     #   all obj paths that are shorter than max-length are filled up with the last
     #   postion/orientation of the object (static paths)
     path_lengths = [len(obj._position) for obj in obj_list]
-    m = max(path_lengths)
+    max_path_len = max(path_lengths)
 
     # objects to tile up and reset below
-    mask_reset = [m!=pl for pl in path_lengths]
-    reset_obj = [obj for obj,mask in zip(obj_list,mask_reset) if mask]
-    reset_obj_m0 = [pl for pl,mask in zip(path_lengths,mask_reset) if mask]
+    mask_reset = [max_path_len != pl for pl in path_lengths]
+    reset_obj = [obj for obj, mask in zip(obj_list, mask_reset) if mask]
+    reset_obj_m0 = [pl for pl, mask in zip(path_lengths, mask_reset) if mask]
 
-    if m>1:
-        for obj,m0 in zip(reset_obj, reset_obj_m0):
+    if max_path_len > 1:
+        for obj, m0 in zip(reset_obj, reset_obj_m0):
             # length to be tiled
-            m_tile = m-m0
+            m_tile = max_path_len - m0
             # tile up position
-            tile_pos = np.tile(obj._position[-1], (m_tile,1))
+            tile_pos = np.tile(obj._position[-1], (m_tile, 1))
             obj._position = np.concatenate((obj._position, tile_pos))
             # tile up orientation
-            tile_orient = np.tile(obj._orientation.as_quat()[-1], (m_tile,1))
+            tile_orient = np.tile(obj._orientation.as_quat()[-1], (m_tile, 1))
+            # FUTURE use Rotation.concatenate() requires scipy>=1.8 and python 3.8
             tile_orient = np.concatenate((obj._orientation.as_quat(), tile_orient))
             obj._orientation = R.from_quat(tile_orient)
 
     # combine information form all sensors to generate pos_obs with-------------
     #   shape (m * concat all sens flat pixel, 3)
     #   allows sensors with different pixel shapes <- relevant?
-    poso =[[r.apply(sens.pixel.reshape(-1,3)) + p
-            for r,p in zip(sens._orientation, sens._position)]
-           for sens in sensors]
-    poso = np.concatenate(poso,axis=1).reshape(-1,3)
+    poso = [
+        [
+            r.apply(sens.pixel.reshape(-1, 3)) + p
+            for r, p in zip(sens._orientation, sens._position)
+        ]
+        for sens in sensors
+    ]
+    poso = np.concatenate(poso, axis=1).reshape(-1, 3)
     n_pp = len(poso)
-    n_pix = int(n_pp/m)
+    n_pix = int(n_pp / max_path_len)
 
     # group similar source types----------------------------------------------
     groups = {}
-    for ind,src in enumerate(src_list):
-        if src._object_type=='CustomSource':
-            group_key = src.field_B_lambda if kwargs['field']=='B' else src.field_H_lambda
+    for ind, src in enumerate(src_list):
+        if src._object_type == "CustomSource":
+            group_key = (
+                src.field_B_lambda if kwargs["field"] == "B" else src.field_H_lambda
+            )
         else:
             group_key = src._object_type
         if group_key not in groups:
-            groups[group_key] = {'sources':[], 'order':[], 'source_type': src._object_type}
-        groups[group_key]['sources'].append(src)
-        groups[group_key]['order'].append(ind)
+            groups[group_key] = {
+                "sources": [],
+                "order": [],
+                "source_type": src._object_type,
+            }
+        groups[group_key]["sources"].append(src)
+        groups[group_key]["order"].append(ind)
 
     # evaluate each group in one vectorized step -------------------------------
-    B = np.empty((l,m,n_pix,3))                         # allocate B
+    B = np.empty((num_of_src_list, max_path_len, n_pix, 3))  # allocate B
     for group in groups.values():
-        lg = len(group['sources'])
-        gr = group['sources']
+        lg = len(group["sources"])
+        gr = group["sources"]
         src_dict = get_src_dict(gr, n_pix, n_pp, poso)  # compute array dict for level1
-        B_group = getBH_level1(field=kwargs['field'], **src_dict)       # compute field
-        B_group = B_group.reshape((lg,m,n_pix,3))       # reshape (2% slower for large arrays)
-        for i in range(lg):                             # put into dedicated positions in B
-            B[group['order'][i]] = B_group[i]
+        B_group = getBH_level1(field=kwargs["field"], **src_dict)  # compute field
+        B_group = B_group.reshape(
+            (lg, max_path_len, n_pix, 3)
+        )  # reshape (2% slower for large arrays)
+        for gr_ind in range(lg):  # put into dedicated positions in B
+            B[group["order"][gr_ind]] = B_group[gr_ind]
 
     # reshape output ----------------------------------------------------------------
     # rearrange B when there is at least one Collection with more than one source
-    if l > l0:
-        for i,src in enumerate(sources):
-            if src._object_type == 'Collection':
+    if num_of_src_list > num_of_sources:
+        for src_ind, src in enumerate(sources):
+            if src._object_type == "Collection":
                 col_len = len(format_obj_input(src, allow="sources"))
-                B[i] = np.sum(B[i:i+col_len],axis=0)    # set B[i] to sum of slice
-                B = np.delete(B,np.s_[i+1:i+col_len],0) # delete remaining part of slice
+                # set B[i] to sum of slice
+                B[src_ind] = np.sum(B[src_ind : src_ind + col_len], axis=0)
+                B = np.delete(
+                    B, np.s_[src_ind + 1 : src_ind + col_len], 0
+                )  # delete remaining part of slice
 
     # apply sensor rotations (after summation over collections to reduce rot.apply operations)
     #   note: replace by math.prod with python 3.8 or later
-    k_pixel = int(np.product(pix_shape[:-1])) # total number of pixel positions
-    for i,sens in enumerate(sensors):         # cycle through all sensors
-        if not unrotated_sensors[i]:          # apply operations only to rotated sensors
-            if static_sensor_rot[i]:          # special case: same rotation along path
-                # select part where rot is applied
-                Bpart = B[:,:,i*k_pixel:(i+1)*k_pixel]
-                # change shape from (l0,m,k_pixel,3) to (P,3) for rot package
-                Bpart_flat = np.reshape(Bpart, (k_pixel*l0*m,3))
-                # apply sensor rotation
-                Bpart_flat_rot = sens._orientation[0].inv().apply(Bpart_flat)
-                # overwrite Bpart in B
-                B[:,:,i*k_pixel:(i+1)*k_pixel] = np.reshape(Bpart_flat_rot, (l0,m,k_pixel,3))
-            else:                         # general case: different rotations along path
-                for j in range(m): # THIS LOOP IS EXTREMELY SLOW !!!! github issue #283
-                    Bpart = B[:,j,i*k_pixel:(i+1)*k_pixel]           # select part
-                    Bpart_flat = np.reshape(Bpart, (k_pixel*l0,3))   # flatten for rot package
-                    Bpart_flat_rot = sens._orientation[j].inv().apply(Bpart_flat)  # apply rotat
-                    B[:,j,i*k_pixel:(i+1)*k_pixel] = np.reshape(Bpart_flat_rot, (l0,k_pixel,3))
+    for sens_ind, sens in enumerate(sensors):  # cycle through all sensors
+        if not unrotated_sensors[sens_ind]:  # apply operations only to rotated sensors
+            # select part where rot is applied
+            Bpart = B[:, :, pix_inds[sens_ind] : pix_inds[sens_ind + 1]]
+            # change shape to (P,3) for rot package
+            Bpart_orig_shape = Bpart.shape
+            Bpart_flat = np.reshape(Bpart, (-1, 3))
+            # apply sensor rotation
+            if static_sensor_rot[sens_ind]:  # special case: same rotation along path
+                sens_orient = sens._orientation[0]
+            else:
+                sens_orient = R.from_quat(
+                    np.tile( # tile for each source from list
+                        np.repeat( # same orientation path index for all indices
+                            sens._orientation.as_quat(), pix_nums[sens_ind], axis=0
+                        ),
+                        (num_of_sources, 1),
+                    )
+                )
+            Bpart_flat_rot = sens_orient.inv().apply(Bpart_flat)
+            # overwrite Bpart in B
+            B[:, :, pix_inds[sens_ind] : pix_inds[sens_ind + 1]] = np.reshape(
+                Bpart_flat_rot, Bpart_orig_shape
+            )
 
     # rearrange sensor-pixel shape
-    sens_px_shape = (k,) + pix_shape
-    B = B.reshape((l0,m)+sens_px_shape)
-
-    #
-    if kwargs['sumup']:
+    if pix_all_same:
+        B = B.reshape((num_of_sources, max_path_len, num_of_sensors, *pix_shapes[0]))
+        # aggregate pixel values
+        if pixel_agg is not None:
+            B = pixel_agg_func(B, axis=tuple(range(3 - B.ndim, -1)))
+    else:  # pixel_agg is not None when pix_all_same, checked with
+        Bsplit = np.split(B, pix_inds[1:-1], axis=2)
+        Bagg = [np.expand_dims(pixel_agg_func(b, axis=2), axis=2) for b in Bsplit]
+        B = np.concatenate(Bagg, axis=2)
+
+    # sumup over sources
+    if kwargs["sumup"]:
         B = np.sum(B, axis=0, keepdims=True)
 
     # reduce all size-1 levels
-    if kwargs['squeeze']:
+    if kwargs["squeeze"]:
         B = np.squeeze(B)
+    elif pixel_agg is not None:
+        # add missing dimension since `pixel_agg` reduces pixel
+        # dimensions to zero. Only needed if `squeeze is False``
+        B = np.expand_dims(B, axis=-2)
 
     # reset tiled objects
-    for obj,m0 in zip(reset_obj, reset_obj_m0):
+    for obj, m0 in zip(reset_obj, reset_obj_m0):
         obj._position = obj._position[:m0]
         obj._orientation = obj._orientation[:m0]
 
diff --git a/magpylib/_src/fields/field_wrap_BH_level3.py b/magpylib/_src/fields/field_wrap_BH_level3.py
index dce4ac868..c1b9e2cfa 100644
--- a/magpylib/_src/fields/field_wrap_BH_level3.py
+++ b/magpylib/_src/fields/field_wrap_BH_level3.py
@@ -1,7 +1,9 @@
 from magpylib._src.fields.field_wrap_BH_level2 import getBH_level2
 
 
-def getB(sources=None, observers=None, sumup=False, squeeze=True, **kwargs):
+def getB(
+    sources=None, observers=None, sumup=False, squeeze=True, pixel_agg=None, **kwargs
+):
     """Compute B-field in [mT] for given sources and observers.
 
     Field implementations can be directly accessed (avoiding the object oriented
@@ -34,6 +36,11 @@ def getB(sources=None, observers=None, sumup=False, squeeze=True, **kwargs):
         If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g. only
         a single sensor or only a single source) are eliminated.
 
+    pixel_agg: str, default=`None`
+        Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+        which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+        With this option, observer inputs with different (pixel) shapes are allowed.
+
     Other Parameters (Direct interface)
     -----------------------------------
     position: array_like, shape (3,) or (n,3), default=`(0,0,0)`
@@ -138,10 +145,20 @@ def getB(sources=None, observers=None, sumup=False, squeeze=True, **kwargs):
      [ 1.14713551  2.29427102 -0.22065346]
      [-2.48213467 -2.48213467 -0.79683487]]
     """
-    return getBH_level2(sources, observers, sumup=sumup, squeeze=squeeze, field='B', **kwargs)
-
-
-def getH(sources=None, observers=None, sumup=False, squeeze=True, **kwargs):
+    return getBH_level2(
+        sources,
+        observers,
+        sumup=sumup,
+        squeeze=squeeze,
+        pixel_agg=pixel_agg,
+        field="B",
+        **kwargs
+    )
+
+
+def getH(
+    sources=None, observers=None, sumup=False, squeeze=True, pixel_agg=None, **kwargs
+):
     """Compute H-field in [kA/m] for given sources and observers.
 
     Field implementations can be directly accessed (avoiding the object oriented
@@ -174,6 +191,11 @@ def getH(sources=None, observers=None, sumup=False, squeeze=True, **kwargs):
         If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g. only
         a single sensor or only a single source) are eliminated.
 
+    pixel_agg: str, default=`None`
+        Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+        which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+        With this option, observer inputs with different (pixel) shapes are allowed.
+
     Other Parameters (Direct interface)
     -----------------------------------
     position: array_like, shape (3,) or (n,3), default=`(0,0,0)`
@@ -278,4 +300,12 @@ def getH(sources=None, observers=None, sumup=False, squeeze=True, **kwargs):
      [ 0.91286143  1.82572286 -0.17559045]
      [-1.97522001 -1.97522001 -0.63410104]]
     """
-    return getBH_level2(sources, observers, sumup=sumup, squeeze=squeeze, field='H', **kwargs)
+    return getBH_level2(
+        sources,
+        observers,
+        sumup=sumup,
+        squeeze=squeeze,
+        pixel_agg=pixel_agg,
+        field="H",
+        **kwargs
+    )
diff --git a/magpylib/_src/input_checks.py b/magpylib/_src/input_checks.py
index 57dae9623..14bacc55e 100644
--- a/magpylib/_src/input_checks.py
+++ b/magpylib/_src/input_checks.py
@@ -9,7 +9,12 @@
 )
 from magpylib._src.defaults.defaults_classes import default_settings
 from magpylib import _src
-from magpylib._src.utility import format_obj_input, wrong_obj_msg, LIBRARY_SOURCES, LIBRARY_SENSORS
+from magpylib._src.utility import (
+    format_obj_input,
+    wrong_obj_msg,
+    LIBRARY_SOURCES,
+    LIBRARY_SENSORS,
+)
 
 
 #################################################################
@@ -22,7 +27,7 @@ def all_same(lst: list) -> bool:
     return lst[1:] == lst[:-1]
 
 
-def is_array_like(inp, msg:str):
+def is_array_like(inp, msg: str):
     """ test if inp is array_like: type list, tuple or ndarray
     inp: test object
     msg: str, error msg
@@ -30,7 +35,8 @@ def is_array_like(inp, msg:str):
     if not isinstance(inp, (list, tuple, np.ndarray)):
         raise MagpylibBadUserInput(msg)
 
-def make_float_array(inp, msg:str):
+
+def make_float_array(inp, msg: str):
     """transform inp to array with dtype=float, throw error with bad input
     inp: test object
     msg: str, error msg
@@ -42,7 +48,7 @@ def make_float_array(inp, msg:str):
     return inp_array
 
 
-def check_array_shape(inp: np.ndarray, dims:tuple, shape_m1:int, msg:str):
+def check_array_shape(inp: np.ndarray, dims: tuple, shape_m1: int, msg: str):
     """check if inp shape is allowed
     inp: test object
     dims: list, list of allowed dims
@@ -50,9 +56,9 @@ def check_array_shape(inp: np.ndarray, dims:tuple, shape_m1:int, msg:str):
     msg: str, error msg
     """
     if inp.ndim in dims:
-        if inp.shape[-1]==shape_m1:
+        if inp.shape[-1] == shape_m1:
             return None
-        if shape_m1 == 'any':
+        if shape_m1 == "any":
             return None
     raise MagpylibBadUserInput(msg)
 
@@ -64,21 +70,22 @@ def check_input_zoom(inp):
             "Input parameter `zoom` must be a number `zoom>=0`.\n"
             f"Instead received {inp}."
         )
-    if inp<0:
+    if inp < 0:
         raise MagpylibBadUserInput(
             "Input parameter `zoom` must be a number `zoom>=0`.\n"
             f"Instead received {inp}."
         )
 
+
 def check_input_animation(inp):
     """check show animation input"""
     ERR_MSG = (
         "Input parameter `animation` must be boolean or a positive number.\n"
         f"Instead received {inp}."
-        )
+    )
     if not isinstance(inp, numbers.Number):
         raise MagpylibBadUserInput(ERR_MSG)
-    if inp<0:
+    if inp < 0:
         raise MagpylibBadUserInput(ERR_MSG)
 
 
@@ -89,7 +96,9 @@ def check_input_animation(inp):
 
 def check_start_type(inp):
     """start input must be int or str"""
-    if not (isinstance(inp, (int, np.integer)) or (isinstance(inp, str) and inp == 'auto')):
+    if not (
+        isinstance(inp, (int, np.integer)) or (isinstance(inp, str) and inp == "auto")
+    ):
         raise MagpylibBadUserInput(
             f"Input parameter `start` must be integer value or 'auto'.\n"
             f"Instead received {repr(inp)}."
@@ -105,13 +114,12 @@ def check_degree_type(inp):
         )
 
 
-
 def check_field_input(inp, origin):
     """check field input"""
     if isinstance(inp, str):
-        if inp == 'B':
+        if inp == "B":
             return True
-        if inp == 'H':
+        if inp == "H":
             return False
     raise MagpylibBadUserInput(
         f"{origin} input can only be `field='B'` or `field='H'`.\n"
@@ -132,14 +140,15 @@ def validate_field_lambda(val, bh):
 
         out = val(np.array([[1, 2, 3], [4, 5, 6]]))
         out_shape = np.array(out).shape
-        case2 = out_shape!=(2, 3)
+        case2 = out_shape != (2, 3)
 
         if case2:
             raise MagpylibBadUserInput(
                 f"Input parameter `field_{bh}_lambda` must be a callable function"
                 " and return a field ndarray of shape (n,3) when its `observer`"
                 " input is of shape (n,3).\n"
-                f"Instead received shape {out_shape}.")
+                f"Instead received shape {out_shape}."
+            )
     return val
 
 
@@ -147,6 +156,7 @@ def validate_field_lambda(val, bh):
 #################################################################
 # CHECK - FORMAT
 
+
 def check_format_input_orientation(inp, init_format=False):
     """checks orientation input returns in formatted form
     - inp must be None or Rotation object
@@ -163,16 +173,17 @@ def check_format_input_orientation(inp, init_format=False):
     if not isinstance(inp, (Rotation, type(None))):
         raise MagpylibBadUserInput(
             f"Input parameter `orientation` must be `None` or scipy `Rotation` object.\n"
-            f"Instead received type {type(inp)}.")
+            f"Instead received type {type(inp)}."
+        )
     # handle None input and compute inpQ
     if inp is None:
-        inpQ = np.array((0,0,0,1))
-        inp=Rotation.from_quat(inpQ)
+        inpQ = np.array((0, 0, 0, 1))
+        inp = Rotation.from_quat(inpQ)
     else:
         inpQ = inp.as_quat()
     # return
     if init_format:
-        return np.reshape(inpQ, (-1,4))
+        return np.reshape(inpQ, (-1, 4))
     return inp, inpQ
 
 
@@ -183,12 +194,14 @@ def check_format_input_anchor(inp):
     if isinstance(inp, numbers.Number) and inp == 0:
         return np.array((0.0, 0.0, 0.0))
 
-    return check_format_input_vector(inp,
-        dims=(1,2),
+    return check_format_input_vector(
+        inp,
+        dims=(1, 2),
         shape_m1=3,
-        sig_name='anchor',
-        sig_type='`None` or `0` or array_like (list, tuple, ndarray) with shape (3,)',
-        allow_None=True)
+        sig_name="anchor",
+        sig_type="`None` or `0` or array_like (list, tuple, ndarray) with shape (3,)",
+        allow_None=True,
+    )
 
 
 def check_format_input_axis(inp):
@@ -201,25 +214,27 @@ def check_format_input_axis(inp):
     - return as ndarray shape (3,)
     """
     if isinstance(inp, str):
-        if inp == 'x':
-            return np.array((1,0,0))
-        if inp == 'y':
-            return np.array((0,1,0))
-        if inp == 'z':
-            return np.array((0,0,1))
+        if inp == "x":
+            return np.array((1, 0, 0))
+        if inp == "y":
+            return np.array((0, 1, 0))
+        if inp == "z":
+            return np.array((0, 0, 1))
         raise MagpylibBadUserInput(
             "Input parameter `axis` must be array_like shape (3,) or one of ['x', 'y', 'z'].\n"
-            f"Instead received string {inp}.\n")
+            f"Instead received string {inp}.\n"
+        )
 
-    inp = check_format_input_vector(inp,
+    inp = check_format_input_vector(
+        inp,
         dims=(1,),
         shape_m1=3,
-        sig_name='axis',
-        sig_type="array_like (list, tuple, ndarray) with shape (3,) or one of ['x', 'y', 'z']")
+        sig_name="axis",
+        sig_type="array_like (list, tuple, ndarray) with shape (3,) or one of ['x', 'y', 'z']",
+    )
 
-    if np.all(inp==0):
-        raise MagpylibBadUserInput(
-            "Input parameter `axis` must not be (0,0,0).\n")
+    if np.all(inp == 0):
+        raise MagpylibBadUserInput("Input parameter `axis` must not be (0,0,0).\n")
     return inp
 
 
@@ -236,19 +251,18 @@ def check_format_input_angle(inp):
     if isinstance(inp, numbers.Number):
         return float(inp)
 
-    return check_format_input_vector(inp,
+    return check_format_input_vector(
+        inp,
         dims=(1,),
-        shape_m1='any',
-        sig_name='angle',
-        sig_type='int, float or array_like (list, tuple, ndarray) with shape (n,)')
+        shape_m1="any",
+        sig_name="angle",
+        sig_type="int, float or array_like (list, tuple, ndarray) with shape (n,)",
+    )
 
 
 def check_format_input_scalar(
-    inp,
-    sig_name,
-    sig_type,
-    allow_None=False,
-    forbid_negative=False):
+    inp, sig_name, sig_type, allow_None=False, forbid_negative=False
+):
     """check sclar input and return in formatted form
     - must be scalar or None (if allowed)
     - must be float compatible
@@ -260,7 +274,8 @@ def check_format_input_scalar(
 
     ERR_MSG = (
         f"Input parameter `{sig_name}` must be {sig_type}.\n"
-        f"Instead received {repr(inp)}.")
+        f"Instead received {repr(inp)}."
+    )
 
     if not isinstance(inp, numbers.Number):
         raise MagpylibBadUserInput(ERR_MSG)
@@ -268,12 +283,13 @@ def check_format_input_scalar(
     inp = float(inp)
 
     if forbid_negative:
-        if inp<0:
+        if inp < 0:
             raise MagpylibBadUserInput(ERR_MSG)
     return inp
 
 
-def check_format_input_vector(inp,
+def check_format_input_vector(
+    inp,
     dims,
     shape_m1,
     sig_name,
@@ -281,7 +297,7 @@ def check_format_input_vector(inp,
     reshape=False,
     allow_None=False,
     forbid_negative0=False,
-    ):
+):
     """checks vector input and returns in formatted form
     - inp must be array_like
     - convert inp to ndarray with dtype float
@@ -295,21 +311,28 @@ def check_format_input_vector(inp,
         if inp is None:
             return None
 
-    is_array_like(inp,
+    is_array_like(
+        inp,
         f"Input parameter `{sig_name}` must be {sig_type}.\n"
-        f"Instead received type {type(inp)}."
+        f"Instead received type {type(inp)}.",
     )
-    inp = make_float_array(inp,
-        f"Input parameter `{sig_name}` must contain only float compatible entries.\n"
+    inp = make_float_array(
+        inp,
+        f"Input parameter `{sig_name}` must contain only float compatible entries.\n",
     )
-    check_array_shape(inp, dims=dims, shape_m1=shape_m1, msg=(
-        f"Input parameter `{sig_name}` must be {sig_type}.\n"
-        f"Instead received array_like with shape {inp.shape}.")
+    check_array_shape(
+        inp,
+        dims=dims,
+        shape_m1=shape_m1,
+        msg=(
+            f"Input parameter `{sig_name}` must be {sig_type}.\n"
+            f"Instead received array_like with shape {inp.shape}."
+        ),
     )
     if reshape:
-        return np.reshape(inp, (-1,3))
+        return np.reshape(inp, (-1, 3))
     if forbid_negative0:
-        if np.any(inp<=0):
+        if np.any(inp <= 0):
             raise MagpylibBadUserInput(
                 f"Input parameter `{sig_name}` cannot have values <= 0."
             )
@@ -320,15 +343,17 @@ def check_format_input_vertices(inp):
     """checks vertices input and returns in formatted form
     - vector check with dim = (n,3) but n must be >=2
     """
-    inp = check_format_input_vector(inp,
+    inp = check_format_input_vector(
+        inp,
         dims=(2,),
         shape_m1=3,
-        sig_name='vertices',
-        sig_type='`None` or array_like (list, tuple, ndarray) with shape (n,3)',
-        allow_None=True)
+        sig_name="vertices",
+        sig_type="`None` or array_like (list, tuple, ndarray) with shape (n,3)",
+        allow_None=True,
+    )
 
     if inp is not None:
-        if inp.shape[0]<2:
+        if inp.shape[0] < 2:
             raise MagpylibBadUserInput(
                 "Input parameter `vertices` must have more than one vertex."
             )
@@ -342,23 +367,26 @@ def check_format_input_cylinder_segment(inp):
     - check if phi2-phi1 > 360
     - return error msg
     """
-    inp = check_format_input_vector(inp,
+    inp = check_format_input_vector(
+        inp,
         dims=(1,),
         shape_m1=5,
-        sig_name='CylinderSegment.dimension',
+        sig_name="CylinderSegment.dimension",
         sig_type=(
-            'array_like of the form (r1, r2, h, phi1, phi2) with r1<r2,'
-            'phi1<phi2 and phi2-phi1<=360'),
-        allow_None=True)
+            "array_like of the form (r1, r2, h, phi1, phi2) with r1<r2,"
+            "phi1<phi2 and phi2-phi1<=360"
+        ),
+        allow_None=True,
+    )
 
     if inp is None:
         return None
 
     r1, r2, h, phi1, phi2 = inp
-    case2 = r1>r2
-    case3 = phi1>phi2
-    case4 = (phi2 - phi1)>360
-    case5 = (r1<0) | (r2<=0) | (h<=0)
+    case2 = r1 > r2
+    case3 = phi1 > phi2
+    case4 = (phi2 - phi1) > 360
+    case5 = (r1 < 0) | (r2 <= 0) | (h <= 0)
     if case2 | case3 | case4 | case5:
         raise MagpylibBadUserInput(
             f"Input parameter `CylinderSegment.dimension` must be array_like of the form"
@@ -372,14 +400,15 @@ def check_format_input_backend(inp):
     """checks show-backend input and returns Non if bad input value"""
     if inp is None:
         inp = default_settings.display.backend
-    if inp in ('matplotlib', 'plotly'):
+    if inp in ("matplotlib", "plotly"):
         return inp
     raise MagpylibBadUserInput(
         "Input parameter `backend` must be one of `('matplotlib', 'plotly', None)`.\n"
-        f"Instead received {inp}.")
+        f"Instead received {inp}."
+    )
 
 
-def check_format_input_observers(inp):
+def check_format_input_observers(inp, pixel_agg=None):
     """
     checks observer input and returns a list of sensor objects
     """
@@ -402,34 +431,38 @@ def check_format_input_observers(inp):
 
     # now inp can still be [pos_vec, sens, coll] or just a pos_vec
 
-    try: # try if input is just a pos_vec
+    try:  # try if input is just a pos_vec
         inp = np.array(inp, dtype=float)
-        return [_src.obj_classes.Sensor(pixel=inp)]
-    except (TypeError, ValueError): # if not, it must be [pos_vec, sens, coll]
-        sensors=[]
+        pix_shapes = [(1, 3) if inp.shape == (3,) else inp.shape]
+        return [_src.obj_classes.Sensor(pixel=inp)], pix_shapes
+    except (TypeError, ValueError):  # if not, it must be [pos_vec, sens, coll]
+        sensors = []
         for obj in inp:
             if getattr(obj, "_object_type", "") == "Sensor":
                 sensors.append(obj)
             elif getattr(obj, "_object_type", "") == "Collection":
-                child_sensors = format_obj_input(obj, allow='sensors')
+                child_sensors = format_obj_input(obj, allow="sensors")
                 if not child_sensors:
                     raise MagpylibBadUserInput(wrong_obj_msg(obj, allow="observers"))
                 sensors.extend(child_sensors)
-            else: # if its not a Sensor or a Collection it can only be a pos_vec
+            else:  # if its not a Sensor or a Collection it can only be a pos_vec
                 try:
                     obj = np.array(obj, dtype=float)
                     sensors.append(_src.obj_classes.Sensor(pixel=obj))
-                except Exception: # or some unwanted crap
+                except Exception:  # or some unwanted crap
                     raise MagpylibBadUserInput(wrong_obj_msg(obj, allow="observers"))
 
         # all pixel shapes must be the same
-        pix_shapes = [s._pixel.shape for s in sensors]
-        if not all_same(pix_shapes):
+        pix_shapes = [
+            (1, 3) if s.pixel.shape == (3,) else s.pixel.shape for s in sensors
+        ]
+        if pixel_agg is None and not all_same(pix_shapes):
             raise MagpylibBadUserInput(
-                'Different observer input detected.'
-                ' All sensor pixel and position vector inputs must'
-                ' be of similar shape !')
-        return sensors
+                "Different observer input shape detected."
+                " All observer inputs must be of similar shape, unless a"
+                " numpy pixel aggregator is provided, e.g. `pixel_agg='mean'`!"
+            )
+        return sensors, pix_shapes
 
 
 def check_format_input_obj(
@@ -497,18 +530,19 @@ def check_format_input_obj(
 ############################################################################################
 # SHOW AND GETB CHECKS
 
+
 def check_dimensions(sources):
     """check if all sources have dimension (or similar) initialized"""
     # pylint: disable=protected-access
     for s in sources:
-        obj_type = getattr(s, '_object_type', None)
-        if obj_type in ('Cuboid', 'Cylinder', 'CylinderSegment'):
+        obj_type = getattr(s, "_object_type", None)
+        if obj_type in ("Cuboid", "Cylinder", "CylinderSegment"):
             if s.dimension is None:
                 raise MagpylibMissingInput(f"Parameter `dimension` of {s} must be set.")
-        elif obj_type in ('Sphere', 'Loop'):
+        elif obj_type in ("Sphere", "Loop"):
             if s.diameter is None:
                 raise MagpylibMissingInput(f"Parameter `diameter` of {s} must be set.")
-        elif obj_type == 'Line':
+        elif obj_type == "Line":
             if s.vertices is None:
                 raise MagpylibMissingInput(f"Parameter `vertices` of {s} must be set.")
 
@@ -517,13 +551,44 @@ def check_excitations(sources):
     """check if all sources have exitation initialized"""
     # pylint: disable=protected-access
     for s in sources:
-        obj_type = getattr(s, '_object_type', None)
-        if obj_type in ('Cuboid', 'Cylinder', 'Sphere', 'CylinderSegment'):
+        obj_type = getattr(s, "_object_type", None)
+        if obj_type in ("Cuboid", "Cylinder", "Sphere", "CylinderSegment"):
             if s.magnetization is None:
-                raise MagpylibMissingInput(f"Parameter `magnetization` of {s} must be set.")
-        elif obj_type in ('Loop', 'Line'):
+                raise MagpylibMissingInput(
+                    f"Parameter `magnetization` of {s} must be set."
+                )
+        elif obj_type in ("Loop", "Line"):
             if s.current is None:
                 raise MagpylibMissingInput(f"Parameter `current` of {s} must be set.")
-        elif obj_type == 'Dipole':
+        elif obj_type == "Dipole":
             if s.moment is None:
                 raise MagpylibMissingInput(f"Parameter `moment` of {s} must be set.")
+
+
+def check_format_pixel_agg(pixel_agg):
+    """
+    check if pixel_agg input is acceptable
+    return the respective numpy function
+    """
+
+    PIXEL_AGG_ERR_MSG = (
+        "Input `pixel_agg` must be a reference to a numpy callable that reduces" +
+        " an array shape like 'mean', 'std', 'median', 'min', ...\n" +
+        f"Instead received {pixel_agg}."
+    )
+
+    if pixel_agg is None:
+        return None
+
+    # test numpy reference
+    try:
+        pixel_agg_func = getattr(np, pixel_agg)
+    except AttributeError as err:
+        raise AttributeError(PIXEL_AGG_ERR_MSG) from err
+
+    # test pixel agg function reduce
+    x = np.array([[[(1,2,3)]*2]*3]*4)
+    if not isinstance(pixel_agg_func(x), numbers.Number):
+        raise AttributeError(PIXEL_AGG_ERR_MSG)
+
+    return pixel_agg_func
diff --git a/magpylib/_src/obj_classes/class_BaseGetBH.py b/magpylib/_src/obj_classes/class_BaseGetBH.py
index 148970e0c..d091a0cf9 100644
--- a/magpylib/_src/obj_classes/class_BaseGetBH.py
+++ b/magpylib/_src/obj_classes/class_BaseGetBH.py
@@ -7,10 +7,9 @@
 
 
 class BaseGetBH:
-    """provides getB and getH methods for source objects
-    """
+    """provides getB and getH methods for source objects"""
 
-    def getB(self, *observers, squeeze=True):
+    def getB(self, *observers, squeeze=True, pixel_agg=None):
         """Compute the B-field in units of [mT] generated by the source.
 
         Parameters
@@ -25,6 +24,11 @@ def getB(self, *observers, squeeze=True):
             If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g.
             only a single source) are eliminated.
 
+        pixel_agg: str, default=`None`
+            Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+            which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+            With this option, observer inputs with different (pixel) shapes are allowed.
+
         Returns
         -------
         B-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3)
@@ -57,10 +61,16 @@ def getB(self, *observers, squeeze=True):
         >>>   [ -0.77638847   0.          -5.15004352]]]
         """
         observers = format_star_input(observers)
-        return getBH_level2(self, observers, sumup=False, squeeze=squeeze, field='B')
-
-
-    def getH(self, *observers, squeeze=True):
+        return getBH_level2(
+            self,
+            observers,
+            sumup=False,
+            squeeze=squeeze,
+            pixel_agg=pixel_agg,
+            field="B",
+        )
+
+    def getH(self, *observers, squeeze=True, pixel_agg=None):
         """Compute the H-field in units of [kA/m] generated by the source.
 
         Parameters
@@ -75,6 +85,11 @@ def getH(self, *observers, squeeze=True):
             If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g.
             only a single source) are eliminated.
 
+        pixel_agg: str, default=`None`
+            Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+            which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+            With this option, observer inputs with different (pixel) shapes are allowed.
+
         Returns
         -------
         H-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3)
@@ -109,4 +124,11 @@ def getH(self, *observers, squeeze=True):
           [ -0.61783031   0.          -4.09827441]]]
         """
         observers = format_star_input(observers)
-        return getBH_level2(self, observers, sumup=False, squeeze=squeeze, field='H')
+        return getBH_level2(
+            self,
+            observers,
+            sumup=False,
+            squeeze=squeeze,
+            pixel_agg=pixel_agg,
+            field="H",
+        )
diff --git a/magpylib/_src/obj_classes/class_Collection.py b/magpylib/_src/obj_classes/class_Collection.py
index 03eb2b953..90f0e4efd 100644
--- a/magpylib/_src/obj_classes/class_Collection.py
+++ b/magpylib/_src/obj_classes/class_Collection.py
@@ -484,7 +484,7 @@ def _validate_getBH_inputs(self, *inputs):
             sources, sensors = self, inputs
         return sources, sensors
 
-    def getB(self, *inputs, squeeze=True):
+    def getB(self, *inputs, squeeze=True, pixel_agg=None):
         """Compute B-field in [mT] for given sources and observer inputs.
 
         Parameters
@@ -499,6 +499,11 @@ def getB(self, *inputs, squeeze=True):
             If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g.
             only a single source) are eliminated.
 
+        pixel_agg: str, default=`None`
+            Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+            which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+            With this option, observer inputs with different (pixel) shapes are allowed.
+
         Returns
         -------
         B-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3)
@@ -531,9 +536,17 @@ def getB(self, *inputs, squeeze=True):
         """
 
         sources, sensors = self._validate_getBH_inputs(*inputs)
-        return getBH_level2(sources, sensors, sumup=False, squeeze=squeeze, field="B")
 
-    def getH(self, *inputs, squeeze=True):
+        return getBH_level2(
+            sources,
+            sensors,
+            sumup=False,
+            squeeze=squeeze,
+            pixel_agg=pixel_agg,
+            field="B",
+        )
+
+    def getH(self, *inputs, squeeze=True, pixel_agg=None):
         """Compute H-field in [kA/m] for given sources and observer inputs.
 
         Parameters
@@ -548,6 +561,11 @@ def getH(self, *inputs, squeeze=True):
             If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g.
             only a single source) are eliminated.
 
+        pixel_agg: str, default=`None`
+            Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+            which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+            With this option, observer inputs with different (pixel) shapes are allowed.
+
         Returns
         -------
         H-field: ndarray, shape squeeze(m, k, n1, n2, ..., 3)
@@ -581,7 +599,14 @@ def getH(self, *inputs, squeeze=True):
 
         sources, sensors = self._validate_getBH_inputs(*inputs)
 
-        return getBH_level2(sources, sensors, sumup=False, squeeze=squeeze, field="H")
+        return getBH_level2(
+            sources,
+            sensors,
+            sumup=False,
+            squeeze=squeeze,
+            pixel_agg=pixel_agg,
+            field="H",
+        )
 
 
 class Collection(BaseGeo, BaseCollection):
diff --git a/magpylib/_src/obj_classes/class_Sensor.py b/magpylib/_src/obj_classes/class_Sensor.py
index f32e24b67..9e3f19958 100644
--- a/magpylib/_src/obj_classes/class_Sensor.py
+++ b/magpylib/_src/obj_classes/class_Sensor.py
@@ -10,7 +10,7 @@
 
 
 class Sensor(BaseGeo, BaseDisplayRepr):
-    """ Magnetic field sensor.
+    """Magnetic field sensor.
 
     Can be used as `observers` input for magnetic field computation.
 
@@ -85,7 +85,7 @@ def __init__(
 
         # instance attributes
         self.pixel = pixel
-        self._object_type = 'Sensor'
+        self._object_type = "Sensor"
 
         # init inheritance
         BaseGeo.__init__(self, position, orientation, style=style, **kwargs)
@@ -101,19 +101,18 @@ def pixel(self):
 
     @pixel.setter
     def pixel(self, pix):
-        """ Set sensor pixel positions in the local sensor coordinates.
+        """Set sensor pixel positions in the local sensor coordinates.
         Must be an array_like, float compatible with shape (..., 3)
         """
         self._pixel = check_format_input_vector(
             pix,
-            dims=range(1,20),
+            dims=range(1, 20),
             shape_m1=3,
-            sig_name='pixel',
-            sig_type='array_like (list, tuple, ndarray) with shape (n1, n2, ..., 3)',
+            sig_name="pixel",
+            sig_type="array_like (list, tuple, ndarray) with shape (n1, n2, ..., 3)",
         )
 
-
-    def getB(self, *sources, sumup=False, squeeze=True):
+    def getB(self, *sources, sumup=False, squeeze=True, pixel_agg=None):
         """Compute the B-field in units of [mT] as seen by the sensor.
 
         Parameters
@@ -129,6 +128,11 @@ def getB(self, *sources, sumup=False, squeeze=True):
             If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g. only
             a single sensor or only a single source) are eliminated.
 
+        pixel_agg: str, default=`None`
+            Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+            which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+            With this option, observer inputs with different (pixel) shapes are allowed.
+
         Returns
         -------
         B-field: ndarray, shape squeeze(l, m, n1, n2, ..., 3)
@@ -165,10 +169,11 @@ def getB(self, *sources, sumup=False, squeeze=True):
          [0.         1.01415383 1.01415383]]
         """
         sources = format_star_input(sources)
-        return getBH_level2(sources, self, sumup=sumup, squeeze=squeeze, field='B')
-
+        return getBH_level2(
+            sources, self, sumup=sumup, squeeze=squeeze, pixel_agg=pixel_agg, field="B"
+        )
 
-    def getH(self, *sources, sumup=False, squeeze=True):
+    def getH(self, *sources, sumup=False, squeeze=True, pixel_agg=None):
         """Compute the H-field in units of [kA/m] as seen by the sensor.
 
         Parameters
@@ -184,6 +189,11 @@ def getH(self, *sources, sumup=False, squeeze=True):
             If `True`, the output is squeezed, i.e. all axes of length 1 in the output (e.g. only
             a single sensor or only a single source) are eliminated.
 
+        pixel_agg: str, default=`None`
+            Reference to a compatible numpy aggregator function like `'min'` or `'mean'`,
+            which is applied to observer output values, e.g. mean of all sensor pixel outputs.
+            With this option, observer inputs with different (pixel) shapes are allowed.
+
         Returns
         -------
         H-field: ndarray, shape squeeze(l, m, n1, n2, ..., 3)
@@ -220,4 +230,6 @@ def getH(self, *sources, sumup=False, squeeze=True):
          [0.         0.80703798 0.80703798]]
         """
         sources = format_star_input(sources)
-        return getBH_level2(sources, self, sumup=sumup, squeeze=squeeze, field='H')
+        return getBH_level2(
+            sources, self, sumup=sumup, squeeze=squeeze, pixel_agg=pixel_agg, field="H"
+        )


Test patch:
diff --git a/tests/test_getBH_level2.py b/tests/test_getBH_level2.py
index 649af1f7c..0517a29c3 100644
--- a/tests/test_getBH_level2.py
+++ b/tests/test_getBH_level2.py
@@ -1,5 +1,7 @@
 import numpy as np
+import pytest
 import magpylib as magpy
+from magpylib._src.exceptions import MagpylibBadUserInput
 
 
 def test_getB_level2_input_simple():
@@ -295,3 +297,138 @@ def test_squeeze_sumup():
     B2 = magpy.getB(ss, s, squeeze=False, sumup=True)
 
     assert B1.shape == B2.shape
+
+
+def test_pixel_agg():
+    """test pixel aggregator"""
+    src1 = magpy.magnet.Cuboid((0,0,1000),(1,1,1)).move([[1,0,0]])
+    sens1 = magpy.Sensor(position=(0,0,1), pixel=np.zeros((4,5,3)), style_label='sens1 pixel(4,5)')
+    sens2 = sens1.copy(position=(0,0,2), style_label='sens2 pixel(4,5)')
+    sens3 = sens1.copy(position=(0,0,3), style_label='sens3 pixel(4,5)')
+    sens_col = magpy.Collection(sens1, sens2, sens3)
+
+    B1 = magpy.getB(src1, sens_col, squeeze=False, pixel_agg=None)
+    np.testing.assert_array_equal(B1.shape, (1, 2, 3, 4, 5, 3))
+
+    B2 = magpy.getB(src1, sens_col, squeeze=False, pixel_agg='mean')
+    np.testing.assert_array_equal(B2.shape, (1, 2, 3, 1, 3))
+
+    B3 = magpy.getB(src1, sens_col, squeeze=True, pixel_agg=None)
+    np.testing.assert_array_equal(B3.shape, (2, 3, 4, 5, 3))
+
+    B4 = magpy.getB(src1, sens_col, squeeze=True, pixel_agg='mean')
+    np.testing.assert_array_equal(B4.shape, (2, 3, 3))
+
+
+def test_pixel_agg_heterogeneous_pixel_shapes():
+    """test pixel aggregator with heterogeneous pixel shapes"""
+    src1 = magpy.magnet.Cuboid((0,0,1000),(1,1,1))
+    src2 = magpy.magnet.Sphere((0,0,1000),1, position=(2,0,0))
+    sens1 = magpy.Sensor(position=(0,0,1), pixel=[0,0,0], style_label='sens1, pixel.shape = (3,)')
+    sens2 = sens1.copy(position=(0,0,2), pixel=[1,1,1], style_label='sens2,  pixel.shape = (3,)')
+    sens3 = sens1.copy(position=(0,0,3), pixel=[2,2,2],  style_label='sens3,  pixel.shape = (3,)')
+    sens4 = sens1.copy(style_label='sens4,  pixel.shape = (3,)')
+    sens5 = sens2.copy(pixel=np.zeros((4,5,3))+1, style_label='sens5,  pixel.shape = (3,)')
+    sens6 = sens3.copy(pixel=np.zeros((4,5,1,3))+2, style_label='sens6,  pixel.shape = (4,5,1,3)')
+    src_col = magpy.Collection(src1, src2)
+    sens_col1 = magpy.Collection(sens1, sens2, sens3)
+    sens_col2 = magpy.Collection(sens4, sens5, sens6)
+    sens_col1.rotate_from_angax([45], 'z', anchor = (5,0,0))
+    sens_col2.rotate_from_angax([45], 'z', anchor = (5,0,0))
+
+    # different pixel shapes withoug pixel_agg should raise an error
+    with pytest.raises(MagpylibBadUserInput):
+        magpy.getB(src1, sens_col2, pixel_agg=None)
+
+    # bad pixel_agg numpy reference
+    with pytest.raises(AttributeError):
+        magpy.getB(src1, sens_col2, pixel_agg='bad_aggregator')
+
+    # good pixel_agg numpy reference, but non-reducing function
+    with pytest.raises(AttributeError):
+        magpy.getB(src1, sens_col2, pixel_agg='array')
+
+    B1 = magpy.getB(src1, sens_col1, squeeze=False, pixel_agg='max')
+    np.testing.assert_array_equal(B1.shape, (1, 2, 3, 1, 3))
+
+    B2 = magpy.getB(src1, sens_col2, squeeze=False, pixel_agg='max')
+    np.testing.assert_array_equal(B2.shape, (1, 2, 3, 1, 3))
+
+    B3 = magpy.getB(src1, sens_col1, squeeze=True)
+    np.testing.assert_array_equal(B3.shape, (2, 3, 3))
+
+    B4 = magpy.getB(src1, sens_col2, squeeze=True, pixel_agg='mean')
+    np.testing.assert_array_equal(B4.shape, (2, 3, 3))
+
+    # B3 and B4 should deliver the same results since pixel all have the same
+    # positions respectively for each sensor, so mean equals single value
+    np.testing.assert_allclose(B3, B4)
+
+    # Testing autmatic vs manual aggregation (mean) with different pixel shapes
+    B_by_sens_agg_1 = magpy.getB(src_col, sens_col2, squeeze=False, pixel_agg='mean')
+    B_by_sens_agg_2 = []
+    for sens in sens_col2:
+        B = magpy.getB(src_col, sens, squeeze=False)
+        B = B.mean(axis=tuple(range(3 - B.ndim, -1)))
+        B = np.expand_dims(B, axis=-2)
+        B_by_sens_agg_2.append(B)
+    B_by_sens_agg_2 = np.concatenate(B_by_sens_agg_2, axis=2)
+
+    np.testing.assert_allclose(B_by_sens_agg_1, B_by_sens_agg_2)
+
+
+def test_pixel_agg3():
+    """ test for various inputs"""
+    B1 = np.array([0.03122074, 0.03122074, 0.03122074])
+
+    e0 = np.array((1,1,1))
+    e1 = [(1,1,1)]
+    e2 = [(1,1,1)]*2
+    e3 = [(1,1,1)]*3
+
+    s0 = magpy.magnet.Cuboid(e0, e0)
+    c0 = magpy.Collection(s0)
+    s1 = magpy.magnet.Cuboid( e0, e0)
+    s2 = magpy.magnet.Cuboid(-e0, e0)
+    c1 = magpy.Collection(c0, s1, s2)
+
+    x0 = magpy.Sensor(pixel=e0)
+    x1 = magpy.Sensor(pixel=e1)
+    x2 = magpy.Sensor(pixel=e2)
+    x3 = magpy.Sensor(pixel=e3)
+
+    c2 = x0 + x1 + x2 + x3
+
+    for src, src_sh in zip(
+        [s0, c0, [s0, c0], c1, [s0, c0, c1, s1]],
+        [1, 1, 2, 1, 4]
+        ):
+        for obs, obs_sh in zip(
+            [e0, e1, e2, e3, x0, x1, x2, x3, c2, [x0, x2, x3]],
+            [1]*8 + [4,3]
+            ):
+            for px_agg in ['mean', 'average', 'min']:
+                np.testing.assert_allclose(
+                    magpy.getB(src, obs, pixel_agg=px_agg),
+                    np.squeeze(np.tile(B1, (src_sh,obs_sh,1))),
+                    rtol=1e-5,
+                    atol=1e-8
+                )
+
+    # same check with a path
+    s0.position = [(0,0,0)]*5
+    for src, src_sh in zip(
+        [s0, c0, [s0, c0], c1, [s0, c0, c1, s1]],
+        [1, 1, 2, 1, 4]
+        ):
+        for obs, obs_sh in zip(
+            [e0, e1, e2, e3, x0, x1, x2, x3, c2, [x0, x2, x3]],
+            [1]*8 + [4,3]
+            ):
+            for px_agg in ['mean', 'average', 'min']:
+                np.testing.assert_allclose(
+                    magpy.getB(src, obs, pixel_agg=px_agg),
+                    np.squeeze(np.tile(B1, (src_sh, 5, obs_sh, 1))),
+                    rtol=1e-5,
+                    atol=1e-8
+                )