<commit_before>#!/usr/bin/env python3

import sys
import os
import os.path
import math
import czifile
import numpy as np
from scipy.misc import toimage
import json
import re
from lxml import etree

class LazyCziConverter (object):

    def __init__(self, czifilename, renormalize=False):
        self._fo = czifile.CziFile(czifilename)

        # sanity check dimensions
        self._Z = self._fo.axes.find('Z')
        if self._Z < 0:
            self._Z = None
        self._Y = self._fo.axes.find('Y')
        self._X = self._fo.axes.find('X')
        assert self._Y >= 0
        assert self._X == self._Y + 1
        self._slc = slice(self._Y,self._X+1)

        # check interleaved pixel config
        assert (self._slc.stop+1) == len(self._fo.axes), (self._slc.stop, len(self._fo.axes))
        assert self._fo.shape[self._slc.stop] in [1, 3], self._fo.shape
        assert self._fo.axes[self._slc] == 'YX', self._fo.axes[self._slc]

        # check multi-channel config
        self._C = self._fo.axes.find('C')
        assert self._C >= 0
        channel = self._fo.start[self._C]
        assert self._fo.shape[self._C] == 1 or self._fo.shape[self._slc.stop] == 1, \
            "do not understand multi-channel interleaved shape %s" % self._fo.shape

        # sort out segments into channel zoom tiers
        self._channel_tiers = [ dict() for c in range(self._fo.shape[self._C]) ]
        self._channel_tier_maps = [ dict() for c in range(self._fo.shape[self._C]) ]

        v0 = None
        v1 = None
        tile_size = (0, 0)

        for entry in self._fo.subblock_directory:
            assert len(entry.shape) == len(self._fo.shape)
            assert entry.axes == self._fo.axes

            channel = entry.start[self._C]
            
            # HACK: discard all but middle plane if there is a Z stack
            if self._Z is not None:
                if entry.start[self._Z] != (entry.shape[self._Z]//2):
                    continue
        
            # infer zoom from canvas size to tile size ratio
            zoom = [c//t for c,t in zip(entry.shape[self._slc], entry.stored_shape[self._slc])]
            assert zoom[0] == zoom[1]
            zoom = zoom[0]
        
            if zoom not in self._channel_tiers[channel]:
                self._channel_tiers[channel][zoom] = []

            # tile bounding box in canvas coordinates
            bbox = (entry.start[self._slc], tuple(map(lambda s, l: s+l, entry.start[self._slc], entry.shape[self._slc])))
            
            self._channel_tiers[channel][zoom].append((bbox, entry))

            # Figure total shape too...
            if v0 is None:
                v0 = bbox[0]
                v1 = bbox[1]
            else:
                v0 = [min(v, b) for v,b in zip(v0, bbox[0])]
                v1 = [max(v, b) for v,b in zip(v1, bbox[1])]
            tile_size = tuple(map(lambda a, b: max(a, b), tile_size, entry.stored_shape[self._slc]))

        # build up tile bbox maps for each zoom tier for efficient intersection tests
        for channel in range(self._fo.shape[self._C]):
            for zoom, bbox_entries in self._channel_tiers[channel].items():
                # pack as (v0.y, v0.x, v1.y, v1.x) for each entry
                bboxes = np.zeros( (len(bbox_entries), 4), dtype=np.int32 )
                for i in range(bboxes.shape[0]):
                    bboxes[i,0:2] = np.array(bbox_entries[i][0][0], dtype=np.int32)
                    bboxes[i,2:4] = np.array(bbox_entries[i][0][1], dtype=np.int32)
                # for fancy-indexing by boolean arrays, repeat entry indexes into array
                boxmap = np.array( range(len(bbox_entries)), dtype=np.int32 )
                self._channel_tier_maps[channel][zoom] = (bboxes, boxmap)

        cm = etree.fromstring(self._fo.metadata.encode('utf-8'))
        channels = cm.findall('Metadata/DisplaySetting/Channels/Channel')
        assert channels, 'found no Metadata/DisplaySetting/Channels/Channel elements in CZI metadata'
        self._channel_names_long = [ c.get('Name') for c in channels ]
        self._channel_names = [ c.find('ShortName').text for c in channels ]
        self._channel_colors = [ c.find('Color').text if c.find('Color') is not None else None for c in channels ]
        
        self._bbox_native = (tuple(v0), tuple(v1))
        self._bbox_zeroed = ((0, 0), (v1[0]-v0[0], v1[1]-v0[1]))
        self._tile_size = tile_size
        self._zoom_levels = self._channel_tiers[0].keys()
        self._zoom_levels.sort()
        
        sys.stderr.write('CZI %s tile-size %s %s\n  channels: %s\n  bounding-box: %s native or shape %s\n  zoom levels: %s\n' % (
            ' '.join(map(lambda d, s: '%s=%d' % (d, s), self._fo.axes, self._fo.shape)),
            'x'.join(map(str, self._tile_size)), self._fo.dtype,
            ', '.join([
                '%s (%s %s)' % (self._channel_names[i], self._channel_names_long[i], self._channel_colors[i])
                for i in range(self._fo.shape[self._C])
            ]),
            self._bbox_native, 'x'.join(map(str, self._bbox_zeroed[1])), self._zoom_levels
        ))

        if renormalize and self._fo.dtype not in [np.uint8]:
            sys.stderr.write('Finding per-channel value ranges for dynamic normalization...\n')
            self._channel_ranges = [ self._value_range(channel) for channel in range(self._fo.shape[self._C]) ]
        else:
            self._channel_ranges = None

        # HACK: try to configure a tile cache for row-major tile grid traversals
        # 1. assume we will output tiles no larger than current tile size
        # 2. assume that up to 3 tile rows might be active due to overlapping source tiles
        try:
            overlap_factor = max([float(e.text) for e in cm.findall('Metadata/Experiment/ExperimentBlocks/AcquisitionBlock/SubDimensionSetups/RegionsSetup/SampleHolder/Overlap') ])
        except:
            overlap_factor = max([float(e.text) for e in cm.findall('Metadata/Experiment/ExperimentBlocks/AcquisitionBlock/TilesSetup/SampleHolder/Overlap') ])

        row_tile_count = math.ceil(self._bbox_zeroed[1][1] / (self._tile_size[1] - self._tile_size[1] * overlap_factor))
        self._tile_cache_size = row_tile_count * 3 + 1
        # use dict as cache... key: (entry, dtype), value: [data, clock]
        self._tile_cache = dict()
        self._tile_cache_clock = 0

        sys.stderr.write('Estimating %d tiles per row at zoom level 1 or %d tile cache size\n' % (row_tile_count, self._tile_cache_size))

        # get per-dimension distances and turn meter value into micrometer
        self.mpps = dict([
            (dx.get('Id'), float(dx.find('Value').text) * 1E6)
            for dx in cm.findall('Metadata/Scaling/Items/Distance')
        ])
        sys.stderr.write('Image reported microns per pixel: %s\n' % self.mpps)

    def canvas_size(self):
        return self._bbox_zeroed[1]

    def num_channels(self):
        return len(self._channel_tiers)
            
    def _value_range(self, channelno):
        v0, v1 = None, None
        zooms = self._channel_tiers[channelno].keys()
        zooms.sort()
        for bbox, entry in self._channel_tiers[channelno][1]:
            data = self._entry_asarray(entry)
            if v0 is None:
                v0 = data.min()
                v1 = data.max()
            else:
                v0 = min(v0, data.min())
                v1 = max(v1, data.max())
        sys.stderr.write('Channel %d: %s .. %s (%s)\n' % (channelno, v0, v1, self._channel_names[channelno]))
        return (v0, v1)

    def _get_intersecting_bbox_entries(self, channelno, zoom, bbox_native):
        """Find CZI tiles that intersect bbox in native canvas coordinates."""

        bboxes, boxmap = self._channel_tier_maps[channelno][zoom]

        # non-intersects are either above or below the desired bbox
        y_nonintersects = (bboxes[:,2] <= bbox_native[0]) + (bboxes[:,0] >= bbox_native[2])
        x_nonintersects = (bboxes[:,3] <= bbox_native[1]) + (bboxes[:,1] >= bbox_native[3])
        nonintersects = y_nonintersects + x_nonintersects

        # use fancy-indexing with intersection boolean map to get sparse list of entry indices<commit_msg>Fixed sort() function for a dictionary keys.<commit_after>#!/usr/bin/env python3

import sys
import os
import os.path
import math
import czifile
import numpy as np
from scipy.misc import toimage
import json
import re
from lxml import etree

class LazyCziConverter (object):

    def __init__(self, czifilename, renormalize=False):
        self._fo = czifile.CziFile(czifilename)

        # sanity check dimensions
        self._Z = self._fo.axes.find('Z')
        if self._Z < 0:
            self._Z = None
        self._Y = self._fo.axes.find('Y')
        self._X = self._fo.axes.find('X')
        assert self._Y >= 0
        assert self._X == self._Y + 1
        self._slc = slice(self._Y,self._X+1)

        # check interleaved pixel config
        assert (self._slc.stop+1) == len(self._fo.axes), (self._slc.stop, len(self._fo.axes))
        assert self._fo.shape[self._slc.stop] in [1, 3], self._fo.shape
        assert self._fo.axes[self._slc] == 'YX', self._fo.axes[self._slc]

        # check multi-channel config
        self._C = self._fo.axes.find('C')
        assert self._C >= 0
        channel = self._fo.start[self._C]
        assert self._fo.shape[self._C] == 1 or self._fo.shape[self._slc.stop] == 1, \
            "do not understand multi-channel interleaved shape %s" % self._fo.shape

        # sort out segments into channel zoom tiers
        self._channel_tiers = [ dict() for c in range(self._fo.shape[self._C]) ]
        self._channel_tier_maps = [ dict() for c in range(self._fo.shape[self._C]) ]

        v0 = None
        v1 = None
        tile_size = (0, 0)

        for entry in self._fo.subblock_directory:
            assert len(entry.shape) == len(self._fo.shape)
            assert entry.axes == self._fo.axes

            channel = entry.start[self._C]
            
            # HACK: discard all but middle plane if there is a Z stack
            if self._Z is not None:
                if entry.start[self._Z] != (entry.shape[self._Z]//2):
                    continue
        
            # infer zoom from canvas size to tile size ratio
            zoom = [c//t for c,t in zip(entry.shape[self._slc], entry.stored_shape[self._slc])]
            assert zoom[0] == zoom[1]
            zoom = zoom[0]
        
            if zoom not in self._channel_tiers[channel]:
                self._channel_tiers[channel][zoom] = []

            # tile bounding box in canvas coordinates
            bbox = (entry.start[self._slc], tuple(map(lambda s, l: s+l, entry.start[self._slc], entry.shape[self._slc])))
            
            self._channel_tiers[channel][zoom].append((bbox, entry))

            # Figure total shape too...
            if v0 is None:
                v0 = bbox[0]
                v1 = bbox[1]
            else:
                v0 = [min(v, b) for v,b in zip(v0, bbox[0])]
                v1 = [max(v, b) for v,b in zip(v1, bbox[1])]
            tile_size = tuple(map(lambda a, b: max(a, b), tile_size, entry.stored_shape[self._slc]))

        # build up tile bbox maps for each zoom tier for efficient intersection tests
        for channel in range(self._fo.shape[self._C]):
            for zoom, bbox_entries in self._channel_tiers[channel].items():
                # pack as (v0.y, v0.x, v1.y, v1.x) for each entry
                bboxes = np.zeros( (len(bbox_entries), 4), dtype=np.int32 )
                for i in range(bboxes.shape[0]):
                    bboxes[i,0:2] = np.array(bbox_entries[i][0][0], dtype=np.int32)
                    bboxes[i,2:4] = np.array(bbox_entries[i][0][1], dtype=np.int32)
                # for fancy-indexing by boolean arrays, repeat entry indexes into array
                boxmap = np.array( range(len(bbox_entries)), dtype=np.int32 )
                self._channel_tier_maps[channel][zoom] = (bboxes, boxmap)

        cm = etree.fromstring(self._fo.metadata.encode('utf-8'))
        channels = cm.findall('Metadata/DisplaySetting/Channels/Channel')
        assert channels, 'found no Metadata/DisplaySetting/Channels/Channel elements in CZI metadata'
        self._channel_names_long = [ c.get('Name') for c in channels ]
        self._channel_names = [ c.find('ShortName').text for c in channels ]
        self._channel_colors = [ c.find('Color').text if c.find('Color') is not None else None for c in channels ]
        
        self._bbox_native = (tuple(v0), tuple(v1))
        self._bbox_zeroed = ((0, 0), (v1[0]-v0[0], v1[1]-v0[1]))
        self._tile_size = tile_size
        self._zoom_levels = self._channel_tiers[0].keys()
        self._zoom_levels = sorted(_zoom_levels)
        
        sys.stderr.write('CZI %s tile-size %s %s\n  channels: %s\n  bounding-box: %s native or shape %s\n  zoom levels: %s\n' % (
            ' '.join(map(lambda d, s: '%s=%d' % (d, s), self._fo.axes, self._fo.shape)),
            'x'.join(map(str, self._tile_size)), self._fo.dtype,
            ', '.join([
                '%s (%s %s)' % (self._channel_names[i], self._channel_names_long[i], self._channel_colors[i])
                for i in range(self._fo.shape[self._C])
            ]),
            self._bbox_native, 'x'.join(map(str, self._bbox_zeroed[1])), self._zoom_levels
        ))

        if renormalize and self._fo.dtype not in [np.uint8]:
            sys.stderr.write('Finding per-channel value ranges for dynamic normalization...\n')
            self._channel_ranges = [ self._value_range(channel) for channel in range(self._fo.shape[self._C]) ]
        else:
            self._channel_ranges = None

        # HACK: try to configure a tile cache for row-major tile grid traversals
        # 1. assume we will output tiles no larger than current tile size
        # 2. assume that up to 3 tile rows might be active due to overlapping source tiles
        try:
            overlap_factor = max([float(e.text) for e in cm.findall('Metadata/Experiment/ExperimentBlocks/AcquisitionBlock/SubDimensionSetups/RegionsSetup/SampleHolder/Overlap') ])
        except:
            overlap_factor = max([float(e.text) for e in cm.findall('Metadata/Experiment/ExperimentBlocks/AcquisitionBlock/TilesSetup/SampleHolder/Overlap') ])

        row_tile_count = math.ceil(self._bbox_zeroed[1][1] / (self._tile_size[1] - self._tile_size[1] * overlap_factor))
        self._tile_cache_size = row_tile_count * 3 + 1
        # use dict as cache... key: (entry, dtype), value: [data, clock]
        self._tile_cache = dict()
        self._tile_cache_clock = 0

        sys.stderr.write('Estimating %d tiles per row at zoom level 1 or %d tile cache size\n' % (row_tile_count, self._tile_cache_size))

        # get per-dimension distances and turn meter value into micrometer
        self.mpps = dict([
            (dx.get('Id'), float(dx.find('Value').text) * 1E6)
            for dx in cm.findall('Metadata/Scaling/Items/Distance')
        ])
        sys.stderr.write('Image reported microns per pixel: %s\n' % self.mpps)

    def canvas_size(self):
        return self._bbox_zeroed[1]

    def num_channels(self):
        return len(self._channel_tiers)
            
    def _value_range(self, channelno):
        v0, v1 = None, None
        zooms = self._channel_tiers[channelno].keys()
        zooms = sorted(zooms)
        for bbox, entry in self._channel_tiers[channelno][1]:
            data = self._entry_asarray(entry)
            if v0 is None:
                v0 = data.min()
                v1 = data.max()
            else:
                v0 = min(v0, data.min())
                v1 = max(v1, data.max())
        sys.stderr.write('Channel %d: %s .. %s (%s)\n' % (channelno, v0, v1, self._channel_names[channelno]))
        return (v0, v1)

    def _get_intersecting_bbox_entries(self, channelno, zoom, bbox_native):
        """Find CZI tiles that intersect bbox in native canvas coordinates."""

        bboxes, boxmap = self._channel_tier_maps[channelno][zoom]

        # non-intersects are either above or below the desired bbox
        y_nonintersects = (bboxes[:,2] <= bbox_native[0]) + (bboxes[:,0] >= bbox_native[2])
        x_nonintersects = (bboxes[:,3] <= bbox_native[1]) + (bboxes[:,1] >= bbox_native[3])
        nonintersects = y_nonintersects + x_nonintersects

        # use fancy-indexing with intersection boolean map to get sparse list of entry indices
{'fraction_non_alphanumeric': 0.12592637589600292, 'fraction_numerical': 0.012635159761875835, 'mean_word_length': 2.258709422011085, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1849372', 'n_tokens_mistral': 5428, 'n_tokens_neox': 5200, 'n_words': 1295}