Largeimages (#31)

* Improve handling of very large images (>4Gpix).

* Further performance improvements.

* Update JSDoc for compatibility reasons.

src/visiomatic/client/package.json

@@ -36,7 +36,7 @@
 	"devDependencies": {
 		"esbuild": "^0.14.53",
 		"eslint": "^8.23.1",
-		"jsdoc": "^4.0"
+		"jsdoc": "^4.0.4"
 	},
 	"main": "dist/visiomatic.js",
 	"scripts": {

src/visiomatic/server/image.py

@@ -66,7 +66,7 @@ def __init__(
         # Small hack to translate compressed FITS headers
         self.header = fits.Header.fromstring(hdu.header.tostring())
         shape = hdu.data.shape
-        self.data = hdu.data.astype(np.float32).reshape(-1,shape[-2], shape[-1])
+        self.data = hdu.data.reshape(-1,shape[-2], shape[-1])
         self.shape = self.data.shape
         self.bitpix = self.header["BITPIX"]
         self.bitdepth = 32
@@ -213,7 +213,7 @@ def compute_background(self, skip : int = 15) -> tuple[np.ndarray, np.ndarray]:
         """
         # NumPy version
         # Speed up ~x8 by using only a fraction of the lines
-        x = self.data[:, ::(skip + 1), :].reshape(self.data.shape[0],-1).copy()
+        x = self.data[:, ::(skip + 1), :].reshape(self.data.shape[0],-1).astype(np.float32)
         # First, normalize to avoid overflows
         norm = abs(np.nanmean(x))
         if norm > 0.:

src/visiomatic/server/tiled.py

@@ -188,7 +188,7 @@ def __init__(
             settings["thread_count"] if 'sphinx' not in modules else 4
         )
         try:
-            hdus = fits.open(self.filename)
+            hdus = fits.open(self.filename, mode='denywrite')
         except:
             raise(FileNotFoundError(f"Cannot open {filename}")) from None
             return
@@ -347,7 +347,7 @@ def make_mosaic(self, images : List[Image]) -> None:
             self.data_filename = get_data_filename(self.filename)
             self.data = np.memmap(
                 self.data_filename,
-                dtype=images[0].data.dtype,
+                dtype=np.float32,
                 mode='w+',
                 shape=shape
             )
@@ -570,43 +570,68 @@ def make_tiles(self) -> None:
         self.tiles_end = np.cumsum(self.counts, dtype=np.int32)
         self.tiles_start[1:] = self.tiles_end[:-1]
         self.tiles_filename = get_tiles_filename(self.filename)
+        tshape = self.tile_shape
         self.tiles: np.ndarray = np.memmap(
             self.tiles_filename,
             dtype=np.float32,
             mode='w+',
-            shape=(
-                self.ntiles,
-                self.tile_shape[0],
-                self.tile_shape[1],
-                self.tile_shape[2]
-            )
+            shape=(self.ntiles, *tshape)
         )
         ima = np.flip(self.data, axis=1)
-        for r in range(self.nlevels):
-            tiler = Tiler(
-                data_shape = ima.shape,
-                tile_shape = self.tile_shape,
-                mode='constant',
-                channel_dimension=0
-            )
-            self.tiles[self.tiles_start[r]:self.tiles_end[r]] = \
-                tiler.get_all_tiles(ima, copy_data=False)
-            # Pure NumPy approach if in multichannel mode
-            # else use OpenCV (faster but does not work with multiplanar data)
-            ima = ima[
-                :,
-                :(-1 if ima.shape[1]%2 else None),
-                :(-1 if ima.shape[2]%2 else None)
-            ].reshape(
-                ima.shape[0], ima.shape[1]//2, 2, -1, 2
-            ).mean(axis=2).mean(axis=3) if self.nchannels > 1 else \
-            cv2.resize(
-                ima[0].squeeze(),
-                fx=0.5,
-                fy=0.5,
-                dsize=(ima.shape[1]//2, ima.shape[0]//2),
-                interpolation=cv2.INTER_AREA
-            )[None,:,:]
+        tiler = Tiler(
+            data_shape = ima.shape,
+            tile_shape = tshape,
+            mode='constant',
+            channel_dimension=0
+        )
+        # We don't use get_all_tiles() as it does not work inplace
+        np.stack(
+            [
+                tiler.get_tile(ima, x, copy_data=False) \
+                    for x in range(tiler.n_tiles)
+            ],
+            axis=0,
+            out=self.tiles[self.tiles_start[0]:self.tiles_end[0]]
+        )
+        # Prepare array of 4 multidimensional numpy slices
+        slices = [
+            np.s_[:, :tshape[1], :tshape[2]],
+            np.s_[:, :tshape[1], tshape[2]:],
+            np.s_[:, tshape[1]:, :tshape[2]],
+            np.s_[:, tshape[1]:, tshape[2]:]
+        ]
+        # Initialize empty tile twice the size of a regular one, for rebinning
+        tile = np.zeros((tshape[0], tshape[1] * 2, tshape[2] * 2))
+        xdither = np.array((0, 1, 0, 1), dtype=np.int32)
+        ydither = np.array((0, 0, 1, 1), dtype=np.int32)
+        # Bin tiles down on subsequent levels
+        for r in range(1, self.nlevels):
+            # Previous resolution level
+            rp = r - 1
+            nt = self.counts[r]
+            ntp = self.counts[rp]
+            w = self.shapes[r, 2]
+            wp = self.shapes[rp, 2]
+            hp = self.shapes[rp, 1]
+            for n in range(nt):
+                xp = (n % w) * 2
+                yp = (n // w) * 2
+                xps = xp + xdither
+                yps = yp + ydither
+                nps = yps * wp + xps
+                nps[np.logical_or(xps >= wp, yps >= hp)] = -1
+                for i in range(4):
+                    tile[slices[i]] = self.tiles[self.tiles_start[rp] + nps[i]] \
+                        if nps[i] >= 0 else 0.
+                tile_out = self.tiles[self.tiles_start[r] + n]
+                for c in range(tshape[0]):
+                    tile_out[c] = cv2.resize(
+                    tile[c],
+                    fx=0.5,
+                    fy=0.5,
+                    dsize=(tshape[2], tshape[1]),
+                    interpolation=cv2.INTER_AREA
+                )
         del ima, tiler
         # Make sure that memory has been completely mapped before exiting
         self.tiles.flush()