diff --git a/.gitignore b/.gitignore index 32a2b189..3922124d 100644 --- a/.gitignore +++ b/.gitignore @@ -28,3 +28,9 @@ docs/_build /subprojects/* !/subprojects/*.wrap .junie + +# NetCDF files written to the working directory by the test suite +/angstrom_unit.nc +/no_physical_sizes.nc +/no_unit.nc +/parallel_save_test.nc diff --git a/SurfaceTopography/ChangeLog.md b/SurfaceTopography/ChangeLog.md index e6bda451..925c5a3f 100644 --- a/SurfaceTopography/ChangeLog.md +++ b/SurfaceTopography/ChangeLog.md @@ -1,6 +1,41 @@ Change log for SurfaceTopography ================================ +v1.23.0 (not yet released) +-------------------------- + +- BUG: Transposed and downsampled topographies report correct `pixel_size`; affected derivatives and PSDs for anisotropic pixels +- BUG: Masked data is normalized by the number of defined points; changes `rms_height` and moments of files with undefined data +- BUG: Power spectrum no longer halves the `q=0` entry when folding the +q and -q branches (#282) +- BUG: Fixed reshape order for non-square scans in the DI, EZD, MI and PS readers +- BUG: MI reader converts lateral sizes from meters to the requested unit +- BUG: BCR reader treats `voidpixels` as a count and masks against format markers +- BUG: Fixed channel selection in the NMM and OIR readers +- BUG: Fixed EZD magic check (bytes vs. str, inverted logic) +- BUG: ZAG container reader re-opens the ZIP file for each read +- BUG: Rewrote the C++ moment kernels in exact polynomial form (sign errors, missing prefactors, division by zero) +- BUG: `Bicubic.__call__` handles non-contiguous input arrays +- BUG: Fixed bounds of the periodic 2D bearing area +- BUG: `to_uniform()` interpolates at the correct positions for scans not starting at zero +- BUG: Slope detrending uses the length-weighted slope +- BUG: Hann window in the nonuniform power spectrum is origin invariant +- BUG: `fourier_synthesis` zeroes the DC mode instead of setting it to `C(q=1)` +- BUG: `make_sphere(..., periodic=True)` passes the periodic flag on +- BUG: Fixed `direction` semantics of `scan_line_align` +- BUG: Nonuniform `polyfit` solves in scan-centered coordinates; detrending no longer depends on the x origin +- BUG: Fixed Nyquist handling in Fourier interpolation and for anisotropic pixels in `power_spectrum_from_area` +- BUG: Registering an analysis function on a subclass no longer modifies the base class +- BUG: Fixed scale-dependent curvature on log-spaced and reliability-trimmed grids +- ENH: Bicubic interpolation computes spline coefficients on demand, instead of 128 bytes per pixel up front +- ENH: Container PSD integration precomputes bandwidth intervals in one pass, instead of re-reading every file per topography +- ENH: The IBW, MI and MNT readers no longer keep file data in memory after construction +- ENH: Nonuniform height-height autocorrelation moved to a C++ kernel +- ENH: Rigid-sphere scan vectorized over blocks of bounded memory +- ENH: Gaussian process regression reuses a single Cholesky factorization +- ENH: Flood-fill stack in the patch finder is allocated once per call, not once per patch +- TST: Added property-based invariant tests for transposition, translation, scaling and origin invariance +- MAINT: Removed broken `pytest-flake8` from the test dependencies; renamed `tiffile` to `tifffile` + v1.22.0 (08May26) ----------------- diff --git a/SurfaceTopography/Container/Averaging.py b/SurfaceTopography/Container/Averaging.py index 8d96ad10..1a753d65 100644 --- a/SurfaceTopography/Container/Averaging.py +++ b/SurfaceTopography/Container/Averaging.py @@ -34,7 +34,7 @@ from ..Support.Regression import resample from .SurfaceContainer import SurfaceContainer -_log = logging.Logger(__name__) +_log = logging.getLogger(__name__) def log_average(self, function_name, unit, nb_points_per_decade=10, reliable=True, progress_callback=None, **kwargs): diff --git a/SurfaceTopography/Container/IO/CE.py b/SurfaceTopography/Container/IO/CE.py index 3546e15c..370e89f8 100644 --- a/SurfaceTopography/Container/IO/CE.py +++ b/SurfaceTopography/Container/IO/CE.py @@ -77,7 +77,18 @@ def __init__(self, zipname, filename): self._filename = filename def __call__(self): - return ZipFile(self._zipname, mode="r").open(self._filename, mode="r") + zipfile = ZipFile(self._zipname, mode="r") + stream = zipfile.open(self._filename, mode="r") + # Close the containing ZipFile together with the member stream; + # otherwise its file descriptor leaks until garbage collection + original_close = stream.close + + def close(): + original_close() + zipfile.close() + + stream.close = close + return stream class CEReader(ContainerReaderBase): diff --git a/SurfaceTopography/Container/IO/ZAG.py b/SurfaceTopography/Container/IO/ZAG.py index 17564889..b790aabf 100644 --- a/SurfaceTopography/Container/IO/ZAG.py +++ b/SurfaceTopography/Container/IO/ZAG.py @@ -57,7 +57,46 @@ from ..SurfaceContainer import LazySurfaceContainer from .Reader import ContainerReaderBase -_log = logging.Logger(__name__) +_log = logging.getLogger(__name__) + + +class ZAGFileOpener(object): + """ + Callable that (re)opens a member of the ZIP archive embedded in a ZAG + file. Readers constructed from ZAG containers hold on to this callable + instead of an open stream, so that reading topographies still works + after `ZAGReader.close` has been called (e.g. when the reader was used + as a context manager by `read_container`). + """ + + def __init__(self, fobj, filename): + self._fobj = fobj # file name or stream + self._filename = filename + + def __call__(self): + # `ZipFile` automatically skips the ZAG header (magic + BMP + # thumbnail) preceding the ZIP archive, since ZIP archives are + # located through their central directory at the end of the file. + if hasattr(self._fobj, "read"): + zipfile = ZipFile(self._fobj, "r") + close_files = [zipfile] + else: + rawfile = open(self._fobj, "rb") + zipfile = ZipFile(rawfile, "r") + close_files = [zipfile, rawfile] + stream = zipfile.open(self._filename, "r") + # Close the containing ZipFile (and raw file) together with the + # member stream; otherwise file descriptors leak until garbage + # collection + original_close = stream.close + + def close(): + original_close() + for f in close_files: + f.close() + + stream.close = close + return stream class ZAGReader(ContainerReaderBase): @@ -113,6 +152,7 @@ def __init__(self, fobj): from ...IO.common import OpenFromAny # Open if a file name is given + self._fobj = fobj if not hasattr(fobj, "read"): # This is a string self._fstream = open(fobj, "rb") @@ -155,10 +195,16 @@ def __init__(self, fobj): # Lazy import to avoid circular dependency during package initialization from ...IO import ZONReader data_path = data.find(self._PATH_TAG).text + # We pass a callable that reopens the ZIP member + # rather than an open stream: the stream `f` + # is closed when this reader is closed, but the + # lazy container must be able to read + # topographies after that (see `read_container`). readers += [ ZONReader( - z.open( - f"{data_uuid}/{data_path}/{self._ZON_UUID}", "r" + ZAGFileOpener( + self._fobj, + f"{data_uuid}/{data_path}/{self._ZON_UUID}", ) ).topography ] diff --git a/SurfaceTopography/Container/IO/__init__.py b/SurfaceTopography/Container/IO/__init__.py index 4e4ab3f0..5562632b 100644 --- a/SurfaceTopography/Container/IO/__init__.py +++ b/SurfaceTopography/Container/IO/__init__.py @@ -68,7 +68,7 @@ def detect_format(fobj): raise CannotDetectFileFormat(msg) -def open_container(fobj, format=None): +def open_container(fobj, format=None, **kwargs): r""" Returns a container reader object for the file `fobj`. @@ -79,6 +79,9 @@ def open_container(fobj, format=None): format : str, optional Specify in which format the file should be interpreted. (Default: None, which means autodetect file format) + **kwargs : dict + Additional keyword arguments passed to the reader's constructor + (e.g. `datafile_keys` or `ignore_filters` of `CEReader`). Returns ------- @@ -87,13 +90,13 @@ def open_container(fobj, format=None): """ # noqa: E501 if not hasattr(fobj, 'read'): # fobj is a path if not os.path.isfile(fobj): - raise FileExistsError("file {} not found".format(fobj)) + raise FileNotFoundError("file {} not found".format(fobj)) if format is None: msg = "" for reader in readers: try: - return reader(fobj) + return reader(fobj, **kwargs) except Exception as err: msg += "tried {}: \n {}\n\n".format(reader.__name__, err) finally: @@ -107,7 +110,7 @@ def open_container(fobj, format=None): raise UnknownFileFormat( "{} not in registered container file formats {}".format( fobj, lookup_reader_by_format.keys())) - return lookup_reader_by_format[format](fobj) + return lookup_reader_by_format[format](fobj, **kwargs) def read_container(fn, format=None, **kwargs): @@ -116,8 +119,14 @@ def read_container(fn, format=None, **kwargs): Parameters ---------- - fobj : str or filelike object + fn : str or filelike object Path of the file or file-like object. + format : str, optional + Specify in which format the file should be interpreted. + (Default: None, which means autodetect file format) + **kwargs : dict + Additional keyword arguments passed to the reader's constructor + (e.g. `datafile_keys` or `ignore_filters` of `CEReader`). Returns ------- @@ -125,9 +134,9 @@ def read_container(fn, format=None, **kwargs): A list of container objects. """ containers = [] - with open_container(fn, format=format) as reader: + with open_container(fn, format=format, **kwargs) as reader: for i in range(reader.nb_containers): - containers += [reader.container(index=i, **kwargs)] + containers += [reader.container(index=i)] return containers @@ -147,12 +156,21 @@ def read_published_container(publication_url, **request_args): container : SurfaceContainer Surface container object read from URL. """ + # Default to a finite timeout so a stale server cannot hang the caller + # indefinitely; can be overridden through `request_args` + request_args.setdefault('timeout', 60) + # If we send json as a request header, then contact.engineering will response with a JSON dictionary - response = requests.get(publication_url, headers={'Accept': 'application/json'}) + custom_headers = request_args.pop('headers', {}) + response = requests.get(publication_url, + headers={'Accept': 'application/json', **custom_headers}, + **request_args) + response.raise_for_status() data = response.json() download_url = data['download_url'] # Then download and read container - container_response = requests.get(download_url, **request_args) + container_response = requests.get(download_url, headers=custom_headers, **request_args) + container_response.raise_for_status() container_file = io.BytesIO(container_response.content) return read_container(container_file) diff --git a/SurfaceTopography/Container/Integration.py b/SurfaceTopography/Container/Integration.py index f01d6ab0..8f8eaa7a 100644 --- a/SurfaceTopography/Container/Integration.py +++ b/SurfaceTopography/Container/Integration.py @@ -28,16 +28,18 @@ from ..Generic.Moments import compute_1d_moment, compute_iso_moment -def _bandwidth_count_from_profile(self, qx, unit, reliable=True): +def _bandwidth_intervals_from_profile(self, unit, reliable=True): r""" - Return number of topographies that include qx in their bandwidth. + Return the wavevector intervals covered by the topographies of this + container as a pair of sorted arrays (lower bounds, upper bounds). + This requires a single pass over the container; counting queries can + then be answered without touching the topographies again (lazy + containers read the data file on every element access). Parameters: ----------- self : SurfaceContainer Collection of Height containers - qx: np.ndarrau of floats - wavevector unit : str Unit of lengths in which the wavevector is defined. reliable : bool, optional @@ -45,11 +47,13 @@ def _bandwidth_count_from_profile(self, qx, unit, reliable=True): Returns ------- - number: np.ndarray - number of topographies having qx in their bandwidth + lower : np.ndarray + Sorted lower bandwidth bounds of the individual topographies. + upper : np.ndarray + Sorted upper bandwidth bounds of the individual topographies. """ - qx = np.abs(qx) - factor = np.zeros_like(qx) + lower = [] + upper = [] for t in self: t = t.to_unit(unit) @@ -58,14 +62,55 @@ def _bandwidth_count_from_profile(self, qx, unit, reliable=True): if short_cutoff is not None: qxmax = min(2 * np.pi / short_cutoff, qxmax) - factor += np.logical_and( - qx >= 2 * np.pi / t.physical_sizes[0], - qx <= qxmax, - ) + qxmin = 2 * np.pi / t.physical_sizes[0] + if qxmax >= qxmin: + lower += [qxmin] + upper += [qxmax] + # else: empty bandwidth interval (reliability cutoff longer than the + # scan); such a topography contains no wavevector and must not enter + # the counting arrays, where the subtraction in `_count_in_intervals` + # would tally it as -1 between its inverted bounds + return np.sort(lower), np.sort(upper) + + +def _count_in_intervals(intervals, nb_topographies, qx): + r""" + Return the number of bandwidth intervals (from + `_bandwidth_intervals_from_profile`) that contain each wavevector in + `qx`. + """ + lower, upper = intervals + qx = np.abs(qx) + # Number of intervals with lower <= qx minus number of intervals with + # upper < qx + count = np.searchsorted(lower, qx, side='right') - np.searchsorted(upper, qx, side='left') # All topographies have qx == 0 wavevector - factor = np.where(qx == 0, len(self), factor) - # TODO: This will be very slow on Topobank - return factor + return np.where(qx == 0, nb_topographies, count) + + +def _bandwidth_count_from_profile(self, qx, unit, reliable=True): + r""" + Return number of topographies that include qx in their bandwidth. + + Parameters: + ----------- + self : SurfaceContainer + Collection of Height containers + qx: np.ndarrau of floats + wavevector + unit : str + Unit of lengths in which the wavevector is defined. + reliable : bool, optional + Only incorporate data deemed reliable. (Default: True) + + Returns + ------- + number: np.ndarray + number of topographies having qx in their bandwidth + """ + return _count_in_intervals( + _bandwidth_intervals_from_profile(self, unit, reliable), len(self), qx + ) def integrate_psd_from_profile(self, factor, unit, window=None, reliable=True): @@ -84,7 +129,7 @@ def integrate_psd_from_profile(self, factor, unit, window=None, reliable=True): .. math:: - \frac{1}{2 \pi} \int_0^\infty dq_x factor(q_x) C^{1D}(q_x) + \frac{1}{2 \pi} \int_{-\infty}^\infty dq_x factor(q_x) C^{1D}(q_x) Discrete @@ -118,13 +163,16 @@ def integrate_psd_from_profile(self, factor, unit, window=None, reliable=True): """ integ = 0 - # TODO: faster: precompute _bandwidth_count_from_profile (defining a piecewise constant function, ) - # a table with bins and number of topographies inside the bins - # this will make us only loop 2N times through the topographies instead of N^2 - # with N the total number of topographies + # Precompute the bandwidth intervals in a single pass over the + # container; the `average` callback below is invoked once per + # topography, and looping over the container inside it would read + # every data file N times (lazy containers construct topographies + # from the file on each element access). + intervals = _bandwidth_intervals_from_profile(self, unit, reliable) + nb_topographies = len(self) def average(qx): - count = _bandwidth_count_from_profile(self, qx, unit, reliable) + count = _count_in_intervals(intervals, nb_topographies, qx) return np.where(count > 0, factor(qx) / count, 0) for t in self: diff --git a/SurfaceTopography/Container/ScaleDependentStatistics.py b/SurfaceTopography/Container/ScaleDependentStatistics.py index 9635e576..fa568981 100644 --- a/SurfaceTopography/Container/ScaleDependentStatistics.py +++ b/SurfaceTopography/Container/ScaleDependentStatistics.py @@ -32,7 +32,7 @@ from ..Exceptions import NoReliableDataError from .SurfaceContainer import SurfaceContainer -_log = logging.Logger(__name__) +_log = logging.getLogger(__name__) def scale_dependent_statistical_property( @@ -116,8 +116,8 @@ def scale_dependent_statistical_property( able to accept one argument (for line scans) and two arguments (for topographies). - >>> s = c.scale_dependent_statistical_property(lambda x, y=None: np.var(x), n=1, distance=[0.1, 1.0, 10], unit='um') - """ + >>> s = c.scale_dependent_statistical_property(lambda x, y=None: np.var(x), n=1, distances=[0.1, 1.0, 10], unit='um') + """ # noqa: E501 results = defaultdict(list) empty = None for i, topography in enumerate(self): @@ -156,8 +156,9 @@ def scale_dependent_statistical_property( len(existing_distances) ) else: - existing_distances = np.array(distances) - unique_distance_index = np.arange(len(distances)) + # `distances` may be a scalar + existing_distances = np.atleast_1d(np.asarray(distances, dtype=float)) + unique_distance_index = np.arange(len(existing_distances)) # For the factor n see 10.1016/j.apsadv.2021.100190 m = np.logical_and(existing_distances > n * lower, existing_distances < upper) existing_distances = existing_distances[m] @@ -172,11 +173,12 @@ def scale_dependent_statistical_property( distance=existing_distances, **kwargs, ) - # Append results to our return values - for i, e, s in zip(unique_distance_index, existing_distances, stat): + # Append results to our return values. (Note: do not reuse the + # outer loop variable `i` here.) + for j, e, s in zip(unique_distance_index, existing_distances, stat): if empty is None: empty = np.zeros_like(s) * np.nan - results[i] += [(e, s)] + results[j] += [(e, s)] else: _log.warning(f"Topography {topography} contributes no data to average.") @@ -189,6 +191,7 @@ def scale_dependent_statistical_property( if distances is not None: # If distances are specified by the user, we return exactly those distances; whether data actually exists is # indicated through the mask of a masked array + distances = np.atleast_1d(np.asarray(distances, dtype=float)) data = np.array( [ ( diff --git a/SurfaceTopography/Container/SurfaceContainer.py b/SurfaceTopography/Container/SurfaceContainer.py index fd6be761..db5dc5ad 100644 --- a/SurfaceTopography/Container/SurfaceContainer.py +++ b/SurfaceTopography/Container/SurfaceContainer.py @@ -26,6 +26,8 @@ import abc from functools import update_wrapper +from ..Support.Deprecation import deprecated as deprecation_warning + class SurfaceContainer(metaclass=abc.ABCMeta): """A list of topographies""" @@ -40,31 +42,48 @@ def __len__(self): def __getitem__(self, item): raise NotImplementedError + def _function_registry(self): + """ + Return the merged dictionary of functions registered on this class + and its bases. See `AbstractTopography._function_registry`. + """ + functions = {} + for klass in reversed(type(self).__mro__): + functions.update(klass.__dict__.get('_functions', {})) + return functions + def apply(self, name, *args, **kwargs): - self._functions[name](self, *args, **kwargs) + return self._function_registry()[name](self, *args, **kwargs) def __getattr__(self, name): - if name in self._functions: + functions = self._function_registry() + if name in functions: def func(*args, **kwargs): - return self._functions[name](self, *args, **kwargs) + return functions[name](self, *args, **kwargs) - update_wrapper(func, self._functions[name]) + update_wrapper(func, functions[name]) return func else: raise AttributeError( "Unkown attribute '{}' and no analysis or pipeline function of this name registered" "(class {}). Available functions: {}".format( - name, self.__class__.__name__, ", ".join(self._functions.keys()) + name, self.__class__.__name__, ", ".join(functions.keys()) ) ) def __dir__(self): - return sorted(super().__dir__() + [*self._functions]) + return sorted(super().__dir__() + [*self._function_registry()]) @classmethod - def register_function(cls, name, function): - cls._functions.update({name: function}) + def register_function(cls, name, function, deprecated=False): + if deprecated: + function = deprecation_warning()(function) + if '_functions' not in cls.__dict__: + # Copy on write: registering on a subclass must not leak into + # the shared base-class registry + cls._functions = {} + cls._functions[name] = function class InMemorySurfaceContainer(SurfaceContainer): diff --git a/SurfaceTopography/FFTTricks.py b/SurfaceTopography/FFTTricks.py index f876b68b..cc8c63a7 100644 --- a/SurfaceTopography/FFTTricks.py +++ b/SurfaceTopography/FFTTricks.py @@ -57,13 +57,18 @@ def make_fft(topography, communicator=None): RuntimeError If the muGrid FFTEngine object's domain decomposition does not match the topography's domain decomposition. """ - # We only initialize this once and attach it to the topography object - if hasattr(topography, '_mufft'): + if communicator is None: + communicator = topography.communicator + + # We only initialize this once and attach it to the topography object. + # Note: The cache is only valid for the communicator it was created + # with; a subsequent call with a different communicator must not return + # the cached engine. + if hasattr(topography, '_mufft') and topography._mufft_communicator is communicator: return topography._mufft if topography.is_domain_decomposed: - fft = muGrid.FFTEngine(topography.nb_grid_pts, - communicator=topography.communicator if communicator is None else communicator) + fft = muGrid.FFTEngine(topography.nb_grid_pts, communicator=communicator) if fft.subdomain_locations != topography.subdomain_locations or \ fft.nb_subdomain_grid_pts != topography.nb_subdomain_grid_pts: raise RuntimeError('muGrid suggested a domain decomposition that ' @@ -71,6 +76,7 @@ def make_fft(topography, communicator=None): else: fft = muGrid.FFTEngine(topography.nb_grid_pts) topography._mufft = fft + topography._mufft_communicator = communicator return fft diff --git a/SurfaceTopography/Generation.py b/SurfaceTopography/Generation.py index 52836275..fcbfa0f5 100644 --- a/SurfaceTopography/Generation.py +++ b/SurfaceTopography/Generation.py @@ -138,11 +138,21 @@ def self_affine_prefactor(nb_grid_pts, physical_sizes, Hurst, rms_height=None, area = np.prod(physical_sizes) if rms_height is not None: + if Hurst <= 0: + # The prefactor expression below has a 0/0 limit at Hurst = 0 + # (the rms height integral diverges logarithmically); evaluating + # it would silently produce an all-NaN topography + raise ValueError( + 'Scaling to a target rms height requires a positive Hurst exponent.') # Assuming no rolloff region fac = 2 * rms_height / np.sqrt(q_min ** (-2 * Hurst) - q_max ** (-2 * Hurst)) * np.sqrt( Hurst * np.pi) elif rms_slope is not None: + if Hurst >= 1: + # Same 0/0 limit at Hurst = 1 for the rms slope integral + raise ValueError( + 'Scaling to a target rms slope requires a Hurst exponent below unity.') fac = 2 * rms_slope / np.sqrt(q_max ** (2 - 2 * Hurst) - q_min ** (2 - 2 * Hurst)) * np.sqrt( (1 - Hurst) * np.pi) @@ -306,6 +316,16 @@ def fourier_synthesis(nb_grid_pts, physical_sizes, karr[x, mask] = rolloff * ran[mask] * q_min ** (-(1 + hurst)) else: karr[mask] = rolloff * ran[mask] * q_min ** (-(0.5 + hurst)) + if x == 0 and psd is None: + # The q=0 mode only controls the mean height, which we fix to + # zero. (Setting q_sq[0] to 1 above merely avoids a + # divide-by-zero; without zeroing this mode the surface would + # acquire a random mean offset of magnitude C(q=1), which is + # unit-dependent and typically dwarfs the rms height.) + if len(nb_grid_pts) == 2: + karr[0, 0] = 0. + else: + karr[0] = 0. if len(nb_grid_pts) == 2: for iy in [0, -1] if ny % 2 == 0 else [0]: # Enforce symmetry diff --git a/SurfaceTopography/Generic/Curvature.py b/SurfaceTopography/Generic/Curvature.py index e0f183a7..c9357dea 100644 --- a/SurfaceTopography/Generic/Curvature.py +++ b/SurfaceTopography/Generic/Curvature.py @@ -32,6 +32,34 @@ from ..HeightContainer import NonuniformLineScanInterface, UniformTopographyInterface +def _curvature_stencil(r, A): + r""" + Evaluate :math:`B(\lambda) = 8A(\lambda) - 2A(2\lambda)` from a sampled + autocorrelation function. + + The autocorrelation may be sampled on an arbitrary (e.g. log-spaced or + reliability-trimmed) grid, so :math:`A` is evaluated by linear + interpolation. This is exact where :math:`2\lambda` falls onto a grid + point, in particular everywhere on a linearly spaced grid that starts + at zero. + """ + valid = np.isfinite(r) & np.isfinite(A) + r = r[valid] + A = A[valid] + # We need A(2 lambda), i.e. we can only evaluate B up to half the + # largest sampled distance; distance zero is excluded because the + # curvature expression divides by lambda^2 + mask = np.logical_and(r > 0, 2 * r <= r[-1]) + rm = r[mask] + B = 8 * np.interp(rm, r, A) - 2 * np.interp(2 * rm, r, A) + # Truncate at the first negative value + nz = np.nonzero(B < 0)[0] + if len(nz) > 0: + rm = rm[:nz[0]] + B = B[:nz[0]] + return rm, B + + def scale_dependent_curvature_from_profile(topography, **kwargs): r""" Compute the one-dimensional scale-dependent curvature. @@ -62,14 +90,8 @@ def scale_dependent_curvature_from_profile(topography, **kwargs): Curvature. (Units: 1/length) """ # noqa: E501 r, A = topography.autocorrelation_from_profile(**kwargs) - n = (len(r) + 1) // 2 - r = r[1:n] - B = 8 * A[1:n] - 2 * A[2::2] - nz = np.nonzero(B < 0)[0] - if len(nz) > 0: - n = nz[0] - # Important: The following expression relies on the fact that r is equally spaced! - return r[:n], np.sqrt(B[:n]) / r[:n] ** 2 + r, B = _curvature_stencil(r, A) + return r, np.sqrt(B) / r ** 2 def scale_dependent_curvature_from_area(topography, **kwargs): @@ -102,14 +124,8 @@ def scale_dependent_curvature_from_area(topography, **kwargs): Curvature. (Units: 1/length) """ # noqa: E501 r, A = topography.autocorrelation_from_area(**kwargs) - n = (len(r) + 1) // 2 - r = r[1:n] - B = 8 * A[1:n] - 2 * A[2::2] - nz = np.nonzero(B < 0)[0] - if len(nz) > 0: - n = nz[0] - # Important: The following expression relies on the fact that r is equally spaced! - return 2 * r[:n], np.sqrt(B[:n]) / r[:n] ** 2 + r, B = _curvature_stencil(r, A) + return 2 * r, np.sqrt(B) / r ** 2 # Register analysis functions from this module diff --git a/SurfaceTopography/Generic/Moments.py b/SurfaceTopography/Generic/Moments.py index f1b63a52..1a328532 100644 --- a/SurfaceTopography/Generic/Moments.py +++ b/SurfaceTopography/Generic/Moments.py @@ -27,11 +27,14 @@ The example below computes the PSD of a surface topography and integrates it to obtain the RMS height. +>>> import numpy as np >>> from SurfaceTopography.Generation import fourier_synthesis ->>> t = fourier_synthesis((256, 256), (256,256), hurst=0.8, rms_height=1, short_cutoff=4, long_cutoff=64) +>>> from SurfaceTopography.Generic.Moments import compute_1d_moment +>>> t = fourier_synthesis((256, 256), (256, 256), hurst=0.8, rms_height=1, short_cutoff=4, long_cutoff=64) >>> hrms = t.rms_height_from_profile() ->>> hrms_from_psd = np.sqrt(compute_1d_moment(,) ->>> assert abs(hrms_from_psd / hrms - 1) < 0.1 +>>> q, C = t.power_spectrum_from_profile(reliable=False) +>>> hrms_from_psd = np.sqrt(compute_1d_moment(q, C, order=0)) +>>> assert abs(hrms_from_psd / hrms - 1) < 0.1 However, because the averaging of the PSD before integration introduces errors, it is better to sum directly the raw spectum, @@ -49,9 +52,12 @@ def compute_1d_moment(x, y, order=1, cumulative=False): .. math:: - m_\alpha = \int dx y x^{\alpha} + m_\alpha = \frac{1}{\pi} \int dx\, y x^{\alpha} - using trapezoidal interpolation of the integrand + using trapezoidal interpolation of the integrand. (The prefactor + :math:`1/\pi` corresponds to the moment of a two-sided function + :math:`y(|x|)` with the :math:`1/(2\pi)` Fourier-transform convention + used for power spectra in this package.) Parameters: ----------- @@ -89,9 +95,12 @@ def compute_iso_moment(x, y, order=1, cumulative=False): .. math:: - m_\alpha = \int dx_1 dx_2 y_{2D}(x_1, x_2) |\vec x|^{\alpha} = 2\pi \int dx y(x) x^{\alpha + 1} + m_\alpha = \frac{1}{(2\pi)^2} \int dx_1 dx_2\, y_{2D}(x_1, x_2) |\vec x|^{\alpha} + = \frac{1}{2\pi} \int dx\, y(x) x^{\alpha + 1} - using trapezoidal interpolation of the integrand + using trapezoidal interpolation of the integrand. (The prefactor + corresponds to the :math:`1/(2\pi)^2` Fourier-transform convention used + for power spectra in this package.) Parameters: ----------- diff --git a/SurfaceTopography/Generic/ScaleDependentStatistics.py b/SurfaceTopography/Generic/ScaleDependentStatistics.py index 64c2c66b..29b6e05f 100644 --- a/SurfaceTopography/Generic/ScaleDependentStatistics.py +++ b/SurfaceTopography/Generic/ScaleDependentStatistics.py @@ -128,7 +128,13 @@ def scale_dependent_statistical_property(self, func, n=1, scale_factor=None, dis d = self.derivative(n=n, scale_factor=scale_factor, distance=distance, interpolation=interpolation, progress_callback=progress_callback) if distance is None: - distance = scale_factor * np.mean(self.physical_sizes) + # The scale factor scales the stencil of the derivative, i.e. the + # physical distance is scale_factor * n * pixel_size (see + # `Uniform.Derivative.derivative`, which converts distances to scale + # factors via scale_factor = distance / (n * px)). Note that scale + # factors are only supported for uniform topographies, so + # `pixel_size` is guaranteed to exist here. + distance = np.asarray(scale_factor) * n * np.mean(self.pixel_size) if self.dim == 1: try: if scale_factor is not None: @@ -159,6 +165,7 @@ def scale_dependent_statistical_property(self, func, n=1, scale_factor=None, dis retvals = np.array([func(dx, dy) for dx, dy in zip(*d)]) except TypeError: # If there is a single distance, the reliability analysis is skipped + d = list(d) # may be a tuple, which does not support item assignment for i in range(0, len(d)): if len(d[i]) < threshold: d[i] = np.nan @@ -166,6 +173,7 @@ def scale_dependent_statistical_property(self, func, n=1, scale_factor=None, dis short_cutoff = self.short_reliability_cutoff() if reliable else None if short_cutoff is not None: + distance = np.asarray(distance) mask = distance > short_cutoff * n / 2 if mask.sum() == 0: raise NoReliableDataError('Dataset contains no reliable data.') diff --git a/SurfaceTopography/Generic/ScanningProbe.py b/SurfaceTopography/Generic/ScanningProbe.py index d2a3b0e2..192686b1 100644 --- a/SurfaceTopography/Generic/ScanningProbe.py +++ b/SurfaceTopography/Generic/ScanningProbe.py @@ -64,8 +64,9 @@ def scanning_probe_reliability_cutoff(self, tip_radius, safety_factor=1 / 2, xto used to compute the reliability cutoff. """ lower, upper = self.bandwidth() - # We need to normalize the bracket search to avoid numerical issues - fac = np.exp((np.log(lower) + np.log(upper) / 2)) + # We need to normalize the bracket search to avoid numerical issues; + # this is the geometric mean of the two bandwidth limits + fac = np.exp((np.log(lower) + np.log(upper)) / 2) target_curvature = fac * safety_factor / tip_radius diff --git a/SurfaceTopography/Generic/Slope.py b/SurfaceTopography/Generic/Slope.py index 0f792a7b..46438e6d 100644 --- a/SurfaceTopography/Generic/Slope.py +++ b/SurfaceTopography/Generic/Slope.py @@ -100,8 +100,9 @@ def scale_dependent_slope_from_area(topography, **kwargs): return r[1:], np.sqrt(2 * A[1:]) / r[1:] -# Register analysis functions from this module +# Register analysis functions from this module. Note: the areal variant is +# not registered for nonuniform line scans, which have no +# `autocorrelation_from_area`. UniformTopographyInterface.register_function('scale_dependent_slope_from_profile', scale_dependent_slope_from_profile) NonuniformLineScanInterface.register_function('scale_dependent_slope_from_profile', scale_dependent_slope_from_profile) UniformTopographyInterface.register_function('scale_dependent_slope_from_area', scale_dependent_slope_from_area) -NonuniformLineScanInterface.register_function('scale_dependent_slope_from_area', scale_dependent_slope_from_area) diff --git a/SurfaceTopography/HeightContainer.py b/SurfaceTopography/HeightContainer.py index d8fc038e..ff2262e1 100644 --- a/SurfaceTopography/HeightContainer.py +++ b/SurfaceTopography/HeightContainer.py @@ -37,6 +37,7 @@ from .Metadata import InfoModel from .Support import doi +from .Support.Deprecation import deprecated as deprecation_warning class AbstractTopography(object): @@ -64,27 +65,46 @@ def __init__(self, unit=None, info={}, communicator=MPI.COMM_WORLD): self._info = InfoModel(**info) self._communicator = communicator + def _function_registry(self): + """ + Return the dictionary of registered analysis and pipeline functions + for this object. Functions registered on any class in the MRO are + visible; functions registered on a subclass are confined to that + subclass (e.g. functions registered on `Topography` do not appear + on line scans). An instance-level `_functions` attribute overrides + the class registry (used by converter classes that redirect + dispatch). + """ + try: + return self.__dict__['_functions'] + except KeyError: + functions = {} + for klass in reversed(type(self).__mro__): + functions.update(klass.__dict__.get('_functions', {})) + return functions + def apply(self, name, *args, **kwargs): - self._functions[name](self, *args, **kwargs) + return self._function_registry()[name](self, *args, **kwargs) def __getattr__(self, name): - if name in self._functions: + functions = self._function_registry() + if name in functions: def func(*args, **kwargs): - return self._functions[name](self, *args, **kwargs) + return functions[name](self, *args, **kwargs) - update_wrapper(func, self._functions[name]) + update_wrapper(func, functions[name]) return func else: raise AttributeError( "Unkown attribute '{}' and no analysis or pipeline function of this name registered (class {}). " "Available functions: {}".format( - name, self.__class__.__name__, ", ".join(self._functions.keys()) + name, self.__class__.__name__, ", ".join(functions.keys()) ) ) def __dir__(self): - return sorted(super().__dir__() + [*self._functions]) + return sorted(super().__dir__() + [*self._function_registry()]) def __getstate__(self): """ @@ -183,10 +203,13 @@ def __setstate__(self, state): @property def info(self) -> dict: - """Return info dictionary""" - return self.parent_topography._info.model_copy(update=self._info).model_dump( - exclude_none=True - ) + """ + Return the info dictionary of the parent topography, updated with + the entries of this decorator. + """ + info = self.parent_topography.info + info.update(self._info.model_dump(exclude_none=True)) + return info @property def nb_subdomain_grid_pts(self): @@ -199,13 +222,33 @@ def pipeline(self): class TopographyInterface(object): @classmethod def register_function(cls, name, function, deprecated=False): # noqa: N805 - # FIXME! Wrap in warnings.deprecated decorator, will be introduced in Python 3.13 - if function.__name__ != "func_with_doi": + if deprecated: + function = deprecation_warning()(function) + if not getattr(function, '__has_doi__', False): # We want the `dois` argument for all pipeline functions. If no # doi has been specified, we simply wrap it in an empty decorator. - cls._functions.update({name: doi()(function)}) - else: - cls._functions.update({name: function}) + # (Note: the doi decorator marks its wrappers with `__has_doi__`; + # checking the function name would be defeated by + # functools.wraps.) + function = doi()(function) + if '_functions' not in cls.__dict__: + # Copy on write: registering a function on a subclass (e.g. a + # 2D-only analysis on `Topography`) must not leak into the + # shared registry of the base class, where it would become + # visible on line scans as well. + cls._functions = {} + cls._functions[name] = function + + @classmethod + def _all_functions(cls): # noqa: N805 + """ + Return the merged registry of analysis and pipeline functions + registered on this class and its bases. + """ + functions = {} + for klass in reversed(cls.__mro__): + functions.update(klass.__dict__.get('_functions', {})) + return functions class UniformTopographyInterface(TopographyInterface, metaclass=abc.ABCMeta): @@ -283,6 +326,10 @@ def positions_and_heights(self, **kwargs): return p, h def __eq__(self, other): + if not isinstance(other, UniformTopographyInterface): + return NotImplemented + if self.nb_grid_pts != other.nb_grid_pts: + return False return Reduction(self._communicator).all( self.unit == other.unit and self.info == other.info @@ -291,6 +338,12 @@ def __eq__(self, other): and np.allclose(self.heights(), other.heights()) ) + # Height containers compare by value but are mutable in principle; + # identity-based hashing nevertheless allows them to be used in sets + # and as dictionary keys. (Defining `__eq__` without `__hash__` would + # make them unhashable.) + __hash__ = object.__hash__ + def __getitem__(self, i): return self.heights()[i] @@ -304,7 +357,11 @@ def is_uniform(self): @property def is_reentrant(self): - return np.min(np.diff(self.positions())) <= 0 + positions = self.positions() + if len(positions) < 2: + # A line scan with less than two points cannot be reentrant + return False + return np.min(np.diff(positions)) <= 0 @property @abc.abstractmethod @@ -346,6 +403,10 @@ def has_undefined_data(self): return False def __eq__(self, other): + if not isinstance(other, NonuniformLineScanInterface): + return NotImplemented + if self.nb_grid_pts != other.nb_grid_pts: + return False return Reduction(self._communicator).all( self.unit == other.unit and self.info == other.info @@ -353,5 +414,8 @@ def __eq__(self, other): and np.allclose(self.positions_and_heights(), other.positions_and_heights()) ) + # See note on `UniformTopographyInterface.__hash__` + __hash__ = object.__hash__ + def __getitem__(self, i): return self.positions()[i], self.heights()[i] diff --git a/SurfaceTopography/IO/AL3D.py b/SurfaceTopography/IO/AL3D.py index 11417712..0a7b07c5 100644 --- a/SurfaceTopography/IO/AL3D.py +++ b/SurfaceTopography/IO/AL3D.py @@ -96,13 +96,21 @@ def read_height_data(self, f): invalid_pixel_value = float(self._header['InvalidPixelValue']) dtype = np.single nx, ny = int(self._header['Cols']), int(self._header['Rows']) + # The row stride is already in bytes (rows are padded to multiples + # of 8 bytes); do not multiply by the item size again, which would + # read four times too much data (into the texture image and beyond) rowstride = (nx * np.dtype(dtype).itemsize + 7) // 8 * 8 - buffer = f.read(rowstride * ny * np.dtype(dtype).itemsize) + buffer = f.read(rowstride * ny) data = as_strided(np.frombuffer(buffer, dtype=dtype), shape=(ny, nx), strides=(rowstride, np.dtype(dtype).itemsize)) mask = np.isnan(data) if not np.isnan(invalid_pixel_value): - mask = np.logical_or(mask, np.abs(data - invalid_pixel_value) < self._INVALID_RELTOL * invalid_pixel_value) + # Note: the tolerance must be scaled with the magnitude of the + # marker; without abs() a negative marker would make the mask + # unconditionally false + mask = np.logical_or( + mask, + np.abs(data - invalid_pixel_value) < self._INVALID_RELTOL * np.abs(invalid_pixel_value)) return np.ma.masked_array(data.T, mask=mask.T) @property diff --git a/SurfaceTopography/IO/BCR.py b/SurfaceTopography/IO/BCR.py index 5bdfe184..f606e9b1 100644 --- a/SurfaceTopography/IO/BCR.py +++ b/SurfaceTopography/IO/BCR.py @@ -214,9 +214,24 @@ def topography(self, channel_index=None, physical_sizes=None, # or pcolormesh(t.heights().T) for origin in lower left and # with inverted y axis (cartesian coordinate system) - invalid_pixel_value = float(self._metadata['voidpixels']) + # `voidpixels` holds the *number* of void (undefined) pixels in the + # file, not the marker value; if the key is missing, the file + # contains no void pixels. Void pixels are marked by the maximum + # value of the respective data type: 32767 for 16-bit integer + # (bcrstm) files and single-precision FLT_MAX (~3.4028235e38) for + # floating-point (bcrf) files. + try: + nb_void_pixels = float(self._metadata.get('voidpixels', 0)) + except ValueError: + nb_void_pixels = 0 + if nb_void_pixels > 0: + if self._file_type == 'bcrstm': + mask = data == 32767 + else: + mask = data >= 3.4028e38 + data = np.ma.masked_array(data, mask=mask) topography = Topography( - np.ma.masked_array(data, mask=data == invalid_pixel_value), + data, physical_sizes=(sx, sy), unit=channel.unit, info=_info, diff --git a/SurfaceTopography/IO/DATX.py b/SurfaceTopography/IO/DATX.py index f15e6223..88814b3d 100644 --- a/SurfaceTopography/IO/DATX.py +++ b/SurfaceTopography/IO/DATX.py @@ -116,7 +116,7 @@ def recursive_replace(raw_metadata, context): category = x_converter[0] unit = x_converter[1] values = x_converter[-1] - if category != b'LateralCat' and unit != b'Pixels': + if category != b'LateralCat' or unit != b'Pixels': raise UnsupportedFormatFeature('DATX reader only supports `LateralCat` with `Pixels` unit for ' 'X converter.') physical_sizes_x = self._nb_grid_pts[0] * values[1] # in units of meters! @@ -124,7 +124,7 @@ def recursive_replace(raw_metadata, context): category = y_converter[0] unit = y_converter[1] values = y_converter[-1] - if category != b'LateralCat' and unit != b'Pixels': + if category != b'LateralCat' or unit != b'Pixels': raise UnsupportedFormatFeature('DATX reader only supports `LateralCat` with `Pixels` unit for ' 'Y converter.') physical_sizes_y = self._nb_grid_pts[1] * values[1] # in units of meters! @@ -133,7 +133,7 @@ def recursive_replace(raw_metadata, context): height_unit = z_converter[1] values = z_converter[-1] if category != b'HeightCat': - raise UnsupportedFormatFeature('DATX reader only supports `LateralCat` for Z converter.') + raise UnsupportedFormatFeature('DATX reader only supports `HeightCat` for Z converter.') height_unit = mangle_length_unit_utf8(height_unit.decode('ascii')) if height_unit != self._unit: diff --git a/SurfaceTopography/IO/DI.py b/SurfaceTopography/IO/DI.py index 7b64c3f7..45ce4ec4 100644 --- a/SurfaceTopography/IO/DI.py +++ b/SurfaceTopography/IO/DI.py @@ -305,8 +305,13 @@ def topography( fobj.seek(offset) rawdata = fobj.read(nx * ny * dtype.itemsize) - unscaleddata = np.frombuffer(rawdata, count=nx * ny, dtype=dtype).reshape( - nx, ny + # The data is stored line by line, i.e. the buffer has C-order + # shape (ny, nx). Transposing yields the (nx, ny) array with the + # x index first that `Topography` expects. + unscaleddata = ( + np.frombuffer(rawdata, count=nx * ny, dtype=dtype) + .reshape(ny, nx) + .T ) # internal information from file @@ -323,7 +328,7 @@ def topography( # with inverted y axis (cartesian coordinate system) surface = Topography( - np.fliplr(unscaleddata.T), + np.fliplr(unscaleddata), physical_sizes=(sx, sy), unit=channel.unit, info=_info, diff --git a/SurfaceTopography/IO/DZI.py b/SurfaceTopography/IO/DZI.py index c9eb508c..f4e6fffc 100644 --- a/SurfaceTopography/IO/DZI.py +++ b/SurfaceTopography/IO/DZI.py @@ -151,7 +151,12 @@ def write_data(fn, subdata, pixel_physical_sizes, physical_offsets): # Get heights and rescale to interval 0, 1 mx, mn = data.max(), data.min() - data = (data - mn) / (mx - mn) + if mx > mn: + data = (data - mn) / (mx - mn) + else: + # A perfectly flat topography would otherwise divide by zero and + # produce an all-NaN image pyramid + data = np.zeros_like(data) # Write configuration file if meta_format == "xml": diff --git a/SurfaceTopography/IO/EZD.py b/SurfaceTopography/IO/EZD.py index 296a3338..84be41cc 100644 --- a/SurfaceTopography/IO/EZD.py +++ b/SurfaceTopography/IO/EZD.py @@ -68,7 +68,7 @@ def __init__(self, file_path): with OpenFromAny(self._file_path, "rb") as fobj: metadata = {} line = fobj.readline() - if line == self._MAGIC: + if line != self._MAGIC.encode("latin-1"): raise FileFormatMismatch("This is not a NanoSurf easyScan data file") section_name = "DataSet" @@ -218,8 +218,13 @@ def topography( fobj.seek(self._start_of_data + offset) rawdata = fobj.read(nx * ny * dtype.itemsize) - unscaleddata = np.frombuffer(rawdata, count=nx * ny, dtype=dtype).reshape( - nx, ny + # The data is stored line by line, i.e. the buffer has C-order + # shape (ny, nx). Transposing yields the (nx, ny) array with the + # x index first that `Topography` expects. + unscaleddata = ( + np.frombuffer(rawdata, count=nx * ny, dtype=dtype) + .reshape(ny, nx) + .T ) # internal information from file @@ -236,7 +241,7 @@ def topography( # with inverted y axis (cartesian coordinate system) topography = Topography( - np.fliplr(unscaleddata.T), + np.fliplr(unscaleddata), physical_sizes=(sx, sy), unit=channel.unit, info=_info, diff --git a/SurfaceTopography/IO/FRT.py b/SurfaceTopography/IO/FRT.py index 12c79405..9d0252ca 100644 --- a/SurfaceTopography/IO/FRT.py +++ b/SurfaceTopography/IO/FRT.py @@ -457,7 +457,7 @@ class FRTReader(ReaderBase): ('filter', 'I'), ('reference_type', 'I'), ('layer_stack_id', 'I'), - ('reference_material_id', 'gint32'), + ('reference_material_id', 'i'), ('reference_constant', 'd'), ('material_thickness', 'd'), ], @@ -513,7 +513,7 @@ def __init__(self, file_path): block_id, block_size = unpack('= 1000: # Only substantial blocks - # Extract block prefix information - element_offset = 0 - elements_per_block = 0 - if i >= 16: - prefix = compressed_blocks_raw[i - 16:i] - element_offset = struct.unpack("= 1000: # Only substantial blocks + # Extract block prefix information + element_offset = 0 + elements_per_block = 0 + if i >= 16: + prefix = compressed_blocks_raw[i - 16:i] + element_offset = struct.unpack(" MagicMatch: if len(buffer) < 4: @@ -121,7 +112,7 @@ def can_read(cls, buffer: bytes) -> MagicMatch: ("part_thickness", ">f"), ("sw_llc", ">h"), ("target_range", ">f"), - ("rad_crv_veasure_eeq", "h"), ("min_mod", ">i"), ("min_mod_count", ">i"), ("phase_res", ">h"), @@ -159,84 +150,84 @@ def can_read(cls, buffer: bytes) -> MagicMatch: ("fda_res", ">h"), ("scan_descr", "20s"), ("n_fiducials_a", ">h"), - ("fiducials_a", "14f"), + ("fiducials_a", ">14f"), ("pixel_width", ">f"), ("pixel_height", ">f"), ("exit_pupil_diam", ">f"), ("light_level_pct", ">f"), - ("coords_state", "i"), + ("coords_x_pos", ">f"), + ("coords_y_pos", ">f"), + ("coords_z_pos", ">f"), + ("coords_x_rot", ">f"), + ("coords_y_rot", ">f"), + ("coords_z_rot", ">f"), + ("coherence_mode", ">h"), + ("surface_filter", ">h"), ("sys_err_file_name", "28s"), ("zoom_descr", "8s"), - ("alpha_part", "f"), + ("beta_part", ">f"), + ("dist_part", ">f"), + ("cam_split_loc_x", ">h"), + ("cam_split_loc_y", ">h"), + ("cam_split_trans_x", ">h"), + ("cam_split_trans_y", ">h"), ("material_a", "24s"), ("material_b", "24s"), - ("cam_split_unused", "h"), (None, "2b"), - ("dmi_ctr_x", "f"), + ("dmi_ctr_y", ">f"), + ("sph_dist_corr", ">h"), (None, "2b"), - ("sph_dist_part_na", "f"), + ("sph_dist_part_radius", ">f"), + ("sph_dist_cal_na", ">f"), + ("sph_dist_cal_radius", ">f"), + ("surface_type", ">h"), + ("ac_surface_type", ">h"), + ("z_position", ">f"), + ("power_multiplier", ">f"), + ("focus_multiplier", ">f"), + ("rad_crv_vocus_sal_lactor", ">f"), + ("rad_crv_vower_ral_lactor", ">f"), + ("ftp_left_pos", ">f"), + ("ftp_right_pos", ">f"), + ("ftp_pitch_pos", ">f"), + ("ftp_roll_pos", ">f"), + ("min_mod_pct", ">f"), + ("max_inten", ">i"), + ("ring_of_fire", ">h"), (None, "1b"), ("rc_orientation", "B"), - ("rc_distance", "f"), + ("rc_angle", ">f"), + ("rc_diameter", ">f"), ("rem_fringes_mode", ">h"), (None, "1b"), ("ftpsi_phase_res", "B"), - ("frames_acquired", "h"), + ("cavity_type", ">h"), + ("cam_frame_rate", ">f"), + ("tune_range", ">f"), + ("cal_pix_loc_x", ">h"), + ("cal_pix_loc_y", ">h"), + ("n_tst_cal_pts", ">h"), + ("n_ref_cal_pts", ">h"), + ("tst_cal_pts", ">4f"), + ("ref_cal_pts", ">4f"), + ("tst_cal_pix_opd", ">f"), + ("ref_cal_pix_opd", ">f"), + ("sys_serial2", ">i"), + ("flash_phase_dc_mask", ">f"), + ("flash_phase_alias_mask", ">f"), + ("flash_phase_filter", ">f"), ("scan_direction", "B"), (None, "1b"), ("pre_fda_filter", ">h"), (None, "4b"), - ("ftpsi_res_factor", "i"), (None, "8b"), ], name="header1", @@ -274,21 +265,21 @@ def can_read(cls, buffer: bytes) -> MagicMatch: ("asphere_att8", ">f"), ("asphere_aperture_pct", ">f"), ("asphere_optimized_r0", ">f"), - ("iff_state", "i"), ("iff_idr_filename", "42s"), ("iff_ise_filename", "42s"), (None, "2b"), ("asphere_eqn_r0", ">f"), ("asphere_eqn_k", ">f"), ("asphere_eqn_coef", ">21f"), - ("awm_enable", "i"), + ("awm_vacuum_wavelength_nm", ">f"), + ("awm_air_wavelength_nm", ">f"), + ("awm_air_temperature_degc", ">f"), + ("awm_air_pressure_mmhg", ">f"), + ("awm_air_rel_humidity_pct", ">f"), + ("awm_air_quality", ">f"), + ("awm_input_power_mw", ">f"), ("asphere_optimizations", ">i"), ("asphere_optimization_mode", ">i"), ("asphere_optimized_k", ">f"), @@ -301,7 +292,7 @@ def can_read(cls, buffer: bytes) -> MagicMatch: (None, "2b"), ("n_fiducials_d", ">h"), ("fiducials_d", ">14f"), - ("gpi_enc_zoom_mag", "f"), ("asphere_max_distortion", ">f"), ("asphere_distortion_uncert", ">f"), ("field_stop_name", "12s"), diff --git a/SurfaceTopography/IO/Mitutoyo.py b/SurfaceTopography/IO/Mitutoyo.py index 3b305bd1..722c1fa7 100644 --- a/SurfaceTopography/IO/Mitutoyo.py +++ b/SurfaceTopography/IO/Mitutoyo.py @@ -26,8 +26,8 @@ import logging import re -from datetime import datetime +import dateutil.parser import numpy as np import openpyxl @@ -98,7 +98,7 @@ def __init__(self, fobj): # check if positions are distributed uniformly _diff = np.diff(_x) - if np.all(np.isclose(_diff, _diff[0])): + if len(_diff) == 0 or np.all(np.isclose(_diff, _diff[0])): self._uniform = True else: self._uniform = False @@ -135,9 +135,6 @@ def __init__(self, fobj): cell.value, ) - # try to infer heights unit from roughness metrics - _h_unit = _roughness_metrics_list[0]["unit"] - # get creation date _date_string = _metadata["E2"].value @@ -149,6 +146,14 @@ def __init__(self, fobj): _cut_off_dict = cut_off_regex.match(_cut_off_string).groupdict() _x_unit = _cut_off_dict["unit"] + # try to infer heights unit from roughness metrics; fall back to + # the lateral unit if no roughness metrics are present in the + # spreadsheet + if len(_roughness_metrics_list) > 0: + _h_unit = _roughness_metrics_list[0]["unit"] + else: + _h_unit = _x_unit + # convert x unit to h unit self._x = _x * get_unit_conversion_factor(_x_unit, _h_unit) self._unit = _h_unit @@ -160,7 +165,10 @@ def __init__(self, fobj): self._physical_sizes = np.max(self._x) self._info = { - "acquisition_time": datetime.strptime(_date_string, "%d-%b-%Y"), + # Note: `dateutil` parses month names independent of the + # current locale (datetime.strptime with '%b' fails on + # non-English locales for the same file) + "acquisition_time": dateutil.parser.parse(_date_string, dayfirst=True), "instrument": {"vendor": "Mitutoyo"}, "raw_metadata": { "roughness_metrics": _roughness_metrics_list, @@ -182,6 +190,7 @@ def __init__(self, fobj): info=self._info, ) ] + wb.close() # Return list of channels (here only a single one) @property diff --git a/SurfaceTopography/IO/NC.py b/SurfaceTopography/IO/NC.py index 72106603..6794199c 100644 --- a/SurfaceTopography/IO/NC.py +++ b/SurfaceTopography/IO/NC.py @@ -285,14 +285,16 @@ def channels(self): ny = self._y_dim nb_grid_pts = (nx, ny) else: - # This is a nonuniform line scan + # This is a nonuniform line scan; the number of points is stored + # in the 'n' dimension (there is no 'x' dimension here, which is + # why we ended up in this branch) uniform = False try: # netCDF4 - nx = self._x_dim.size + nx = self._n_dim.size except AttributeError: # scipy.io.netcdf_file - nx = self._x_dim + nx = self._n_dim nb_grid_pts = (nx,) return [ ChannelInfo( @@ -449,7 +451,10 @@ def write_nc_uniform(topography, fobj, format="NETCDF3_64BIT_OFFSET"): Dataset = _SpecialNetCDFFile kwargs = dict(version=format_to_scipy_version[format], maskandscale=True) var_kwargs = {} - if not topography.is_domain_decomposed and topography.communicator.rank > 1: + if not topography.is_domain_decomposed and topography.communicator.rank > 0: + # Only the root rank writes the file if the topography is not + # decomposed; with `rank > 1` both ranks 0 and 1 would write the + # same file concurrently and corrupt it return with Dataset(fobj, "w", **kwargs) as nc: # Serialize info dictionary as JSON and write to NetCDF file diff --git a/SurfaceTopography/IO/NMM.py b/SurfaceTopography/IO/NMM.py index ce5121ee..8bfd9daf 100644 --- a/SurfaceTopography/IO/NMM.py +++ b/SurfaceTopography/IO/NMM.py @@ -255,8 +255,8 @@ def channels(self): return [ ChannelInfo( self, - 0, # Channel index - name=f"Scan {i+1}", # There is only a single channel + i, # Channel index; each scan is exposed as its own channel + name=f"Scan {i+1}", dim=2, nb_grid_pts=self._nb_grid_pts, physical_sizes=self._physical_sizes, diff --git a/SurfaceTopography/IO/OIR.py b/SurfaceTopography/IO/OIR.py index c251dad8..718cf7b1 100644 --- a/SurfaceTopography/IO/OIR.py +++ b/SurfaceTopography/IO/OIR.py @@ -520,7 +520,11 @@ def _validate_metadata(self): tags={ # The suffix _0 is probably the frame number, but I have # never seen files with multiple frames. - "reader": lambda stream_obj: np.frombuffer( + # Note: The loop variables are bound as default + # arguments; a plain closure would capture them + # by reference and every channel would end up + # reading the *last* channel's data. + "reader": lambda stream_obj, prefix=prefix, uuid=uuid, nx=nx, ny=ny, dtype=dtype: np.frombuffer( # noqa: E501 data[f"{prefix}_{uuid}_0"](stream_obj), dtype ) .reshape((ny, nx)) @@ -566,6 +570,11 @@ def channels(self): for fn, r in self._readers: for c in r.channels: c.reader = self + # Renumber the channel to its position in the concatenated + # list; it otherwise retains the index within its own file + # and `topography()` would return the data of a different + # channel. + c.index = len(channels) c.tags["fn"] = fn channels += [c] return channels diff --git a/SurfaceTopography/IO/OPD.py b/SurfaceTopography/IO/OPD.py index c316cad2..6a40548b 100644 --- a/SurfaceTopography/IO/OPD.py +++ b/SurfaceTopography/IO/OPD.py @@ -26,7 +26,7 @@ import numpy as np -from ..Exceptions import MetadataAlreadyFixedByFile +from ..Exceptions import CorruptFile, MetadataAlreadyFixedByFile from ..UniformLineScanAndTopography import Topography from .common import OpenFromAny from .Reader import ChannelInfo, ReaderBase @@ -102,9 +102,16 @@ def read_block(f): pixel_size = 1.0 aspect = 1.0 mult = 1.0 + wavelength = None for n, t, L, a in blocks: if L <= 0: continue + # Remember the block start; each block occupies exactly the + # length L declared in the directory, independent of how + # many bytes we actually interpret. (Consuming fixed byte + # counts would silently shift all subsequent block offsets + # if L differs.) + block_start = f.tell() if n == "RAW DATA" or n == "RAW_DATA" or n == "OPD" or n == "Raw": nx, ny, elsize = unpack("= len(self._channels): raise RuntimeError( - f"There is only a single channel. Channel index must be {self._default_channel_index}." + f"Channel index is {channel_index} but must be between 0 and " + f"{len(self._channels) - 1}." ) if physical_sizes is not None: diff --git a/SurfaceTopography/IO/PS.py b/SurfaceTopography/IO/PS.py index 0ec4f6e1..d87e300b 100644 --- a/SurfaceTopography/IO/PS.py +++ b/SurfaceTopography/IO/PS.py @@ -29,7 +29,7 @@ import io import numpy as np -from tiffile import TiffFile, TiffFileError +from tifffile import TiffFile, TiffFileError from ..Exceptions import CorruptFile, FileFormatMismatch, MetadataAlreadyFixedByFile from ..Support.UnitConversion import get_unit_conversion_factor @@ -207,8 +207,12 @@ def topography(self, channel_index=None, physical_sizes=None, t = self._header['data_type'] dtype = np.int16 if t == self._DATA_TYPE_INT16 else \ np.int32 if t == self._DATA_TYPE_INT32 else np.single + # The data is stored line by line, i.e. the buffer has + # C-order shape (ny, nx). Transposing yields the (nx, ny) + # array with the x index first that `Topography` expects. + nx, ny = self._nb_grid_pts height_data = \ - np.frombuffer(raw_data, dtype=dtype, count=np.prod(self._nb_grid_pts)).reshape(self._nb_grid_pts).T + np.frombuffer(raw_data, dtype=dtype, count=nx * ny).reshape((ny, nx)).T _info = self._info.copy() _info.update(info) diff --git a/SurfaceTopography/IO/Reader.py b/SurfaceTopography/IO/Reader.py index a9b575f5..1cd843f2 100644 --- a/SurfaceTopography/IO/Reader.py +++ b/SurfaceTopography/IO/Reader.py @@ -238,6 +238,15 @@ def index(self): """Unique integer channel index.""" return self._index + @index.setter + def index(self, value): + """ + Set the channel index. This is used by container readers that + aggregate the channels of multiple file readers into a single list + and need to renumber them. + """ + self._index = value + @property def name(self): """ @@ -1295,6 +1304,8 @@ def metadata(self): def topography( self, channel_index=None, + channel_id=None, + height_channel_index=None, physical_sizes=None, height_scale_factor=None, unit=None, @@ -1306,17 +1317,11 @@ def topography( if subdomain_locations is not None or nb_subdomain_grid_pts is not None: raise RuntimeError("This reader does not support MPI parallelization.") - if channel_index is None: - channel_index = self._default_channel_index - - channels = self.channels - if channel_index < 0 or channel_index >= len(channels): - raise RuntimeError( - f"Channel index is {channel_index} but must be between 0 and {len(channels) - 1}." - ) - - # Get channel information - channel = channels[channel_index] + # Get channel information (this also supports the `channel_id` and + # `height_channel_index` selection methods of the base class) + channel, channel_index = self._resolve_channel( + channel_index, channel_id, height_channel_index + ) if physical_sizes is None: physical_sizes = channel.physical_sizes @@ -1346,4 +1351,8 @@ def topography( periodic=False if periodic is None else periodic, info=_info, ) + if height_scale_factor is None: + # A declarative reader is not required to provide height-scale + # metadata; without it, the heights are returned unscaled + return topo return topo.scale(height_scale_factor) diff --git a/SurfaceTopography/IO/SDF.py b/SurfaceTopography/IO/SDF.py index c5390efe..8975feed 100644 --- a/SurfaceTopography/IO/SDF.py +++ b/SurfaceTopography/IO/SDF.py @@ -102,14 +102,21 @@ def _parse_ascii_header(content): return header -def _read_ascii_data(data_section, num_points, num_profiles, z_scale): +def _read_ascii_data(data_section, num_points, num_profiles, z_scale, data_type): """Parse the ASCII data section.""" # Replace BAD markers with NaN data_section = data_section.replace("BAD", "NAN") # Parse data - data = np.fromstring(data_section, sep=" ", dtype=np.float64) + data = np.array(data_section.split(), dtype=np.float64) data = data.reshape(num_profiles, num_points) + + # ASCII files with an integer DataType may also mark invalid points + # with the same numeric marker as the binary variant + invalid_value = INVALID_VALUE_MAP[data_type] + if not np.isnan(invalid_value): + data[data == invalid_value] = np.nan + data *= z_scale return data @@ -309,7 +316,7 @@ def topography( data *= z_scale else: # Parse ASCII data - data = _read_ascii_data(self._data_section, nx, ny, z_scale) + data = _read_ascii_data(self._data_section, nx, ny, z_scale, data_type) # Transpose to get (nx, ny) ordering data = data.T diff --git a/SurfaceTopography/IO/Text.py b/SurfaceTopography/IO/Text.py index f7389975..89e97dd5 100644 --- a/SurfaceTopography/IO/Text.py +++ b/SurfaceTopography/IO/Text.py @@ -267,7 +267,9 @@ def parse_metadata(self, line): def __init__(self, file_path): # Open file and parse self._channel_names = [] - self._metadata = defaultdict(None) + # Note: this is a plain dictionary; `defaultdict(None)` behaves + # identically (None is not callable and cannot construct defaults) + self._metadata = {} self._data = defaultdict(list) self._dim = 2 with OpenFromAny(file_path, "r") as fobj: @@ -464,7 +466,7 @@ def topography( if channel_index is None: channel_index = self._default_channel_index - if channel_index < 0 or channel_index > len(self._channel_names): + if channel_index < 0 or channel_index >= len(self._channel_names): raise RuntimeError( f"There are only {len(self._channel_names)} channels, but channel " f"index is {channel_index}." diff --git a/SurfaceTopography/IO/VK.py b/SurfaceTopography/IO/VK.py index 792e4b59..77537067 100644 --- a/SurfaceTopography/IO/VK.py +++ b/SurfaceTopography/IO/VK.py @@ -255,9 +255,13 @@ def read_vk34_header(self, f): ) ) + # Note: the physical size is the number of pixels times the pixel + # size (pixel convention, like everywhere else in this library and + # in other implementations of this file format), not the distance + # between the first and last pixel centers self._physical_sizes = ( - float(self._data["width"] - 1) * self._header["x_length_per_pixel"], - float(self._data["height"] - 1) * self._header["y_length_per_pixel"], + float(self._data["width"]) * self._header["x_length_per_pixel"], + float(self._data["height"]) * self._header["y_length_per_pixel"], ) self._unit = "pm" diff --git a/SurfaceTopography/IO/WSXM.py b/SurfaceTopography/IO/WSXM.py index 8ea77582..247b9b88 100644 --- a/SurfaceTopography/IO/WSXM.py +++ b/SurfaceTopography/IO/WSXM.py @@ -150,7 +150,12 @@ def __init__(self, file_path): z_amplitude = float(z_amplitude) min_value = float(metadata["Miscellaneous"]["Minimum"]) max_value = float(metadata["Miscellaneous"]["Maximum"]) - height_scale_factor = z_amplitude / (max_value - min_value) + if max_value > min_value: + height_scale_factor = z_amplitude / (max_value - min_value) + else: + # A perfectly flat scan has zero data range; any scale + # factor reproduces it + height_scale_factor = 1 self._channels = [ ChannelInfo( @@ -234,7 +239,6 @@ def topography( periodic=False if periodic is None else periodic, info=_info, ) - print(height_data.min(), height_data.max()) return topo.scale(channel.height_scale_factor) channels.__doc__ = ReaderBase.channels.__doc__ diff --git a/SurfaceTopography/IO/X3P.py b/SurfaceTopography/IO/X3P.py index e43e9acd..a956152b 100644 --- a/SurfaceTopography/IO/X3P.py +++ b/SurfaceTopography/IO/X3P.py @@ -147,6 +147,11 @@ def __init__(self, file_path): else: self._height_scale_factor = None + # The z-axis may carry an offset: h = raw * Increment + Offset. + # This is typically present for integer data types. + offset = cz.find("Offset") + self._height_offset = float(offset.text) if offset is not None else 0.0 + # Parse record3 matrix_dimension = record3.find("MatrixDimension") nb_grid_pts_x = int(matrix_dimension.find("SizeX").text) @@ -167,6 +172,14 @@ def __init__(self, file_path): data_link = record3.find("DataLink") self._name_of_binary_file = data_link.find("PointDataLink").text + # Optional file marking invalid (undefined) data points, + # one bit per point (ISO 5436-2). This is the only way + # invalid points can be encoded for integer data types. + valid_points_link = data_link.find("ValidPointsLink") + self._name_of_valid_points_file = ( + valid_points_link.text if valid_points_link is not None else None + ) + # Check if binary file exists and has a reasonable size binary_info = x3p.getinfo(self._name_of_binary_file) expected_file_size = ( @@ -322,6 +335,27 @@ def topography( .T ) + if self._name_of_valid_points_file is not None: + # One bit per point, least-significant bit first, in the + # same point order as the height data + validdata = x3p.open(self._name_of_valid_points_file).read() + valid = ( + np.unpackbits( + np.frombuffer(validdata, dtype=np.uint8), bitorder="little" + )[: nx * ny] + .reshape(ny, nx) + .T.astype(bool) + ) + if not valid.all(): + height_data = np.ma.masked_array(height_data, mask=~valid) + + if self._height_offset != 0: + # h = raw * Increment + Offset; we fold the offset into the raw + # data so that the height scale factor can still be applied + # through the pipeline + fac = self._height_scale_factor if self._height_scale_factor is not None else 1 + height_data = height_data + self._height_offset / fac + topo = Topography( height_data, self._physical_sizes, @@ -408,12 +442,13 @@ def write_x3p( dx = sx * scale_to_meters / nx dy = sy * scale_to_meters / ny - # Get height data and convert to meters - heights = self.heights() - if np.ma.isMaskedArray(heights): - # X3P doesn't directly support masked data - use NaN for undefined - heights = np.ma.filled(heights, np.nan) - heights = heights * scale_to_meters + # Get height data and convert to meters. Undefined (masked) data points + # are encoded as NaN for floating-point data types; for integer data + # types they are stored as the minimum raw value and marked as invalid + # in the validity file (written below). + heights = np.ma.filled(self.heights(), np.nan) * scale_to_meters + invalid_mask = np.isnan(heights) + has_invalid = invalid_mask.any() # For integer types, we need to scale the data if dtype in ("I", "L"): @@ -422,6 +457,11 @@ def write_x3p( height_range = height_max - height_min if height_range == 0: height_range = 1.0 + if has_invalid: + # NaN cannot be cast to an integer; the actual value stored for + # invalid points is arbitrary since they are flagged in the + # validity file + heights = np.where(invalid_mask, height_min, heights) if dtype == "I": scale_factor = height_range / 65535 heights = ((heights - height_min) / scale_factor).astype(np_dtype) @@ -522,6 +562,17 @@ def write_x3p( ElementTree.SubElement(data_link, "PointDataLink").text = "bindata/data.bin" ElementTree.SubElement(data_link, "MD5ChecksumPointData").text = md5_hash + # Validity data: one bit per point, least-significant bit first, in the + # same point order as the height data + valid_data = None + if has_invalid: + valid_data = np.packbits( + ~invalid_mask.T.reshape(-1), bitorder="little" + ).tobytes() + ElementTree.SubElement(data_link, "ValidPointsLink").text = "bindata/valid.bin" + ElementTree.SubElement(data_link, "MD5ChecksumValidPoints").text = \ + hashlib.md5(valid_data).hexdigest().upper() + # Record4: Checksum file reference (optional, but included for completeness) record4 = ElementTree.SubElement(root, "Record4") ElementTree.SubElement(record4, "ChecksumFile").text = "md5checksum.hex" @@ -566,18 +617,13 @@ def serialize_element(elem, indent=0): # Create checksum file content checksum_content = f"{md5_hash} *bindata/data.bin\n" - # Write ZIP file - if isinstance(fobj, str): - with ZipFile(fobj, "w") as x3p: - x3p.writestr("main.xml", xml_string.encode("utf-8")) - x3p.writestr("bindata/data.bin", binary_data) - x3p.writestr("md5checksum.hex", checksum_content.encode("utf-8")) - else: - # File-like object - with ZipFile(fobj, "w") as x3p: - x3p.writestr("main.xml", xml_string.encode("utf-8")) - x3p.writestr("bindata/data.bin", binary_data) - x3p.writestr("md5checksum.hex", checksum_content.encode("utf-8")) + # Write ZIP file (`ZipFile` accepts both file names and file-like objects) + with ZipFile(fobj, "w") as x3p: + x3p.writestr("main.xml", xml_string.encode("utf-8")) + x3p.writestr("bindata/data.bin", binary_data) + if valid_data is not None: + x3p.writestr("bindata/valid.bin", valid_data) + x3p.writestr("md5checksum.hex", checksum_content.encode("utf-8")) UniformTopographyInterface.register_function("to_x3p", write_x3p) diff --git a/SurfaceTopography/IO/XYZ.py b/SurfaceTopography/IO/XYZ.py index 18797846..995bbb5e 100644 --- a/SurfaceTopography/IO/XYZ.py +++ b/SurfaceTopography/IO/XYZ.py @@ -24,6 +24,7 @@ import logging import re +from array import array from collections import defaultdict import numpy as np @@ -44,7 +45,7 @@ from .common import OpenFromAny from .Reader import ChannelInfo, ReaderBase -_log = logging.Logger(__name__) +_log = logging.getLogger(__name__) def read_text_header_hfm(fobj, unit, height_scale_factor): @@ -223,7 +224,9 @@ def read_csv(fobj, sep=None, usecols=None, skiprows=0): for i in range(skiprows): fobj.readline() line = fobj.readline() - data = defaultdict(list) + # Accumulate raw doubles (8 bytes per value); lists of Python floats + # would need roughly four times the memory during parsing + data = defaultdict(lambda: array('d')) min_cols = None nb_cols = None if usecols is not None: @@ -251,11 +254,11 @@ def read_csv(fobj, sep=None, usecols=None, skiprows=0): if usecols is None: # If no columns are given, we return all columns for key, value in enumerate(line): - data[key] += [float(value)] + data[key].append(float(value)) else: # If columns are given by the user, only collect the data from those colums for i, key in enumerate(usecols): - data[key] += [float(line[i])] + data[key].append(float(line[i])) line = fobj.readline() nb_lines += 1 if usecols is None: diff --git a/SurfaceTopography/IO/ZON.py b/SurfaceTopography/IO/ZON.py index 8fadba68..32fbd203 100644 --- a/SurfaceTopography/IO/ZON.py +++ b/SurfaceTopography/IO/ZON.py @@ -37,7 +37,7 @@ import zstandard from numpy.lib.stride_tricks import as_strided -from ..Exceptions import FileFormatMismatch, MetadataAlreadyFixedByFile +from ..Exceptions import FileFormatMismatch, MetadataAlreadyFixedByFile, UnsupportedFormatFeature from ..UniformLineScanAndTopography import Topography from .binary import decode from .common import OpenFromAny @@ -80,10 +80,13 @@ def _read_array(f, dtype=np.dtype(" qL), - np.where(q < qr, - self.cr, - self.cr * (q / qr) ** (-2 - 2 * self.hurst_exponent)), - 0) + # Note: only evaluate the power law where it applies; evaluating it + # at q = 0 (e.g. through np.where) would produce division warnings + # and infinities + q = np.asarray(q, dtype=float) + C = np.zeros_like(q) + inside = np.logical_and(q < qs, q > qL) + plateau = np.logical_and(inside, q < qr) + decay = np.logical_and(inside, q >= qr) + C[plateau] = self.cr + C[decay] = self.cr * (q[decay] / qr) ** (-2 - 2 * self.hurst_exponent) + return C def power_spectrum_profile(self, q): # TODO: I think this is only the simplified formula. @@ -215,7 +220,11 @@ def variance_derivative(self, order, shortcut_wavelength=None, longcut_wavelengt if shortcut_wavevector > rolloff_wavevector and longcut_wavevector < shortcut_wavevector: # self-affine region if self.hurst_exponent == order: - self_affine = c0 / (2 * np.pi) * np.log(shortcut_wavevector / rolloff_wavevector) + # The self-affine region spans from the larger of the + # rolloff and longcut wavevectors to the shortcut + # wavevector, like in the general-order branch below + self_affine = c0 / (2 * np.pi) * np.log( + shortcut_wavevector / max(rolloff_wavevector, longcut_wavevector)) else: self_affine = (c0 / (2 * np.pi) * (shortcut_wavevector ** (2 * (order - self.hurst_exponent)) @@ -253,14 +262,21 @@ def generate_roughness(self, dx = pixel_size sx = nx * dx - c0 = self.cr * (2 * np.pi / self.rolloff_wavelength) ** (2 + 2 * self.hurst_exponent) + # Build a model with the (possibly overridden) cutoffs and hand its + # PSD to the Fourier synthesis. This honors shortcut, rolloff *and* + # longcut; passing `c0`/`long_cutoff` instead would only describe + # the rolloff knee and silently ignore `longcut_wavelength`. + model = SelfAffine(cr=self.cr, + rolloff_wavelength=self.rolloff_wavelength, + hurst_exponent=self.hurst_exponent, + longcut_wavelength=longcut_wavelength, + shortcut_wavelength=shortcut_wavelength, + unit=self.unit) np.random.seed(seed) roughness = fourier_synthesis((n_pixels, n_pixels), (sx, sx), - hurst=self.hurst_exponent, - c0=c0, - short_cutoff=shortcut_wavelength, - long_cutoff=self.rolloff_wavelength, + psd=model.power_spectrum_isotropic, + short_cutoff=shortcut_wavelength if shortcut_wavelength > 0 else None, unit=self.unit, **kwargs ) diff --git a/SurfaceTopography/Nonuniform/Autocorrelation.py b/SurfaceTopography/Nonuniform/Autocorrelation.py index b720f90d..3b324787 100644 --- a/SurfaceTopography/Nonuniform/Autocorrelation.py +++ b/SurfaceTopography/Nonuniform/Autocorrelation.py @@ -71,37 +71,10 @@ def height_height_autocorrelation(self, distances=None): distances = np.linspace(0, size, res) else: distances = np.asarray(distances, dtype=float) - A = np.zeros_like(distances) x, h = self.positions_and_heights() - s = self.derivative(1) - # FIXME!!! This is slow - for i in range(len(x) - 1): - for j in range(len(x) - 1): - # Determine lower and upper distance between segment i, i+1 and - # segment j, j+1 - x1 = x[i] - x2 = x[j] - h1 = h[i] - h2 = h[j] - s1 = s[i] - s2 = s[j] - b1 = np.maximum(x1, x2 - distances) - b2 = np.minimum(x[i + 1], x[j + 1] - distances) - b = (b1 + b2) / 2 - db = (b2 - b1) / 2 - m = db > 0 - if m.sum() > 0: - b = b[m] - db = db[m] - # f1[x_] := (h1 + s1*(x - x1)) - # f2[x_] := (h2 + s2*(x - x2)) - # FullSimplify[Integrate[f1[x]*f2[x + d], - # {x, b - db, b + db}]] - # = 2 * f1[b] * f2[b + d] * db + 2 * s1 * s2 * db ** 3 / 3 - A[m] += 2 * (h1 + s1 * (b - x1)) * ( - h2 + s2 * (b + distances[m] - x2)) * db + 2 * ( - s1 * s2 * db ** 3) / 3 + A = _SurfaceTopography.nonuniform_height_height_autocorrelation( + np.asarray(x, dtype=float), np.asarray(h, dtype=float), distances) return distances, A diff --git a/SurfaceTopography/Nonuniform/Converters.py b/SurfaceTopography/Nonuniform/Converters.py index 4c57aabc..9b495095 100644 --- a/SurfaceTopography/Nonuniform/Converters.py +++ b/SurfaceTopography/Nonuniform/Converters.py @@ -56,7 +56,9 @@ def __init__(self, topography, nb_points=None, padding=0, nb_interpolate=None, p automatically computed from the mean grid spacing and `nb_interpolate` if set to None. (Default: None) padding : int, optional - Number of padding grid points, zeros appended to the data. + Number of padding grid points appended to the data; they hold + the value of the final data point (np.interp clamps to the + boundary values). (Default: 0) nb_interpolate : int, optional Number of grid points to between closest points on surface. @@ -73,7 +75,7 @@ def __init__(self, topography, nb_points=None, padding=0, nb_interpolate=None, p # This is populated with functions from the nonuniform topography, but # this is a uniform topography - self._functions = UniformLineScan._functions + self._functions = UniformLineScan._all_functions() def _update_nb_points_and_pixel_size(self): """Automatically compute `nb_points` and `pixel_size` if it is None""" @@ -161,8 +163,13 @@ def positions(self): def heights(self): """ Computes the rescaled profile. """ - x = self.positions() - return np.interp(x, *self.parent_topography.positions_and_heights()) + parent_x, parent_h = self.parent_topography.positions_and_heights() + # The uniform grid starts at zero, but the parent topography's + # x-coordinates may start anywhere; shift the interpolation points + # by the parent's origin. (np.interp clamps outside the data range, + # so without the shift most of the grid would be filled with the + # first or last height value.) + return np.interp(self.positions() + parent_x[0], parent_x, parent_h) # Register pipeline functions from this module diff --git a/SurfaceTopography/Nonuniform/Detrending.py b/SurfaceTopography/Nonuniform/Detrending.py index 61582d4f..3faac687 100644 --- a/SurfaceTopography/Nonuniform/Detrending.py +++ b/SurfaceTopography/Nonuniform/Detrending.py @@ -28,6 +28,7 @@ """ import numpy as np +from numpy.polynomial import polynomial as npoly from SurfaceTopography.HeightContainer import NonuniformLineScanInterface from SurfaceTopography.NonuniformLineScan import DecoratedNonuniformTopography @@ -75,17 +76,49 @@ def polyfit(self, deg): Array with coefficients :math:`a_k`. """ # noqa: E501 x, h = self.positions_and_heights() - dx = np.diff(x) + # Solve the normal equations in coordinates centered on the scan: + # raw monomials x^k make the system ill-conditioned when the scan + # is short compared to its distance from x = 0. The coefficients + # are transformed back to monomials of the absolute position below. + x0 = (x[0] + x[-1]) / 2 + xs = x - x0 + dx = np.diff(xs) k = np.arange(deg + 1).reshape(-1, 1) - b = np.sum(((2 * h[:-1] + h[1:]) * x[:-1] ** k + - (2 * h[1:] + h[:-1]) * x[1:] ** k) * dx, + b = np.sum(((2 * h[:-1] + h[1:]) * xs[:-1] ** k + + (2 * h[1:] + h[:-1]) * xs[1:] ** k) * dx, axis=1) L = k.reshape(1, -1, 1) k = k.reshape(-1, 1, 1) - A = np.sum((2 * x[:-1] ** (k + L) + 2 * x[1:] ** (k + L) + - x[:-1] ** k * x[1:] ** L + x[1:] ** k * x[:-1] ** L) * dx, + A = np.sum((2 * xs[:-1] ** (k + L) + 2 * xs[1:] ** (k + L) + + xs[:-1] ** k * xs[1:] ** L + xs[1:] ** k * xs[:-1] ** L) * dx, axis=2) - return np.linalg.solve(A, b) + a = np.linalg.solve(A, b) + # Expand p(x) = sum_k a_k (x - x0)^k into monomial coefficients of x. + # (numpy's polypow trims trailing zeros, hence the copy into a + # fixed-size output array.) + p = np.zeros(deg + 1) + for n, a_n in enumerate(a): + term = a_n * npoly.polypow([-x0, 1.0], n) + p[:len(term)] += term + return p + + +def _slope_detrend_coeffs(topography): + """ + Compute the coefficients [a0, a1] that minimize the rms slope of the + detrended profile and center it around zero. + + The constant slope that minimizes the rms slope is the length-weighted + mean of the derivative, i.e. (h(x_max) - h(x_min)) / (x_max - x_min). + (A plain mean over the per-segment slopes would weight each segment + equally, irrespective of its length.) The offset a0 removes the mean of + the tilt-corrected profile; note that the trend a0 + a1 x is evaluated + at the absolute positions x, which do not necessarily start at zero. + """ + x, h = topography.positions_and_heights() + a1 = (h[-1] - h[0]) / (x[-1] - x[0]) + a0 = topography.mean() - a1 * (x[0] + x[-1]) / 2 + return [a0, a1] class DetrendedNonuniformTopography(DecoratedNonuniformTopography): @@ -104,7 +137,7 @@ class DetrendedNonuniformTopography(DecoratedNonuniformTopography): # same as 'rms-tilt', deprecate 'height' in the future 'height': lambda self: self.parent_topography.polyfit(1), 'mad-tilt': lambda self: self.parent_topography.mad_polyfit(1), - 'slope': lambda self: [self.parent_topography.mean(), self.parent_topography.derivative(1).mean()], + 'slope': lambda self: _slope_detrend_coeffs(self.parent_topography), 'rms-curvature': lambda self: self.parent_topography.polyfit(2), # same as 'rms-curvature', deprecate 'curvature' in the future 'curvature': lambda self: self.parent_topography.polyfit(2), diff --git a/SurfaceTopography/Nonuniform/Interpolation.py b/SurfaceTopography/Nonuniform/Interpolation.py index a593fa4c..5073e518 100644 --- a/SurfaceTopography/Nonuniform/Interpolation.py +++ b/SurfaceTopography/Nonuniform/Interpolation.py @@ -42,7 +42,7 @@ def interpolate_linear(self): def interpolate_cubic(self): r""" - Returns a linear interpolation function based on the topography's heights. + Returns a cubic interpolation function based on the topography's heights. """ if self.is_reentrant: raise ReentrantDataError('This topography is reentrant (i.e. it contains overhangs). Interpolation is not ' diff --git a/SurfaceTopography/Nonuniform/PowerSpectrum.py b/SurfaceTopography/Nonuniform/PowerSpectrum.py index bc457bca..edfb489c 100644 --- a/SurfaceTopography/Nonuniform/PowerSpectrum.py +++ b/SurfaceTopography/Nonuniform/PowerSpectrum.py @@ -40,8 +40,13 @@ def sinc(x): def dsinc(x): """Derivative of the sinc function, d/dx [sin(x)/x].""" - tol = 1e-6 - x = np.asarray(x) + # Note: cos(x) - sinc(x) is formed by cancellation of two O(1) terms + # and loses accuracy for small x; the Taylor branch is accurate to + # machine precision up to about x = 1e-2, so the crossover sits where + # both branches are accurate. (A threshold of 1e-6 would leave a band + # with relative errors up to 1e-3 just above it.) + tol = 1e-2 + x = np.asarray(x, dtype=float) small_values = np.abs(x) < tol if small_values.sum() > 0: ret = np.zeros_like(x) @@ -106,7 +111,9 @@ def apply_window(x, y, window=None): """ if window == 'hann': length = x.max() - x.min() - return (2 / 3) ** (1 / 2) * (1 - np.cos(2 * np.pi * x / length)) * y + # The window argument must be relative to the start of the scan; + # the x-coordinates do not necessarily start at zero. + return (2 / 3) ** (1 / 2) * (1 - np.cos(2 * np.pi * (x - x.min()) / length)) * y elif window is None or window == 'None': return y else: @@ -206,8 +213,6 @@ def power_spectrum(self, reliable=True, algorithm='fft', wavevectors=None, nb_in y = apply_window(x, y, window=window) L = x[-1] - x[0] - if wavevectors is None: - wavevectors = 2 * np.pi * np.arange(int(L / np.diff(x).min())) / L y_q = np.zeros_like(wavevectors, dtype=complex) for x1, x2, y1, y2 in zip(x[:-1], x[1:], y[:-1], y[1:]): dx = x2 - x1 diff --git a/SurfaceTopography/Nonuniform/ScalarParameters.py b/SurfaceTopography/Nonuniform/ScalarParameters.py index 0bea5bf3..f0bd3c4c 100644 --- a/SurfaceTopography/Nonuniform/ScalarParameters.py +++ b/SurfaceTopography/Nonuniform/ScalarParameters.py @@ -63,15 +63,22 @@ def moment(topography, alpha): Returns ------- - moment : float or array - Root-mean square height. + moment : float + Moment of order `alpha` of the heights. """ # noqa: E501 x, h = topography.positions_and_heights() dx = np.diff(x) if len(x) <= 1: return 0.0 L = x[-1] - x[0] - return 1 / (alpha + 1) * np.sum(dx * (h[1:] ** (alpha + 1) - h[:-1] ** (alpha + 1)) / (h[1:] - h[:-1])) / L + # The quotient (h2^(alpha+1) - h1^(alpha+1)) / (h2 - h1) given in the + # docstring equals the complete homogeneous symmetric polynomial of + # degree alpha in h1 and h2. The polynomial form is exact for equal + # heights (where the quotient is 0/0, e.g. for quantized instrument + # data) and avoids cancellation for nearly equal heights. + h1, h2 = h[:-1], h[1:] + s = sum(h1 ** k * h2 ** (alpha - k) for k in range(alpha + 1)) + return np.sum(dx * s) / ((alpha + 1) * L) def rms_height(self): @@ -144,7 +151,7 @@ def rms_slope(self): def rms_curvature(self): r""" - Computes root-mean square slope fluctuation of the line scan: + Computes root-mean square curvature fluctuation of the line scan: Parameters ---------- @@ -153,8 +160,8 @@ def rms_curvature(self): Returns ------- - rms_slope : float - Root-mean square slope. + rms_curvature : float + Root-mean square curvature. """ x = self.positions() d2 = self.derivative(n=2) @@ -168,9 +175,14 @@ def rms_curvature(self): NonuniformLineScanInterface.register_function( 'mean', lambda self: _SurfaceTopography.nonuniform_mean(*self.positions_and_heights())) NonuniformLineScanInterface.register_function('moment', moment) -NonuniformLineScanInterface.register_function('rms_height_from_profile', rms_height, deprecated=True) +# Note: The `rms_*_from_profile` names are not marked as deprecated here +# (although they used to carry an inactive `deprecated` flag): they are the +# primary names of the corresponding uniform analysis functions, so +# deprecating only the nonuniform variants would make generic code warn +# depending on the data type it happens to operate on. +NonuniformLineScanInterface.register_function('rms_height_from_profile', rms_height) NonuniformLineScanInterface.register_function('Rq', rms_height) -NonuniformLineScanInterface.register_function('rms_slope_from_profile', rms_slope, deprecated=True) +NonuniformLineScanInterface.register_function('rms_slope_from_profile', rms_slope) NonuniformLineScanInterface.register_function('Rdq', rms_slope) -NonuniformLineScanInterface.register_function('rms_curvature_from_profile', rms_curvature, deprecated=True) +NonuniformLineScanInterface.register_function('rms_curvature_from_profile', rms_curvature) NonuniformLineScanInterface.register_function('Rddq', rms_curvature) diff --git a/SurfaceTopography/Nonuniform/VariableBandwidth.py b/SurfaceTopography/Nonuniform/VariableBandwidth.py index 8e9fc39c..54a79ce8 100644 --- a/SurfaceTopography/Nonuniform/VariableBandwidth.py +++ b/SurfaceTopography/Nonuniform/VariableBandwidth.py @@ -53,8 +53,10 @@ def checkerboard_detrend_profile(self, subdivisions, tol=1e-6): subdivisions : int Number of subdivisions. tol : float - Tolerance for searching for existing data points at domain boundaries. - (Default: 1e-6) + Tolerance for searching for existing data points at domain + boundaries, relative to the width of a subdivision. (An absolute + tolerance would silently misbehave depending on the length unit of + the position data.) (Default: 1e-6) Returns ------- @@ -66,6 +68,9 @@ def checkerboard_detrend_profile(self, subdivisions, tol=1e-6): x, y = self.positions_and_heights() + # Convert the relative tolerance into an absolute one + atol = tol * (x[-1] - x[0]) / subdivisions + subdivided_line_scans = [] for i in range(subdivisions): # Subdivide interval @@ -82,7 +87,7 @@ def checkerboard_detrend_profile(self, subdivisions, tol=1e-6): # Put additional data points on the left and right boundaries, if there # is none already in the data set at exactly those points - if sub_ileft != 0 and sub_xleft < x[sub_ileft] - tol: + if sub_ileft != 0 and sub_xleft < x[sub_ileft] - atol: # Linear interpolation to boundary point sub_yleft = y[sub_ileft - 1] + (sub_xleft - x[sub_ileft - 1]) / ( x[sub_ileft] - x[sub_ileft - 1]) * ( @@ -91,7 +96,7 @@ def checkerboard_detrend_profile(self, subdivisions, tol=1e-6): sub_x = np.append([sub_xleft], sub_x) sub_y = np.append([sub_yleft], sub_y) - if sub_iright != len(x) and sub_xright > x[sub_iright - 1] + tol: + if sub_iright != len(x) and sub_xright > x[sub_iright - 1] + atol: # Linear interpolation to boundary point sub_yright = y[sub_iright - 1] + ( sub_xright - x[sub_iright - 1]) / ( @@ -151,14 +156,19 @@ def variable_bandwidth_from_profile(self, quantities='bh', reliable=True, resamp raise ValueError('`variable_bandwidth_from_profile` does not support resampling.') magnification = 1 - min_nb_grid_pts, = self.nb_grid_pts magnifications = [] bandwidths = [] rms_heights = [] - while min_nb_grid_pts >= nb_grid_pts_cutoff: + while True: subdivided_line_scans = self.checkerboard_detrend_profile(magnification) min_nb_grid_pts = min( [line.nb_grid_pts[0] for line in subdivided_line_scans]) + # Check the actual subdivision *before* recording it; otherwise the + # last recorded magnification can contain segments with fewer than + # `nb_grid_pts_cutoff` points (a detrended two-point segment has an + # rms height of exactly zero), biasing the smallest-bandwidth datum. + if min_nb_grid_pts < nb_grid_pts_cutoff: + break magnifications += [magnification] bandwidths += [subdivided_line_scans[0].physical_sizes[0]] rms_heights += [ diff --git a/SurfaceTopography/Nonuniform/__init__.py b/SurfaceTopography/Nonuniform/__init__.py index 4308ed94..429255cc 100644 --- a/SurfaceTopography/Nonuniform/__init__.py +++ b/SurfaceTopography/Nonuniform/__init__.py @@ -24,5 +24,5 @@ # """ -Module containing all functions operating on uniform topographies +Module containing all functions operating on nonuniform topographies """ diff --git a/SurfaceTopography/NonuniformLineScan.py b/SurfaceTopography/NonuniformLineScan.py index e3a2d3c5..d35451c9 100644 --- a/SurfaceTopography/NonuniformLineScan.py +++ b/SurfaceTopography/NonuniformLineScan.py @@ -139,11 +139,8 @@ def unit(self): else: return self._unit - @property - def info(self) -> dict: - info = self.parent_topography.info - info.update(self._info.model_dump(exclude_none=True)) - return info + # Note: `info` is merged with the parent topography's dictionary in the + # `DecoratedTopography` base class @property def physical_sizes(self): diff --git a/SurfaceTopography/Pipeline.py b/SurfaceTopography/Pipeline.py index 763fe7a7..b6a5a53e 100644 --- a/SurfaceTopography/Pipeline.py +++ b/SurfaceTopography/Pipeline.py @@ -78,10 +78,19 @@ def __setstate__(self, state): super().__setstate__(superstate) def __getattr__(self, name): - if name in self._kwargs: - return self._kwargs[name] - else: - return getattr(self.parent_topography, name) + # Expose keyword arguments as attributes; everything else + # goes through the regular dispatch of the parent class. + # (Delegating to `getattr(self.parent_topography, name)` + # here would silently apply chained pipeline functions to + # the *parent*, discarding this transformation.) + if name != '_kwargs': + try: + kwargs = object.__getattribute__(self, '_kwargs') + except AttributeError: + kwargs = {} + if name in kwargs: + return kwargs[name] + return super().__getattr__(name) def heights(self): return func(self.parent_topography, *self._args, **self._kwargs) diff --git a/SurfaceTopography/ScanningProbe/RigidScan.py b/SurfaceTopography/ScanningProbe/RigidScan.py index 40d3ae8b..6a85e78e 100644 --- a/SurfaceTopography/ScanningProbe/RigidScan.py +++ b/SurfaceTopography/ScanningProbe/RigidScan.py @@ -25,6 +25,7 @@ import numpy as np +from ..HeightContainer import NonuniformLineScanInterface, UniformTopographyInterface from ..NonuniformLineScan import NonuniformLineScan from ..UniformLineScanAndTopography import UniformLineScan @@ -50,26 +51,43 @@ def scan_with_rigid_sphere(topography, radius): raise ValueError("Only one-dimensional scans are supported at present.") positions, heights = topography.positions_and_heights() - scanned_heights = [] - for x in positions: - left = np.searchsorted(positions, x - radius) - right = np.searchsorted(positions, x + radius) - - # import matplotlib.pyplot as plt - # plt.figure() - # plt.plot(positions[left:right], heights[left:right], 'k-') - # plt.plot(positions[left:right], - np.sqrt(radius ** 2 - (positions[left:right] - x) ** 2), 'k--') - # plt.plot(positions[left:right], - # - np.sqrt(radius ** 2 - (positions[left:right] - x) ** 2) - heights[left:right], 'r-') - # plt.show() - - scanned_heights += [ - np.max( - heights[left:right] - + np.sqrt(radius**2 - (positions[left:right] - x) ** 2) - - radius - ) - ] + nb_pts = len(positions) + + # For each scan position, the tip contacts the data points within + # [x - radius, x + radius]; the scanned height is the maximum over that + # window of the height plus the local tip profile. + lefts = np.searchsorted(positions, positions - radius) + rights = np.searchsorted(positions, positions + radius) + widths = rights - lefts + + # The windows are evaluated block-wise: a block of scan positions is + # padded to its widest window and reduced with a masked max. The block + # size is chosen such that the temporary array stays below a fixed + # element count, so memory use is bounded irrespective of scan size + # and tip radius (a single scan-sized block could otherwise allocate + # nb_pts * max_window elements). + target_nb_elements = 2 ** 22 # 32 MB of doubles + scanned_heights = np.empty(nb_pts) + start = 0 + while start < nb_pts: + nb_block = min(nb_pts - start, + max(1, target_nb_elements // max(1, widths[start]))) + max_width = widths[start:start + nb_block].max() + # Shrinking the block can only shrink its widest window, so after + # this step nb_block * max_width <= target_nb_elements holds + nb_block = min(nb_block, max(1, target_nb_elements // max_width)) + max_width = widths[start:start + nb_block].max() + + block = slice(start, start + nb_block) + col = lefts[block].reshape(-1, 1) + np.arange(max_width).reshape(1, -1) + valid = col < rights[block].reshape(-1, 1) + col = np.minimum(col, nb_pts - 1) + distance = positions[col] - positions[block].reshape(-1, 1) + tip_heights = heights[col] + np.sqrt( + np.maximum(radius ** 2 - distance * distance, 0)) - radius + scanned_heights[block] = np.max( + np.where(valid, tip_heights, -np.inf), axis=1) + start += nb_block return scanned_heights @@ -90,6 +108,10 @@ def pipeline_scan_with_rigid_sphere(self, radius): topography : :obj:`SurfaceTopography.UniformLineScan` or :obj:`SurfaceTopography.NonuniformLineScan` Topography with scannned heights on the same grid as the topography. """ + if self.dim != 1: + raise ValueError( + "Scanning with a rigid sphere is only supported for line scans." + ) info_dict = dict( instrument=dict( name="Scanning rigid sphere simulation", @@ -109,7 +131,9 @@ def pipeline_scan_with_rigid_sphere(self, radius): unit=self.unit, info=info_dict, ) - elif isinstance(self, UniformLineScan): + elif self.is_uniform: + # Note: structural check rather than isinstance, so that decorated + # line scans (detrended, scaled, ...) can be scanned as well scanned_heights = scan_with_rigid_sphere(self, radius) return UniformLineScan( scanned_heights, @@ -118,18 +142,19 @@ def pipeline_scan_with_rigid_sphere(self, radius): unit=self.unit, info=info_dict, ) - elif isinstance(self, NonuniformLineScan): + else: scanned_heights = scan_with_rigid_sphere(self, radius) return NonuniformLineScan( self.positions(), scanned_heights, unit=self.unit, info=info_dict ) - else: - raise ValueError("Unexpected topography instance", type(self)) -UniformLineScan.register_function( +# Register on the interfaces so that decorated topographies (detrended, +# scaled, etc.) can also be scanned; the function itself checks that the +# data is one-dimensional +UniformTopographyInterface.register_function( "scan_with_rigid_sphere", pipeline_scan_with_rigid_sphere ) -NonuniformLineScan.register_function( +NonuniformLineScanInterface.register_function( "scan_with_rigid_sphere", pipeline_scan_with_rigid_sphere ) diff --git a/SurfaceTopography/Special.py b/SurfaceTopography/Special.py index bc1c749d..380d1952 100644 --- a/SurfaceTopography/Special.py +++ b/SurfaceTopography/Special.py @@ -93,20 +93,29 @@ def make_topography_from_function(fun, physical_sizes, # serial code nb_subdomain_grid_pts = nb_grid_pts nb_grid_pts = None - topography = Topography(heights=np.zeros(nb_subdomain_grid_pts), - physical_sizes=physical_sizes, - subdomain_locations=subdomain_locations, - nb_grid_pts=nb_grid_pts, - decomposition="subdomain", - **kwargs, - ) + # We need a topography object to compute the positions, but the heights + # are only known afterwards; construct a second topography with the + # actual heights so that the constructor's validation (array shape, + # conversion to float, masking of NaNs as undefined data) is applied. + # (Assigning to `_heights` directly would bypass all of this.) + dummy = Topography(heights=np.zeros(nb_subdomain_grid_pts), + physical_sizes=physical_sizes, + subdomain_locations=subdomain_locations, + nb_grid_pts=nb_grid_pts, + decomposition="subdomain", + **kwargs, + ) cx, cy = centre - x, y = topography.positions() + x, y = dummy.positions() - topography._heights = fun(x - cx, y - cy) - - return topography + return Topography(heights=fun(x - cx, y - cy), + physical_sizes=physical_sizes, + subdomain_locations=subdomain_locations, + nb_grid_pts=nb_grid_pts, + decomposition="subdomain", + **kwargs, + ) def make_sphere(radius, nb_grid_pts, physical_sizes, centre=None, @@ -214,18 +223,21 @@ def get_r(res, size, centre, subd_loc, subd_res): "Should be 'sphere' or 'paraboloid'".format(kind))) if dim == 1: - ret_top = UniformLineScan(h + offset, physical_sizes) + ret_top = UniformLineScan(h + offset, physical_sizes, + periodic=periodic) else: ret_top = Topography(h + offset, physical_sizes, + periodic=periodic, decomposition='subdomain', nb_grid_pts=nb_grid_pts, subdomain_locations=subdomain_locations, communicator=communicator) - if standoff == "undefined": - return ret_top - else: - return ret_top.fill_undefined_data(standoff_val) + # For `standoff="undefined"` the region outside the sphere contains NaNs, + # which the constructors above turn into masked (undefined) data points. + # For a numeric standoff the heights outside the sphere were already set + # to the (finite) standoff value above; no filling is necessary. + return ret_top class PlasticTopography(DecoratedUniformTopography): diff --git a/SurfaceTopography/Support/Bibliography.py b/SurfaceTopography/Support/Bibliography.py index 4f0d83c9..fd8d30d5 100644 --- a/SurfaceTopography/Support/Bibliography.py +++ b/SurfaceTopography/Support/Bibliography.py @@ -76,13 +76,22 @@ def func_with_doi(*args, **kwargs): doi._n = 0 doi.dois.update(self._add_these_dois) doi._n += 1 - retvals = func(*args, **kwargs) - doi._n -= 1 - if doi._n == 0: - # We have reached the point where the original `dois` argument - # was passed. Create a new set such that subsequent calls don't - # contaminate the bibliography. - doi.dois = set() - return retvals + try: + # The bookkeeping must be unwound even if the wrapped + # function raises; otherwise the class-level state stays + # permanently corrupted and all subsequent analysis calls + # silently dump DOIs into the caller's set. + return func(*args, **kwargs) + finally: + doi._n -= 1 + if doi._n == 0: + # We have reached the point where the original `dois` + # argument was passed. Create a new set such that + # subsequent calls don't contaminate the bibliography. + doi.dois = set() + # Mark the wrapper so that `register_function` can detect that a + # function already records DOIs (checking the function name does + # not work because functools.wraps restores the original name) + func_with_doi.__has_doi__ = True return func_with_doi diff --git a/SurfaceTopography/Support/Deprecation.py b/SurfaceTopography/Support/Deprecation.py index 58ce09e6..38af87a6 100644 --- a/SurfaceTopography/Support/Deprecation.py +++ b/SurfaceTopography/Support/Deprecation.py @@ -44,7 +44,10 @@ def _get_warn_str(version=version, alternative=alternative): @wraps(func) def deprecated_func(*args, **kwargs): - warnings.warn(_get_warn_str(version, alternative), DeprecationWarning) + # stacklevel=2 attributes the warning to the caller rather than + # to this wrapper + warnings.warn(_get_warn_str(version, alternative), DeprecationWarning, + stacklevel=2) return func(*args, **kwargs) docstring = deprecated_func.__doc__ or "" diff --git a/SurfaceTopography/Support/JSON.py b/SurfaceTopography/Support/JSON.py index f237c256..676ab27f 100644 --- a/SurfaceTopography/Support/JSON.py +++ b/SurfaceTopography/Support/JSON.py @@ -18,13 +18,25 @@ def nan_to_none(obj): elif isinstance(obj, list) or isinstance(obj, set): return [nan_to_none(v) for v in obj] elif isinstance(obj, MaskedArray): - return [None if m else nan_to_none(v) for v, m in zip(obj.data, obj.mask)] + # Note: `getmaskarray` always returns a full boolean array, while + # `obj.mask` can be the scalar `nomask` (which cannot be iterated) + mask = np.ma.getmaskarray(obj) + if obj.ndim == 0: + return None if bool(mask) else nan_to_none(obj.item()) + return [ + # Recurse row by row for multidimensional arrays + nan_to_none(np.ma.masked_array(v, mask=m)) if np.ndim(v) > 0 + else (None if m else nan_to_none(v)) + for v, m in zip(obj.data, mask) + ] elif isinstance(obj, np.ndarray) or isinstance(obj, ArrayImpl): if obj.ndim == 0: return nan_to_none(obj.item()) else: return [nan_to_none(v) for v in obj] - elif isinstance(obj, float) and np.isnan(obj): + elif isinstance(obj, float) and not np.isfinite(obj): + # NaN as well as +/-Inf have no JSON representation ('NaN' and + # 'Infinity' produced by the stdlib encoder are not valid JSON) return None return obj @@ -123,3 +135,10 @@ def encode(self, obj, *args, **kwargs): # https://stackoverflow.com/questions/28639953/python-json-encoder-convert-nans-to-null-instead obj = nan_to_none(obj) return super().encode(obj, *args, **kwargs) + + def iterencode(self, obj, *args, **kwargs): + # Note: `json.dump` calls `iterencode` directly and never goes + # through `encode`; without this override, `json.dump` with this + # encoder would emit invalid JSON (bare NaN/Infinity tokens) + obj = nan_to_none(obj) + return super().iterencode(obj, *args, **kwargs) diff --git a/SurfaceTopography/Support/Regression.py b/SurfaceTopography/Support/Regression.py index 9efafbde..4212f2b7 100644 --- a/SurfaceTopography/Support/Regression.py +++ b/SurfaceTopography/Support/Regression.py @@ -31,10 +31,11 @@ import logging import numpy as np +from scipy.linalg import cho_factor, cho_solve from ..Exceptions import NoReliableDataError -_log = logging.Logger(__name__) +_log = logging.getLogger(__name__) def make_grid(collocation, min_value, max_value, nb_points=None, nb_points_per_decade=10, dectol=0.01): @@ -57,7 +58,10 @@ def make_grid(collocation, min_value, max_value, nb_points=None, nb_points_per_d Maximum value. nb_points : int, optional Number of bins for averaging. Bins are automatically determined if set - to None. (Default: None) + to None. Caution: for 'log' collocation, this is the number of bin + *edges*, i.e. `nb_points` - 1 collocation points are returned, while + for 'linear' and 'quadratic' collocation it is the number of + collocation points (`nb_points` + 1 edges). (Default: None) nb_points_per_decade : int, optional Number of points per decade for log-spaced collocation points. (Default: None) @@ -274,11 +278,13 @@ def gaussian_process_regression(output_x, x, values, kernel=gaussian_kernel, noi # Add noise to observation covariance matrix obs_cov += noise_variance * np.identity(len(x)) - # Compute kernel coefficients - coeff = np.linalg.solve(obs_cov, values) + # The observation covariance matrix is symmetric positive definite; + # factorize it once and reuse the factorization for both the predictive + # mean and the predictive variance below + obs_cov_factor = cho_factor(obs_cov) - # Covariance between test outputs - test_cov = kernel(output_x.reshape(-1, 1), output_x.reshape(1, -1)) + # Compute kernel coefficients + coeff = cho_solve(obs_cov_factor, values) # Covariance between observation and test outputs obs_test_cov = kernel(x.reshape(-1, 1), output_x.reshape(1, -1)) @@ -286,11 +292,15 @@ def gaussian_process_regression(output_x, x, values, kernel=gaussian_kernel, noi # Compute predictive mean pred_mean = coeff.dot(obs_test_cov) - # Compute predictive covariance - pred_cov = test_cov - obs_test_cov.T.dot(np.linalg.solve(obs_cov, obs_test_cov)) + # Compute predictive variance. Only the diagonal of the predictive + # covariance is returned; computing the full matrix would additionally + # cost O(nb_output_points^2) memory. The kernels are elementwise + # functions, so kernel(output_x, output_x) is that diagonal. + pred_var = kernel(output_x, output_x) - \ + np.sum(obs_test_cov * cho_solve(obs_cov_factor, obs_test_cov), axis=0) # Return mean and variance - return pred_mean, pred_cov.diagonal() + return pred_mean, pred_var def resample(x, values, collocation='log', nb_points=None, min_value=None, max_value=None, nb_points_per_decade=10, diff --git a/SurfaceTopography/Support/UnitConversion.py b/SurfaceTopography/Support/UnitConversion.py index d48e6353..080df370 100644 --- a/SurfaceTopography/Support/UnitConversion.py +++ b/SurfaceTopography/Support/UnitConversion.py @@ -22,10 +22,12 @@ # SOFTWARE. # +import re + import numpy as np length_units = {'Gm': 1e9, 'Mm': 1e6, 'km': 1000.0, 'm': 1.0, 'mm': 1e-3, 'µm': 1e-6, 'um': 1e-6, 'nm': 1e-9, - 'Å': 1e-10, 'pm': 1e-12, 'fm': 1e-15} + 'Å': 1e-10, 'A': 1e-10, 'pm': 1e-12, 'fm': 1e-15} voltage_units = {'GV': 1e9, 'MV': 1e6, 'kV': 1000.0, 'V': 1.0, 'mV': 1e-3, 'µV': 1e-6, 'nV': 1e-9, 'pV': 1e-12, 'fV': 1e-15} @@ -178,8 +180,17 @@ def suggest_length_unit(scale, lower_in_meters, upper_in_meters): """ if scale == 'linear': v = max(abs(lower_in_meters), abs(upper_in_meters)) + if v == 0: + # All-zero data can be represented in any unit + return 'm' m10 = 3 * int(np.floor(np.log10(v) / 3)) elif scale == 'log': + if upper_in_meters <= 0: + # Nothing can be displayed on a log axis anyway + return 'm' + if lower_in_meters <= 0: + # Base the suggestion on the upper bound only + lower_in_meters = upper_in_meters u10 = int(np.ceil(np.log10(upper_in_meters))) l10 = int(np.floor(np.log10(lower_in_meters))) m10 = 3 * int(np.ceil((l10 + u10) / 6) - 1) @@ -241,7 +252,12 @@ def suggest_length_unit_for_data(scale, data, unit): def find_length_unit_in_string(s): """Check the string `s` contains any length information""" - for unit, normalized_unit in length_units_to_utf8.items(): - if s.find(unit) >= 0: - return normalized_unit + # Match whole tokens only. Substring matching would produce false + # positives, e.g. the alias 'A' (Angstrom) would match 'X Axis' and + # 'nm' embedded in a longer word; a false unit silently rescales data. + for token in re.split(r'[\s()\[\]{},;:=/-]+', s): + if token in length_units: + return mangle_length_unit_utf8(token) + if token in length_units_to_utf8: + return length_units_to_utf8[token] return None diff --git a/SurfaceTopography/Support/__init__.py b/SurfaceTopography/Support/__init__.py index f7707736..74124d6c 100644 --- a/SurfaceTopography/Support/__init__.py +++ b/SurfaceTopography/Support/__init__.py @@ -49,7 +49,9 @@ def fold_fft_half(arr, n): """ result = arr[:n // 2, ...] result[1:n // 2, ...] += arr[n - 1:(n + 1) // 2:-1, ...] - result /= 2 + # The entry at zero frequency appears just once in the FFT output and + # received no mirrored contribution above; it must not be halved. + result[1:n // 2, ...] /= 2 return result diff --git a/SurfaceTopography/Uniform/Autocorrelation.py b/SurfaceTopography/Uniform/Autocorrelation.py index dba9881e..afdf58e9 100644 --- a/SurfaceTopography/Uniform/Autocorrelation.py +++ b/SurfaceTopography/Uniform/Autocorrelation.py @@ -222,9 +222,13 @@ def autocorrelation_from_area(self, reliable=True, collocation='log', nb_points= nx, ny = self.nb_grid_pts sx, sy = self.physical_sizes - # The factor of two comes from the fact that the short cutoff is estimated - # from the curvature but the ACF is the slope, see 10.1016/j.apsadv.2021.100190 - short_cutoff = self.short_reliability_cutoff(np.mean(self.pixel_size)) if reliable else np.mean(self.pixel_size) + # The factor of two (below, in the `min_radius` arguments) comes from the + # fact that the short cutoff is estimated from the curvature but the ACF + # is the slope, see 10.1016/j.apsadv.2021.100190. The lower bound of two + # pixels mirrors the profile ACF (`autocorrelation_from_profile`), which + # uses the same convention. + short_cutoff = self.short_reliability_cutoff(2 * np.mean(self.pixel_size)) if reliable \ + else 2 * np.mean(self.pixel_size) # Compute FFT and normalize if self.is_periodic: @@ -239,7 +243,7 @@ def autocorrelation_from_area(self, reliable=True, collocation='log', nb_points= # Radial average r_val, r_edges, A_val, _ = resample_radial(A_xy, physical_sizes=(sx, sy), nb_points=nb_points, nb_points_per_decade=nb_points_per_decade, collocation=collocation, - min_radius=short_cutoff, max_radius=(sx + sy) / 4, + min_radius=short_cutoff / 2, max_radius=(sx + sy) / 4, method=resampling_method) else: p = self.heights() @@ -263,7 +267,7 @@ def autocorrelation_from_area(self, reliable=True, collocation='log', nb_points= # Radial average r_val, r_edges, A_val, _ = resample_radial(A_xy, physical_sizes=(sx, sy), collocation=collocation, nb_points=nb_points, nb_points_per_decade=nb_points_per_decade, - min_radius=short_cutoff, max_radius=(sx + sy) / 2, full=False, + min_radius=short_cutoff / 2, max_radius=(sx + sy) / 2, full=False, method=resampling_method) if return_map: diff --git a/SurfaceTopography/Uniform/BearingArea.py b/SurfaceTopography/Uniform/BearingArea.py index 5e22bf9c..10e8c9f9 100644 --- a/SurfaceTopography/Uniform/BearingArea.py +++ b/SurfaceTopography/Uniform/BearingArea.py @@ -123,13 +123,19 @@ def __init__(self, h, is_periodic): self._is_periodic = is_periodic if self._is_periodic: + # Note: `np.roll` needs explicit axes; with a tuple shift and no + # axis it would flatten the array and roll by the sum of the + # shifts. + h10 = np.roll(h, -1, axis=0) + h01 = np.roll(h, -1, axis=1) + h11 = np.roll(h, (-1, -1), axis=(0, 1)) self._el_min_heights = np.sort(np.ma.compressed([ - np.minimum(np.minimum(h, np.roll(h, (-1, 0))), np.roll(h, (0, -1))), - np.minimum(np.minimum(np.roll(h, (-1, -1)), np.roll(h, (-1, 0))), np.roll(h, (0, -1))) + np.minimum(np.minimum(h, h10), h01), + np.minimum(np.minimum(h11, h10), h01) ])) self._el_max_heights = np.sort(np.ma.compressed([ - np.maximum(np.maximum(h, np.roll(h, (-1, 0))), np.roll(h, (0, -1))), - np.maximum(np.maximum(np.roll(h, (-1, -1)), np.roll(h, (-1, 0))), np.roll(h, (0, -1))) + np.maximum(np.maximum(h, h10), h01), + np.maximum(np.maximum(h11, h10), h01) ])) else: self._el_min_heights = np.sort(np.ma.compressed([ diff --git a/SurfaceTopography/Uniform/Converters.py b/SurfaceTopography/Uniform/Converters.py index a013b153..a19073b1 100644 --- a/SurfaceTopography/Uniform/Converters.py +++ b/SurfaceTopography/Uniform/Converters.py @@ -57,7 +57,7 @@ def __init__(self, topography, info={}): # This is populated with functions from the nonuniform topography, but # this is a uniform topography - self._functions = NonuniformLineScan._functions + self._functions = NonuniformLineScan._all_functions() # Implement abstract methods of AbstractHeightContainer diff --git a/SurfaceTopography/Uniform/Derivative.py b/SurfaceTopography/Uniform/Derivative.py index 44b1c3bb..3f298bd8 100644 --- a/SurfaceTopography/Uniform/Derivative.py +++ b/SurfaceTopography/Uniform/Derivative.py @@ -31,7 +31,6 @@ from ..Exceptions import UndefinedDataError from ..HeightContainer import UniformTopographyInterface from ..Support import toiter -from ..UniformLineScanAndTopography import Topography class FourierDerivative: @@ -137,6 +136,29 @@ def fourier(self, phase): return self._operator.fourier(phase) +def _operator_direction(op): + """ + Return the Cartesian axis an axis-aligned derivative operator acts + along, or None if the direction cannot be determined (e.g. for mixed + stencils). + """ + direction = getattr(op, 'direction', None) + if direction is not None: + return direction + try: + stencil = np.asarray(op.stencil) + except AttributeError: + return None + # An axis-aligned operator has nonzero coefficients that vary along a + # single axis only. (The stencil array itself may be zero-padded, so + # its shape is not a reliable indicator.) + nonzero = np.nonzero(stencil) + axes = [axis for axis, indices in enumerate(nonzero) if len(set(indices)) > 1] + if len(axes) == 1: + return axes[0] + return None + + # # Stencils for first and second derivatives # @@ -330,8 +352,9 @@ def derivative( unity (line scan), then an array of the same shape as the topography is returned. Otherwise, the first array index contains the direction of the derivative. If the topgography is nonperiodic, - then all returning array with have shape one less than the input - arrays. + then the returned arrays are trimmed at the boundary where the + derivative stencil sticks out of the data region; the number of + trimmed points is the stencil extent times the scale factor. """ if self.physical_sizes is None: raise ValueError( @@ -446,7 +469,7 @@ def derivative( interpolation_required = ( np.any(s - s.astype(int) != 0) or interpolation == "fourier" ) - if interpolation_required and interpolation == "disabled": + if interpolation_required and interpolation == "disable": raise ValueError( "Interpolation is required to compute derivative at the " "desired scale but is explicitly disabled through the " @@ -513,8 +536,16 @@ def derivative( if not is_periodic: _der = trim_nonperiodic(_der, s, op) - # We need to divide by the grid spacing to make this a derivative - _der /= scaled_pixel_size[i] ** n + # We need to divide by the grid spacing to make this a + # derivative. The grid spacing is that of the direction the + # operator acts along; the operator index only coincides with + # the direction for the default one-operator-per-direction case + # (e.g. a single y-direction operator must not be normalized by + # the x grid spacing). + direction = _operator_direction(op) + if direction is None: + direction = i + _der /= scaled_pixel_size[direction] ** n # Mask array if it has NaNs if np.sum(~np.isfinite(_der)) > 0: @@ -592,8 +623,9 @@ def fourier_derivative(self, scale_factor=None, distance=None, mask_function=Non unity (line scan), then an array of the same shape as the topography is returned. Otherwise, the first array index contains the direction of the derivative. If the topgography is nonperiodic, - then all returning array with have shape one less than the input - arrays. + then the returned arrays are trimmed at the boundary where the + derivative stencil sticks out of the data region; the number of + trimmed points is the stencil extent times the scale factor. """ dim = self.dim return self.derivative( @@ -606,6 +638,8 @@ def fourier_derivative(self, scale_factor=None, distance=None, mask_function=Non ) -# Register analysis functions from this module +# Register analysis functions from this module. Note: `fourier_derivative` +# is registered on the interface (it is dimension-agnostic); registering on +# `Topography` would make it unavailable on decorated topographies. UniformTopographyInterface.register_function("derivative", derivative) -Topography.register_function("fourier_derivative", fourier_derivative) +UniformTopographyInterface.register_function("fourier_derivative", fourier_derivative) diff --git a/SurfaceTopography/Uniform/Detrending.py b/SurfaceTopography/Uniform/Detrending.py index e2df268c..09662e7d 100644 --- a/SurfaceTopography/Uniform/Detrending.py +++ b/SurfaceTopography/Uniform/Detrending.py @@ -46,6 +46,12 @@ def polyfit_line_scan(self, deg): """ x, h = self.positions_and_heights() x /= self.physical_sizes[0] + # Exclude undefined (masked) data points from the fit; np.polyfit on a + # masked array would silently fit the raw values under the mask + mask = np.ma.getmaskarray(h) + if mask.any(): + x = x[~mask] + h = np.asarray(h[~mask]) coeffs = np.polyfit(x, h, deg) return np.array(coeffs)[::-1] @@ -129,7 +135,7 @@ def polyfit2_topography(self, full_output=False): coeffs ordered as follows - {5} + {0} x + {1} y + {2} x^2 + {3} y^2 + {4} xy + {0} + {1} x + {2} y + {3} x^2 + {4} y^2 + {5} xy """ arr = self.heights() diff --git a/SurfaceTopography/Uniform/Filtering.py b/SurfaceTopography/Uniform/Filtering.py index 1f2e244d..03284697 100644 --- a/SurfaceTopography/Uniform/Filtering.py +++ b/SurfaceTopography/Uniform/Filtering.py @@ -81,6 +81,21 @@ def __setstate__(self, state): def nb_grid_pts(self): return tuple(n // f for n, f in zip(self.parent_topography.nb_grid_pts, self._factor)) + @property + def physical_sizes(self): + # The downsampled grid has a pixel size of `factor` times the parent + # pixel size. If the downsampling factor does not divide the number + # of grid points, trailing points are discarded and the physical size + # shrinks accordingly. (The pixel size reported by the base class is + # physical_sizes / nb_grid_pts, which then remains consistent with + # the actual sample spacing.) + return tuple( + n * f * p + for n, f, p in zip( + self.nb_grid_pts, self._factor, self.parent_topography.pixel_size + ) + ) + def heights(self): heights = self.parent_topography.heights() fx, fy = self._factor @@ -226,7 +241,7 @@ class FourierFilteredUniformTopography(DecoratedUniformTopography): name = 'filtered_topography' def __init__(self, topography, - filter_function=lambda qx, qy: (np.abs(qx) <= 1) * np.abs(qy) <= 1, + filter_function=lambda q: q <= 1, isotropic=True, info={}): @@ -271,9 +286,9 @@ def filter_function(self, *args): if self.dim == 2 and not self.is_filter_isotropic \ and len(args) != 2: - raise ("ValueError: qx, qy expected") + raise ValueError("qx, qy expected") elif self.dim == 1 and len(args) != 1: - raise ("ValueError: q expected") + raise ValueError("q expected") return self._filter_function(*args) diff --git a/SurfaceTopography/Uniform/GeometryAnalysis.py b/SurfaceTopography/Uniform/GeometryAnalysis.py index 50dd5a2f..c82150b3 100644 --- a/SurfaceTopography/Uniform/GeometryAnalysis.py +++ b/SurfaceTopography/Uniform/GeometryAnalysis.py @@ -110,8 +110,11 @@ def assign_patch_numbers_profile(mask, periodic): # Patches are odd numbers patch_ids += 1 if periodic and mask[-1]: - # Assign same patch id to first and last patch - patch_ids[patch_ids == patch_ids[-1]] = patch_ids[0] + # Assign same patch id to first and last patch. Note that the + # patch containing the first pixel always has (pre-transform) + # id 1 here; `patch_ids[0]` would be the id of the *second* + # pixel, which is not necessarily part of the first patch. + patch_ids[patch_ids == patch_ids[-1]] = 1 # Patches are odd numbers, set even numbers to zero patch_ids += 1 patch_ids[(patch_ids & 0x1).astype(bool)] = 0 diff --git a/SurfaceTopography/Uniform/Imputation.py b/SurfaceTopography/Uniform/Imputation.py index fbf7d373..80719f2c 100644 --- a/SurfaceTopography/Uniform/Imputation.py +++ b/SurfaceTopography/Uniform/Imputation.py @@ -101,37 +101,42 @@ def heights(self): pixel_index[patch_mask] = np.arange(nb_patch) pixel_index[perimeter_mask] = np.arange(nb_patch, nb_pixels) - # Assemble Laplace matrix; diagonal terms - i0 = np.arange(nb_patch) - j0 = np.arange(nb_patch) - - # Off-diagonal terms - i1 = pixel_index[patch_mask] - j1 = np.roll(pixel_index, 1, 0)[patch_mask] - i2 = pixel_index[patch_mask] - j2 = np.roll(pixel_index, -1, 0)[patch_mask] - + # Assemble Laplace matrix. Each patch pixel couples to its + # nearest neighbors; for nonperiodic topographies, stencil + # legs that cross the domain boundary must be dropped (and + # the diagonal reduced correspondingly, yielding a natural + # boundary condition). An unconditional `np.roll` would wrap + # around the boundary and couple edge pixels to unrelated + # pixels (whose `pixel_index` entry is zero, aliasing patch + # unknown #0). + legs = [(1, 0), (-1, 0)] if dim == 2: - i3 = pixel_index[patch_mask] - j3 = np.roll(pixel_index, 1, 1)[patch_mask] - i4 = pixel_index[patch_mask] - j4 = np.roll(pixel_index, -1, 1)[patch_mask] - - # Laplace matrix from coordinates - laplace = scipy.sparse.coo_matrix( - (np.concatenate((-4 * np.ones(nb_patch), np.ones(nb_patch), np.ones(nb_patch), - np.ones(nb_patch), np.ones(nb_patch), np.ones(nb_perimeter))), - (np.concatenate((i0, i1, i2, i3, i4, np.arange(nb_patch, nb_pixels))), - np.concatenate((j0, j1, j2, j3, j4, np.arange(nb_patch, nb_pixels))))), - shape=(nb_pixels, nb_pixels)) - else: - # Laplace matrix from coordinates - laplace = scipy.sparse.coo_matrix( - (np.concatenate((-2 * np.ones(nb_patch), np.ones(nb_patch), np.ones(nb_patch), - np.ones(nb_perimeter))), - (np.concatenate((i0, i1, i2, np.arange(nb_patch, nb_pixels))), - np.concatenate((j0, j1, j2, np.arange(nb_patch, nb_pixels))))), - shape=(nb_pixels, nb_pixels)) + legs += [(1, 1), (-1, 1)] + + diagonal = np.zeros(nb_patch) + rows = [] + cols = [] + for shift, axis in legs: + neighbor = np.roll(pixel_index, shift, axis) + valid = np.ones_like(patch_mask) + if not self.is_periodic: + # `np.roll(a, 1)[p]` is the neighbor a[p - 1], so + # for shift == 1 the wrapped (invalid) entries are + # at the low edge of the axis, and at the high edge + # for shift == -1 + edge = [slice(None)] * dim + edge[axis] = 0 if shift == 1 else -1 + valid[tuple(edge)] = False + sel = np.logical_and(patch_mask, valid) + rows += [pixel_index[sel]] + cols += [neighbor[sel]] + diagonal[pixel_index[sel]] -= 1 + + laplace = scipy.sparse.coo_matrix( + (np.concatenate([diagonal] + [np.ones(len(r)) for r in rows] + [np.ones(nb_perimeter)]), + (np.concatenate([np.arange(nb_patch)] + rows + [np.arange(nb_patch, nb_pixels)]), + np.concatenate([np.arange(nb_patch)] + cols + [np.arange(nb_patch, nb_pixels)]))), + shape=(nb_pixels, nb_pixels)) # Dirichlet boundary conditions (heights on perimeter) rhs = np.zeros(nb_pixels) diff --git a/SurfaceTopography/Uniform/Integration.py b/SurfaceTopography/Uniform/Integration.py index 35e5daa6..1e88e7e2 100644 --- a/SurfaceTopography/Uniform/Integration.py +++ b/SurfaceTopography/Uniform/Integration.py @@ -21,6 +21,8 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE. # +import inspect + import numpy as np from ..Exceptions import NoReliableDataError, UndefinedDataError @@ -78,10 +80,17 @@ def integrate_psd(self, factor=lambda q: 1, window=None, reliable=True, ): raise NoReliableDataError('Dataset contains no reliable data.') C_raw = C_raw * mask - qvec = self.fftfreq() + # Dispatch on the number of arguments the factor accepts rather than + # catching TypeError, which would silently mask genuine TypeErrors + # raised inside a user-supplied factor function try: + nb_factor_args = len(inspect.signature(factor).parameters) + except (TypeError, ValueError): + nb_factor_args = 1 + if self.dim == 2 and nb_factor_args >= 2: + qvec = self.fftfreq() return np.sum(C_raw * factor(*qvec)) / np.prod(self.physical_sizes) - except TypeError: + else: return np.sum(C_raw * factor(q)) / np.prod(self.physical_sizes) @@ -94,7 +103,7 @@ def integrate_psd_from_profile(self, factor=lambda qx: 1, window=None, reliable= .. math:: - \frac{1}{2 \pi} \int_0^\infty dq_x factor(q_x) C^{1D}(q_x) + \frac{1}{2 \pi} \int_{-\infty}^\infty dq_x factor(q_x) C^{1D}(q_x) Discrete @@ -173,7 +182,7 @@ def moment_power_spectrum(self, order=0, window=None, reliable=True, ): .. math:: - \frac{1}{2 \pi} \int_0^\infty dq_x |q|^{\alpha} C^{1D}(q_x) + \frac{1}{2 \pi} \int_{-\infty}^\infty dq_x |q|^{\alpha} C^{1D}(q_x) Discrete diff --git a/SurfaceTopography/Uniform/Interpolation.py b/SurfaceTopography/Uniform/Interpolation.py index 28413429..285ac870 100644 --- a/SurfaceTopography/Uniform/Interpolation.py +++ b/SurfaceTopography/Uniform/Interpolation.py @@ -50,7 +50,10 @@ def linear_interpolator_line_scan(x, periodic=None): if not periodic: if np.any(scaled_x < 0) or np.any(scaled_x > nx-1): raise ValueError('Cannot interpolate outside of physical domain for nonperiodic line scans.') - int_x = np.array(scaled_x, dtype=int) + # Note: `floor` (not truncation towards zero) is required so that + # negative positions on periodic line scans interpolate rather than + # extrapolate + int_x = np.array(np.floor(scaled_x), dtype=int) frac_x = scaled_x - int_x return (1 - frac_x) * heights[int_x % nx] + frac_x * heights[(int_x + 1) % nx] @@ -62,8 +65,11 @@ def linear_interpolator_topography(x, y, periodic=None): if not periodic: if np.any(scaled_x < 0) or np.any(scaled_x > nx-1) or np.any(scaled_y < 0) or np.any(scaled_y > ny-1): raise ValueError('Cannot interpolate outside of physical domain for nonperiodic topographies.') - int_x = np.array(scaled_x, dtype=int) - int_y = np.array(scaled_y, dtype=int) + # Note: `floor` (not truncation towards zero) is required so that + # negative positions on periodic topographies interpolate rather + # than extrapolate + int_x = np.array(np.floor(scaled_x), dtype=int) + int_y = np.array(np.floor(scaled_y), dtype=int) frac_x = scaled_x - int_x frac_y = scaled_y - int_y @@ -164,7 +170,7 @@ def bicubic_interpolator_topography(x, y, derivative=0): interp_derxx = interp_derxx / dx ** 2 interp_deryy = interp_deryy / dy ** 2 - interp_derxy = interp_derxx / dx / dy + interp_derxy = interp_derxy / dx / dy return interp_field, interp_derx, interp_dery, interp_derxx, interp_deryy, interp_derxy @@ -201,13 +207,17 @@ def interpolate_fourier(self, nb_grid_pts): # the entries at the nyquist frequency are the superposition of the # positive and negative frequency. When we increase the fourier domain, we # will split that again between positive and negative, therefore we divide - # the value by two and copy it to the corresponding negative vector - if ny % 2 == 0: + # the value by two and copy it to the corresponding negative vector. + # Note: this splitting must only happen in directions that are actually + # enlarged; if the grid size in a direction is unchanged, the entry + # remains the (folded) Nyquist entry of the output spectrum and halving + # it would corrupt the data. + if ny % 2 == 0 and nb_grid_pts[1] > ny: # the nyquist frequency also contains the symmetric and will be # twice too big when the symmetric will be included smallspectrum[:, -1] = smallspectrum[:, -1] / 2 - if nx % 2 == 0: + if nx % 2 == 0 and nb_grid_pts[0] > nx: # the nyquist frequency also contains the symmetric and will be # twice too big when the symmetric will be included smallspectrum[int(nx / 2), :] = smallspectrum[int(nx / 2), :] / 2 @@ -223,7 +233,7 @@ def interpolate_fourier(self, nb_grid_pts): # To ensure the same behaviour in the x direction # (see test_fourier_interpolate_transpose_symmetry`)` # we have to mirror the entries at the niquist frequency by hand - if snx % 2 == 0: + if snx % 2 == 0 and nb_grid_pts[0] > nx: bigspectrum[i, :sny] = smallspectrum[i, :sny] return Topography(np.fft.irfft2(bigspectrum, s=nb_grid_pts) @@ -240,6 +250,9 @@ class MirrorStichedTopography(DecoratedUniformTopography): """ def __init__(self, parent_topography, info={}): + if parent_topography.dim != 2: + raise ValueError('Mirror stitching is only implemented for ' + 'topographies (not line scans).') if parent_topography.communicator.Get_size() > 1: raise (NotImplementedError("MirrorStichedTopography " "not domain decomposable")) @@ -257,6 +270,16 @@ def physical_sizes(self): def nb_grid_pts(self): return [2 * s for s in self.parent_topography.nb_grid_pts] + @property + def nb_subdomain_grid_pts(self): + # This class is serial only (see constructor); the subdomain is the + # full (mirrored) domain, not the parent's subdomain + return self.nb_grid_pts + + @property + def subdomain_locations(self): + return (0, 0) + def heights(self): h = self.parent_topography.heights() return np.block([[h[:, :], h[:, ::-1]], @@ -264,15 +287,17 @@ def heights(self): def positions(self): nx, ny = self.nb_grid_pts - lnx, lny = self.nb_subdomain_grid_pts sx, sy = self.physical_sizes - return np.meshgrid( - (self.subdomain_locations[0] + np.arange(lnx)) * sx / nx, - (self.subdomain_locations[1] + np.arange(lny)) * sy / ny, - indexing='ij') - - -Topography.register_function("mirror_stitch", MirrorStichedTopography) -Topography.register_function("interpolate_linear", interpolate_linear) -Topography.register_function("interpolate_bicubic", interpolate_bicubic) + return np.meshgrid(np.arange(nx) * sx / nx, + np.arange(ny) * sy / ny, + indexing='ij') + + +# Note: these are registered on the interface (rather than on `Topography`) +# because decorated topographies do not derive from `Topography`; the +# 2D-only functions carry explicit dimension checks with clear error +# messages for line scans +UniformTopographyInterface.register_function("mirror_stitch", MirrorStichedTopography) +UniformTopographyInterface.register_function("interpolate_linear", interpolate_linear) +UniformTopographyInterface.register_function("interpolate_bicubic", interpolate_bicubic) UniformTopographyInterface.register_function('interpolate_fourier', interpolate_fourier) diff --git a/SurfaceTopography/Uniform/PowerSpectrum.py b/SurfaceTopography/Uniform/PowerSpectrum.py index c3f8ebd4..2e84e1dd 100644 --- a/SurfaceTopography/Uniform/PowerSpectrum.py +++ b/SurfaceTopography/Uniform/PowerSpectrum.py @@ -126,7 +126,12 @@ def power_spectrum_from_profile(self, window=None, reliable=True, resampling_met if self.dim == 2: q = np.resize(q, (C_all.shape[1], q.shape[0])).T.ravel() C_all = np.ravel(C_all) - q, _, C, _ = resample(q, C_all, min_value=q[1], max_value=short_cutoff, collocation=collocation, + # The smallest positive wavevector limits the resampling range; note + # that q = 0 has either been stripped above (log collocation) or + # carries the (uninteresting) mean of the topography. Using q[1] + # here would discard the fundamental mode when q = 0 has already + # been stripped. + q, _, C, _ = resample(q, C_all, min_value=np.min(q[q > 0]), max_value=short_cutoff, collocation=collocation, nb_points=nb_points, nb_points_per_decade=nb_points_per_decade, method=resampling_method) return q, C @@ -183,7 +188,10 @@ def power_spectrum_from_area(self, window=None, reliable=True, collocation='log' nx, ny = self.nb_grid_pts sx, sy = self.physical_sizes - qmax = 2 * np.pi * nx / (2 * sx) + # Radial averaging is meaningful up to the smaller of the two Nyquist + # frequencies; using only the x-Nyquist would make the result depend on + # the orientation of the map for anisotropic pixels + qmax = np.pi * min(nx / sx, ny / sy) if reliable: # Update qmax diff --git a/SurfaceTopography/Uniform/ScalarParameters.py b/SurfaceTopography/Uniform/ScalarParameters.py index 138fcd18..9bca8741 100644 --- a/SurfaceTopography/Uniform/ScalarParameters.py +++ b/SurfaceTopography/Uniform/ScalarParameters.py @@ -57,14 +57,11 @@ def Rq(topography): if topography.is_domain_decomposed: raise NotImplementedError("`Rq` does not support MPI-decomposed topographies.") - n = np.prod(topography.nb_grid_pts) - reduction = Reduction(topography._communicator) + # Note: all means and sums must be normalized by the number of *defined* + # (unmasked) data points; masked points do not contribute to the sums profile = topography.heights() return np.sqrt( - reduction.sum( - (profile - reduction.sum(profile, axis=0) / topography.nb_grid_pts[0]) ** 2 - ) - / n + np.ma.sum((profile - np.ma.mean(profile, axis=0)) ** 2) / np.ma.count(profile) ) @@ -89,10 +86,11 @@ def Sq(topography): "Areal rms height can only be computed for topographies, not line scans." ) elif topography.dim == 2: - n = np.prod(topography.nb_grid_pts) - reduction = Reduction(topography._communicator) + # Normalize by the number of *defined* (unmasked) data points profile = topography.heights() - return np.sqrt(reduction.sum((profile - reduction.sum(profile) / n) ** 2) / n) + return np.sqrt( + np.ma.sum((profile - np.ma.mean(profile)) ** 2) / np.ma.count(profile) + ) else: raise ValueError(f"Cannot handle topographies of dimension {topography.dim}") diff --git a/SurfaceTopography/Uniform/ScanLineAlignment.py b/SurfaceTopography/Uniform/ScanLineAlignment.py index 009dfba9..eba72442 100644 --- a/SurfaceTopography/Uniform/ScanLineAlignment.py +++ b/SurfaceTopography/Uniform/ScanLineAlignment.py @@ -63,9 +63,12 @@ def __init__(self, topography, direction='x', mode='median', degree=1, topography : Topography 2D topography to align. direction : str, optional - Scan line direction: 'x' (rows, default) or 'y' (columns). - 'x' means each row is a scan line (fast scan in x-direction). - For AFM, 'x' is typically the fast scan direction. + Direction along which the scan lines run: 'x' (default) or 'y'. + 'x' means each scan line runs along the x direction (fast scan + in x), i.e. each line has constant y. For AFM, 'x' is typically + the fast scan direction. Note that the first index of the height + array is the x index in this library, so scan lines along x are + *columns* of the height array. mode : str, optional Method for computing line centers during alignment: - 'median': Use median for robust offset alignment (default). @@ -85,6 +88,8 @@ def __init__(self, topography, direction='x', mode='median', degree=1, "2D topographies") if degree < 0: raise ValueError("Polynomial degree must be non-negative") + if direction not in ('x', 'y'): + raise ValueError("Direction must be 'x' or 'y'") super().__init__(topography, info=info) self._direction = direction self._mode = mode @@ -167,8 +172,11 @@ def _compute_alignment(self, heights): heights : ndarray 2D array of height values. """ - # Transpose if aligning columns instead of rows - if self._direction == 'y': + # The first array index is the x index in this library. A scan line + # running along the x direction therefore corresponds to a *column* + # heights[:, j] of the height array; transpose so that the loop below + # (which iterates over the first index) walks over scan lines. + if self._direction == 'x': heights = heights.T nx, ny = heights.shape @@ -237,8 +245,9 @@ def heights(self): if self._line_coeffs is None: self._compute_alignment(heights) - # Transpose if aligning columns instead of rows - if self._direction == 'y': + # See note in `_compute_alignment`: scan lines along x are columns + # of the height array + if self._direction == 'x': heights = heights.T nx, ny = heights.shape @@ -253,7 +262,7 @@ def heights(self): heights[i, :] += self._offsets[i] # Transpose back if we transposed earlier - if self._direction == 'y': + if self._direction == 'x': heights = heights.T return heights diff --git a/SurfaceTopography/Uniform/VariableBandwidth.py b/SurfaceTopography/Uniform/VariableBandwidth.py index e6347200..51d883e8 100644 --- a/SurfaceTopography/Uniform/VariableBandwidth.py +++ b/SurfaceTopography/Uniform/VariableBandwidth.py @@ -85,10 +85,19 @@ def checkerboard_detrend_profile(self, subdivisions, order=1, return_plane=False x = x.reshape(-1) region_index = region_index.reshape(-1) - b = np.array([np.bincount(region_index, h * (x ** i)) for i in range(order + 1)]) - C = np.array([[np.bincount(region_index, x ** (k + i)) for i in range(order + 1)] for k in range(order + 1)]) + # Undefined (masked) data points must not contribute to the fit; they + # enter both the right-hand side and the normal-equation matrix with + # zero weight. (`np.bincount` silently strips masks and would otherwise + # ingest the raw values stored under the mask.) + w = 1.0 - np.ma.getmaskarray(h) + hf = np.ma.filled(h, 0) + + b = np.array([np.bincount(region_index, w * hf * (x ** i)) for i in range(order + 1)]) + C = np.array([[np.bincount(region_index, w * x ** (k + i)) for i in range(order + 1)] for k in range(order + 1)]) a = np.linalg.solve(C.T, b.T.reshape(b.T.shape + (1,))).T[0] + # Subtracting from the original (possibly masked) heights preserves the + # mask on undefined data points detrended_h = h - np.sum([a[i, region_index] * x ** i for i in range(order + 1)], axis=0) detrended_h.shape = shape @@ -166,8 +175,13 @@ def checkerboard_detrend_area(self, subdivisions, order=1, return_plane=False): y = y.reshape(-1) region_index = region_index.reshape(-1) - b = np.array([np.bincount(region_index, h * (x ** i) * (y ** j)) for i, j in ij]) - C = np.array([[np.bincount(region_index, x ** (k + i) * y ** (l + j)) for i, j in ij] for k, l in ij]) + # See `checkerboard_detrend_profile`: masked data points enter the fit + # with zero weight + w = 1.0 - np.ma.getmaskarray(h) + hf = np.ma.filled(h, 0) + + b = np.array([np.bincount(region_index, w * hf * (x ** i) * (y ** j)) for i, j in ij]) + C = np.array([[np.bincount(region_index, w * x ** (k + i) * y ** (l + j)) for i, j in ij] for k, l in ij]) a = np.linalg.solve(C.T, b.T.reshape(b.T.shape + (1,))).T[0] detrended_h = h - np.sum([a[k, region_index] * (x ** i) * (y ** j) for k, (i, j) in enumerate(ij)], axis=0) @@ -196,7 +210,7 @@ def variable_bandwidth_from_profile(self, quantities='bh', reliable=True, resamp - 'm': Magnification (Unit: dimensionless) - 'b': Bandwidth (Unit: length) - 'h': Statistical property (Unit: default length, see func) - - 's': RMS detrending slope (Unit: dimensionless) + - 'g': RMS detrending slope/gradient (Unit: dimensionless) For example, 'mbh' return a tuple with the three entries magnification, bandwidth, rms height. (Default: 'bh') @@ -317,7 +331,7 @@ def variable_bandwidth_from_area(self, quantities='bh', reliable=True, resamplin corresponding to the respective magnification. """ if resampling_method is not None: - raise ValueError('`variable_bandwidth_from_profile` does not support resampling.') + raise ValueError('`variable_bandwidth_from_area` does not support resampling.') magnification = 1 physical_sizes = np.array(self.physical_sizes) diff --git a/SurfaceTopography/UniformLineScanAndTopography.py b/SurfaceTopography/UniformLineScanAndTopography.py index 2aaaaf7a..88875d4a 100644 --- a/SurfaceTopography/UniformLineScanAndTopography.py +++ b/SurfaceTopography/UniformLineScanAndTopography.py @@ -478,7 +478,14 @@ def dim(self): @property def pixel_size(self): - return self.parent_topography.pixel_size + # Derive the pixel size from the physical sizes and number of grid + # points of the decorated topography itself, not of the parent. + # Decorators that change the geometry (e.g. transposition or + # downsampling) override `physical_sizes` and `nb_grid_pts`, and the + # pixel size must remain consistent with those. + return tuple( + s / n for s, n in zip(self.physical_sizes, self.nb_grid_pts) + ) @property def unit(self): @@ -487,11 +494,8 @@ def unit(self): else: return self._unit - @property - def info(self): - info = self.parent_topography.info - info.update(self._info.model_dump(exclude_none=True)) - return info + # Note: `info` is merged with the parent topography's dictionary in the + # `DecoratedTopography` base class @property def physical_sizes(self): @@ -518,7 +522,7 @@ def subdomain_slices(self): @property def area_per_pt(self): - return self.parent_topography.area_per_pt + return np.prod(self.pixel_size) def positions(self, **kwargs): return self.parent_topography.positions(**kwargs) @@ -575,21 +579,11 @@ def height_scale_factor(self): def position_scale_factor(self): return get_unit_conversion_factor(self.parent_topography.unit, self.unit) - @property - def pixel_size(self): - """Compute rescaled pixel sizes.""" - return tuple(self.position_scale_factor * s for s in super().pixel_size) - @property def physical_sizes(self): """Compute rescaled physical sizes.""" return tuple(self.position_scale_factor * s for s in super().physical_sizes) - @property - def area_per_pt(self): - """Compute rescaled physical sizes.""" - return np.prod(self.pixel_size) - def positions(self, **kwargs): """Compute the rescaled positions.""" if self.dim == 1: @@ -700,11 +694,29 @@ def physical_sizes(self): sx, sy = self.parent_topography.physical_sizes return sy, sx + @property + def nb_subdomain_grid_pts(self): + if self.dim == 1: + return self.parent_topography.nb_subdomain_grid_pts + else: + nx, ny = self.parent_topography.nb_subdomain_grid_pts + return ny, nx + + @property + def subdomain_locations(self): + if self.dim == 1: + return self.parent_topography.subdomain_locations + else: + ix, iy = self.parent_topography.subdomain_locations + return iy, ix + def heights(self): """Computes the rescaled profile.""" return self.parent_topography.heights().T def positions(self, **kwargs): + if self.dim == 1: + return self.parent_topography.positions(**kwargs) X, Y = self.parent_topography.positions(**kwargs) return Y.T, X.T @@ -736,14 +748,19 @@ def offset(self): return self._offset @offset.setter - def offset(self, offset, offsety=None): - if offsety is None: - self._offset = offset - else: - self._offset = (offset, offsety) + def offset(self, offset): + # Note: property setters receive exactly one value; a second + # parameter would be unreachable + self._offset = offset def heights(self): """Computes the translated profile.""" + if self.dim == 1: + # Line scans have a single offset; accept scalars as well as + # (possibly longer, from the default value) tuples + offset = self.offset + offsetx = offset[0] if np.ndim(offset) > 0 else offset + return np.roll(self.parent_topography.heights(), offsetx, axis=0) offsetx, offsety = self.offset return np.roll( np.roll(self.parent_topography.heights(), offsetx, axis=0), offsety, axis=1 @@ -764,7 +781,7 @@ def __init__(self, topography_a, topography_b, info={}): super().__init__(topography_a, info=info) - def combined_val(prop_a, prop_b, propname): + def check_combined_val(prop_a, prop_b, propname): """ topographies can have a fixed or dynamic, adaptive nb_grid_pts (or other attributes). This function assures that -- if this function @@ -775,20 +792,19 @@ def combined_val(prop_a, prop_b, propname): prop_b -- field of other topography propname -- field identifier (for error messages only) """ - if prop_a is None: - return prop_b - else: - if prop_b is not None: - assert prop_a == prop_b, "{} incompatible:{} <-> {}".format( - propname, prop_a, prop_b - ) - return prop_a + # Note: raise a proper exception; an `assert` would vanish under + # `python -O` and allow silent superposition of incompatible + # grids + if prop_a is not None and prop_b is not None and prop_a != prop_b: + raise ValueError( + "{} incompatible: {} <-> {}".format(propname, prop_a, prop_b) + ) - self._dim = combined_val(topography_a.dim, topography_b.dim, "dim") - self._nb_grid_pts = combined_val( + check_combined_val(topography_a.dim, topography_b.dim, "dim") + check_combined_val( topography_a.nb_grid_pts, topography_b.nb_grid_pts, "nb_grid_pts" ) - self._size = combined_val( + check_combined_val( topography_a.physical_sizes, topography_b.physical_sizes, "physical_sizes" ) self.parent_topography_a = topography_a diff --git a/cpp/autocorrelation.cpp b/cpp/autocorrelation.cpp index f6ee1cbd..cc7d59af 100644 --- a/cpp/autocorrelation.cpp +++ b/cpp/autocorrelation.cpp @@ -38,10 +38,10 @@ SOFTWARE. std::tuple nonuniform_autocorrelation( - Eigen::Ref x, - Eigen::Ref h, + Eigen::Ref x, + Eigen::Ref h, double physical_size, - std::optional> distances_opt) + std::optional> distances_opt) { const auto nb_grid_pts = x.size(); @@ -53,6 +53,13 @@ std::tuple nonuniform_autocorrelation( Eigen::ArrayXd distances; if (distances_opt) { distances = *distances_opt; + // The normalization below divides by (physical_size - distance); + // distances at or beyond the physical size have no data and would + // silently produce NaNs or wrong-signed values + if ((distances >= physical_size).any()) { + throw std::invalid_argument( + "All distances must be smaller than the physical size of the line scan."); + } } else { distances = Eigen::ArrayXd::LinSpaced(nb_grid_pts, 0.0, physical_size * (nb_grid_pts - 1) / nb_grid_pts); } @@ -75,10 +82,17 @@ std::tuple nonuniform_autocorrelation( double b = (b1 + b2) / 2; double db = (b2 - b1) / 2; if (db > 0) { - // f1[x_] := (h1 + s1*(x - x1)) - // f2[x_] := (h2 + s2*(x - x2)) - // FullSimplify[Integrate[f1[x]*f2[x + d], {x, b - db, b + db}]] - // = 2 * f1[b] * f2[b + d] * db + 2 * s1 * s2 * db ** 3 / 3 + // This accumulates the height-difference autocorrelation + // A(d) = (1/2) < [h(x+d) - h(x)]^2 >, + // i.e. the conventional factor 1/2 is folded into the + // expression below. With + // f1[x_] := (h1 + s1*(x - x1)) + // f2[x_] := (h2 + s2*(x - x2)) + // the difference on the overlap is + // f2[x + d] - f1[x] = z - (s1 - s2)*(x - b) + // where z = f2[b + d] - f1[b], and + // (1/2) Integrate[(f2[x + d] - f1[x])^2, {x, b - db, b + db}] + // = db*z^2 + (s1 - s2)^2 * db^3 / 3 double z = h2 - s2 * x2 + (b + distances(k)) * s2 - h1 + s1 * x1 - b * s1; double ds = s1 - s2; acf(k) += (db * (3 * z * z + ds * ds * db * db)) / 3; @@ -93,3 +107,53 @@ std::tuple nonuniform_autocorrelation( return {distances, acf}; } + + +Eigen::ArrayXd nonuniform_height_height_autocorrelation( + Eigen::Ref x, + Eigen::Ref h, + Eigen::Ref distances) +{ + const auto nb_grid_pts = x.size(); + + if (h.size() != nb_grid_pts) { + throw std::runtime_error("x- and h-arrays must contain identical number of data points."); + } + + const auto nb_distance_pts = distances.size(); + Eigen::ArrayXd acf = Eigen::ArrayXd::Zero(nb_distance_pts); + + /* This is the product (height-height) autocorrelation + * A(d) = Integrate[h(x) h(x + d)] + * over the overlap of the two piecewise-linear segments, without the + * factor 1/2 and without normalization by the overlap length. It mirrors + * the pure Python loop it replaces (Nonuniform/Autocorrelation.py) but + * runs in constant additional memory: + * f1[x_] := (h1 + s1*(x - x1)) + * f2[x_] := (h2 + s2*(x - x2)) + * Integrate[f1[x]*f2[x + d], {x, b - db, b + db}] + * = 2 * f1[b] * f2[b + d] * db + 2 * s1 * s2 * db^3 / 3 + */ + for (Eigen::Index i = 0; i < nb_grid_pts - 1; ++i) { + double x1 = x(i); + double h1 = h(i); + double s1 = (h(i + 1) - h1) / (x(i + 1) - x1); + for (Eigen::Index j = 0; j < nb_grid_pts - 1; ++j) { + double x2 = x(j); + double h2 = h(j); + double s2 = (h(j + 1) - h2) / (x(j + 1) - x2); + for (Eigen::Index k = 0; k < nb_distance_pts; ++k) { + double b1 = std::max(x1, x2 - distances(k)); + double b2 = std::min(x(i + 1), x(j + 1) - distances(k)); + double b = (b1 + b2) / 2; + double db = (b2 - b1) / 2; + if (db > 0) { + acf(k) += 2 * (h1 + s1 * (b - x1)) * (h2 + s2 * (b + distances(k) - x2)) * db + + 2 * s1 * s2 * db * db * db / 3; + } + } + } + } + + return acf; +} diff --git a/cpp/autocorrelation.h b/cpp/autocorrelation.h index 1effaf4d..2d623998 100644 --- a/cpp/autocorrelation.h +++ b/cpp/autocorrelation.h @@ -39,9 +39,14 @@ SOFTWARE. #include std::tuple nonuniform_autocorrelation( - Eigen::Ref x, - Eigen::Ref h, + Eigen::Ref x, + Eigen::Ref h, double physical_size, - std::optional> distances = std::nullopt); + std::optional> distances = std::nullopt); + +Eigen::ArrayXd nonuniform_height_height_autocorrelation( + Eigen::Ref x, + Eigen::Ref h, + Eigen::Ref distances); #endif diff --git a/cpp/bearing_area.cpp b/cpp/bearing_area.cpp index 8e489b24..e36a8588 100644 --- a/cpp/bearing_area.cpp +++ b/cpp/bearing_area.cpp @@ -21,13 +21,12 @@ SOFTWARE. */ #include -#include #include "bearing_area.h" -Eigen::ArrayXd nonuniform_bearing_area(Eigen::Ref x, Eigen::Ref h, - Eigen::Ref el_sort_by_max, Eigen::Ref heights) { +Eigen::ArrayXd nonuniform_bearing_area(Eigen::Ref x, Eigen::Ref h, + Eigen::Ref el_sort_by_max, Eigen::Ref heights) { if (x.size() != h.size()) { throw std::runtime_error("`x` and `h` must have the same size"); } @@ -73,8 +72,8 @@ Eigen::ArrayXd nonuniform_bearing_area(Eigen::Ref x, Eigen::Ref< } -Eigen::ArrayXd uniform1d_bearing_area(Eigen::Ref topography_h, bool periodic, - Eigen::Ref heights) { +Eigen::ArrayXd uniform1d_bearing_area(Eigen::Ref topography_h, bool periodic, + Eigen::Ref heights) { /* Bearing area values for each input height */ Eigen::ArrayXd fractional_bearing_areas(heights.size()); @@ -125,8 +124,8 @@ double _triangle(double h1_in, double h2_in, double h3_in, double h) { } -Eigen::ArrayXd uniform2d_bearing_area(Eigen::Ref topography_h, bool periodic, - Eigen::Ref heights) { +Eigen::ArrayXd uniform2d_bearing_area(Eigen::Ref topography_h, bool periodic, + Eigen::Ref heights) { /* Number of grid points for looping */ const auto nx{periodic ? topography_h.rows() : topography_h.rows()-1}; const auto ny{periodic ? topography_h.cols() : topography_h.cols()-1}; @@ -135,14 +134,16 @@ Eigen::ArrayXd uniform2d_bearing_area(Eigen::Ref topography_h, bool Eigen::ArrayXd fractional_bearing_areas(heights.size()); /* Compute bearing areas */ - for (int j{0}; j < heights.size(); j++) { + for (Eigen::Index j{0}; j < heights.size(); j++) { double bearing_area{0}; - int projected_area{0}; + /* Note: 64-bit counter; an `int` would overflow at 2^31 triangles */ + std::int64_t projected_area{0}; - /* This is assuming column-major storage */ - for (int x{0}; x < nx; x++) { + /* The inner loop runs over the last (fast) index of the row-major + storage */ + for (Eigen::Index x{0}; x < nx; x++) { const auto x1{x < topography_h.rows()-1 ? x+1 : 0}; - for (int y{0}; y < ny; y++) { + for (Eigen::Index y{0}; y < ny; y++) { const auto y1{y < topography_h.cols()-1 ? y+1 : 0}; const double h00{topography_h(x, y)}; const double h10{topography_h(x1, y)}; diff --git a/cpp/bearing_area.h b/cpp/bearing_area.h index dd3579a5..14fd6433 100644 --- a/cpp/bearing_area.h +++ b/cpp/bearing_area.h @@ -22,9 +22,9 @@ SOFTWARE. #include "eigen_helper.h" -Eigen::ArrayXd nonuniform_bearing_area(Eigen::Ref x, Eigen::Ref h, - Eigen::Ref el_sort_by_max, Eigen::Ref heights); -Eigen::ArrayXd uniform1d_bearing_area(Eigen::Ref topography_h, bool periodic, - Eigen::Ref heights); -Eigen::ArrayXd uniform2d_bearing_area(Eigen::Ref topography_h, bool periodic, - Eigen::Ref heights); +Eigen::ArrayXd nonuniform_bearing_area(Eigen::Ref x, Eigen::Ref h, + Eigen::Ref el_sort_by_max, Eigen::Ref heights); +Eigen::ArrayXd uniform1d_bearing_area(Eigen::Ref topography_h, bool periodic, + Eigen::Ref heights); +Eigen::ArrayXd uniform2d_bearing_area(Eigen::Ref topography_h, bool periodic, + Eigen::Ref heights); diff --git a/cpp/bicubic.cpp b/cpp/bicubic.cpp index cc4a215c..f434e79f 100644 --- a/cpp/bicubic.cpp +++ b/cpp/bicubic.cpp @@ -49,7 +49,7 @@ Bicubic::Bicubic(const Eigen::Ref& values, has_derivativey_{derivativey_opt.has_value()}, derivativex_(has_derivativex_ ? n1_ * n2_ : 0), derivativey_(has_derivativey_ ? n1_ * n2_ : 0), - coeff_(NPARA, n1_ * n2_) + cached_cell_{-1} { // Copy values to internal storage (row-major to linear) for (int i = 0; i < n1_; ++i) { @@ -137,15 +137,9 @@ Bicubic::Bicubic(const Eigen::Ref& values, this->A_ = this->A_.inverse(); /* - * Compute all spline coefficients and store them. + * Spline coefficients are computed on demand in + * get_spline_coefficients(); see the comment on cached_cell_. */ - for (int i1 = 0; i1 < this->n1_; i1++) { - for (int i2 = 0; i2 < this->n2_; i2++) { - this->coeff_.col(_row_major(i1, i2, this->n1_, this->n2_)) = - compute_spline_coefficients(i1, i2, this->values_, this->has_derivativex_, this->derivativex_, - this->has_derivativey_, this->derivativey_); - } - } } @@ -353,14 +347,27 @@ Bicubic::call(py::object py_x, py::object py_y, int derivative) double x = py_x.cast(); double y = py_y.cast(); double v, dx_out, dy_out; - eval(x, y, v, dx_out, dy_out); - return py::float_(v); + if (derivative == 0) { + eval(x, y, v); + return py::float_(v); + } else if (derivative == 1) { + eval(x, y, v, dx_out, dy_out); + return py::make_tuple(v, dx_out, dy_out); + } else { + double d2x, d2y, d2xy; + eval(x, y, v, dx_out, dy_out, d2x, d2y, d2xy); + return py::make_tuple(v, dx_out, dy_out, d2x, d2y, d2xy); + } } } - // Array inputs - convert to numpy arrays - py::array_t x_arr = py::array_t::ensure(py_x); - py::array_t y_arr = py::array_t::ensure(py_y); + // Array inputs - convert to numpy arrays. Note: `c_style` forces a + // C-contiguous copy for strided views; the evaluation loop below walks + // the buffers linearly and would silently read wrong (or out-of-bounds) + // elements from non-contiguous arrays. + using carray_t = py::array_t; + carray_t x_arr = carray_t::ensure(py_x); + carray_t y_arr = carray_t::ensure(py_y); if (!x_arr || !y_arr) { throw std::invalid_argument("Could not convert inputs to arrays."); diff --git a/cpp/bicubic.h b/cpp/bicubic.h index dfd11647..a225e5f3 100644 --- a/cpp/bicubic.h +++ b/cpp/bicubic.h @@ -95,15 +95,27 @@ class Bicubic { bool has_derivativex_, has_derivativey_; Eigen::ArrayXd derivativex_, derivativey_; - /* spline coefficients */ - Eigen::Array coeff_; + /* spline coefficients of the most recently used cell. Coefficients are + computed on demand: storing them for all cells requires 16 doubles + (128 bytes) per pixel, more than 2 GB for a 4096 x 4096 map. Query + points typically arrive ordered, so consecutive evaluations hit the + same cell and the single-cell cache amortizes the recomputation. */ + ptrdiff_t cached_cell_; + Eigen::Matrix cached_coeff_; /* lhs matrix */ Eigen::Matrix A_; - Eigen::Array + const Eigen::Matrix & get_spline_coefficients(int i1, int i2) { - return this->coeff_.col(_row_major(i1, i2, this->n1_, this->n2_)); + ptrdiff_t cell = _row_major(i1, i2, this->n1_, this->n2_); + if (cell != this->cached_cell_) { + this->cached_coeff_ = compute_spline_coefficients( + i1, i2, this->values_, this->has_derivativex_, this->derivativex_, + this->has_derivativey_, this->derivativey_); + this->cached_cell_ = cell; + } + return this->cached_coeff_; } Eigen::Matrix diff --git a/cpp/eigen_helper.h b/cpp/eigen_helper.h index 9a7b9cc8..0309e1ee 100644 --- a/cpp/eigen_helper.h +++ b/cpp/eigen_helper.h @@ -25,9 +25,15 @@ SOFTWARE. #ifndef __EIGEN_HELPER_H #define __EIGEN_HELPER_H +#include + #include -using ArrayXl = Eigen::Array; +/* Note: This must be a fixed-width 64-bit integer to match the int64 arrays + (e.g. from np.argsort) passed from the Python side; `long` is 32 bits on + Windows (LLP64), where a plain-long array type would make every binding + taking an ArrayXl reject its arguments. */ +using ArrayXl = Eigen::Array; using RowMajorXXd = Eigen::Array; using RowMajorXXi = Eigen::Array; using RowMajorXXb = Eigen::Array; diff --git a/cpp/module.cpp b/cpp/module.cpp index 5c366434..6be27d4f 100644 --- a/cpp/module.cpp +++ b/cpp/module.cpp @@ -43,29 +43,44 @@ PYBIND11_MODULE(_SurfaceTopography, mod) { mod.def("nonuniform_autocorrelation", &nonuniform_autocorrelation, "Height-difference autocorrelation of nonuniform line scans", py::arg("x"), py::arg("h"), py::arg("physical_size"), - py::arg("distances") = std::nullopt); + py::arg("distances") = std::nullopt, + py::call_guard()); + mod.def("nonuniform_height_height_autocorrelation", &nonuniform_height_height_autocorrelation, + "Height-height (product) autocorrelation of nonuniform line scans; " + "returns the unnormalized integral for each distance", + py::arg("x"), py::arg("h"), py::arg("distances"), + py::call_guard()); // === Patch/geometry analysis (from former C module) === mod.def("assign_patch_numbers", &assign_patch_numbers, "Assign unique numbers to connected patches", - py::arg("map"), py::arg("periodic"), py::arg("stencil") = std::nullopt); + py::arg("map"), py::arg("periodic"), py::arg("stencil") = std::nullopt, + py::call_guard()); mod.def("assign_segment_numbers", &assign_segment_numbers, - "Assign unique numbers to connected 1D segments", + "Assign unique numbers to connected 1D segments; the map is " + "always treated as periodic", py::arg("map")); mod.def("distance_map", &distance_map, - "Compute distance from each point to nearest marked point", - py::arg("map")); + "Compute distance from each point to nearest marked point; " + "distances are always measured with periodic boundary conditions", + py::arg("map"), py::call_guard()); mod.def("closest_patch_map", &closest_patch_map, - "Compute the tag of the closest patch for each point", - py::arg("map")); + "Compute the tag of the closest patch for each point; distances " + "are always measured with periodic boundary conditions", + py::arg("map"), py::call_guard()); mod.def("shortest_distance", &shortest_distance, - "Compute shortest distance between patches", - py::arg("fromc"), py::arg("fromp"), py::arg("to"), py::arg("maxd") = -1); + "Compute shortest distance between patches; distances are always " + "measured with periodic boundary conditions", + py::arg("fromc"), py::arg("fromp"), py::arg("to"), py::arg("maxd") = -1, + py::call_guard()); mod.def("correlation_function", &correlation_function, - "Compute real-space correlation function between two maps", - py::arg("map1"), py::arg("map2"), py::arg("max_dist")); + "Compute real-space correlation function between two maps; the " + "maps are always treated as periodic", + py::arg("map1"), py::arg("map2"), py::arg("max_dist"), + py::call_guard()); mod.def("perimeter_length", &perimeter_length, - "Compute total perimeter length of marked regions", + "Compute total perimeter length of marked regions; the map is " + "always treated as periodic", py::arg("map")); // === Bicubic interpolation (from former C module) === @@ -81,9 +96,15 @@ PYBIND11_MODULE(_SurfaceTopography, mod) { py::arg("x"), py::arg("y"), py::arg("derivative") = 0); // === Bearing area (from former C++ module) === - mod.def("nonuniform_bearing_area", &nonuniform_bearing_area, "Bearing area of a nonuniform line scan"); - mod.def("uniform1d_bearing_area", &uniform1d_bearing_area, "Bearing area of a uniform line scan"); - mod.def("uniform2d_bearing_area", &uniform2d_bearing_area, "Bearing area of a topography map"); + mod.def("nonuniform_bearing_area", &nonuniform_bearing_area, "Bearing area of a nonuniform line scan", + py::arg("x"), py::arg("h"), py::arg("el_sort_by_max"), py::arg("heights"), + py::call_guard()); + mod.def("uniform1d_bearing_area", &uniform1d_bearing_area, "Bearing area of a uniform line scan", + py::arg("h"), py::arg("periodic"), py::arg("heights"), + py::call_guard()); + mod.def("uniform2d_bearing_area", &uniform2d_bearing_area, "Bearing area of a topography map", + py::arg("h"), py::arg("periodic"), py::arg("heights"), + py::call_guard()); // === Moments (from former C++ module) === mod.def("nonuniform_mean", &nonuniform_moment<1>, "Mean of a nonuniform line scan", diff --git a/cpp/moments.h b/cpp/moments.h index 3d2a21a1..fb355240 100644 --- a/cpp/moments.h +++ b/cpp/moments.h @@ -32,11 +32,18 @@ template class _LineScanMoment { public: static double eval(double h1, double h2) { - if (std::abs(h2 - h1) < 1e-12) { - return 0; + /* (order+1) times the mean of h^order over a segment with linearly + interpolated end point heights h1 and h2. Analytically this is + (h2^(order+1) - h1^(order+1)) / (h2 - h1), which equals the + complete homogeneous symmetric polynomial of degree `order` in h1 + and h2. The polynomial form is exact for equal heights (where the + quotient is 0/0) and does not suffer from cancellation for nearly + equal heights. */ + double sum{0}; + for (int k{0}; k <= order; k++) { + sum += std::pow(h1, k) * std::pow(h2, order - k); } - // This is the generic expression, but it has numerical issues when h1 and h2 are close to each other - return (std::pow(h2, order+1) - std::pow(h1, order+1)) / (h2 - h1); + return sum; } }; @@ -67,7 +74,7 @@ class _LineScanMoment<3> { template -double nonuniform_moment(Eigen::Ref topography_x, Eigen::Ref topography_h, +double nonuniform_moment(Eigen::Ref topography_x, Eigen::Ref topography_h, double ref_h) { if (topography_x.size() != topography_h.size()) { throw std::runtime_error("`topography_x` and `topography_h` must have the same size"); @@ -80,7 +87,7 @@ double nonuniform_moment(Eigen::Ref topography_x, Eigen::Ref::eval(hi, hi1); @@ -91,14 +98,15 @@ double nonuniform_moment(Eigen::Ref topography_x, Eigen::Ref -double uniform1d_moment(Eigen::Ref topography_h, bool periodic, double ref_h) { +double uniform1d_moment(Eigen::Ref topography_h, bool periodic, double ref_h) { /* Accumulator for moment */ double moment{0}; - int physical_size{0}; + /* Note: 64-bit counter; an `int` would overflow for large grids */ + std::int64_t physical_size{0}; /* Compute moment */ const auto maxi{periodic ? topography_h.size() : topography_h.size()-1}; - for (int i{0}; i < maxi; i++) { + for (Eigen::Index i{0}; i < maxi; i++) { const auto i1{i < topography_h.size()-1 ? i+1 : 0}; const double hi{topography_h(i) - ref_h}, hi1{topography_h(i1) - ref_h}; /* Check for NaNs and only add if there are no NaNs */ @@ -114,81 +122,46 @@ double uniform1d_moment(Eigen::Ref topography_h, bool periodic, template class _TriangleMoment { -public: - static double eval(double h1_in, double h2_in, double h3_in) { - double h1{h1_in}, h2{h2_in}, h3{h3_in}; - - /* Sort h1, h2, h3 in ascending order */ - if (h1 > h2) std::swap(h1, h2); - if (h2 > h3) std::swap(h2, h3); - if (h1 > h2) std::swap(h1, h2); - - /* Compute moment */ - return ((std::pow(h2, order+2) - std::pow(h1, order+2)) / (h2 - h1) + - (std::pow(h3, order+2) - std::pow(h2, order+2)) / (h3 - h2)) / (order + 2) - - (h1*(std::pow(h2, order+1) - std::pow(h1, order+1)) / (h2 - h1) + - h3*(std::pow(h3, order+1) - std::pow(h2, order+1)) / (h3 - h2)) / (order + 1); - } -}; - -template <> -class _TriangleMoment<1> { public: static double eval(double h1, double h2, double h3) { - /* Sort h1, h2, h3 in ascending order */ - if (h1 > h2) std::swap(h1, h2); - if (h2 > h3) std::swap(h2, h3); - if (h1 > h2) std::swap(h1, h2); - - /* Compute moment */ - return (4*h2*h2 - h1*h1 - h3*h3 - h1*h2 - h2*h3) / 6; - } -}; - -template <> -class _TriangleMoment<2> { -public: - static double eval(double h1, double h2, double h3) { - /* Sort h1, h2, h3 in ascending order */ - if (h1 > h2) std::swap(h1, h2); - if (h2 > h3) std::swap(h2, h3); - if (h1 > h2) std::swap(h1, h2); - - /* Compute moment */ - return (6*h2*h2*h2 - h1*h1*h1 - h3*h3*h3 - h1*h1*h2 - h1*h2*h2 - h2*h2*h3 - h2*h3*h3) / 12; - } -}; - -template <> -class _TriangleMoment<3> { -public: - static double eval(double h1, double h2, double h3) { - /* Sort h1, h2, h3 in ascending order */ - if (h1 > h2) std::swap(h1, h2); - if (h2 > h3) std::swap(h2, h3); - if (h1 > h2) std::swap(h1, h2); - - /* Compute moment */ - return (8*h2*h2*h2*h2 - h1*h1*h1*h1 - h3*h3*h3*h3 - h1*h1*h1*h2 - h1*h1*h2*h2 - h1*h2*h2*h2 - h2*h2*h2*h3 - - h2*h2*h3*h3 - h2*h3*h3*h3) / 20; + /* Mean of h^order over a triangle whose corner heights are h1, h2 + and h3 (linear interpolation). Using the barycentric integral + formula + int_T l1^i l2^j l3^k dA = 2 A i! j! k! / (i+j+k+2)! + the mean evaluates to + = 2 / ((order+1) (order+2)) * H_order(h1, h2, h3) + where H_order is the complete homogeneous symmetric polynomial of + degree `order`. This expression is exact, requires no sorting of + the corner heights and remains valid for degenerate (equal + height) corners. */ + double sum{0}; + for (int i{0}; i <= order; i++) { + for (int j{0}; j <= order - i; j++) { + sum += std::pow(h1, i) * std::pow(h2, j) * std::pow(h3, order - i - j); + } + } + return 2 * sum / ((order + 1) * (order + 2)); } }; template -double uniform2d_moment(Eigen::Ref topography_h, bool periodic, double ref_h) { +double uniform2d_moment(Eigen::Ref topography_h, bool periodic, double ref_h) { /* Number of grid points for looping */ const auto nx{periodic ? topography_h.rows() : topography_h.rows()-1}; const auto ny{periodic ? topography_h.cols() : topography_h.cols()-1}; /* Accumulator for moment */ double moment{0}; - int projected_area{0}; + /* Note: 64-bit counter; an `int` would overflow at 2^31 triangles, + i.e. maps larger than 32768 x 32768 pixels */ + std::int64_t projected_area{0}; - /* Compute moment. Loop assumes column-major storage */ - for (int x{0}; x < nx; x++) { + /* Compute moment. The inner loop runs over the last (fast) index of the + row-major storage. */ + for (Eigen::Index x{0}; x < nx; x++) { const auto x1{x < topography_h.rows()-1 ? x+1 : 0}; - for (int y{0}; y < ny; y++) { + for (Eigen::Index y{0}; y < ny; y++) { const auto y1{y < topography_h.cols()-1 ? y+1 : 0}; const double h00{topography_h(x, y)}; const double h10{topography_h(x1, y)}; diff --git a/cpp/patchfinder.cpp b/cpp/patchfinder.cpp index fcad4554..1a54684e 100644 --- a/cpp/patchfinder.cpp +++ b/cpp/patchfinder.cpp @@ -49,10 +49,12 @@ static const std::vector> default_stencil = { void fill_patch(Eigen::Index nx, Eigen::Index ny, const RowMajorXXb &map, std::ptrdiff_t i0, std::ptrdiff_t j0, int p, bool periodic, - const std::vector> &stencil, RowMajorXXi &id) + const std::vector> &stencil, RowMajorXXi &id, + Stack &stack) { - Stack stack(DEFAULT_STACK_SIZE); - + /* The caller passes a shared stack that is reused across patches; + it is empty on entry and drained again before this function + returns, so no per-patch allocation is necessary. */ stack.push(i0, j0); id(i0, j0) = p; @@ -61,11 +63,13 @@ void fill_patch(Eigen::Index nx, Eigen::Index ny, const RowMajorXXb &map, stack.pop_bottom(i, j); for (const auto &[di, dj] : stencil) { - // Periodic boundary conditions + // Periodic boundary conditions. Note: user-provided stencils + // can contain offsets larger than the grid dimensions; a single + // `if` would wrap only one period and index out of bounds. std::ptrdiff_t jj = j + dj; if (periodic) { - if (jj < 0) jj += ny; - if (jj > ny - 1) jj -= ny; + while (jj < 0) jj += ny; + while (jj > ny - 1) jj -= ny; } else { if (jj < 0) continue; if (jj > ny - 1) continue; @@ -74,8 +78,8 @@ void fill_patch(Eigen::Index nx, Eigen::Index ny, const RowMajorXXb &map, // Periodic boundary conditions std::ptrdiff_t ii = i + di; if (periodic) { - if (ii < 0) ii += nx; - if (ii > nx - 1) ii -= nx; + while (ii < 0) ii += nx; + while (ii > nx - 1) ii -= nx; } else { if (ii < 0) continue; if (ii > nx - 1) continue; @@ -115,11 +119,15 @@ std::tuple assign_patch_numbers( RowMajorXXi id = RowMajorXXi::Zero(nx, ny); int p = 0; + /* Allocate the flood-fill stack once; allocating it inside fill_patch + would malloc/free the buffer for every individual patch. */ + Stack stack(DEFAULT_STACK_SIZE); + for (Eigen::Index i = 0; i < nx; ++i) { for (Eigen::Index j = 0; j < ny; ++j) { if (map(i, j) && id(i, j) == 0) { p++; - fill_patch(nx, ny, map, i, j, p, periodic, stencil, id); + fill_patch(nx, ny, map, i, j, p, periodic, stencil, id, stack); } } } @@ -177,7 +185,7 @@ std::tuple assign_segment_numbers(Eigen::Ref map) void track_distance(Eigen::Index nx, Eigen::Index ny, const RowMajorXXb &map, - RowMajorXXd &dist, RowMajorXXi &next) + RowMajorXXd &dist) { Stack stack(DEFAULT_STACK_SIZE); @@ -210,7 +218,6 @@ void track_distance(Eigen::Index nx, Eigen::Index ny, const RowMajorXXb &map, // Is i0, j0 closer than what is currently stored? if (d < dist(i, j)) { dist(i, j) = d; - next(i, j) = i0 * ny + j0; // Loop over all neighbors for (int joff = -1; joff <= 1; ++joff) { @@ -247,11 +254,8 @@ RowMajorXXd distance_map(Eigen::Ref map) // This stores the distance to the closest point on the contour RowMajorXXd dist = RowMajorXXd::Constant(nx, ny, nx * ny); - // This stores the index of the closest point - RowMajorXXi next = RowMajorXXi::Constant(nx, ny, nx * ny); - // Track distances from contact edge - track_distance(nx, ny, map, dist, next); + track_distance(nx, ny, map, dist); return dist; } diff --git a/cpp/stack.h b/cpp/stack.h index 784c1a22..d5d19a3e 100644 --- a/cpp/stack.h +++ b/cpp/stack.h @@ -39,6 +39,7 @@ SOFTWARE. #include #include +#include #define DEBUG_STACK_MAGIC_START std::size_t(0xDEADBEEF) #define DEBUG_STACK_MAGIC_END std::size_t(0xBEEFDEAD) @@ -52,11 +53,18 @@ class Stack { bp_ = 0; is_empty_ = true; data_ = malloc(buffer_size_); + if (!data_) { + throw std::bad_alloc(); + } } ~Stack() { free(data_); } + /* The class manages a raw buffer; copying would lead to a double free */ + Stack(const Stack &) = delete; + Stack &operator=(const Stack &) = delete; + bool is_empty() { return is_empty_; } @@ -270,14 +278,15 @@ class Stack { void *data_; /* Buffer containing the actual data. */ void expand(size_t new_size) { - printf("Expanding stack size to %3.2f MB.\n", - ((double) new_size)/(1024*1024)); void *new_data = malloc(new_size); #ifdef DEBUG_STACK_PRINT std::cout << "tp_ = " << tp_ << ", bp_ = " << bp_ << ", top_ = " << top_ << ", buffer_size_ = " << buffer_size_ << ", is_empty = " << is_empty_ << std::endl; #endif if (!new_data) { - printf("Failed to allocate new stack!\n"); + /* Note: propagate the failure; dereferencing the null pointer in the + memcpy below would crash the interpreter instead of raising + MemoryError */ + throw std::bad_alloc(); } if (tp_ > bp_) { assert(top_ == 0); diff --git a/meson.build b/meson.build index 69bded19..560392cb 100644 --- a/meson.build +++ b/meson.build @@ -3,6 +3,13 @@ project( 'SurfaceTopography', # Project name 'cpp', # Project type. We need a C++ compiler. default_options : ['cpp_std=c++17'], # Yes, we need C++17, at least for constexpr + # This must be 'python', not 'python3': builds run inside a virtual + # environment (pip/mesonpy build isolation), which provides 'python' on + # all platforms but no 'python3.exe' on Windows -- there, 'python3' + # resolves to some other interpreter on PATH that does not have the + # DiscoverVersion build requirement installed. When building with bare + # meson outside a venv (e.g. on Debian without python-is-python3), + # create a venv first. version: run_command('python', '-m', 'DiscoverVersion', check: true).stdout().strip(), ) diff --git a/pyproject.toml b/pyproject.toml index 4d585d28..e0d19d3a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,7 +58,7 @@ dependencies = [ "python-dateutil", "pyyaml", "requests", - "tiffile", + "tifffile", "xmltodict", "olefile", "pandas", @@ -68,10 +68,8 @@ dependencies = [ [project.optional-dependencies] test = [ - "flake8<8", "pytest", "pytest-cov", - "pytest-flake8", "runtests" ] diff --git a/test/IO/test_asc.py b/test/IO/test_asc.py index 4c43d692..44973f66 100644 --- a/test/IO/test_asc.py +++ b/test/IO/test_asc.py @@ -141,7 +141,9 @@ def test_wyko_matrix8(file_format_examples, filename="matrix-8.txt"): t = r.topography() assert t.unit == "nm" np.testing.assert_allclose(t.physical_sizes, (950400, 1267200)) - np.testing.assert_allclose(t.rms_height_from_area(), 74424.357775) + # Note: about half of this file's data points are undefined; sums and + # means are normalized by the number of *defined* points + np.testing.assert_allclose(t.rms_height_from_area(), 7804.797789741755) def test_single_column(file_format_examples, filename="single_column.txt"): diff --git a/test/IO/test_datx.py b/test/IO/test_datx.py index ca8ec56d..8a52665b 100644 --- a/test/IO/test_datx.py +++ b/test/IO/test_datx.py @@ -80,4 +80,4 @@ def test_datx2_metadata(file_format_examples): np.testing.assert_allclose(t.max(), 11.037771, rtol=1e-6) np.testing.assert_allclose(t.min(), -116.920823, rtol=1e-6) - np.testing.assert_allclose(t.rms_height_from_area(), 1.905869, rtol=1e-6) + np.testing.assert_allclose(t.rms_height_from_area(), 1.9089580434424194, rtol=1e-6) diff --git a/test/IO/test_frt.py b/test/IO/test_frt.py index 9527f09b..ea427166 100644 --- a/test/IO/test_frt.py +++ b/test/IO/test_frt.py @@ -72,12 +72,12 @@ def test_frt1_metadata(file_format_examples): assert t.unit == 'm' assert t.info['instrument']['vendor'] == 'FRT' - np.testing.assert_allclose(t.rms_height_from_area(), 2.047476e-05, rtol=1e-6) - np.testing.assert_allclose(t.rms_height_from_profile(), 1.23256e-05, rtol=1e-6) + np.testing.assert_allclose(t.rms_height_from_area(), 4.915785774480653e-06, rtol=1e-6) + np.testing.assert_allclose(t.rms_height_from_profile(), 1.951522762933369e-06, rtol=1e-6) t = t.detrend('curvature') - np.testing.assert_allclose(t.rms_height_from_area(), 3.463934e-06, rtol=1e-4) - np.testing.assert_allclose(t.rms_height_from_profile(), 1.248258e-06, rtol=1e-4) + np.testing.assert_allclose(t.rms_height_from_area(), 4.157523910393923e-06, rtol=1e-4) + np.testing.assert_allclose(t.rms_height_from_profile(), 1.0968828309450876e-06, rtol=1e-4) assert t.has_undefined_data @@ -103,11 +103,11 @@ def test_frt2_metadata(file_format_examples): assert t.unit == 'm' assert t.info['instrument']['vendor'] == 'FRT' - np.testing.assert_allclose(t.rms_height_from_area(), 1.853335e-05, rtol=1e-6) - np.testing.assert_allclose(t.rms_height_from_profile(), 1.439208e-05, rtol=1e-6) + np.testing.assert_allclose(t.rms_height_from_area(), 8.826763181113768e-06, rtol=1e-6) + np.testing.assert_allclose(t.rms_height_from_profile(), 7.663174319074244e-06, rtol=1e-6) t = t.detrend('curvature') - np.testing.assert_allclose(t.rms_height_from_area(), 7.405055e-06, rtol=1e-4) - np.testing.assert_allclose(t.rms_height_from_profile(), 7.332406e-06, rtol=1e-4) + np.testing.assert_allclose(t.rms_height_from_area(), 8.688677385181484e-06, rtol=1e-4) + np.testing.assert_allclose(t.rms_height_from_profile(), 7.627057775705951e-06, rtol=1e-4) assert t.has_undefined_data diff --git a/test/IO/test_gwy.py b/test/IO/test_gwy.py index e000c2fd..a7fb91ef 100644 --- a/test/IO/test_gwy.py +++ b/test/IO/test_gwy.py @@ -71,12 +71,14 @@ def test_gwyddion_metadata(file_format_examples): assert t.unit == "m" - np.testing.assert_allclose(t.rms_height_from_area(), 8.355506e-09, rtol=1e-6) - np.testing.assert_allclose(t.rms_height_from_profile(), 7.14371e-09, rtol=1e-6) + # Note: this file contains undefined data points; sums and means are + # normalized by the number of *defined* points + np.testing.assert_allclose(t.rms_height_from_area(), 8.356451162720974e-09, rtol=1e-6) + np.testing.assert_allclose(t.rms_height_from_profile(), 7.115125208610971e-09, rtol=1e-6) t = t.detrend("curvature") - np.testing.assert_allclose(t.rms_height_from_area(), 6.499288e-09, rtol=1e-6) - np.testing.assert_allclose(t.rms_height_from_profile(), 6.429779e-09, rtol=1e-6) + np.testing.assert_allclose(t.rms_height_from_area(), 6.500515512793928e-09, rtol=1e-6) + np.testing.assert_allclose(t.rms_height_from_profile(), 6.43099193553545e-09, rtol=1e-6) def test_gwyddion_undefined(file_format_examples): diff --git a/test/IO/test_ibw.py b/test/IO/test_ibw.py index a4445239..0a3625ae 100644 --- a/test/IO/test_ibw.py +++ b/test/IO/test_ibw.py @@ -87,7 +87,9 @@ def test_init(self): ['HeightRetrace', 'AmplitudeRetrace', 'PhaseRetrace', 'ZSensorRetrace']) self.assertEqual(reader._default_channel, 0) - self.assertEqual(reader.data['wave_header']['next'], 114425520) + # The reader must not pin the wave data in memory; it is loaded on + # demand in `topography()` + self.assertFalse(hasattr(reader, 'data')) def test_channels(self): reader = IBWReader(self.file_path) diff --git a/test/IO/test_io.py b/test/IO/test_io.py index 3f890205..4546e5cf 100644 --- a/test/IO/test_io.py +++ b/test/IO/test_io.py @@ -179,7 +179,9 @@ def _convert_filelist(filelist): ) explicit_physical_sizes = _convert_filelist( - ["matrix-5.txt", "mat-1.mat", "example-2d.npy"] + # Note: SRTM (.hgt) files do not carry physical sizes; the user needs + # to provide them + ["matrix-5.txt", "mat-1.mat", "example-2d.npy", "N46E013.hgt"] ) text_example_memory_list = [ diff --git a/test/IO/test_mi.py b/test/IO/test_mi.py index 4a918a37..b565aa5f 100644 --- a/test/IO/test_mi.py +++ b/test/IO/test_mi.py @@ -61,7 +61,11 @@ def test_read_header(): # Check if metadata has been read in correctly assert loader.channels[0].dim == 2 assert loader.channels[0].nb_grid_pts == (256, 256) - assert loader.channels[0].physical_sizes == (2e-05, 2e-05) + # `xLength`/`yLength` are stored in meters (2e-05 m) but the channel unit + # is µm, so the physical sizes must be reported as 20 µm + np.testing.assert_allclose(loader.channels[0].physical_sizes, (20.0, 20.0)) + # Non-length channels (V) keep the lateral sizes in meters + np.testing.assert_allclose(loader.channels[1].physical_sizes, (2e-05, 2e-05)) assert ( loader.channels[0].info["raw_metadata"]["DisplayOffset"] == "8.8577270507812517e-004" diff --git a/test/IO/test_plux.py b/test/IO/test_plux.py index 43676a38..bf2fc236 100644 --- a/test/IO/test_plux.py +++ b/test/IO/test_plux.py @@ -71,9 +71,9 @@ def test_plux_metadata(file_format_examples): assert t.info['instrument']['vendor'] == 'Sensofar' assert t.info['instrument']['name'] == 'S neox' - np.testing.assert_allclose(t.rms_height_from_area(), 2.15209, rtol=1e-6) - np.testing.assert_allclose(t.rms_height_from_profile(), 1.582635, rtol=1e-6) + np.testing.assert_allclose(t.rms_height_from_area(), 2.162500740959119, rtol=1e-6) + np.testing.assert_allclose(t.rms_height_from_profile(), 1.589988838291435, rtol=1e-6) t = t.detrend('curvature') - np.testing.assert_allclose(t.rms_height_from_area(), 1.435114, rtol=1e-4) - np.testing.assert_allclose(t.rms_height_from_profile(), 1.202886, rtol=1e-4) + np.testing.assert_allclose(t.rms_height_from_area(), 1.4420593851629364, rtol=1e-4) + np.testing.assert_allclose(t.rms_height_from_profile(), 1.2085860010079734, rtol=1e-4) diff --git a/test/IO/test_vk.py b/test/IO/test_vk.py index eb74979b..09c6f272 100644 --- a/test/IO/test_vk.py +++ b/test/IO/test_vk.py @@ -62,8 +62,9 @@ def test_vk3_metadata(file_format_examples): assert ny == 768 sx, sy = t.physical_sizes - np.testing.assert_allclose(sx, 704847000, rtol=1e-6) - np.testing.assert_allclose(sy, 528463000, rtol=1e-6) + # Physical size follows the pixel convention (nb_pixels * pixel_size) + np.testing.assert_allclose(sx, 705536000, rtol=1e-6) + np.testing.assert_allclose(sy, 529152000, rtol=1e-6) assert t.unit == 'pm' assert t.info['instrument']['vendor'] == 'Keyence' @@ -84,8 +85,8 @@ def test_vk4_metadata(file_format_examples): assert ny == 768 sx, sy = t.physical_sizes - np.testing.assert_allclose(sx, 1396330551, rtol=1e-6) - np.testing.assert_allclose(sy, 1046906679, rtol=1e-6) + np.testing.assert_allclose(sx, 1397695488, rtol=1e-6) + np.testing.assert_allclose(sy, 1048271616, rtol=1e-6) assert t.unit == 'pm' assert t.info['instrument']['vendor'] == 'Keyence' @@ -106,8 +107,8 @@ def test_vk6_metadata(file_format_examples): assert ny == 1536 sx, sy = t.physical_sizes - np.testing.assert_allclose(sx, 97169043, rtol=1e-6) - np.testing.assert_allclose(sy, 72864915, rtol=1e-6) + np.testing.assert_allclose(sx, 97216512, rtol=1e-6) + np.testing.assert_allclose(sy, 72912384, rtol=1e-6) assert t.unit == 'pm' assert t.info['instrument']['vendor'] == 'Keyence' diff --git a/test/IO/test_x3p.py b/test/IO/test_x3p.py index 9f7eaa5a..57d23111 100644 --- a/test/IO/test_x3p.py +++ b/test/IO/test_x3p.py @@ -48,8 +48,13 @@ def test_read(file_format_examples): assert surface.unit == 'm' assert surface.is_uniform assert surface.has_undefined_data - np.testing.assert_allclose(surface.rms_height_from_area(), 9.528212249587946e-05, rtol=1e-6) - np.testing.assert_allclose(surface.interpolate_undefined_data().rms_gradient(), 0.15300265543799388, rtol=1e-6) + # Note: this file contains 1.7% undefined data points. The reference + # value used to be 9.53e-05, which was an artifact of the mean height + # being normalized by the total instead of the defined point count (the + # heights sit at about -0.00572, so the resulting systematic offset + # dwarfed the true roughness of this surface). + np.testing.assert_allclose(surface.rms_height_from_area(), 2.3756839548083504e-07, rtol=1e-6) + np.testing.assert_allclose(surface.interpolate_undefined_data().rms_gradient(), 0.15300264662900961, rtol=1e-6) assert surface.info['instrument']['name'] == 'Mountains Map Technology Software (DIGITAL SURF, version 6.2)' assert surface.info['instrument']['vendor'] == 'DIGITAL SURF' diff --git a/test/IO/test_zag.py b/test/IO/test_zag.py index 4c633fae..8a62d1e1 100644 --- a/test/IO/test_zag.py +++ b/test/IO/test_zag.py @@ -23,23 +23,69 @@ # SOFTWARE. # +import io +import os +import struct +import zipfile + +import numpy as np import pytest -from SurfaceTopography.Container.IO import ZAGReader +from SurfaceTopography.Container.IO import ZAGReader, read_container +from SurfaceTopography.IO import ZONReader + + +@pytest.fixture +def synthetic_zag(file_format_examples, tmp_path): + """ + Build a minimal ZAG container (header + BMP thumbnail + ZIP archive) + that wraps the `zon-1.zon` example file, mirroring the layout parsed + by `ZAGReader`. + """ + with open(os.path.join(file_format_examples, "zon-1.zon"), "rb") as f: + zon = f.read() + + buf = io.BytesIO() + fake_bmp = b"BM" + b"\x00" * 62 # dummy thumbnail + buf.write(b"KPK0" + struct.pack("abc/item.xml" + "", + ) + z.writestr( + "abc/item.xml", + "data0" + "", + ) + z.writestr(f"abc/data0/{ZAGReader._ZON_UUID}", zon) + fn = tmp_path / "zag-1.zag" + fn.write_bytes(buf.getvalue()) + return str(fn) -@pytest.mark.skip -def test_zag(): - file_path = "/home/pastewka/Downloads/zag-1.zag" - with ZAGReader(file_path) as r: +def test_zag(synthetic_zag, file_format_examples): + with ZAGReader(synthetic_zag) as r: c = r.container(0) + assert len(c) == 1 + t = c[0] + assert t.dim == 2 - import matplotlib.pyplot as plt + # Heights must match a direct read of the wrapped ZON file + t_ref = ZONReader( + os.path.join(file_format_examples, "zon-1.zon") + ).topography() + np.testing.assert_allclose(c[0].heights(), t_ref.heights()) - for t in c: - plt.figure() - t.plot() - print(t.info) - plt.show() +def test_zag_read_container_outlives_reader(synthetic_zag, file_format_examples): + # `read_container` closes the reader (and its stream) before returning + # the lazy container; element access must still work afterwards + (c,) = read_container(synthetic_zag) + t = c[0] + t_ref = ZONReader( + os.path.join(file_format_examples, "zon-1.zon") + ).topography() + np.testing.assert_allclose(t.heights(), t_ref.heights()) diff --git a/test/file_format_examples/workflowtest.npy b/test/file_format_examples/workflowtest.npy new file mode 100644 index 00000000..9ff7bead Binary files /dev/null and b/test/file_format_examples/workflowtest.npy differ diff --git a/test/test_downsample.py b/test/test_downsample.py index 6aa89b8b..30bec29a 100644 --- a/test/test_downsample.py +++ b/test/test_downsample.py @@ -48,6 +48,9 @@ def test_downsample_nth(): t_down = t.downsample(2, mode="nth") assert t_down.nb_grid_pts == (5, 4) assert t_down.physical_sizes == (sx, sy) + # Pixel size must grow by the downsampling factor + np.testing.assert_allclose(t_down.pixel_size, (2 * sx / nx, 2 * sy / ny)) + np.testing.assert_allclose(t_down.area_per_pt, 4 * sx * sy / (nx * ny)) expected_heights = heights[::2, ::2] np.testing.assert_allclose(t_down.heights(), expected_heights) @@ -92,6 +95,17 @@ def test_downsample_non_integer_multiple(): assert t_down.nb_grid_pts == (5, 4) expected_heights = heights[:10, :8][::2, ::2] np.testing.assert_allclose(t_down.heights(), expected_heights) + # Pixel size is factor times the parent pixel size; since trailing + # points are dropped, the physical size shrinks accordingly + np.testing.assert_allclose(t_down.pixel_size, (2 * sx / nx, 2 * sy / ny)) + np.testing.assert_allclose( + t_down.physical_sizes, (5 * 2 * sx / nx, 4 * 2 * sy / ny) + ) + # Positions must coincide with the positions of the retained samples + x, y = t_down.positions() + xp, yp = t.positions() + np.testing.assert_allclose(x, xp[:10, :8][::2, ::2]) + np.testing.assert_allclose(y, yp[:10, :8][::2, ::2]) t_down = t.downsample(2, mode="average") assert t_down.nb_grid_pts == (5, 4) diff --git a/test/test_invariants.py b/test/test_invariants.py new file mode 100644 index 00000000..a3525030 --- /dev/null +++ b/test/test_invariants.py @@ -0,0 +1,172 @@ +# +# Copyright 2026 Lars Pastewka +# +# ### MIT license +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# + +""" +Property-based invariant tests: analysis results must be unchanged (or +transform in a known way) under transposition, translation, scaling and +shifts of the coordinate origin, and pipeline decorators must report +consistent geometry. +""" + +import numpy as np +import pytest + +from SurfaceTopography import NonuniformLineScan, Topography, UniformLineScan + + +@pytest.fixture +def anisotropic_topography(): + """Periodic topography with anisotropic pixels (px = 4 py).""" + rng = np.random.RandomState(42) + nx, ny = 64, 32 + heights = rng.uniform(-1, 1, (nx, ny)) + return Topography(heights, (8.0, 1.0), periodic=True) + + +def test_transpose_geometry(anisotropic_topography): + t = anisotropic_topography + tt = t.transpose() + nx, ny = t.nb_grid_pts + sx, sy = t.physical_sizes + px, py = t.pixel_size + assert tt.nb_grid_pts == (ny, nx) + assert tt.physical_sizes == (sy, sx) + np.testing.assert_allclose(tt.pixel_size, (py, px)) + # Geometry consistency of the decorator itself + np.testing.assert_allclose( + tt.pixel_size, np.asarray(tt.physical_sizes) / np.asarray(tt.nb_grid_pts)) + np.testing.assert_allclose(tt.heights(), t.heights().T) + + +def test_transpose_invariance_of_scalar_parameters(anisotropic_topography): + """Scalar roughness parameters that do not single out a scan direction + must be invariant under transposition, also for anisotropic pixels.""" + t = anisotropic_topography + tt = t.transpose() + np.testing.assert_allclose(tt.rms_height_from_area(), t.rms_height_from_area()) + np.testing.assert_allclose(tt.rms_gradient(), t.rms_gradient()) + np.testing.assert_allclose(tt.rms_laplacian(), t.rms_laplacian()) + np.testing.assert_allclose(tt.bandwidth(), t.bandwidth()) + + +def test_translate_invariance(anisotropic_topography): + """Integer-pixel translations of a periodic topography change nothing + but the position of the heights.""" + t = anisotropic_topography + tt = t.translate(offset=(5, 7)) + np.testing.assert_allclose( + tt.heights(), np.roll(np.roll(t.heights(), 5, axis=0), 7, axis=1)) + assert tt.physical_sizes == t.physical_sizes + assert tt.nb_grid_pts == t.nb_grid_pts + np.testing.assert_allclose(tt.rms_height_from_area(), t.rms_height_from_area()) + np.testing.assert_allclose(tt.rms_gradient(), t.rms_gradient()) + + +def test_scale_invariants(anisotropic_topography): + """Linear analysis functions scale linearly with a static height scale.""" + t = anisotropic_topography + factor = 2.5 + ts = t.scale(factor) + np.testing.assert_allclose(ts.rms_height_from_area(), + factor * t.rms_height_from_area()) + np.testing.assert_allclose(ts.rms_gradient(), factor * t.rms_gradient()) + # The bearing area transforms with the height scale + c = 0.3 + np.testing.assert_allclose(ts.bearing_area(factor * c), t.bearing_area(c)) + + +def test_decorator_geometry_consistency(anisotropic_topography): + """All decorators must report pixel_size == physical_sizes / nb_grid_pts.""" + t = anisotropic_topography + for decorated in [t.scale(2.0), + t.detrend('center'), + t.transpose(), + t.translate(offset=(1, 2)), + t.transpose().detrend('center'), + t.scale(3.0).transpose()]: + np.testing.assert_allclose( + decorated.pixel_size, + np.asarray(decorated.physical_sizes) / np.asarray(decorated.nb_grid_pts), + err_msg=f'inconsistent geometry for {decorated}') + + +def test_bearing_area_bounds_and_complement(): + """The bearing area is a CDF: it is 1 below the minimum, 0 above the + maximum, and the bearing area of the inverted surface complements it.""" + rng = np.random.RandomState(0) + heights = rng.uniform(-1, 1, (32, 33)) + for periodic in [False, True]: + t = Topography(heights, (2.0, 1.0), periodic=periodic) + ti = Topography(-heights, (2.0, 1.0), periodic=periodic) + assert t.bearing_area(heights.min() - 0.1) == pytest.approx(1.0) + assert t.bearing_area(heights.max() + 0.1) == pytest.approx(0.0) + for c in [-0.5, -0.1, 0.0, 0.2, 0.7]: + np.testing.assert_allclose(ti.bearing_area(-c), + 1 - t.bearing_area(c), atol=1e-12) + + +def test_uniform_line_scan_scale_and_transpose_roundtrip(): + rng = np.random.RandomState(1) + h = rng.uniform(-1, 1, 128) + t = UniformLineScan(h, 4.0, periodic=True) + np.testing.assert_allclose(t.scale(3.0).rms_height_from_profile(), + 3.0 * t.rms_height_from_profile()) + tt = t.translate(offset=17) + np.testing.assert_allclose(tt.heights(), np.roll(h, 17)) + np.testing.assert_allclose(tt.rms_height_from_profile(), + t.rms_height_from_profile()) + + +def test_nonuniform_detrend_origin_invariance(): + """Detrending must not depend on where the scan sits on the x-axis. + + This is a regression test for the ill-conditioned normal equations in + the nonuniform `polyfit`, which are now solved in centered coordinates. + """ + rng = np.random.RandomState(2) + x = np.sort(rng.uniform(0, 1, 100)) + x[0] = 0 + h = 0.1 * rng.randn(100) + 0.5 * x - 0.8 * x * x + t = NonuniformLineScan(x, h) + for offset in [-1e3, 1e3, 1e5]: + t_shifted = NonuniformLineScan(x + offset, h) + # Evaluating the detrending polynomial a0 + a1 x + a2 x^2 at + # x ~ offset is limited by double precision to an absolute error + # of order eps * offset^2, even for exact coefficients + atol = max(1e-9, 100 * np.finfo(float).eps * offset ** 2) + for mode in ['mean', 'median', 'rms-tilt', 'slope', 'rms-curvature']: + np.testing.assert_allclose( + t_shifted.detrend(mode).heights(), + t.detrend(mode).heights(), atol=atol, + err_msg=f'detrend mode {mode} is not origin invariant ' + f'for offset {offset}') + + +def test_nonuniform_polyfit_recovers_exact_polynomial(): + rng = np.random.RandomState(3) + x = np.sort(rng.uniform(0, 1, 50)) + x[0] = 0 + h = 2.0 + 3.0 * x - 4.0 * x * x + t = NonuniformLineScan(x, h) + np.testing.assert_allclose(t.polyfit(2), [2.0, 3.0, -4.0], atol=1e-9) diff --git a/test/test_moments.py b/test/test_moments.py index 03e78f37..b384a56e 100644 --- a/test/test_moments.py +++ b/test/test_moments.py @@ -137,33 +137,41 @@ def test_with_nan_values(self): # ============================================================================= # Tests for uniform2d_mean # ============================================================================= -# NOTE: uniform2d_mean uses triangle-based integration which computes moments -# differently than simple averaging. For a constant surface at height h with -# ref_h=0, the result is 0 (not h). These functions are not currently used -# in the codebase but are tested for coverage. +# uniform2d_mean uses triangle-based integration: each pixel is split into +# two triangles over which the heights are interpolated linearly. class TestUniform2dMean: """Tests for uniform2d_mean function.""" - def test_constant_topography_returns_zero(self): - """Constant topography with ref_h=0 returns 0 due to triangle integration.""" + def test_constant_topography(self): + """Mean of a constant topography is that constant.""" h = np.full((4, 4), 5.0) - # Triangle-based moment returns 0 for constant surface mean_periodic = cpp.uniform2d_mean(h, True) mean_nonperiodic = cpp.uniform2d_mean(h, False) - assert_allclose(mean_periodic, 0.0, atol=1e-10) - assert_allclose(mean_nonperiodic, 0.0, atol=1e-10) + assert_allclose(mean_periodic, 5.0, rtol=1e-10) + assert_allclose(mean_nonperiodic, 5.0, rtol=1e-10) - def test_varying_topography_nonzero(self): - """Varying topography should return non-zero moment.""" + def test_linear_ramp(self): + """Mean of a linear (nonperiodic) ramp is the midpoint height.""" h = np.array([[0.0, 0.0, 0.0], [1.0, 1.0, 1.0], [2.0, 2.0, 2.0]]) mean = cpp.uniform2d_mean(h, False) - # Should return some non-zero value for varying data - assert mean != 0.0 + assert_allclose(mean, 1.0, rtol=1e-10) + + def test_vs_numpy(self): + """ + For periodic topographies every grid point is shared by exactly six + triangles, so the triangle-based mean equals np.mean exactly. + """ + np.random.seed(42) + h = np.random.randn(64, 64) + + mean = cpp.uniform2d_mean(h, True) + + assert_allclose(mean, np.mean(h), atol=1e-10) # ============================================================================= @@ -236,6 +244,22 @@ def test_constant_zero_variance(self): assert_allclose(var, 0.0, atol=1e-10) + def test_terraced_variance(self): + """ + Variance of a terraced +1/-1 surface (rows alternating) about its + mean. Within the linear interpolation, each pixel is a ramp from -1 + to +1 (or constant), analogous to the 1D alternating case. + """ + h = np.empty((4, 4)) + h[0::2, :] = 1.0 + h[1::2, :] = -1.0 + + var = cpp.uniform2d_variance(h, True, ref_h=0.0) + + # Every triangle has corner heights (+1, -1, +1) up to sign; + # = 2/12 * H_2(1, -1, 1) = 2/12 * (1+1+1-1+1-1) = 1/3 + assert_allclose(var, 1.0 / 3.0, rtol=1e-10) + # ============================================================================= # Tests for 3rd moment functions @@ -303,6 +327,41 @@ def test_uniform1d_moment4_positive(self): m4 = cpp.uniform1d_moment4(h, True, ref_h=0.0) + # Each segment ramps from -1 to +1; = H_4(-1, 1)/5 = 1/5 + assert_allclose(m4, 1.0 / 5.0, rtol=1e-10) + + def test_uniform1d_moment4_flat_segments(self): + """ + 4th moment of a constant profile about a different reference; flat + segments must contribute their exact moment (the analytic quotient + is a removable 0/0 singularity here). + """ + h = np.full(5, 2.0) + + m4 = cpp.uniform1d_moment4(h, False, ref_h=0.0) + + assert_allclose(m4, 16.0, rtol=1e-10) + + def test_nonuniform_moment4_flat_segments(self): + """Same as above for the nonuniform variant.""" + x = np.array([0.0, 1.0, 3.0]) + h = np.array([5.0, 5.0, 5.0]) + + m4 = cpp.nonuniform_moment4(x, h, ref_h=0.0) + + assert_allclose(m4, 625.0, rtol=1e-10) + + def test_uniform2d_moment4_quantized_no_nan(self): + """ + Quantized (rounded) height data produces many triangles with equal + corner heights; the moment must remain finite (no 0/0). + """ + np.random.seed(42) + h = np.round(np.random.randn(16, 16)) + + m4 = cpp.uniform2d_moment4(h, True, ref_h=0.0) + + assert np.isfinite(m4) assert m4 > 0 def test_uniform2d_moment4_varying_positive(self): @@ -351,14 +410,13 @@ def test_uniform1d_variance_vs_numpy(self): # Should be close for smooth data assert_allclose(var_cpp, var_numpy, rtol=0.1) - def test_uniform2d_mean_constant_is_zero(self): - """uniform2d_mean of constant returns 0 due to triangle integration.""" + def test_uniform2d_mean_constant(self): + """uniform2d_mean of a constant surface is that constant.""" h = np.full((10, 10), 7.5) mean_cpp = cpp.uniform2d_mean(h, True) - # Triangle-based integration returns 0 for constant surfaces - assert_allclose(mean_cpp, 0.0, atol=1e-10) + assert_allclose(mean_cpp, 7.5, rtol=1e-10) # ============================================================================= @@ -393,9 +451,9 @@ def test_uniform2d_2x2(self): mean = cpp.uniform2d_mean(h, False) - # Triangle-based integration returns non-zero for varying data - # Note: the exact value depends on triangle integration details - assert np.isfinite(mean) + # Two triangles with corner heights (0, 1, 1) and (1, 2, 1): + # means 2/3 and 4/3, average 1 + assert_allclose(mean, 1.0, rtol=1e-10) def test_all_nan_uniform1d(self): """All NaN values should result in NaN or 0.""" diff --git a/test/test_nonuniform_line_scan.py b/test/test_nonuniform_line_scan.py index f78b3c6c..2d45978c 100644 --- a/test/test_nonuniform_line_scan.py +++ b/test/test_nonuniform_line_scan.py @@ -165,6 +165,43 @@ def test_detrend_slope(file_format_examples): assert not t.detrend("slope").is_periodic +def test_detrend_slope_coefficients(): + # The rms-slope-minimizing constant slope is the *length-weighted* mean + # of the derivative, i.e. (h[-1] - h[0]) / L, and the offset must remove + # the mean of the tilt-corrected profile + t = NonuniformLineScan(x=[0.0, 1.0, 10.0], y=[0.0, 1.0, 1.0]) + d = t.detrend("slope") + a0, a1 = d.coeffs + np.testing.assert_allclose(a1, 0.1) + np.testing.assert_allclose(d.mean(), 0, atol=1e-15) + + +def test_detrend_slope_origin_invariance(): + # Detrending must not depend on where the x-axis starts + x = np.array([0.0, 1.3, 2.4, 3.5, 4.8, 6.0]) + h = np.array([1.0, -0.5, 2.0, 0.7, -1.2, 0.3]) + for mode in ["center", "height", "slope", "curvature"]: + r0 = NonuniformLineScan(x, h).detrend(mode).rms_height_from_profile() + r5 = NonuniformLineScan(x + 5.0, h).detrend(mode).rms_height_from_profile() + np.testing.assert_allclose(r0, r5, err_msg=f"mode {mode}") + + +def test_to_uniform_offset_origin(): + # The uniform grid starts at zero; the parent scan may start anywhere. + # The interpolation must be shifted by the parent's origin, otherwise + # np.interp clamps and fills the grid with the boundary values. + x = np.array([5.0, 6.0, 7.0, 8.0]) + h = np.array([1.0, 2.0, 3.0, 4.0]) + u = NonuniformLineScan(x, h).to_uniform(nb_points=4, padding=0) + np.testing.assert_allclose(u.heights(), h) + # PSD of the interpolated scan must not depend on the x origin + t0 = NonuniformLineScan(x - 5.0, h) + t5 = NonuniformLineScan(x, h) + q0, C0 = t0.power_spectrum_from_profile(reliable=False) + q5, C5 = t5.power_spectrum_from_profile(reliable=False) + np.testing.assert_allclose(C0, C5) + + def test_detrend_curvature(file_format_examples): t = XYZReader(os.path.join(file_format_examples, "xy-1.txt")).topography() assert not t.detrend("curvature").is_periodic diff --git a/test/test_pipeline.py b/test/test_pipeline.py index a98c52f3..469780e3 100644 --- a/test/test_pipeline.py +++ b/test/test_pipeline.py @@ -297,6 +297,17 @@ def test_transposed_topography(): assert sx == sy2 assert sy == sx2 assert (surf.heights() == surf2.heights().T).all() + # Pixel size and subdomain information must be transposed as well + px, py = surf.pixel_size + px2, py2 = surf2.pixel_size + assert px == py2 + assert py == px2 + assert surf2.nb_subdomain_grid_pts == (ny, nx) + assert surf2.subdomain_locations == tuple(reversed(surf.subdomain_locations)) + np.testing.assert_almost_equal(surf2.area_per_pt, surf.area_per_pt) + # Scalar analysis results must be invariant under transposition, also for + # anisotropic pixels + np.testing.assert_almost_equal(surf2.rms_gradient(), surf.rms_gradient()) def test_undefined_data_and_squeeze(): diff --git a/test/test_scalar_parameters.py b/test/test_scalar_parameters.py index 3575b145..04f60525 100644 --- a/test/test_scalar_parameters.py +++ b/test/test_scalar_parameters.py @@ -306,8 +306,9 @@ def test_rms_slope_from_area(): def test_rms_height_with_undefined_data(file_format_examples): t = read_topography(os.path.join(file_format_examples, "opd-3.opd")) assert t.has_undefined_data - np.testing.assert_allclose(t.rms_height_from_profile(), 0.011449207840819613) - np.testing.assert_allclose(t.rms_height_from_area(), 0.011782116700392838) + # Sums and means are normalized by the number of *defined* data points + np.testing.assert_allclose(t.rms_height_from_profile(), 0.011449215992795511) + np.testing.assert_allclose(t.rms_height_from_area(), 0.011782127367582199) @pytest.mark.skipif( @@ -682,8 +683,21 @@ def func(dx, dy=None): # Makes sure there is a significant difference by removing tip artefacts, so we are doing a meaningful test assert hrms_r[2] / hrms_tip_artefacts_removed[2] > 4 - # now we test that we indeed removed the tip artefacts when integrating the PSD - assert abs(hrms_f_reliable[2] / hrms_tip_artefacts_removed[2] - 1) < 0.2 + # The tip-artefact-cleaned scale-dependent curvature must match the + # scale-dependent curvature of the original (unartefacted) topography at + # the same distances + r_clean, k_clean = t_artefacted.scale_dependent_curvature_from_profile( + reliable=True, resampling_method=None + ) + r_ref, k_ref = t.scale_dependent_curvature_from_profile(resampling_method=None) + np.testing.assert_allclose( + np.max(k_clean), np.max(np.interp(r_clean, r_ref, k_ref)), rtol=0.2 + ) + # PSD integration weights small scales differently than the + # second-derivative stencil (whose transfer function rolls off towards + # the stencil's Nyquist frequency), so these two measures of the rms + # curvature over the reliable band agree only roughly + assert abs(hrms_f_reliable[2] / hrms_tip_artefacts_removed[2] - 1) < 0.5 def test_integrate_psd_from_profile_remove_tip_artefacts_areal_scan(): @@ -916,8 +930,12 @@ def func(dx, dy=None): # Makes sure there is a significant difference by removing tip artefacts, so we are doing a meaningful test assert hrms_r[2] / hrms_tip_artefacts_removed[2] > 4 # now we test that we indeed removed the tip artefacts when integrating the PSD, - # i.e. that reliable PSD integration is approx equivalent to reliable SDRPs - assert abs(hrms_f_reliable[2] / hrms_tip_artefacts_removed[2] - 1) < 0.2 + # i.e. that reliable PSD integration is approx equivalent to reliable SDRPs. + # PSD integration weights small scales differently than the + # second-derivative stencil (whose transfer function rolls off towards + # the stencil's Nyquist frequency), so these two measures agree only + # roughly + assert abs(hrms_f_reliable[2] / hrms_tip_artefacts_removed[2] - 1) < 0.5 # Assert the container gives the same results: c_artefacted = InMemorySurfaceContainer( diff --git a/test/test_scan_line_alignment.py b/test/test_scan_line_alignment.py index 650f6c93..e542fd08 100644 --- a/test/test_scan_line_alignment.py +++ b/test/test_scan_line_alignment.py @@ -44,7 +44,7 @@ def test_scan_line_align_removes_offsets(): h[i, :] += offsets[i] topo = Topography(h, (1, 1), unit='um') - aligned = topo.scan_line_align() + aligned = topo.scan_line_align(direction='y') # Check that aligned heights have reduced line-to-line variation aligned_heights = aligned.heights() @@ -65,7 +65,7 @@ def test_scan_line_align_removes_tilt(): h[i, :] = slope * y topo = Topography(h, (1, 1), unit='um') - aligned = topo.scan_line_align() + aligned = topo.scan_line_align(direction='y') # Check that slopes within lines are reduced aligned_heights = aligned.heights() @@ -75,17 +75,22 @@ def test_scan_line_align_removes_tilt(): assert abs(slope) < 0.001 -def test_scan_line_align_y_direction(): - """Test scan line alignment in y-direction.""" +def test_scan_line_align_x_direction(): + """Test scan line alignment in the (default) x-direction. + + The first array index is the x index, so scan lines running along the + x direction are *columns* of the height array; each column has constant + y and is acquired as one fast scan. + """ nx, ny = 32, 32 h = np.zeros((nx, ny)) - # Add offsets to columns instead of rows + # Add a per-scan-line offset; each scan line along x is a column np.random.seed(42) for j in range(ny): h[:, j] += np.random.randn() * 5 topo = Topography(h, (1, 1), unit='um') - aligned = topo.scan_line_align(direction='y') + aligned = topo.scan_line_align(direction='x') aligned_heights = aligned.heights() col_means = [aligned_heights[:, j].mean() for j in range(ny)] @@ -106,7 +111,7 @@ def test_scan_line_align_preserves_features(): h[i, :] += np.random.randn() * 0.1 * np.arange(ny) / ny # Tilt topo = Topography(h, (1, 1), unit='um') - aligned = topo.scan_line_align() + aligned = topo.scan_line_align(direction='y') # The bump should still be visible after alignment # (center should be higher than corners) @@ -132,7 +137,7 @@ def test_scan_line_align_masked_data(): h = np.ma.array(h, mask=mask) topo = Topography(h, (1, 1), unit='um') - aligned = topo.scan_line_align() + aligned = topo.scan_line_align(direction='y') # Should complete without error aligned_heights = aligned.heights() @@ -153,7 +158,7 @@ def test_scan_line_align_pipeline(): topo = Topography(h, (1, 1), unit='um') # Should work in pipeline: first align scan lines, then global detrend - result = topo.scan_line_align().detrend() + result = topo.scan_line_align(direction='y').detrend() assert result.heights().shape == (32, 32) # The pipeline should be preserved @@ -169,7 +174,7 @@ def test_scan_line_align_mode_mean(): h[i, :] += np.random.randn() * 5 topo = Topography(h, (1, 1), unit='um') - aligned = topo.scan_line_align(mode='mean') + aligned = topo.scan_line_align(direction='y', mode='mean') aligned_heights = aligned.heights() line_means = [aligned_heights[i, :].mean() for i in range(nx)] @@ -182,7 +187,7 @@ def test_scan_line_align_periodicity(): topo = Topography(h, (1, 1), unit='um', periodic=True) assert topo.is_periodic - aligned = topo.scan_line_align() + aligned = topo.scan_line_align(direction='y') assert not aligned.is_periodic @@ -195,7 +200,7 @@ def test_scan_line_align_coefficients(): h[i, :] = 2 + 3 * np.arange(ny) / ny # Same tilt for all lines topo = Topography(h, (1, 1), unit='um') - aligned = topo.scan_line_align() + aligned = topo.scan_line_align(direction='y') # Access coefficients coeffs = aligned.line_coeffs @@ -229,7 +234,7 @@ def test_scan_line_align_pickling(): h[i, :] += i * 0.5 topo = Topography(h, (1, 1), unit='um') - aligned = topo.scan_line_align() + aligned = topo.scan_line_align(direction='y') # Trigger computation original_heights = aligned.heights() @@ -259,7 +264,7 @@ def test_scan_line_align_degree_0(): h[i, :] += 5 * t topo = Topography(h, (1, 1), unit='um') - aligned = topo.scan_line_align(degree=0) + aligned = topo.scan_line_align(direction='y', degree=0) # Offsets should be removed aligned_heights = aligned.heights() @@ -293,7 +298,7 @@ def test_scan_line_align_degree_2(): topo = Topography(h, (1, 1), unit='um') # With degree=2, quadratic term should be removed - aligned_deg2 = topo.scan_line_align(degree=2) + aligned_deg2 = topo.scan_line_align(direction='y', degree=2) aligned_heights_2 = aligned_deg2.heights() # Check residual curvature is small curvatures_2 = [np.polyfit(t, aligned_heights_2[i, :], 2)[0] @@ -327,7 +332,7 @@ def test_scan_line_align_degree_3(): h[i, :] = a * t**3 + b * t**2 + c * t + d topo = Topography(h, (1, 1), unit='um') - aligned = topo.scan_line_align(degree=3) + aligned = topo.scan_line_align(direction='y', degree=3) # Check coefficients shape assert aligned.line_coeffs.shape == (nx, 4) @@ -362,7 +367,7 @@ def test_scan_line_align_pickling_with_degree(): h[i, :] = np.random.randn() * t**2 + np.random.randn() * t + i * 0.5 topo = Topography(h, (1, 1), unit='um') - aligned = topo.scan_line_align(degree=2) + aligned = topo.scan_line_align(direction='y', degree=2) # Trigger computation original_heights = aligned.heights() @@ -399,7 +404,7 @@ def test_scan_line_align_scanner_bow_correction(): topo = Topography(h, (10, 20), unit='um') # Correct with degree=2 - aligned = topo.scan_line_align(degree=2) + aligned = topo.scan_line_align(direction='y', degree=2) aligned_heights = aligned.heights() # The parabolic bow should be removed