diff --git a/CHANGELOG.md b/CHANGELOG.md index 0c47866..bfed637 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,30 @@ Observes [Semantic Versioning](https://semver.org/spec/v2.0.0.html) standard and [Keep a Changelog](https://keepachangelog.com/en/1.0.0/) convention. +## [0.8.4] - 2026-08-21 + ++ Fix - `prairie_view_loader.py` route multipage frames by the acquisition XML. Page offsets + were computed from a global channel-by-plane stride and applied to a file list + `get_prairieview_filenames` had already filtered by channel and plane, so the offset was + applied twice. With two channels roughly half the output frames were never assigned, and + because the output array is allocated with `np.empty` they held uninitialised memory rather + than raising. `get_prairieview_file_pages()` now returns the filename and page the XML names + for each frame, and the multipage branch writes exactly those. Note the previous formula was + correct when the file list was not filtered — a single interleaved file holding every channel + and plane — which is why the defect went unnoticed. ++ Fix - `prairie_view_loader.py` carry z positions forward across frames. `PVStateShard` + records state *changes*, so a frame whose z has not moved omits `positionCurrent` entirely + and inherits it; counting the elements in a cycle counted re-declarations rather than planes, + and a 3-plane bidirectional-Z recording failed outright with "Number of z fields does not + match number of depths". Positions are accumulated from the document-level shard through + every preceding frame, and a depth that is declared nowhere now raises instead of entering + `fieldZ` as `None`, which the length check could not detect. ++ Fix - `prairie_view_loader.py` handle an acquisition that stopped part-way through its final + cycle. `num_frames` is floored by `num_planes`, so the early planes are named more frames + than the late ones; the incomplete trailing cycle is now dropped with a warning, keeping + every plane the same length. Fewer frames than expected still raises, since that indicates + missing data rather than a partial cycle. + ## [0.8.3] - 2026-07-31 + Fix - `prairie_view_loader.py` correct plane-to-file mapping for bidirectional Z scans. diff --git a/element_interface/prairie_view_loader.py b/element_interface/prairie_view_loader.py index 1a6f22e..e444590 100644 --- a/element_interface/prairie_view_loader.py +++ b/element_interface/prairie_view_loader.py @@ -113,6 +113,46 @@ def get_prairieview_filenames( fnames = np.unique([f.attrib["filename"] for f in frames]).tolist() return fnames if not return_pln_chn else (fnames, plane_idx, channel) + def get_prairieview_file_pages(self, plane_idx=None, channel=None): + """Ordered `(filename, page)` pairs for one plane and channel. + + The acquisition XML states, per frame, which file and which page inside it + holds that frame. Reading those pairs directly is the only safe way to + assemble a multipage series: computing page offsets from a global + channel-by-plane stride assumes every channel and plane shares one + interleaved page sequence, which is false when PrairieView writes each + channel to its own file. Applying such a stride to an already-filtered + file list applies the offset twice. + + Pages are 1-based in the XML and returned as declared. Order follows the + XML, which is the recording order. + + Returns: + (pairs, plane_idx, channel) where pairs is [(filename, page), ...] + """ + _, plane_idx, channel = self.get_prairieview_filenames( + plane_idx=plane_idx, channel=channel, return_pln_chn=True + ) + channel_search = f"/[@channel='{channel}']" + bidi_map = self.meta.get("_bidi_z_index_map") + multiplane = self.meta["num_planes"] > 1 + + pairs = [] + for sequence in self._xml_root.findall(".//Sequence[@cycle]"): + if multiplane and bidi_map: + cycle_num = int(sequence.attrib.get("cycle")) + target_idx = plane_idx if cycle_num % 2 == 1 else bidi_map[plane_idx] + frames = sequence.findall(f"Frame[@index='{target_idx}']") + elif multiplane: + frames = sequence.findall(f"Frame[@index='{plane_idx}']") + else: + frames = sequence.findall("Frame") + for frame in frames: + for f in frame.findall(f"File{channel_search}"): + page = f.attrib.get("page") + pairs.append((f.attrib["filename"], int(page) if page else 1)) + return pairs, plane_idx, channel + def write_single_bigtiff( self, plane_idx=None, @@ -155,20 +195,38 @@ def write_single_bigtiff( logger.warning( "Ignoring `gb_per_file` argument for multi-page tiff (NotYetImplemented)" ) - # For multi-page tiff - the pages are organized as: - # (channel x slice x frame) - each page is (height x width) - # - TODO: verify this is the case for Bruker multi-page tiff - # This implementation is partially based on the reference code from `scanreader` package - https://github.com/atlab/scanreader - # See: https://github.com/atlab/scanreader/blob/2a021a85fca011c17e553d0e1c776998d3f2b2d8/scanreader/scans.py#L337 - slice_step = self.meta["num_channels"] - frame_step = self.meta["num_channels"] * self.meta["num_planes"] - slice_idx = self.meta["plane_indices"].index(plane_idx) - channel_idx = self.meta["channels"].index(channel) - - page_indices = [ - frame_idx * frame_step + slice_idx * slice_step + channel_idx - for frame_idx in range(self.meta["num_frames"]) - ] + # Each output frame is taken from the file and page the XML names for + # it. An earlier implementation instead computed page offsets from a + # global channel-by-plane stride, which double-counts the offset when + # the file list has already been filtered by channel and plane, and + # silently leaves frames unwritten. + file_pages, plane_idx, channel = self.get_prairieview_file_pages( + plane_idx=plane_idx, channel=channel + ) + # `num_frames` is floored by `num_planes`, so an acquisition that + # stopped part-way through its last cycle names more frames for the + # early planes than for the late ones. Every plane must come out the + # same length, so the incomplete cycle is dropped rather than making + # one plane longer than another. + if len(file_pages) > self.meta["num_frames"]: + logger.warning( + "The XML names %d frames for plane %s channel %s but only %d " + "complete cycles were acquired; dropping the %d trailing " + "frame(s) of the incomplete final cycle.", + len(file_pages), + plane_idx, + channel, + self.meta["num_frames"], + len(file_pages) - self.meta["num_frames"], + ) + file_pages = file_pages[: self.meta["num_frames"]] + elif len(file_pages) < self.meta["num_frames"]: + raise ValueError( + f"The XML names {len(file_pages)} frames for plane {plane_idx} " + f"channel {channel}, fewer than the {self.meta['num_frames']} " + f"the metadata reports. Frames are missing, so the movie is " + f"not being written." + ) combined_data = np.empty( [ @@ -178,30 +236,42 @@ def write_single_bigtiff( ], dtype=np.uint16, # use unsigned int 16 instead of int. int is defined as 32 or 64 bit based on the platform -> this will inflated a 16 bit tiff by 2 to 4 times! ) - start_page = 0 - try: - for input_file in tiff_names: - with tifffile.TiffFile((self.prairieview_dir / input_file).as_posix()) as tffl: - # Get indices in this tiff file and in output array - final_page_in_file = start_page + len(tffl.pages) - is_page_in_file = lambda page: page in range( - start_page, final_page_in_file + written = np.zeros(self.meta["num_frames"], dtype=bool) + + # Group the output positions by source file so each file opens once. + by_file = {} + for out_idx, (fname, page) in enumerate(file_pages): + by_file.setdefault(fname, []).append((out_idx, page)) + + for input_file, entries in by_file.items(): + try: + with tifffile.TiffFile( + (self.prairieview_dir / input_file).as_posix() + ) as tffl: + n_pages = len(tffl.pages) + # XML pages are 1-based; a page outside the file means the + # XML and the file on disk disagree. + for out_idx, page in entries: + if not 1 <= page <= n_pages: + raise ValueError( + f"XML names page {page} of {input_file}, which " + f"holds {n_pages} pages." + ) + keys = [page - 1 for _, page in entries] + combined_data[[out_idx for out_idx, _ in entries]] = ( + tffl.asarray(key=keys) ) - pages_in_file = filter(is_page_in_file, page_indices) - file_indices = [page - start_page for page in pages_in_file] - global_indices = [ - is_page_in_file(page) for page in page_indices - ] - - # Read from this tiff file (if needed) - if len(file_indices) > 0: - # this line looks a bit ugly but is memory efficient. Do not separate - combined_data[global_indices] = tffl.asarray( - key=file_indices - ) - start_page += len(tffl.pages) - except Exception as e: - raise Exception(f"Error in processing tiff file {input_file}: {e}") + written[[out_idx for out_idx, _ in entries]] = True + except Exception as e: + raise Exception(f"Error in processing tiff file {input_file}: {e}") + + if not written.all(): + missing = int((~written).sum()) + raise ValueError( + f"{missing} of {written.size} frames were never written for " + f"plane {plane_idx} channel {channel}. The movie would contain " + f"uninitialised memory, so it is not being saved." + ) output_tiff_fullpath = output_dir / f"{output_tiff_stem}.tif" tifffile.imwrite( @@ -401,27 +471,57 @@ def _extract_prairieview_metadata(xml_filepath: str): "/SubindexedValue/[@subindex='{subindex}']" ) + _z_subpath = ( + "PVStateShard/PVStateValue/[@key='positionCurrent']" + "/SubindexedValues/[@index='ZAxis']" + "/SubindexedValue/[@subindex='{subindex}']" + ) + + def _cycle_z(target_cycle): + """One z position per frame of `target_cycle`, in acquisition order. + + `PVStateShard` records state *changes*, so a frame whose z has not + moved since the previous frame omits `positionCurrent` entirely and + inherits it. Counting the elements in a cycle therefore counts + re-declarations, not planes: a recording that happens not to + re-declare z on its first frame reports fewer positions than it has + planes. Values are accumulated from the document-level shard through + every preceding frame instead. + """ + seed = xml_root.find(_z_subpath.format(subindex=active_subindex)) + current = float(seed.attrib["value"]) if seed is not None else None + positions = [] + for sequence in xml_root.findall(".//Sequence[@cycle]"): + cycle = sequence.attrib.get("cycle") + for frame in sequence.findall("Frame"): + declared = frame.find(_z_subpath.format(subindex=active_subindex)) + if declared is not None: + current = float(declared.attrib["value"]) + if cycle == target_cycle: + if current is None: + # A caller downstream would otherwise receive None as + # a depth, and the length check below cannot see it. + raise ValueError( + f"No z position for subindex {active_subindex} is " + f"declared at the document level or in any frame " + f"preceding cycle {target_cycle}." + ) + positions.append(current) + if cycle == target_cycle: + return positions + return positions + if bidirection_z: # With bidirectional Z, even-numbered cycles scan planes in - # reverse z-order. Extract z-positions from cycle 1 (forward) - # so that fieldZ aligns with the plane_indices ordering. - z_fields = [ - float(z.attrib.get("value")) - for z in xml_root.findall( - _z_xpath.format(cycle="1", subindex=active_subindex) - ) - ] + # reverse z-order. Take z-positions from a forward (odd) cycle so + # that fieldZ aligns with the plane_indices ordering. + z_fields = _cycle_z("1") # Build mapping: plane_idx → Frame[@index] in backward (even) cycles. - # In forward cycles, Frame[@index] matches the plane ordering from - # cycle 1 directly. In backward cycles the z-positions are reversed, - # so a different index is needed to reach the same physical plane. - z_fields_bwd = [ - float(z.attrib.get("value")) - for z in xml_root.findall( - _z_xpath.format(cycle="2", subindex=active_subindex) - ) - ] + # In forward cycles, Frame[@index] matches the forward plane ordering + # directly. In backward cycles the z-positions are reversed, so a + # different index is needed to reach the same physical plane. + z_fields_bwd = _cycle_z("2") fwd_indices = sorted(plane_indices) fwd_z = dict(zip(fwd_indices, z_fields)) @@ -439,17 +539,14 @@ def _extract_prairieview_metadata(xml_filepath: str): f"Backward z-values: {bwd_z}" ) else: - z_fields = [ - float(z.attrib.get("value")) - for z in xml_root.findall( - _z_xpath.format(cycle="2", subindex=active_subindex) - ) - ] + z_fields = _cycle_z("2") bidi_z_index_map = None - assert ( - len(z_fields) == n_depths - ), "Number of z fields does not match number of depths." + assert len(z_fields) == n_depths, ( + f"Recovered {len(z_fields)} z positions for {n_depths} planes. Each " + f"frame of a cycle should yield one position once inherited values " + f"are carried forward." + ) metainfo = dict( num_fields=n_depths, diff --git a/element_interface/version.py b/element_interface/version.py index ca67c3c..56a1710 100644 --- a/element_interface/version.py +++ b/element_interface/version.py @@ -1,3 +1,3 @@ """Package metadata""" -__version__ = "0.8.3" +__version__ = "0.8.4"