Skills › Scientific Agent Skills › pathml › image_qc.py
image_qc.py
327 lines · skills/pathml/scripts/image_qc.py at 5f2bc03
1#!/usr/bin/env python3
2"""Compute bounded synthetic/local image QC and a coarse tissue-like mask."""
4from __future__ import annotations
6import argparse
7import math
8from pathlib import Path
9from typing import Any
11from _common import (
12 CliError,
13 MAX_IMAGE_BYTES,
14 MAX_PIXELS,
15 atomic_write_bytes,
16 checked_input_file,
17 checked_root,
18 emit_json,
19 run_cli,
20)
23IMAGE_SUFFIXES = {".ppm", ".pgm", ".png", ".jpg", ".jpeg", ".tif", ".tiff"}
24WHITESPACE = b" \t\r\n\f\v"
27def _next_token(payload: bytes, offset: int) -> tuple[bytes, int]:
28 length = len(payload)
29 while offset < length:
30 if payload[offset] in WHITESPACE:
31 offset += 1
32 continue
33 if payload[offset] == ord("#"):
34 newline = payload.find(b"\n", offset)
35 if newline < 0:
36 raise CliError("unterminated PNM comment")
37 offset = newline + 1
38 continue
39 break
40 start = offset
41 while offset < length and payload[offset] not in WHITESPACE:
42 if payload[offset] == ord("#"):
43 break
44 offset += 1
45 if start == offset:
46 raise CliError("invalid PNM header")
47 return payload[start:offset], offset
50def _read_pnm(path: Path, max_pixels: int) -> tuple[int, int, bytes]:
51 payload = path.read_bytes()
52 magic, offset = _next_token(payload, 0)
53 width_token, offset = _next_token(payload, offset)
54 height_token, offset = _next_token(payload, offset)
55 maxval_token, offset = _next_token(payload, offset)
56 if magic not in {b"P5", b"P6"}:
57 raise CliError("only binary P5/P6 PNM images are supported without Pillow")
58 try:
59 width = int(width_token)
60 height = int(height_token)
61 maxval = int(maxval_token)
62 except ValueError as exc:
63 raise CliError("PNM dimensions and max value must be integers") from exc
64 if width <= 0 or height <= 0 or width * height > max_pixels:
65 raise CliError(f"PNM pixel count must be in [1, {max_pixels}]")
66 if not 1 <= maxval <= 255:
67 raise CliError("PNM max value must be between 1 and 255")
68 if offset >= len(payload) or payload[offset] not in WHITESPACE:
69 raise CliError("PNM header must end with whitespace")
70 if payload[offset : offset + 2] == b"\r\n":
71 offset += 2
72 else:
73 offset += 1
74 channels = 3 if magic == b"P6" else 1
75 expected = width * height * channels
76 pixels = payload[offset:]
77 if len(pixels) != expected:
78 raise CliError(
79 f"PNM pixel payload is {len(pixels)} bytes; expected {expected}"
80 )
81 if pixels and max(pixels) > maxval:
82 raise CliError("PNM pixel value exceeds the declared maximum")
83 if maxval != 255:
84 pixels = bytes(round(value * 255 / maxval) for value in pixels)
85 if channels == 1:
86 rgb = bytearray(width * height * 3)
87 for index, value in enumerate(pixels):
88 start = index * 3
89 rgb[start : start + 3] = bytes((value, value, value))
90 pixels = bytes(rgb)
91 return width, height, pixels
94def _read_with_pillow(path: Path, max_pixels: int) -> tuple[int, int, bytes]:
95 try:
96 from PIL import Image
97 except ModuleNotFoundError as exc:
98 raise CliError(
99 "Pillow is required for PNG/JPEG/TIFF input; use a P5/P6 PNM "
100 "file for dependency-free inspection"
101 ) from exc
102 Image.MAX_IMAGE_PIXELS = max_pixels
103 try:
104 with Image.open(path) as image:
105 width, height = image.size
106 if width <= 0 or height <= 0 or width * height > max_pixels:
107 raise CliError(f"image pixel count must be in [1, {max_pixels}]")
108 rgb = image.convert("RGB")
109 return int(width), int(height), rgb.tobytes()
110 except CliError:
111 raise
112 except Exception as exc:
113 raise CliError(f"cannot decode bounded image: {exc}") from exc
116def read_image(path: Path, max_pixels: int) -> tuple[int, int, bytes]:
117 """Read a bounded RGB image with PNM support in the standard library."""
119 if path.suffix.lower() in {".ppm", ".pgm"}:
120 return _read_pnm(path, max_pixels)
121 return _read_with_pillow(path, max_pixels)
124def synthetic_image(
125 width: int,
126 height: int,
127 tissue_fraction: float,
128) -> bytes:
129 """Create a deterministic white image with a centered pink tissue region."""
131 if width <= 0 or height <= 0 or width * height > MAX_PIXELS:
132 raise CliError(f"synthetic pixel count must be in [1, {MAX_PIXELS}]")
133 if not math.isfinite(tissue_fraction) or not 0 <= tissue_fraction <= 1:
134 raise CliError("--tissue-fraction must be finite and in [0, 1]")
135 pixels = bytearray(bytes((245, 245, 245)) * (width * height))
136 scale = math.sqrt(tissue_fraction)
137 tissue_width = min(width, round(width * scale))
138 tissue_height = min(height, round(height * scale))
139 left = (width - tissue_width) // 2
140 top = (height - tissue_height) // 2
141 for row in range(top, top + tissue_height):
142 for column in range(left, left + tissue_width):
143 offset = (row * width + column) * 3
144 pixels[offset : offset + 3] = bytes((180, 75, 135))
145 return bytes(pixels)
148def calculate_qc(
149 width: int,
150 height: int,
151 pixels: bytes,
152 *,
153 saturation_threshold: float,
154 brightness_threshold: float,
155) -> tuple[dict[str, Any], bytes]:
156 """Calculate simple bounded QC metrics and an 8-bit binary mask."""
158 if len(pixels) != width * height * 3:
159 raise CliError("RGB payload length does not match dimensions")
160 if not math.isfinite(saturation_threshold) or not (
161 0 <= saturation_threshold <= 1
162 ):
163 raise CliError("--saturation-threshold must be finite and in [0, 1]")
164 if not math.isfinite(brightness_threshold) or not (
165 0 <= brightness_threshold <= 1
166 ):
167 raise CliError("--brightness-threshold must be finite and in [0, 1]")
169 count = width * height
170 sums = [0, 0, 0]
171 tissue = 0
172 white = 0
173 dark = 0
174 saturated_channel = 0
175 mask = bytearray(count)
176 for pixel_index in range(count):
177 offset = pixel_index * 3
178 red, green, blue = pixels[offset : offset + 3]
179 sums[0] += red
180 sums[1] += green
181 sums[2] += blue
182 high = max(red, green, blue)
183 low = min(red, green, blue)
184 brightness = high / 255.0
185 saturation = 0.0 if high == 0 else (high - low) / high
186 is_tissue = (
187 saturation >= saturation_threshold
188 and brightness <= brightness_threshold
189 )
190 if is_tissue:
191 tissue += 1
192 mask[pixel_index] = 255
193 if brightness >= 0.90 and saturation <= 0.10:
194 white += 1
195 if brightness <= 0.15:
196 dark += 1
197 if high == 255:
198 saturated_channel += 1
200 report = {
201 "width": width,
202 "height": height,
203 "pixel_count": count,
204 "mean_rgb": [round(total / count, 6) for total in sums],
205 "tissue_like_fraction": round(tissue / count, 8),
206 "white_like_fraction": round(white / count, 8),
207 "dark_fraction": round(dark / count, 8),
208 "any_channel_clipped_high_fraction": round(saturated_channel / count, 8),
209 "mask_rule": {
210 "minimum_saturation": saturation_threshold,
211 "maximum_brightness": brightness_threshold,
212 },
213 "clinical_use": False,
214 "note": (
215 "This coarse RGB heuristic is for synthetic/pilot QC only; it is not "
216 "PathML TissueDetectionHE and is not a diagnostic quality decision."
217 ),
218 }
219 return report, bytes(mask)
222def _write_mask(
223 mask: bytes,
224 width: int,
225 height: int,
226 *,
227 destination: str,
228 root: str,
229 force: bool,
230) -> None:
231 header = f"P5\n{width} {height}\n255\n".encode("ascii")
232 atomic_write_bytes(
233 Path(destination),
234 header + mask,
235 root=root,
236 suffixes={".pgm"},
237 force=force,
238 )
241def build_parser() -> argparse.ArgumentParser:
242 parser = argparse.ArgumentParser(
243 description=(
244 "Compute bounded synthetic/local RGB QC and a coarse tissue-like "
245 "mask. No network access."
246 )
247 )
248 subparsers = parser.add_subparsers(dest="command", required=True)
250 synthetic = subparsers.add_parser(
251 "synthetic", help="generate an in-memory synthetic image and report QC"
252 )
253 synthetic.add_argument("--width", type=int, default=256)
254 synthetic.add_argument("--height", type=int, default=256)
255 synthetic.add_argument("--tissue-fraction", type=float, default=0.5)
256 synthetic.add_argument("--saturation-threshold", type=float, default=0.15)
257 synthetic.add_argument("--brightness-threshold", type=float, default=0.95)
258 synthetic.add_argument("--mask-output")
259 synthetic.add_argument("--output")
260 synthetic.add_argument("--root", default=".")
261 synthetic.add_argument("--force", action="store_true")
263 inspect = subparsers.add_parser(
264 "inspect", help="inspect a bounded local raster or P5/P6 PNM image"
265 )
266 inspect.add_argument("--image", required=True)
267 inspect.add_argument("--root", default=".")
268 inspect.add_argument("--max-pixels", type=int, default=MAX_PIXELS)
269 inspect.add_argument("--max-image-bytes", type=int, default=MAX_IMAGE_BYTES)
270 inspect.add_argument("--saturation-threshold", type=float, default=0.15)
271 inspect.add_argument("--brightness-threshold", type=float, default=0.95)
272 inspect.add_argument("--mask-output")
273 inspect.add_argument("--output")
274 inspect.add_argument("--force", action="store_true")
275 return parser
278def main() -> None:
279 args = build_parser().parse_args()
280 root = checked_root(args.root)
281 if args.command == "synthetic":
282 width = args.width
283 height = args.height
284 pixels = synthetic_image(width, height, args.tissue_fraction)
285 source = "synthetic"
286 else:
287 if not 1 <= args.max_pixels <= 100_000_000:
288 raise CliError("--max-pixels must be between 1 and 100000000")
289 if not 1 <= args.max_image_bytes <= MAX_IMAGE_BYTES:
290 raise CliError(
291 f"--max-image-bytes must be between 1 and {MAX_IMAGE_BYTES}"
292 )
293 image = checked_input_file(
294 args.image,
295 root=root,
296 suffixes=IMAGE_SUFFIXES,
297 max_bytes=args.max_image_bytes,
298 )
299 width, height, pixels = read_image(image, args.max_pixels)
300 source = "local_file"
302 report, mask = calculate_qc(
303 width,
304 height,
305 pixels,
306 saturation_threshold=args.saturation_threshold,
307 brightness_threshold=args.brightness_threshold,
308 )
309 report["source"] = source
310 report["path_redacted"] = True
311 if args.mask_output:
312 _write_mask(
313 mask,
314 width,
315 height,
316 destination=args.mask_output,
317 root=str(root),
318 force=args.force,
319 )
320 report["mask_written"] = True
321 else:
322 report["mask_written"] = False
323 emit_json(report, output=args.output, root=root, force=args.force)
326if __name__ == "__main__":
327 raise SystemExit(run_cli(main))