Add StarDenoise source code

StarDenoise/README.md

@@ -0,0 +1,300 @@
+# StarDenoise：天文图像 AI 降噪 Demo
+
+这是一个可直接运行的 PyTorch 天文图像降噪 demo。复现对象是 **SetiAstro Cosmic Clarity Denoise AI 3.6** 的网络与官方预训练权重，并将原项目中面向桌面 GUI 的推理逻辑整理成适合工程集成的命令行模块。
+
+本目录已经包含：
+
+- PyTorch 模型结构与官方 AI 3.6 权重；
+- 从 AstroNoiseNet 官方仓库下载的 M51 短曝光/长曝光测试对；
+- JPG、PNG、TIFF、单色或 RGB FITS 的推理支持；
+- 自动 GPU/CPU 选择、重叠分块、边缘融合、降噪强度与模式控制；
+- FITS 线性数据的自动拉伸、网络推理、反拉伸，以及 FITS header 保留；
+- 下载校验脚本、简单的参考图 PSNR/MAE 评估和单元测试。
+- 三张 ESO 1920×1080 在线图片的完整测试、对比图与指标，见 [`README_1080P.md`](README_1080P.md)。
+
+> 注意：“星图降噪”与“去星”“星点缩小”“反卷积”不同。降噪要压制随机/结构噪声，同时尽量保持星点 PSF、亮度、星云微弱结构和色彩；它不能替代 bias/dark/flat 校准、配准和叠加。
+
+## 1. 快速运行
+
+工作目录：
+
+```powershell
+cd D:\Workman\MyProject\ImageStack\AstroColorProcess\StarDenoise
+conda activate study
+```
+
+环境依赖已经在 `study` 中安装过；在新机器上可执行：
+
+```powershell
+# 根据显卡/CUDA 选择 PyTorch 官方安装命令，然后：
+python -m pip install -r requirements.txt
+```
+
+当前已验证环境是 Python 3.10.20、PyTorch 2.13.0+cpu。由于本次网络下载 CUDA wheel 过慢，先安装 CPU 版保证 demo 即开即用；代码本身已支持 CUDA。要启用本机 RTX 4060，可在网络条件合适时按 [PyTorch 官方安装选择器](https://pytorch.org/get-started/locally/)重装 CUDA wheel，例如：
+
+```powershell
+python -m pip install --upgrade --force-reinstall torch torchvision `
+  --index-url https://download.pytorch.org/whl/cu126
+```
+
+运行内置 M51 样例：
+
+```powershell
+python denoise.py samples\M51_short.jpg outputs\M51_short_denoised.jpg `
+  --device auto --strength 0.9 --mode full --tile-size 256 --overlap 64
+```
+
+运行自己的 FITS：
+
+```powershell
+python denoise.py D:\path\input.fits outputs\input_denoised.fits `
+  --device cuda --strength 0.8 --mode luminance
+```
+
+常用参数：
+
+| 参数 | 默认值 | 说明 |
+|---|---:|---|
+| `--device` | `auto` | 自动选 CUDA，否则 CPU；也可指定 `cuda:0` |
+| `--strength` | `0.9` | 神经网络亮度降噪结果与原图的混合比例，范围 0～1 |
+| `--color-strength` | 同 `strength` | `full` 模式下的色度降噪强度 |
+| `--mode` | `full` | `luminance` 只处理亮度；`full` 亮度 AI + 引导色度滤波；`separate` 三通道各自 AI |
+| `--tile-size` | `256` | 官方模型训练/推理块大小；一般无需修改 |
+| `--overlap` | `64` | tile 重叠宽度，减小接缝 |
+| `--batch-size` | `1` | 8 GB 显存建议 1～2；显存足够再增加 |
+| `--precision` | `fp16` | CUDA 上用半精度；CPU 会自动使用 FP32 |
+| `--linear-stretch` | `auto` | 线性天文图自动 MTF 拉伸后推理并反拉伸，可设 `always/never` |
+
+比较短曝光输入、我们的结果与仓库给出的长曝光参考：
+
+```powershell
+python evaluate.py samples\M51_long_reference.jpg samples\M51_short.jpg
+python evaluate.py samples\M51_long_reference.jpg outputs\M51_short_denoised.jpg
+```
+
+长曝光参考并不是严格像素级 ground truth（曝光、配准、亮度都可能存在差异），因此该指标只能用于 demo 的相对观察，不能当成论文级结论。更可靠的验证见“训练与评估建议”。
+
+本机实测结果（CPU、400×400、FP32、strength=0.9、full）：
+
+| 项目 | 相对 M51 长曝光参考 |
+|---|---:|
+| 短曝光输入 | PSNR 30.588 dB；MAE 0.022448 |
+| 本 demo 输出 | **PSNR 35.322 dB；MAE 0.013325** |
+| 推理时间 | 3.414 秒 |
+
+官方权重严格加载成功，4 个单元测试（结构/参数量、颜色空间往返、零强度恒等、FITS channel-first 往返）全部通过。
+
+如需重新下载资产：
+
+```powershell
+python download_assets.py
+```
+
+权重下载后会检查官方公布的 SHA-256：`f860e234dc6a3ab1acacbea399c235179c86cb8b287f2ab9010bc18c8af08671`。
+
+## 2. 天文图像噪声与降噪原理
+
+天文图像的信号极弱，观测值可简化为：
+
+```text
+观测 = 天体光子 + 天空背景 + 暗电流 + 偏置 + 读出噪声 + 固定图样 + 离群点
+```
+
+主要成分如下：
+
+- **光子散粒噪声**：光子到达服从泊松统计，标准差约为 `sqrt(N)`；它是信号相关噪声，不是固定方差高斯噪声。
+- **天空背景噪声**：月光、光污染和大气辉光也产生光子与散粒噪声。
+- **读出噪声**：传感器读出电路在每个子曝光中引入，常近似高斯分布。
+- **暗电流及热像素**：热产生电子；其平均图样可由 dark frame 校正，但暗电流自身仍有随机性。
+- **偏置、PRNU、渐晕和尘埃**：分别通过 bias、flat 等校准处理；校准残差会表现为结构噪声。
+- **宇宙射线、卫星、坏点**：属于稀疏离群污染，通常应在叠加/专门剔除阶段处理，而非交给普通降噪器“猜”。
+
+[Astropy CCD reduction guide](https://www.astropy.org/ccd-reduction-and-photometry-guide/v/dev/notebooks/01-01-astronomical-CCD-image-components.html) 对 bias、dark、read noise、天空背景和 flat correction 有系统说明；NASA 的[传感器噪声与平方根定律说明](https://science.nasa.gov/wp-content/uploads/2023/09/Astrophotography_Guide.pdf)解释了泊松计数的 `sqrt(N)` 关系。2026 年的 [Denoising the Deep Sky](https://arxiv.org/abs/2601.23276) 进一步把光子散粒、PRNU、暗电流、读出效应、宇宙射线和热像素统一到物理 CCD 噪声合成管线中。
+
+传统流程首先做校准、配准和叠加。对独立随机噪声，平均 `N` 帧后噪声标准差大致降低到 `1/sqrt(N)`；深度模型学习的是一个条件先验：根据大量“短曝光/带噪输入 → 长曝光堆栈/低噪目标”样本，预测给定结构更可能对应真实星点和星云，哪些更可能是噪声。
+
+常见网络有两种输出形式：
+
+1. 直接输出干净图 `x_hat = f(y)`；
+2. 预测噪声残差 `n_hat = f(y)`，再计算 `x_hat = y - n_hat`。
+
+天文降噪比普通照片更怕“看起来漂亮但不真实”：很小的点状结构既可能是恒星，也可能是热像素/噪声；过强感知损失或 GAN 还可能生成不存在的结构。因此不仅要看 PSNR/SSIM，还要检查星点检出率、孔径测光偏差、PSF/FWHM、背景 NMAD、假阳性和弱结构保真度。
+
+## 3. 模型调研与选型
+
+调研日期：2026-08-31。这里把“学术主流”和“深空摄影社区实际可用”分开考虑。
+
+| 工作 | 领域匹配 | 开源/权重 | 规模与部署 | 结论 |
+|---|---|---|---|---|
+| [DnCNN](https://github.com/cszn/DnCNN) | 通用 AWGN | 官方代码和 PyTorch 权重；经典论文，社区最成熟 | 约 0.56M 参数，极轻 | 很适合基线，但高斯噪声假设与天文泊松/读噪/固定图样不完全匹配 |
+| [SCUNet](https://github.com/cszn/SCUNet) | 通用真实图盲降噪 | 官方 PyTorch；Swin + Conv，研究认可度高 | 约 17.95M 参数、FP32 约 72 MB，速度/显存成本较高 | 通用盲降噪强，但不是天文数据训练，弱星点有域偏移风险 |
+| [Astro U-net](https://github.com/Sponka/Astro_U-net) | HST WFC3 科学图像 | 论文、代码和补充数据开放 | U-Net；论文报告 SNR 平均提高 1.63×、恒星恢复率 95.9%、平均 flux error 2.26% | 科学依据强，但针对特定 HST 波段；直接用于消费级 RGB 深空摄影不够即插即用 |
+| [AstroNoiseNet](https://github.com/Steffenhir/AstroNoiseNet) | 真实短/长曝光深空摄影 | MIT，代码开放；PRIDNet + 判别器 | 较复杂，仓库没有便于独立部署的现成 PyTorch checkpoint | 数据构造思路很好；本 demo 的 M51 测试对来自该仓库 |
+| [GraXpert](https://github.com/Steffenhir/GraXpert) | 深空摄影 | GPL-3.0 应用开放，发布 ONNX 模型 | 产品成熟、约 300+ GitHub stars、支持 CLI/PixInsight/Siril 工作流 | 实际认可度最高之一，但模型以 ONNX 应用形式分发，独立 PyTorch 复现和再训练不如本方案直接 |
+| [Cosmic Clarity AI 3.6](https://github.com/setiastro/cosmicclarity) | 深空摄影 | **MIT，官方 PyTorch `.pth` + ONNX** | **2,945,283 参数；纯 FP32 参数约 11.24 MiB；官方 checkpoint 33.8 MB** | **本 demo 选择：领域匹配、体积、开放程度和部署成本最均衡** |
+| [Denoising the Deep Sky](https://github.com/ShuhongLL/Denoising-Deep-Sky) | 科研级 raw CCD | 代码、MuSCAT 数据、PyTorch 训练流程开放 | 物理噪声模型 + U-Net/PMN-U-Net，数据约 124 GB | 2026 年的新工作，科学训练路线很强，但还不够“轻量即用”，适合下一阶段自训 |
+
+选择 Cosmic Clarity 的理由：
+
+- 是针对实际深空摄影训练的模型，不是用自然图像 AWGN 权重硬套；
+- 官方直接发布 PyTorch checkpoint 和网络源代码，许可证为 MIT；
+- 约 2.94M 参数，显著小于 SCUNet，256 tile 可在 8 GB 显卡上轻松运行；
+- 结构是全分辨率残差 CNN，含 dilation 和多层 skip concatenation，适合保存点源与小尺度纹理；
+- 上游工具已形成社区使用，并有 Siril 相关集成；
+- 相比只发布 ONNX 的应用模型，更容易审阅、修改、继续训练和嵌入 Python 工程。
+
+局限也要明确：该模型没有正式 model card 公布训练集规模、相机分布和严格 benchmark；MIT 是仓库整体许可证，权重发布页没有更细的单独数据许可声明。它适合“工程 demo / 深空摄影后期”，若目标是可发表的测光或形态学分析，必须在自己的仪器和目标分布上重新评估或训练。
+
+## 4. 本复现与上游的对应关系
+
+网络完全对应官方 AI 3.6 PyTorch 定义：
+
+- 输入/输出：3 通道，范围 `[0,1]`；单色图复制为三通道，取输出第一通道；
+- 编码通道：`3 → 16 → 32 → 64 → 128 → 256`；
+- 64 和 256 通道入口使用 dilation=2；
+- 每级包含两层 3×3 卷积的 residual block；
+- 解码端逐级拼接 encoder 特征，最终 sigmoid 输出；
+- 模型不做 pooling，空间尺寸始终不变。
+
+本 demo 保留上游的 256 tile、64 overlap、MTF 线性拉伸/反拉伸、Y/Cb/Cr 亮度处理和引导色度降噪；工程化改动包括：
+
+- Hann 权重融合取代简单覆盖/平均，减少 tile 接缝；
+- `torch.inference_mode()`、CUDA FP16 和可调 batch；
+- `torch.load(..., weights_only=True)` 限制反序列化范围；
+- 独立、可测试的 I/O，支持 FITS header 与原始动态范围回写；
+- CLI 参数和 SHA-256 校验。
+
+官方来源和许可证副本保存在 `third_party/`。测试样例来自 AstroNoiseNet 的 `examples/M51_short.jpg` 与 `examples/M51_long.jpg`，均由本项目的下载脚本给出可复现 URL。
+
+## 5. 推荐处理顺序
+
+对消费级深空摄影数据，建议：
+
+```text
+raw 子帧
+  → bias/dark/flat 校准
+  → 星点配准
+  → 异常值剔除与叠加
+  → 背景/梯度校正、颜色校准
+  → 在线性数据上温和降噪（或按训练域选择拉伸后降噪）
+  → 反卷积/星点处理
+  → 非线性拉伸与最终调色
+```
+
+若模型明确在 stretched 图上训练（Cosmic Clarity 属于这一类工程工具），本 demo 的 `auto` 会在线性 FITS 中位数很低时先做可逆 MTF 拉伸。强度建议从 `0.6～0.8` 开始，放大 100% 检查弱星与星云边缘；不要只看缩略图。
+
+## 6. 如何获取自己的训练数据
+
+### 路线 A：自己的相机，最推荐
+
+对每个目标保存校准后的单帧，并生成高 SNR reference：
+
+1. 同一目标、相同相机/增益/温度/滤镜拍摄 20～100+ 个 dithered subframes；
+2. 每个子帧完成 bias/dark/flat 校准；
+3. 亚像素配准到同一坐标；
+4. 从其中一帧或少量帧构造 noisy input；其余帧用 sigma-clipping/加权平均生成 clean target；
+5. 避免把 input 本身同时放进 target stack，否则 target 噪声与输入相关，会高估效果；
+6. 不同观测目标划分 train/val/test，不能把同一张大图随机 patch 后分到不同集合；
+7. 保存 32-bit float linear FITS，以及 gain、read noise、曝光、温度、滤镜等 metadata。
+
+如果没有足够干净 target，可用 **Noise2Noise**：把同一视场奇数帧和偶数帧分别叠加成两张独立噪声图作为 pair。两边天体信号一致、随机噪声独立，使用 L1/L2 训练仍能学习期望的干净信号。也可用 Noise2Void/Self2Self 做单图自监督，但在密集星场要谨慎设计 blind-spot，避免把点源当噪声。
+
+### 路线 B：公开科学数据
+
+- [MuSCAT Raw Image dataset](https://huggingface.co/datasets/ToferFish/MuSCAT-RawImage)：约 124 GB，CC BY-NC 4.0；包含 MuSCAT-3/4 的 raw、overscan-corrected `OS`、calibrated、stacked `MEAN`、mask、bias/dark/flat。可构造 `OS → CALIB` 或 `synthetic → MEAN`，是当前最完整的物理 CCD 降噪数据之一。非商业限制必须注意。
+- [Denoising the Deep Sky 官方代码](https://github.com/ShuhongLL/Denoising-Deep-Sky)：包含 FITS→NPY、数据划分、U-Net/PMN-U-Net 训练、PSF 注入和 PSNR/SSIM/NMAD/FLUX-SNR 评估脚本。
+- [MAST](https://archive.stsci.edu/)：HST 等任务的公开观测档案。可按 Astro U-net 方法选择同一目标/滤镜的短曝光与深曝光或 drizzle 产品；务必读取仪器校准与 data-use 说明。
+- [Astro U-net 论文与代码](https://github.com/Sponka/Astro_U-net)：适合复现 HST WFC3 F555W/F606W 的科学图像训练设计。
+- [AstroNoiseNet](https://github.com/Steffenhir/AstroNoiseNet)：仓库展示真实短/长积分 pair 与数据预处理思路，但没有发布完整训练数据；可借鉴采集协议。
+
+公开数据和自己相机的数据不要直接混合后随机切分。不同相机的 PSF、像元、Bayer/mono、增益、读噪、滤镜和后处理差异很大，应先按 instrument 做域划分和归一化，并保留独立的跨仪器测试集。
+
+## 7. 自训练架构建议
+
+### 最低风险、最快落地
+
+直接继续训练本 demo 的 `CosmicClarityDenoiseCNN`：
+
+- 2.95M 参数，8 GB RTX 4060 可用 `256×256` crop、batch 2～4、AMP；
+- 单色科学 FITS 最好把首尾层改为 1 通道，不要无意义复制 RGB；
+- 输入预测 clean image，或者改为 residual noise prediction 并加全局 skip；
+- 在输入中追加 noise-level/gain/read-noise map，提升不同相机和曝光参数的泛化。
+
+### 追求更高上限
+
+使用小型 **Residual U-Net / NAFNet-small**，4 层尺度、base width 32，约 8～20M 参数。下采样扩大感受野，skip 保留星点定位；用 tile inference 控制显存。科研 raw CCD 可参考 PMN-U-Net，把物理噪声参数或固定图样作为条件输入。
+
+不建议第一版就用扩散模型或 GAN：推理���、训练成本高，而且感知上“锐利”的结果可能改变光度或生成假结构。天文图像首先应优化可测量的真实性。
+
+建议训练配置：
+
+| 项目 | 建议 |
+|---|---|
+| 数据域 | 线性、校准后、未锐化 FITS；每个 filter/instrument 记录 metadata |
+| crop | 256×256；至少一部分 patch 强制含星点/星云，另一部分纯背景 |
+| augmentation | 90° 旋转、翻转；Bayer raw 不能任意旋转/翻转后破坏 CFA 排列 |
+| optimizer | AdamW，初始 LR `1e-4`，cosine decay，50～150 epochs |
+| primary loss | Charbonnier 或 L1，比纯 MSE 更抗离群点 |
+| structure loss | `0.1～0.2 × (1-SSIM)` + 小权重 gradient loss |
+| astronomy loss | aperture-flux L1、PSF/FWHM loss 或 source-mask 加权；背景与源区分别统计 |
+| validation | 按目标/夜晚/仪器分组，绝不能 patch 泄漏 |
+| checkpoint | 以 flux bias + 背景 NMAD + SSIM 综合选模，而非只看训练 loss |
+
+物理噪声合成可在干净 stack 上加入：
+
+```text
+electrons = clean_ADU × gain
+shot      ~ Poisson(electrons)
+dark      ~ Poisson(dark_current × exposure)
+read      ~ Normal(0, read_noise_e)
+PRNU      = signal × pixel_response_map
+outliers  = hot pixels + cosmic ray events
+noisy_ADU = (shot + dark + read + PRNU + outliers) / gain
+```
+
+参数应从真实 bias/dark/flat 与相机 specification 估计，不要只用固定 `Gaussian(σ=25)`。合成数据用于扩充覆盖，最终仍应以真实独立短/长曝光 pair 微调和验证。
+
+## 8. 科学验证清单
+
+在把模型用于正式处理前，至少做以下检查：
+
+1. 背景：均值/中位数偏移、标准差、MAD/NMAD、空间相关噪声功率谱；
+2. 星点：检测 completeness、false positive、centroid 偏差、FWHM/椭率变化；
+3. 测光：多个亮度区间的 aperture flux ratio 与 magnitude bias；
+4. 扩展源：表面亮度、颜色、弱结构差值图；
+5. 泛化：未见过的目标、夜晚、相机、滤镜、增益和曝光；
+6. 伪影：星周黑环、彩斑、tile seam、背景塑料感、被抹掉的弱星；
+7. 与传统基线比较：stacking、BM3D/小波、DnCNN/SCUNet，以及不降噪输入。
+
+论文级评估可参考 Astro U-net 的星点恢复率与 flux error，以及 Denoising the Deep Sky 的 PSNR、SSIM、NMAD、source injection 与 FLUX-SNR。保存原始图、输出图、模型版本、权重 hash、参数和环境版本，保证可追溯。
+
+## 9. 目录结构与许可证
+
+```text
+StarDenoise/
+├── denoise.py                 # CLI 与推理管线
+├── model.py                   # AI 3.6 PyTorch 网络
+├── image_io.py                # JPG/PNG/TIFF/FITS I/O
+├── evaluate.py                # 样例全参考指标
+├── download_assets.py         # 官方资产下载与 hash 校验
+├── models/                    # 官方 checkpoint
+├── samples/                   # 官方 M51 短/长曝光样例
+├── outputs/                   # demo 输出
+├── tests/                     # 单元测试
+└── third_party/               # 上游许可证/审阅源码
+```
+
+Cosmic Clarity 和 AstroNoiseNet 均使用 MIT License；上游许可证副本位于 `third_party/`。本复现应保留这些版权与许可声明。测试样例的使用遵循其来源仓库条款；用于商业产品前，仍建议向权重作者确认模型训练数据与权重的商业授权范围。
+
+## 10. 主要来源
+
+- [Cosmic Clarity 官方仓库（MIT、PyTorch 源码）](https://github.com/setiastro/cosmicclarity)
+- [Cosmic Clarity 官方 Release（AI 3.6 PTH/ONNX 与 SHA-256）](https://github.com/setiastro/cosmicclarity/releases/tag/Linux)
+- [AstroNoiseNet 官方仓库与 M51 样例](https://github.com/Steffenhir/AstroNoiseNet)
+- [GraXpert 官方仓库与 CLI](https://github.com/Steffenhir/GraXpert)
+- [Learning to denoise astronomical images with U-nets（MNRAS）](https://academic.oup.com/mnras/article/503/3/3204/5989737)
+- [Denoising the Deep Sky（论文）](https://arxiv.org/abs/2601.23276)
+- [Denoising the Deep Sky（代码与数据说明）](https://github.com/ShuhongLL/Denoising-Deep-Sky)
+- [DnCNN 官方仓库](https://github.com/cszn/DnCNN)
+- [SCUNet 官方仓库](https://github.com/cszn/SCUNet)
+- [Astropy CCD reduction guide](https://www.astropy.org/ccd-reduction-and-photometry-guide/v/dev/notebooks/01-01-astronomical-CCD-image-components.html)

StarDenoise/README_1080P.md

@@ -0,0 +1,111 @@
+# 1080P 在线天文图片测试报告
+
+测试日期：2026-09-01。目标是验证 `StarDenoise` 在真正的 1920×1080 图像上能否完整运行，并观察不同天文内容对降噪强度的敏感性。
+
+## 测试图片与授权
+
+三张图都从欧洲南方天文台（ESO）官方下载，文件经 Pillow 验证为 `1920×1080 RGB`。
+
+| 测试图 | ESO 页面 | 下载文件 | 署名 |
+|---|---|---|---|
+| Eagle / Omega / Sh2-54 三大星云全景 | [vbg_011](https://www.eso.org/public/products/virtualbackgrounds/vbg_011/) | `inputs_1080p/references/vbg_011_nebulae.jpg` | ESO |
+| M17 恒星形成区 | [vbg_012](https://www.eso.org/public/products/virtualbackgrounds/vbg_012/) | `inputs_1080p/references/vbg_012_messier17.jpg` | ESO |
+| 猎户座星云 | [vbg_013](https://www.eso.org/public/products/virtualbackgrounds/vbg_013/) | `inputs_1080p/references/vbg_013_orion.jpg` | ESO/Igor Chekalin |
+
+除非页面另有说明，ESO 公共图片采用 [Creative Commons Attribution 4.0](https://www.eso.org/public/outreach/copyright/)；再发布图片或本报告中的对比图时必须保留上表完整署名。图中的 ESO 水印来自官方原图，本测试没有移除。
+
+## 为什么添加可控噪声
+
+ESO 发布图片已经经过校准、叠加和后期处理，本身接近低噪成片。如果直接对低噪图运行模型，只能看到轻微平滑，无法客观判断去噪正确性。
+
+因此测试保留下载图作为 reference，并用固定随机种子生成真实感更接近传感器的测试输入：
+
+```text
+shot noise  = Poisson(reference × 350) / 350
+read noise  = Normal(0, 0.010)
+row pattern = 每行每通道 Normal(0, 0.003)
+noisy       = clip(shot + read + row pattern, 0, 1)
+```
+
+生成的 `inputs_1080p/noisy/*.png` 仍然是 1920×1080。脚本保存 noisy PNG 后会重新读取该文件，再送入模型；输出 PNG 也重新读取后才计算指标，所以报告数字与交付文件完全一致。
+
+## 统一测试参数
+
+```text
+model           Cosmic Clarity AI 3.6 PyTorch
+device          CPU / torch 2.13.0+cpu
+resolution      1920×1080 RGB
+tile            256×256
+overlap         64
+mode            full（亮度 AI + 引导色度滤波）
+strength        0.85
+color-strength  0.75
+precision       FP32
+linear-stretch  never（输入已经拉伸）
+```
+
+运行：
+
+```powershell
+cd D:\Workman\MyProject\ImageStack\AstroColorProcess\StarDenoise
+conda activate study
+python benchmark_1080p.py
+```
+
+## 严格落盘结果
+
+| 图片 | CPU 耗时 | 输入 PSNR | 输出 PSNR | 输入 SSIM | 输出 SSIM | 输入 MAE | 输出 MAE | 结论 |
+|---|---:|---:|---:|---:|---:|---:|---:|---|
+| 三大星云全景 | 52.980 s | 30.603 | 27.274 | 0.89247 | 0.78516 | 0.022849 | 0.030590 | **0.85 过强，真实密集星点被平滑** |
+| M17 | 57.588 s | 30.724 | 30.781 | 0.75422 | 0.88554 | 0.021921 | 0.017849 | PSNR基本持平，结构和绝对误差改善 |
+| 猎户座星云 | 59.533 s | 28.392 | **34.168** | 0.50262 | **0.91295** | 0.029182 | **0.011197** | **明显有效** |
+
+平均 CPU 时间约 56.7 秒/张；输出尺寸全部保持 1920×1080，未发现 tile seam。
+
+## 密集星野的低强度复测
+
+三大星云全景包含极密集星点和大量细碎高频结构。统一的 `strength=0.85` 会把一部分真实结构识别成噪声，因此使用同一个 `denoise.py` 精确重跑：
+
+```powershell
+python denoise.py `
+  inputs_1080p\noisy\vbg_011_nebulae_noisy.png `
+  outputs_1080p\vbg_011_nebulae_denoised_strength030.png `
+  --device cpu --strength 0.30 --color-strength 0.25 --mode full `
+  --tile-size 256 --overlap 64 --precision fp32 --linear-stretch never
+```
+
+| 三大星云全景 | PSNR | SSIM | MAE | CPU 时间 |
+|---|---:|---:|---:|---:|
+| 带噪输入 | 30.603 | 0.89247 | 0.022849 | — |
+| strength=0.85 | 27.274 | 0.78516 | 0.030590 | 52.980 s |
+| **strength=0.30** | **32.300** | **0.92579** | **0.018252** | 51.126 s |
+
+这说明失败原因不是分辨率或 tile 推理错误，而是降噪强度与内容不匹配。网络计算量基本与 strength 无关，所以低强度不会显著加速；它只是把更多原始细节混回输出。
+
+## 对比图
+
+每张图从左到右都是：带噪输入、模型输出、ESO reference。
+
+### 三大星云全景——调优后 strength=0.30
+
+![Eagle Omega Sh2-54 tuned](outputs_1080p/comparisons/vbg_011_nebulae_strength030_comparison.jpg)
+
+### M17——strength=0.85
+
+![Messier 17](outputs_1080p/comparisons/vbg_012_messier17_comparison.jpg)
+
+### 猎户座星云——strength=0.85
+
+![Orion Nebula](outputs_1080p/comparisons/vbg_013_orion_comparison.jpg)
+
+## 结论
+
+1. 当前代码可以稳定处理 1920×1080 RGB 图片，输出尺寸正确且无明显分块接缝。
+2. 对猎户座这种大面积弥散结构，`0.8～0.9` 强度效果很好。
+3. 对 M17 这类星点与星云混合画面，强降噪会改善背景/结构一致性，但 PSNR不一定大幅增加，建议从 `0.5～0.7` 开始。
+4. 对极密集星野，建议 `0.25～0.40`；本次 `0.30` 明显优于 `0.85`。
+5. 仅以缩略图判断容易忽略被删除的弱星。产品化时应加入星点密度/高频能量检测，自动限制 strength，并提供 100% 局部预览。
+6. 本测试是受控合成噪声 benchmark，不等价于相机 raw/FITS 物理验证。正式使用仍应测试真实短曝光与长曝光 stack pair、星点 flux 和 FWHM。
+
+完整机器可读数据位于 `outputs_1080p/benchmark_results.json`，测试脚本为 `benchmark_1080p.py`。
+

StarDenoise/benchmark_1080p.py

@@ -0,0 +1,191 @@
+"""Reproducible 1080p benchmark using three CC BY 4.0 ESO images.
+
+The ESO publication images are already processed and comparatively clean, so this
+script creates deterministic Poisson + Gaussian sensor-noise inputs.  The original
+downloads remain untouched and act as full-reference targets.
+"""
+
+from __future__ import annotations
+
+import json
+import math
+import time
+from pathlib import Path
+
+import cv2
+import numpy as np
+import torch
+from PIL import Image, ImageDraw
+
+from denoise import denoise_array
+from model import load_model
+
+
+ROOT = Path(__file__).resolve().parent
+REFERENCE_DIR = ROOT / "inputs_1080p" / "references"
+NOISY_DIR = ROOT / "inputs_1080p" / "noisy"
+OUTPUT_DIR = ROOT / "outputs_1080p"
+COMPARISON_DIR = OUTPUT_DIR / "comparisons"
+WEIGHTS = ROOT / "models" / "deep_denoise_cnn_AI3_6.pth"
+
+IMAGES = {
+    "vbg_011_nebulae": {
+        "filename": "vbg_011_nebulae.jpg",
+        "title": "Eagle / Omega / Sh2-54 nebulae",
+        "source": "https://www.eso.org/public/products/virtualbackgrounds/vbg_011/",
+        "credit": "ESO",
+    },
+    "vbg_012_messier17": {
+        "filename": "vbg_012_messier17.jpg",
+        "title": "Messier 17",
+        "source": "https://www.eso.org/public/products/virtualbackgrounds/vbg_012/",
+        "credit": "ESO",
+    },
+    "vbg_013_orion": {
+        "filename": "vbg_013_orion.jpg",
+        "title": "Orion Nebula",
+        "source": "https://www.eso.org/public/products/virtualbackgrounds/vbg_013/",
+        "credit": "ESO/Igor Chekalin",
+    },
+}
+
+
+def load_rgb(path: Path) -> np.ndarray:
+    with Image.open(path) as image:
+        return np.asarray(image.convert("RGB"), dtype=np.float32) / 255.0
+
+
+def save_rgb(path: Path, image: np.ndarray) -> None:
+    path.parent.mkdir(parents=True, exist_ok=True)
+    array = np.clip(image * 255.0 + 0.5, 0, 255).astype(np.uint8)
+    Image.fromarray(array, "RGB").save(path)
+
+
+def add_sensor_noise(clean: np.ndarray, seed: int) -> np.ndarray:
+    """Add signal-dependent shot noise plus channel-wise read noise.
+
+    ``peak_electrons`` controls the Poisson component.  The values are deliberately
+    moderate: visible at 100% but not so strong that the test becomes artificial.
+    """
+    rng = np.random.default_rng(seed)
+    peak_electrons = 350.0
+    shot = rng.poisson(np.clip(clean, 0, 1) * peak_electrons).astype(np.float32) / peak_electrons
+    read = rng.normal(0.0, 0.010, clean.shape).astype(np.float32)
+    # A small row-correlated component approximates residual banding/read pattern.
+    row_pattern = rng.normal(0.0, 0.003, (clean.shape[0], 1, 3)).astype(np.float32)
+    return np.clip(shot + read + row_pattern, 0.0, 1.0)
+
+
+def psnr(reference: np.ndarray, candidate: np.ndarray) -> float:
+    mse = float(np.mean((reference - candidate) ** 2))
+    return math.inf if mse == 0 else 10.0 * math.log10(1.0 / mse)
+
+
+def ssim(reference: np.ndarray, candidate: np.ndarray) -> float:
+    """Standard local SSIM averaged over RGB channels."""
+    c1, c2 = 0.01**2, 0.03**2
+    scores = []
+    for channel in range(3):
+        x, y = reference[..., channel], candidate[..., channel]
+        mu_x = cv2.GaussianBlur(x, (11, 11), 1.5)
+        mu_y = cv2.GaussianBlur(y, (11, 11), 1.5)
+        sigma_x = cv2.GaussianBlur(x * x, (11, 11), 1.5) - mu_x * mu_x
+        sigma_y = cv2.GaussianBlur(y * y, (11, 11), 1.5) - mu_y * mu_y
+        sigma_xy = cv2.GaussianBlur(x * y, (11, 11), 1.5) - mu_x * mu_y
+        numerator = (2 * mu_x * mu_y + c1) * (2 * sigma_xy + c2)
+        denominator = (mu_x * mu_x + mu_y * mu_y + c1) * (sigma_x + sigma_y + c2)
+        scores.append(float(np.mean(numerator / np.maximum(denominator, 1e-12))))
+    return float(np.mean(scores))
+
+
+def mae(reference: np.ndarray, candidate: np.ndarray) -> float:
+    return float(np.mean(np.abs(reference - candidate)))
+
+
+def make_comparison(path: Path, noisy: np.ndarray, denoised: np.ndarray, reference: np.ndarray, title: str) -> None:
+    panel_size = (640, 360)
+    labels = ("NOISY INPUT", "DENOISED", "ESO REFERENCE")
+    images = (noisy, denoised, reference)
+    canvas = Image.new("RGB", (1920, 410), "#111111")
+    draw = ImageDraw.Draw(canvas)
+    draw.text((12, 8), title, fill="white")
+    for index, (label, array) in enumerate(zip(labels, images)):
+        panel = Image.fromarray(np.clip(array * 255.0 + 0.5, 0, 255).astype(np.uint8), "RGB")
+        panel = panel.resize(panel_size, Image.Resampling.LANCZOS)
+        x = index * panel_size[0]
+        canvas.paste(panel, (x, 50))
+        draw.text((x + 12, 30), label, fill="white")
+    path.parent.mkdir(parents=True, exist_ok=True)
+    canvas.save(path, quality=94)
+
+
+def main() -> None:
+    NOISY_DIR.mkdir(parents=True, exist_ok=True)
+    OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
+    COMPARISON_DIR.mkdir(parents=True, exist_ok=True)
+    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
+    model = load_model(WEIGHTS, device)
+    results = []
+
+    for seed, (key, metadata) in enumerate(IMAGES.items(), start=20260901):
+        reference_path = REFERENCE_DIR / metadata["filename"]
+        reference = load_rgb(reference_path)
+        if reference.shape[:2] != (1080, 1920):
+            raise ValueError(f"Expected 1920x1080: {reference_path} -> {reference.shape}")
+        noisy = add_sensor_noise(reference, seed)
+        noisy_path = NOISY_DIR / f"{key}_noisy.png"
+        save_rgb(noisy_path, noisy)
+        # Benchmark the exact 8-bit file that a user can pass to denoise.py.
+        noisy = load_rgb(noisy_path)
+
+        started = time.perf_counter()
+        denoised = denoise_array(
+            noisy,
+            model,
+            device,
+            strength=0.85,
+            color_strength=0.75,
+            mode="full",
+            tile_size=256,
+            overlap=64,
+            batch_size=1,
+            half_precision=device.type == "cuda",
+            linear_stretch="never",
+        )
+        elapsed = time.perf_counter() - started
+        output_path = OUTPUT_DIR / f"{key}_denoised.png"
+        save_rgb(output_path, denoised)
+        denoised = load_rgb(output_path)
+        make_comparison(
+            COMPARISON_DIR / f"{key}_comparison.jpg",
+            noisy,
+            denoised,
+            reference,
+            metadata["title"],
+        )
+
+        result = {
+            "id": key,
+            **metadata,
+            "resolution": "1920x1080",
+            "device": str(device),
+            "elapsed_seconds": round(elapsed, 3),
+            "noisy_psnr_db": round(psnr(reference, noisy), 3),
+            "denoised_psnr_db": round(psnr(reference, denoised), 3),
+            "noisy_ssim": round(ssim(reference, noisy), 5),
+            "denoised_ssim": round(ssim(reference, denoised), 5),
+            "noisy_mae": round(mae(reference, noisy), 6),
+            "denoised_mae": round(mae(reference, denoised), 6),
+            "input": str(noisy_path.relative_to(ROOT)),
+            "output": str(output_path.relative_to(ROOT)),
+        }
+        results.append(result)
+        print(json.dumps(result, ensure_ascii=False))
+
+    (OUTPUT_DIR / "benchmark_results.json").write_text(
+        json.dumps(results, indent=2, ensure_ascii=False), encoding="utf-8"
+    )
+
+
+if __name__ == "__main__":
+    main()

StarDenoise/denoise.py

@@ -0,0 +1,262 @@
+"""Command-line inference demo for Cosmic Clarity AI 3.6."""
+
+from __future__ import annotations
+
+import argparse
+import contextlib
+import math
+import time
+from pathlib import Path
+from typing import Iterator
+
+import cv2
+import numpy as np
+import torch
+
+from image_io import load_image, save_image
+from model import load_model, model_statistics
+
+
+ROOT = Path(__file__).resolve().parent
+DEFAULT_WEIGHTS = ROOT / "models" / "deep_denoise_cnn_AI3_6.pth"
+
+
+def rgb_to_ycbcr(image: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    matrix = np.array(
+        [[0.299, 0.587, 0.114], [-0.168736, -0.331264, 0.5], [0.5, -0.418688, -0.081312]],
+        dtype=np.float32,
+    )
+    converted = image @ matrix.T
+    return converted[..., 0], converted[..., 1] + 0.5, converted[..., 2] + 0.5
+
+
+def ycbcr_to_rgb(y: np.ndarray, cb: np.ndarray, cr: np.ndarray) -> np.ndarray:
+    converted = np.stack((y, cb - 0.5, cr - 0.5), axis=-1)
+    matrix = np.array(
+        [[1.0, 0.0, 1.402], [1.0, -0.344136, -0.714136], [1.0, 1.772, 0.0]],
+        dtype=np.float32,
+    )
+    return np.clip(converted @ matrix.T, 0.0, 1.0)
+
+
+def _mtf_stretch(image: np.ndarray, target_median: float = 0.25) -> tuple[np.ndarray, float, list[float]]:
+    original_min = float(image.min())
+    stretched = image - original_min
+    medians = [float(np.median(stretched[..., channel])) for channel in range(3)]
+    for channel, median in enumerate(medians):
+        if median == 0:
+            continue
+        x = stretched[..., channel]
+        denominator = median * (target_median + x - 1.0) - target_median * x
+        stretched[..., channel] = ((median - 1.0) * target_median * x) / np.where(
+            np.abs(denominator) < 1e-12, 1e-12, denominator
+        )
+    return np.clip(stretched, 0.0, 1.0), original_min, medians
+
+
+def _mtf_unstretch(image: np.ndarray, original_min: float, original_medians: list[float]) -> np.ndarray:
+    result = image.copy()
+    current_medians = [float(np.median(result[..., channel])) for channel in range(3)]
+    for channel, current in enumerate(current_medians):
+        target = original_medians[channel]
+        if current == 0 or target == 0:
+            continue
+        x = result[..., channel]
+        denominator = current * (target + x - 1.0) - target * x
+        result[..., channel] = ((current - 1.0) * target * x) / np.where(
+            np.abs(denominator) < 1e-12, 1e-12, denominator
+        )
+    return np.clip(result + original_min, 0.0, 1.0)
+
+
+def _guided_filter(guide: np.ndarray, source: np.ndarray, radius: int, epsilon: float) -> np.ndarray:
+    kernel = (2 * radius + 1, 2 * radius + 1)
+    kwargs = {"ddepth": -1, "ksize": kernel, "borderType": cv2.BORDER_REFLECT}
+    mean_i = cv2.boxFilter(guide, **kwargs)
+    mean_p = cv2.boxFilter(source, **kwargs)
+    corr_i = cv2.boxFilter(guide * guide, **kwargs)
+    corr_ip = cv2.boxFilter(guide * source, **kwargs)
+    variance_i = corr_i - mean_i * mean_i
+    covariance_ip = corr_ip - mean_i * mean_p
+    a = covariance_ip / (variance_i + epsilon)
+    b = mean_p - a * mean_i
+    return cv2.boxFilter(a, **kwargs) * guide + cv2.boxFilter(b, **kwargs)
+
+
+def denoise_chroma(y: np.ndarray, cb: np.ndarray, cr: np.ndarray, strength: float) -> tuple[np.ndarray, np.ndarray]:
+    effective = float(np.clip(strength * 2.0, 0.0, 1.0))
+    if effective == 0:
+        return cb, cr
+    radius = 2 + int(round(10 * effective))
+    epsilon = (0.001 + 0.05 * effective) ** 2
+    cb_filtered = _guided_filter(y, cb, radius, epsilon)
+    cr_filtered = _guided_filter(y, cr, radius, epsilon)
+    return (
+        (1.0 - effective) * cb + effective * cb_filtered,
+        (1.0 - effective) * cr + effective * cr_filtered,
+    )
+
+
+def _positions(length: int, tile: int, overlap: int) -> list[int]:
+    if length <= tile:
+        return [0]
+    stride = tile - overlap
+    result = list(range(0, length - tile + 1, stride))
+    if result[-1] != length - tile:
+        result.append(length - tile)
+    return result
+
+
+def _batches(items: list[tuple[int, int, np.ndarray]], size: int) -> Iterator[list[tuple[int, int, np.ndarray]]]:
+    for start in range(0, len(items), size):
+        yield items[start : start + size]
+
+
+def infer_channel(
+    channel: np.ndarray,
+    model: torch.nn.Module,
+    device: torch.device,
+    tile_size: int,
+    overlap: int,
+    batch_size: int,
+    half_precision: bool,
+) -> np.ndarray:
+    if overlap >= tile_size:
+        raise ValueError("overlap must be smaller than tile-size")
+    height, width = channel.shape
+    pad_y, pad_x = max(0, tile_size - height), max(0, tile_size - width)
+    padded = np.pad(channel, ((0, pad_y), (0, pad_x)), mode="reflect")
+    ph, pw = padded.shape
+    tiles = [
+        (y, x, padded[y : y + tile_size, x : x + tile_size])
+        for y in _positions(ph, tile_size, overlap)
+        for x in _positions(pw, tile_size, overlap)
+    ]
+    # A Hann window suppresses tile seams; the floor guarantees edge coverage.
+    window_1d = np.hanning(tile_size).astype(np.float32)
+    window = np.maximum(np.outer(window_1d, window_1d), 1e-3)
+    accumulated = np.zeros_like(padded, dtype=np.float32)
+    weights = np.zeros_like(padded, dtype=np.float32)
+
+    use_amp = half_precision and device.type == "cuda"
+    for batch in _batches(tiles, batch_size):
+        array = np.stack([item[2] for item in batch])
+        tensor = torch.from_numpy(array[:, None]).repeat(1, 3, 1, 1).to(device)
+        amp = torch.autocast(device_type="cuda", dtype=torch.float16) if use_amp else contextlib.nullcontext()
+        with torch.inference_mode(), amp:
+            prediction = model(tensor)[:, 0].float().cpu().numpy()
+        for (y, x, _), tile_prediction in zip(batch, prediction):
+            accumulated[y : y + tile_size, x : x + tile_size] += tile_prediction * window
+            weights[y : y + tile_size, x : x + tile_size] += window
+    return (accumulated / np.maximum(weights, 1e-8))[:height, :width]
+
+
+def denoise_array(
+    pixels: np.ndarray,
+    model: torch.nn.Module,
+    device: torch.device,
+    strength: float = 0.9,
+    color_strength: float | None = None,
+    mode: str = "full",
+    tile_size: int = 256,
+    overlap: int = 64,
+    batch_size: int = 1,
+    half_precision: bool = True,
+    linear_stretch: str = "auto",
+) -> np.ndarray:
+    if not 0.0 <= strength <= 1.0:
+        raise ValueError("strength must be in [0, 1]")
+    if strength == 0.0:
+        return pixels.copy()
+    mono = pixels.shape[-1] == 1
+    rgb = np.repeat(pixels, 3, axis=-1) if mono else pixels[..., :3].copy()
+
+    should_stretch = linear_stretch == "always" or (
+        linear_stretch == "auto" and float(np.median(rgb - rgb.min())) < 0.05
+    )
+    stretch_state = None
+    if should_stretch:
+        rgb, minimum, medians = _mtf_stretch(rgb)
+        stretch_state = (minimum, medians)
+
+    def run(channel: np.ndarray) -> np.ndarray:
+        return infer_channel(channel, model, device, tile_size, overlap, batch_size, half_precision)
+
+    if mono:
+        denoised = ((1.0 - strength) * rgb[..., 0] + strength * run(rgb[..., 0]))[..., None]
+        result = np.repeat(denoised, 3, axis=-1)
+    elif mode == "separate":
+        result = np.stack(
+            [(1.0 - strength) * rgb[..., c] + strength * run(rgb[..., c]) for c in range(3)],
+            axis=-1,
+        )
+    elif mode in ("luminance", "full"):
+        y, cb, cr = rgb_to_ycbcr(rgb)
+        y = (1.0 - strength) * y + strength * run(y)
+        if mode == "full":
+            cb, cr = denoise_chroma(y, cb, cr, strength if color_strength is None else color_strength)
+        result = ycbcr_to_rgb(y, cb, cr)
+    else:
+        raise ValueError(f"Unknown mode: {mode}")
+
+    if stretch_state is not None:
+        result = _mtf_unstretch(result, *stretch_state)
+    return result[..., :1] if mono else result
+
+
+def choose_device(requested: str) -> torch.device:
+    if requested == "auto":
+        return torch.device("cuda" if torch.cuda.is_available() else "cpu")
+    device = torch.device(requested)
+    if device.type == "cuda" and not torch.cuda.is_available():
+        raise RuntimeError("CUDA was requested but is not available")
+    return device
+
+
+def build_parser() -> argparse.ArgumentParser:
+    parser = argparse.ArgumentParser(description="Astronomy-specific Cosmic Clarity AI 3.6 denoising demo")
+    parser.add_argument("input", type=Path, help="Input JPG/PNG/TIFF/FITS")
+    parser.add_argument("output", type=Path, help="Output path; extension selects format")
+    parser.add_argument("--weights", type=Path, default=DEFAULT_WEIGHTS)
+    parser.add_argument("--device", default="auto", help="auto, cpu, cuda, or cuda:N")
+    parser.add_argument("--strength", type=float, default=0.9)
+    parser.add_argument("--color-strength", type=float, default=None)
+    parser.add_argument("--mode", choices=("luminance", "full", "separate"), default="full")
+    parser.add_argument("--tile-size", type=int, default=256)
+    parser.add_argument("--overlap", type=int, default=64)
+    parser.add_argument("--batch-size", type=int, default=1)
+    parser.add_argument("--precision", choices=("fp16", "fp32"), default="fp16")
+    parser.add_argument("--linear-stretch", choices=("auto", "always", "never"), default="auto")
+    return parser
+
+
+def main() -> None:
+    args = build_parser().parse_args()
+    device = choose_device(args.device)
+    model = load_model(args.weights, device)
+    record = load_image(args.input)
+    started = time.perf_counter()
+    output = denoise_array(
+        record.pixels,
+        model,
+        device,
+        strength=args.strength,
+        color_strength=args.color_strength,
+        mode=args.mode,
+        tile_size=args.tile_size,
+        overlap=args.overlap,
+        batch_size=args.batch_size,
+        half_precision=args.precision == "fp16",
+        linear_stretch=args.linear_stretch,
+    )
+    save_image(args.output, record, output)
+    elapsed = time.perf_counter() - started
+    statistics = model_statistics(model)
+    print(f"Model: Cosmic Clarity AI 3.6 ({statistics['parameters']:,} parameters)")
+    print(f"Device: {device}; input: {record.pixels.shape}; elapsed: {elapsed:.3f}s")
+    print(f"Saved: {args.output.resolve()}")
+
+
+if __name__ == "__main__":
+    main()
+

StarDenoise/download_assets.py

@@ -0,0 +1,53 @@
+"""Re-download official model/sample assets with integrity verification."""
+
+from __future__ import annotations
+
+import hashlib
+import shutil
+import urllib.request
+from pathlib import Path
+
+
+ROOT = Path(__file__).resolve().parent
+ASSETS = {
+    ROOT / "models" / "deep_denoise_cnn_AI3_6.pth": (
+        "https://github.com/setiastro/cosmicclarity/releases/download/Linux/deep_denoise_cnn_AI3_6.pth",
+        "f860e234dc6a3ab1acacbea399c235179c86cb8b287f2ab9010bc18c8af08671",
+    ),
+    ROOT / "samples" / "M51_short.jpg": (
+        "https://raw.githubusercontent.com/Steffenhir/AstroNoiseNet/main/examples/M51_short.jpg",
+        None,
+    ),
+    ROOT / "samples" / "M51_long_reference.jpg": (
+        "https://raw.githubusercontent.com/Steffenhir/AstroNoiseNet/main/examples/M51_long.jpg",
+        None,
+    ),
+}
+
+
+def sha256(path: Path) -> str:
+    digest = hashlib.sha256()
+    with path.open("rb") as stream:
+        for block in iter(lambda: stream.read(1024 * 1024), b""):
+            digest.update(block)
+    return digest.hexdigest()
+
+
+def main() -> None:
+    for destination, (url, expected_hash) in ASSETS.items():
+        destination.parent.mkdir(parents=True, exist_ok=True)
+        temporary = destination.with_suffix(destination.suffix + ".download")
+        print(f"Downloading {url}")
+        with urllib.request.urlopen(url) as response, temporary.open("wb") as output:
+            shutil.copyfileobj(response, output)
+        actual_hash = sha256(temporary)
+        if expected_hash and actual_hash != expected_hash:
+            temporary.unlink(missing_ok=True)
+            raise RuntimeError(f"SHA-256 mismatch for {destination.name}: {actual_hash}")
+        temporary.replace(destination)
+        print(f"Saved {destination} (sha256={actual_hash})")
+
+
+if __name__ == "__main__":
+    main()
+

StarDenoise/evaluate.py

@@ -0,0 +1,37 @@
+"""Simple full-reference metrics for the bundled short/long M51 pair."""
+
+from __future__ import annotations
+
+import argparse
+import math
+from pathlib import Path
+
+import numpy as np
+from PIL import Image
+
+
+def load_rgb(path: Path) -> np.ndarray:
+    with Image.open(path) as image:
+        return np.asarray(image.convert("RGB"), dtype=np.float32) / 255.0
+
+
+def psnr(reference: np.ndarray, candidate: np.ndarray) -> float:
+    mse = float(np.mean((reference - candidate) ** 2))
+    return math.inf if mse == 0 else 10.0 * math.log10(1.0 / mse)
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser()
+    parser.add_argument("reference", type=Path)
+    parser.add_argument("candidate", type=Path)
+    args = parser.parse_args()
+    reference, candidate = load_rgb(args.reference), load_rgb(args.candidate)
+    if reference.shape != candidate.shape:
+        raise ValueError(f"Shape mismatch: {reference.shape} vs {candidate.shape}")
+    print(f"PSNR: {psnr(reference, candidate):.3f} dB")
+    print(f"MAE:  {np.mean(np.abs(reference - candidate)):.6f}")
+
+
+if __name__ == "__main__":
+    main()
+

StarDenoise/image_io.py

@@ -0,0 +1,121 @@
+"""Image I/O with FITS metadata and dynamic-range preservation."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+import tifffile
+from astropy.io import fits
+from PIL import Image
+
+
+FITS_EXTENSIONS = {".fit", ".fits", ".fts"}
+TIFF_EXTENSIONS = {".tif", ".tiff"}
+
+
+@dataclass
+class ImageRecord:
+    pixels: np.ndarray  # HWC float32, normalized to [0, 1]
+    offset: float
+    scale: float
+    original_dtype: np.dtype
+    original_layout: str
+    header: Any = None
+
+
+def _to_hwc(array: np.ndarray) -> tuple[np.ndarray, str]:
+    if array.ndim == 2:
+        return array[..., None], "hw"
+    if array.ndim != 3:
+        raise ValueError(f"Only mono/RGB 2-D images are supported, got shape {array.shape}")
+    if array.shape[-1] in (1, 3, 4):
+        return array[..., :3], "hwc"
+    if array.shape[0] in (1, 3, 4):
+        return np.moveaxis(array[:3], 0, -1), "chw"
+    raise ValueError(f"Cannot determine channel axis for shape {array.shape}")
+
+
+def _normalize(array: np.ndarray) -> tuple[np.ndarray, float, float]:
+    finite = np.isfinite(array)
+    if not finite.any():
+        raise ValueError("Image contains no finite pixels")
+    clean = np.nan_to_num(array, copy=True)
+    if np.issubdtype(array.dtype, np.integer):
+        info = np.iinfo(array.dtype)
+        offset = float(max(0, info.min))
+        scale = float(info.max - offset)
+    else:
+        minimum = float(np.min(clean[finite]))
+        maximum = float(np.max(clean[finite]))
+        if minimum >= 0.0 and maximum <= 1.0:
+            offset, scale = 0.0, 1.0
+        else:
+            offset, scale = minimum, maximum - minimum
+    if scale <= 0:
+        scale = 1.0
+    normalized = np.clip((clean.astype(np.float32) - offset) / scale, 0.0, 1.0)
+    return normalized, offset, scale
+
+
+def load_image(path: str | Path) -> ImageRecord:
+    path = Path(path)
+    extension = path.suffix.lower()
+    header = None
+    if extension in FITS_EXTENSIONS:
+        with fits.open(path, memmap=False) as hdul:
+            index = next((i for i, hdu in enumerate(hdul) if hdu.data is not None), None)
+            if index is None:
+                raise ValueError(f"No image HDU in {path}")
+            array = np.asarray(hdul[index].data)
+            header = hdul[index].header.copy()
+    elif extension in TIFF_EXTENSIONS:
+        array = np.asarray(tifffile.imread(path))
+    else:
+        with Image.open(path) as image:
+            if image.mode not in ("L", "I;16", "I", "F"):
+                image = image.convert("RGB")
+            array = np.asarray(image)
+
+    original_dtype = array.dtype
+    hwc, layout = _to_hwc(array)
+    pixels, offset, scale = _normalize(hwc)
+    return ImageRecord(pixels, offset, scale, original_dtype, layout, header)
+
+
+def _restore_layout(array: np.ndarray, layout: str) -> np.ndarray:
+    if layout == "hw":
+        return array[..., 0]
+    if layout == "chw":
+        return np.moveaxis(array, -1, 0)
+    return array
+
+
+def _restore_range(record: ImageRecord, pixels: np.ndarray) -> np.ndarray:
+    restored = np.clip(pixels, 0.0, 1.0) * record.scale + record.offset
+    if np.issubdtype(record.original_dtype, np.integer):
+        info = np.iinfo(record.original_dtype)
+        restored = np.clip(np.rint(restored), info.min, info.max).astype(record.original_dtype)
+    else:
+        restored = restored.astype(record.original_dtype)
+    return _restore_layout(restored, record.original_layout)
+
+
+def save_image(path: str | Path, record: ImageRecord, pixels: np.ndarray) -> None:
+    path = Path(path)
+    path.parent.mkdir(parents=True, exist_ok=True)
+    extension = path.suffix.lower()
+    restored = _restore_range(record, pixels)
+    if extension in FITS_EXTENSIONS:
+        fits.PrimaryHDU(restored, header=record.header).writeto(path, overwrite=True)
+    elif extension in TIFF_EXTENSIONS:
+        tifffile.imwrite(path, restored)
+    else:
+        # JPEG/PNG demo output is display-oriented. Preserve mono when applicable.
+        display = np.clip(pixels * 255.0 + 0.5, 0, 255).astype(np.uint8)
+        if record.original_layout == "hw":
+            display = display[..., 0]
+        Image.fromarray(display).save(path, quality=95)
+

StarDenoise/model.py

@@ -0,0 +1,117 @@
+"""PyTorch reproduction of the SetiAstro Cosmic Clarity AI 3.6 denoiser.
+
+The architecture follows SetiAstroCosmicClarity_denoise.py from the upstream
+MIT-licensed project.  It is intentionally kept independent of the original
+GUI so it can be embedded in a larger image-processing pipeline.
+"""
+
+from __future__ import annotations
+
+from pathlib import Path
+from typing import Any
+
+import torch
+from torch import nn
+
+
+class ResidualBlock(nn.Module):
+    def __init__(self, channels: int) -> None:
+        super().__init__()
+        self.conv1 = nn.Conv2d(channels, channels, 3, padding=1)
+        self.relu = nn.ReLU(inplace=False)
+        self.conv2 = nn.Conv2d(channels, channels, 3, padding=1)
+
+    def forward(self, x: torch.Tensor) -> torch.Tensor:
+        return self.relu(self.conv2(self.relu(self.conv1(x))) + x)
+
+
+class CosmicClarityDenoiseCNN(nn.Module):
+    """Full-resolution residual encoder/decoder used by AI 3.6.
+
+    Despite the encoder/decoder names, the network does not spatially
+    downsample. Dilated convolutions and skip concatenations preserve faint
+    stars and small-scale structure while increasing the receptive field.
+    """
+
+    def __init__(self) -> None:
+        super().__init__()
+        self.encoder1 = self._encoder(3, 16)
+        self.encoder2 = self._encoder(16, 32)
+        self.encoder3 = self._encoder(32, 64, dilation=2)
+        self.encoder4 = self._encoder(64, 128)
+        self.encoder5 = self._encoder(128, 256, dilation=2)
+
+        self.decoder5 = self._decoder(256 + 128, 128)
+        self.decoder4 = self._decoder(128 + 64, 64)
+        self.decoder3 = self._decoder(64 + 32, 32)
+        self.decoder2 = self._decoder(32 + 16, 16)
+        self.decoder1 = nn.Sequential(nn.Conv2d(16, 3, 3, padding=1), nn.Sigmoid())
+
+    @staticmethod
+    def _encoder(in_channels: int, out_channels: int, dilation: int = 1) -> nn.Sequential:
+        return nn.Sequential(
+            nn.Conv2d(
+                in_channels,
+                out_channels,
+                3,
+                padding=dilation,
+                dilation=dilation,
+            ),
+            nn.ReLU(inplace=False),
+            ResidualBlock(out_channels),
+        )
+
+    @staticmethod
+    def _decoder(in_channels: int, out_channels: int) -> nn.Sequential:
+        return nn.Sequential(
+            nn.Conv2d(in_channels, out_channels, 3, padding=1),
+            nn.ReLU(inplace=False),
+            ResidualBlock(out_channels),
+        )
+
+    def forward(self, x: torch.Tensor) -> torch.Tensor:
+        e1 = self.encoder1(x)
+        e2 = self.encoder2(e1)
+        e3 = self.encoder3(e2)
+        e4 = self.encoder4(e3)
+        e5 = self.encoder5(e4)
+        d5 = self.decoder5(torch.cat((e5, e4), dim=1))
+        d4 = self.decoder4(torch.cat((d5, e3), dim=1))
+        d3 = self.decoder3(torch.cat((d4, e2), dim=1))
+        d2 = self.decoder2(torch.cat((d3, e1), dim=1))
+        return self.decoder1(d2)
+
+
+def _safe_torch_load(path: Path, device: torch.device) -> Any:
+    """Load a trusted upstream checkpoint while preferring restricted unpickling."""
+    try:
+        return torch.load(path, map_location=device, weights_only=True)
+    except TypeError:  # torch < 2.0 compatibility
+        return torch.load(path, map_location=device)
+
+
+def load_model(weights: str | Path, device: torch.device) -> CosmicClarityDenoiseCNN:
+    weights_path = Path(weights)
+    if not weights_path.is_file():
+        raise FileNotFoundError(
+            f"Model weights not found: {weights_path}. Run download_assets.py first."
+        )
+    checkpoint = _safe_torch_load(weights_path, device)
+    if not isinstance(checkpoint, dict):
+        raise ValueError("Unexpected checkpoint format: expected a mapping")
+    state_dict = checkpoint.get("model_state_dict", checkpoint)
+    model = CosmicClarityDenoiseCNN().to(device)
+    model.load_state_dict(state_dict, strict=True)
+    model.eval()
+    return model
+
+
+def model_statistics(model: nn.Module) -> dict[str, float | int]:
+    parameters = sum(parameter.numel() for parameter in model.parameters())
+    trainable = sum(parameter.numel() for parameter in model.parameters() if parameter.requires_grad)
+    return {
+        "parameters": parameters,
+        "trainable_parameters": trainable,
+        "fp32_size_mib": parameters * 4 / 1024**2,
+    }
+

StarDenoise/requirements.txt

@@ -0,0 +1,7 @@
+numpy>=1.24
+Pillow>=10.0
+astropy>=5.3
+tifffile>=2023.7
+opencv-python>=4.8
+torch>=2.1
+

StarDenoise/tests/__init__.py

@@ -0,0 +1 @@
+

StarDenoise/tests/test_core.py

@@ -0,0 +1,54 @@
+from __future__ import annotations
+
+import sys
+import tempfile
+import unittest
+from pathlib import Path
+
+import numpy as np
+import torch
+
+
+ROOT = Path(__file__).resolve().parents[1]
+sys.path.insert(0, str(ROOT))
+
+from denoise import denoise_array, rgb_to_ycbcr, ycbcr_to_rgb  # noqa: E402
+from image_io import load_image, save_image  # noqa: E402
+from model import CosmicClarityDenoiseCNN, model_statistics  # noqa: E402
+
+
+class CoreTests(unittest.TestCase):
+    def test_model_shape_and_parameter_count(self) -> None:
+        model = CosmicClarityDenoiseCNN()
+        output = model(torch.rand(1, 3, 16, 16))
+        self.assertEqual(tuple(output.shape), (1, 3, 16, 16))
+        self.assertEqual(model_statistics(model)["parameters"], 2_945_283)
+
+    def test_ycbcr_roundtrip(self) -> None:
+        image = np.random.default_rng(3).random((20, 21, 3), dtype=np.float32)
+        reconstructed = ycbcr_to_rgb(*rgb_to_ycbcr(image))
+        np.testing.assert_allclose(image, reconstructed, atol=3e-5)
+
+    def test_zero_strength_is_identity(self) -> None:
+        image = np.random.default_rng(4).random((9, 10, 3), dtype=np.float32)
+        result = denoise_array(image, torch.nn.Identity(), torch.device("cpu"), strength=0.0)
+        np.testing.assert_array_equal(image, result)
+
+    def test_fits_roundtrip_preserves_channel_layout(self) -> None:
+        from astropy.io import fits
+
+        original = np.arange(3 * 8 * 9, dtype=np.uint16).reshape(3, 8, 9)
+        with tempfile.TemporaryDirectory() as directory:
+            source = Path(directory) / "source.fits"
+            output = Path(directory) / "output.fits"
+            fits.PrimaryHDU(original).writeto(source)
+            record = load_image(source)
+            save_image(output, record, record.pixels)
+            with fits.open(output) as hdul:
+                restored = hdul[0].data
+        self.assertEqual(restored.shape, original.shape)
+        np.testing.assert_array_equal(restored, original)
+
+
+if __name__ == "__main__":
+    unittest.main()

StarDenoise/third_party/AstroNoiseNet_LICENSE.txt

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2022 Steffenhir
+
+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.

StarDenoise/third_party/CosmicClarity_LICENSE.txt

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2024 Franklin marek
+
+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.

StarDenoise/third_party/upstream_denoise.py

@@ -0,0 +1,1225 @@
+import warnings
+from xisf import XISF
+import os
+import sys
+import torch
+import numpy as np
+import torch.nn as nn
+import lz4.block
+import zstandard
+import base64
+import ast
+import platform
+import tifffile as tiff
+from astropy.io import fits
+from PIL import Image
+import tkinter as tk
+from tkinter import simpledialog, messagebox
+from tkinter import ttk
+import argparse  # For command-line argument parsing
+import onnxruntime as ort
+import cv2 
+
+sys.stdout.reconfigure(encoding='utf-8')
+
+
+
+
+#torch.cuda.is_available = lambda: False
+
+# Suppress model loading warnings
+warnings.filterwarnings("ignore")
+
+# Detect broken fp16/mixed-precision support
+def has_broken_fp16():
+    try:
+        cc = torch.cuda.get_device_capability()
+        return cc[0] < 8  # For example, <8 means <Turing (RTX 20xx)
+    except:
+        return True
+
+DISABLE_MIXED_PRECISION = has_broken_fp16()
+if DISABLE_MIXED_PRECISION:
+    print("[INFO] Mixed precision disabled due to unsupported GPU capability.")
+
+# Define the ResidualBlock
+class ResidualBlock(nn.Module):
+    def __init__(self, channels):
+        super(ResidualBlock, self).__init__()
+        self.conv1 = nn.Conv2d(channels, channels, kernel_size=3, padding=1)
+        self.relu = nn.ReLU()
+        self.conv2 = nn.Conv2d(channels, channels, kernel_size=3, padding=1)
+    
+    def forward(self, x):
+        residual = x
+        out = self.conv1(x)
+        out = self.relu(out)
+        out = self.conv2(out)
+        out += residual
+        out = self.relu(out)
+        return out
+
+# Updated DenoiseCNN with Residual Blocks
+class DenoiseCNN(nn.Module):
+    def __init__(self):
+        super(DenoiseCNN, self).__init__()
+        
+        # Encoder (down-sampling path)
+        self.encoder1 = nn.Sequential(
+            nn.Conv2d(3, 16, kernel_size=3, padding=1),  # 1st layer (3 -> 16 feature maps)
+            nn.ReLU(),
+            ResidualBlock(16)  # Replaced Conv2d + ReLU with ResidualBlock
+        )
+        self.encoder2 = nn.Sequential(
+            nn.Conv2d(16, 32, kernel_size=3, padding=1),  # 2nd layer (16 -> 32 feature maps)
+            nn.ReLU(),
+            ResidualBlock(32)  # Replaced Conv2d + ReLU with ResidualBlock
+        )
+        self.encoder3 = nn.Sequential(
+            nn.Conv2d(32, 64, kernel_size=3, padding=2, dilation=2),  # 3rd layer (32 -> 64) with dilation
+            nn.ReLU(),
+            ResidualBlock(64)  # Replaced Conv2d + ReLU with ResidualBlock
+        )
+        self.encoder4 = nn.Sequential(
+            nn.Conv2d(64, 128, kernel_size=3, padding=1),  # 4th layer (64 -> 128 feature maps)
+            nn.ReLU(),
+            ResidualBlock(128)  # Replaced Conv2d + ReLU with ResidualBlock
+        )
+        self.encoder5 = nn.Sequential(
+            nn.Conv2d(128, 256, kernel_size=3, padding=2, dilation=2),  # 5th layer (128 -> 256) with dilation
+            nn.ReLU(),
+            ResidualBlock(256)  # Replaced Conv2d + ReLU with ResidualBlock
+        )
+        
+        # Decoder (up-sampling path with skip connections)
+        self.decoder5 = nn.Sequential(
+            nn.Conv2d(256 + 128, 128, kernel_size=3, padding=1),  # 256 + 128 feature maps from encoder4
+            nn.ReLU(),
+            ResidualBlock(128)  # Replaced Conv2d + ReLU with ResidualBlock
+        )
+        self.decoder4 = nn.Sequential(
+            nn.Conv2d(128 + 64, 64, kernel_size=3, padding=1),  # 128 + 64 feature maps from encoder3
+            nn.ReLU(),
+            ResidualBlock(64)  # Replaced Conv2d + ReLU with ResidualBlock
+        )
+        self.decoder3 = nn.Sequential(
+            nn.Conv2d(64 + 32, 32, kernel_size=3, padding=1),  # 64 + 32 feature maps from encoder2
+            nn.ReLU(),
+            ResidualBlock(32)  # Replaced Conv2d + ReLU with ResidualBlock
+        )
+        self.decoder2 = nn.Sequential(
+            nn.Conv2d(32 + 16, 16, kernel_size=3, padding=1),  # 32 + 16 feature maps from encoder1
+            nn.ReLU(),
+            ResidualBlock(16)  # Replaced Conv2d + ReLU with ResidualBlock
+        )
+        self.decoder1 = nn.Sequential(
+            nn.Conv2d(16, 3, kernel_size=3, padding=1),  # Output layer (16 -> 3 channels for RGB output)
+            nn.Sigmoid()  # Ensure output values are between 0 and 1
+        )
+
+    def forward(self, x):
+        # Encoder
+        e1 = self.encoder1(x)  # First encoding block
+        e2 = self.encoder2(e1)  # Second encoding block
+        e3 = self.encoder3(e2)  # Third encoding block
+        e4 = self.encoder4(e3)  # Fourth encoding block
+        e5 = self.encoder5(e4)  # Fifth encoding block
+        
+        # Decoder with skip connections
+        d5 = self.decoder5(torch.cat([e5, e4], dim=1))  # Concatenate with encoder4 output
+        d4 = self.decoder4(torch.cat([d5, e3], dim=1))  # Concatenate with encoder3 output
+        d3 = self.decoder3(torch.cat([d4, e2], dim=1))  # Concatenate with encoder2 output
+        d2 = self.decoder2(torch.cat([d3, e1], dim=1))  # Concatenate with encoder1 output
+        d1 = self.decoder1(d2)  # Final output layer
+
+        return d1
+
+
+# Get the directory of the executable or the script location
+# Get the directory of the executable or the script location
+exe_dir = os.path.dirname(sys.executable) if getattr(sys, 'frozen', False) else os.path.dirname(__file__)
+
+cached_models = None  # Cache to avoid reloading models unnecessarily
+
+def load_models(exe_dir, use_gpu=True):
+    """
+    Load only the mono denoise model.
+    - PyTorch (CUDA/CPU) if available
+    - Otherwise ONNX (DirectML) if available and GPU requested
+    """
+    global cached_models
+    if cached_models:
+        return cached_models
+
+    # Decide backend (PyTorch preferred; ONNX/DML fallback)
+    if torch.cuda.is_available() and use_gpu:
+        device = torch.device("cuda"); is_onnx = False
+    elif "DmlExecutionProvider" in ort.get_available_providers() and use_gpu:
+        device = "DirectML"; is_onnx = True
+    else:
+        device = torch.device("cpu"); is_onnx = False
+
+    # --- Mono model only ---
+    if not is_onnx:
+        mono = DenoiseCNN().to(device)
+        ckpt = torch.load(os.path.join(exe_dir, "deep_denoise_cnn_AI3_6.pth"), map_location=device)
+        mono.load_state_dict(ckpt.get("model_state_dict", ckpt))
+        mono.eval()
+        mono_model = mono
+    else:
+        mono_model = ort.InferenceSession(
+            os.path.join(exe_dir, "deep_denoise_cnn_AI3_6.onnx"),
+            providers=["DmlExecutionProvider"]
+        )
+
+    cached_models = {
+        "device":     device,
+        "is_onnx":    is_onnx,
+        "mono_model": mono_model
+    }
+
+    backend_str = "ONNX/DML" if is_onnx else (device.type if isinstance(device, torch.device) else str(device))
+    print(f"Loaded mono model → backend={backend_str}")
+    return cached_models
+
+
+
+
+# Function to extract luminance (Y channel) directly using a matrix for 32-bit float precision
+def extract_luminance(image):
+    # Ensure the image is a NumPy array in 32-bit float format
+    if isinstance(image, Image.Image):
+        image = np.array(image).astype(np.float32) / 255.0  # Ensure it's in 32-bit float format
+
+    # Check if the image is grayscale (single channel), and convert it to 3-channel RGB if so
+    if image.shape[-1] == 1:
+        image = np.concatenate([image] * 3, axis=-1)  # Duplicate the single channel to create RGB
+
+    # Conversion matrix for RGB to YCbCr (ITU-R BT.601 standard)
+    conversion_matrix = np.array([[0.299, 0.587, 0.114],
+                                  [-0.168736, -0.331264, 0.5],
+                                  [0.5, -0.418688, -0.081312]], dtype=np.float32)
+
+    # Apply the RGB to YCbCr conversion matrix
+    ycbcr_image = np.dot(image, conversion_matrix.T)
+
+    # Split the channels: Y is luminance, Cb and Cr are chrominance
+    y_channel = ycbcr_image[:, :, 0]
+    cb_channel = ycbcr_image[:, :, 1] + 0.5  # Offset back to [0, 1] range
+    cr_channel = ycbcr_image[:, :, 2] + 0.5  # Offset back to [0, 1] range
+
+    # Return the Y (luminance) channel in 32-bit float and the Cb, Cr channels for later use
+    return y_channel, cb_channel, cr_channel
+
+def merge_luminance(y_channel, cb_channel, cr_channel):
+    # Ensure all channels are 32-bit float and in the range [0, 1]
+    y_channel = np.clip(y_channel, 0, 1).astype(np.float32)
+    cb_channel = np.clip(cb_channel, 0, 1).astype(np.float32) - 0.5  # Adjust Cb to [-0.5, 0.5]
+    cr_channel = np.clip(cr_channel, 0, 1).astype(np.float32) - 0.5  # Adjust Cr to [-0.5, 0.5]
+
+    # Convert YCbCr to RGB in 32-bit float format
+    rgb_image = ycbcr_to_rgb(y_channel, cb_channel, cr_channel)
+
+    return rgb_image
+
+
+# Helper function to convert YCbCr (32-bit float) to RGB
+def ycbcr_to_rgb(y_channel, cb_channel, cr_channel):
+    # Combine Y, Cb, Cr channels into one array for matrix multiplication
+    ycbcr_image = np.stack([y_channel, cb_channel, cr_channel], axis=-1)
+
+    # Conversion matrix for YCbCr to RGB (ITU-R BT.601 standard)
+    conversion_matrix = np.array([[1.0,  0.0, 1.402],
+                                  [1.0, -0.344136, -0.714136],
+                                  [1.0,  1.772, 0.0]], dtype=np.float32)
+
+    # Apply the YCbCr to RGB conversion matrix in 32-bit float
+    rgb_image = np.dot(ycbcr_image, conversion_matrix.T)
+
+    # Clip to ensure the RGB values stay within [0, 1] range
+    rgb_image = np.clip(rgb_image, 0.0, 1.0)
+
+    return rgb_image
+
+def _guided_filter(guide: np.ndarray, src: np.ndarray, radius: int, eps: float) -> np.ndarray:
+    """
+    Fast guided filter using boxFilter (edge-preserving, very fast).
+    guide and src are HxW float32 in [0,1].
+    radius is the neighborhood radius; ksize=(2*radius+1).
+    eps is the regularization term.
+    """
+    r = max(1, int(radius))
+    ksize = (2*r + 1, 2*r + 1)
+
+    mean_I  = cv2.boxFilter(guide, ddepth=-1, ksize=ksize, borderType=cv2.BORDER_REFLECT)
+    mean_p  = cv2.boxFilter(src,   ddepth=-1, ksize=ksize, borderType=cv2.BORDER_REFLECT)
+    mean_Ip = cv2.boxFilter(guide * src, ddepth=-1, ksize=ksize, borderType=cv2.BORDER_REFLECT)
+    cov_Ip  = mean_Ip - mean_I * mean_p
+
+    mean_II = cv2.boxFilter(guide * guide, ddepth=-1, ksize=ksize, borderType=cv2.BORDER_REFLECT)
+    var_I   = mean_II - mean_I * mean_I
+
+    a = cov_Ip / (var_I + eps)
+    b = mean_p - a * mean_I
+
+    mean_a = cv2.boxFilter(a, ddepth=-1, ksize=ksize, borderType=cv2.BORDER_REFLECT)
+    mean_b = cv2.boxFilter(b, ddepth=-1, ksize=ksize, borderType=cv2.BORDER_REFLECT)
+
+    q = mean_a * guide + mean_b
+    return q
+
+
+
+def denoise_chroma(cb: np.ndarray,
+                   cr: np.ndarray,
+                   strength: float,
+                   method: str = "guided",
+                   strength_scale: float = 2.0,
+                   guide_y: np.ndarray | None = None):
+    """
+    Fast chroma-only denoise for Cb/Cr in [0,1] float32.
+    method: 'guided' (default), 'gaussian', 'bilateral'
+    strength_scale: lets chroma smoothing go up to ~2× your slider.
+    guide_y: optional luminance guide (Y in [0,1]); required for 'guided' to be best.
+    """
+    eff = float(np.clip(strength * strength_scale, 0.0, 1.0))
+    if eff <= 0.0:
+        return cb, cr
+
+    cb = cb.astype(np.float32, copy=False)
+    cr = cr.astype(np.float32, copy=False)
+
+    if method == "guided":
+        # Need a guide; if not provided, fall back to Gaussian
+        if guide_y is not None:
+            # radius & eps scale with strength; tuned for strong chroma smoothing but edge-safe
+            radius = 2 + int(round(10 * eff))         # ~2..12  (ksize ~5..25)
+            eps    = (0.001 + 0.05 * eff) ** 2        # small regularization
+            cb_f   = _guided_filter(guide_y, cb, radius, eps)
+            cr_f   = _guided_filter(guide_y, cr, radius, eps)
+        else:
+            method = "gaussian"  # no guide provided → fast fallback
+
+    if method == "gaussian":
+        k     = 1 + 2 * int(round(8 * eff))           # 1,3,5,..,17
+        sigma = max(0.15, 2.4 * eff)
+        cb_f  = cv2.GaussianBlur(cb, (k, k), sigmaX=sigma, sigmaY=sigma, borderType=cv2.BORDER_REFLECT)
+        cr_f  = cv2.GaussianBlur(cr, (k, k), sigmaX=sigma, sigmaY=sigma, borderType=cv2.BORDER_REFLECT)
+
+    if method == "bilateral":
+        # Bilateral is decent but slower than Gaussian; guided is preferred for speed/quality.
+        d      = 5 + 2 * int(round(6 * eff))          # 5..17
+        sigmaC = 25.0 * (0.5 + 3.0 * eff)             # ~12.5..100
+        sigmaS = 3.0  * (0.5 + 6.0 * eff)             # ~1.5..21
+        cb_f = cv2.bilateralFilter(cb, d=d, sigmaColor=sigmaC, sigmaSpace=sigmaS)
+        cr_f = cv2.bilateralFilter(cr, d=d, sigmaColor=sigmaC, sigmaSpace=sigmaS)
+
+    # Blend (maskless)
+    w = eff
+    cb_out = (1.0 - w) * cb + w * cb_f
+    cr_out = (1.0 - w) * cr + w * cr_f
+    return cb_out, cr_out
+
+
+# Function to split an image into chunks with overlap
+def split_image_into_chunks_with_overlap(image, chunk_size, overlap):
+    height, width = image.shape[:2]
+    chunks = []
+    step_size = chunk_size - overlap  # Define how much to step over (overlapping area)
+
+    for i in range(0, height, step_size):
+        for j in range(0, width, step_size):
+            end_i = min(i + chunk_size, height)
+            end_j = min(j + chunk_size, width)
+            chunk = image[i:end_i, j:end_j]
+            chunks.append((chunk, i, j))  # Return chunk and its position
+    return chunks
+
+# Soft blending weights for the overlap area
+def generate_blend_weights(chunk_size, overlap):
+    ramp = np.linspace(0, 1, overlap)
+    flat = np.ones(chunk_size - 2 * overlap)
+    blend_vector = np.concatenate([ramp, flat, ramp[::-1]])  # Create smooth transition
+    blend_matrix = np.outer(blend_vector, blend_vector)  # 2D blending weights
+    return blend_matrix
+
+# Function to stitch overlapping chunks back together with soft blending while ignoring borders for all chunks
+def stitch_chunks_ignore_border(chunks, image_shape, chunk_size, overlap, border_size=16):
+    print(f"Number of chunks: {len(chunks)}")
+    print(f"First few chunks: {[type(c) if not isinstance(c, tuple) else len(c) for c in chunks[:3]]}")
+
+    stitched_image = np.zeros(image_shape, dtype=np.float32)
+    weight_map = np.zeros(image_shape, dtype=np.float32)  # Track blending weights
+    
+    for chunk, i, j in chunks:
+        actual_chunk_h, actual_chunk_w = chunk.shape[:2]
+
+        # Calculate the boundaries for the current chunk, ignoring the border
+        border_h = min(border_size, actual_chunk_h // 2)
+        border_w = min(border_size, actual_chunk_w // 2)
+
+        # Always ignore the border for all chunks, including edges
+        inner_chunk = chunk[border_h:actual_chunk_h-border_h, border_w:actual_chunk_w-border_w]
+        stitched_image[i+border_h:i+actual_chunk_h-border_h, j+border_w:j+actual_chunk_w-border_w] += inner_chunk
+        weight_map[i+border_h:i+actual_chunk_h-border_h, j+border_w:j+actual_chunk_w-border_w] += 1
+
+    # Normalize by weight_map to blend overlapping areas
+    stitched_image /= np.maximum(weight_map, 1)  # Avoid division by zero
+    return stitched_image
+
+# Function to blend two images (before and after)
+def blend_images(before, after, amount):
+    return (1 - amount) * before + amount * after
+
+# Function to get user input for denoise strength (Interactive)
+def get_user_input():
+    # Define global variables to store the user input
+    global use_gpu, denoise_strength, denoise_mode, color_denoise_strength
+
+    def on_submit():
+        # Update globals from UI
+        global use_gpu, denoise_strength, denoise_mode, separate_channels, color_denoise_strength
+        use_gpu = gpu_var.get() == "Yes"
+        denoise_strength = denoise_strength_slider.get()
+        # if the user didn't move the color slider, it still has a value; we take it as-is
+        color_denoise_strength = color_strength_slider.get()
+        denoise_mode = denoise_mode_var.get().lower()
+        separate_channels = sep_var.get()
+        root.quit()
+
+    root = tk.Tk()
+    root.title("Cosmic Clarity Denoise Tool")
+
+    # GPU selection
+    ttk.Label(root, text="Use GPU Acceleration:").pack(pady=5)
+    gpu_var = tk.StringVar(value="Yes")
+    ttk.OptionMenu(root, gpu_var, "Yes", "Yes", "No").pack()
+
+    # Luminance denoise strength
+    ttk.Label(root, text="Luminance Denoise Strength (0–1):").pack(pady=5)
+    denoise_strength_slider = tk.Scale(root, from_=0, to=1, resolution=0.01, orient="horizontal")
+    denoise_strength_slider.set(0.9)
+    denoise_strength_slider.pack()
+
+    # Color denoise strength (separate)
+    ttk.Label(root, text="Color (Chroma) Denoise Strength (0–1):").pack(pady=5)
+    color_strength_slider = tk.Scale(root, from_=0, to=1, resolution=0.01, orient="horizontal")
+    color_strength_slider.set(0.9)  # default to same as luminance
+    color_strength_slider.pack()
+
+    # Denoise mode
+    ttk.Label(root, text="Denoise Mode:").pack(pady=5)
+    denoise_mode_var = tk.StringVar(value="full")
+    ttk.OptionMenu(root, denoise_mode_var, "full", "luminance", "full", "separate").pack()
+
+    # Separate channels
+    ttk.Label(root, text="Process each RGB channel separately:").pack(pady=5)
+    sep_var = tk.BooleanVar(value=False)
+    ttk.Checkbutton(root, variable=sep_var, text="Separate channels").pack()
+
+    # Submit
+    ttk.Button(root, text="Submit", command=on_submit).pack(pady=20)
+
+    root.mainloop()
+    root.destroy()
+
+    return use_gpu, denoise_strength, denoise_mode, separate_channels, color_denoise_strength
+
+
+# Function to show progress during chunk processing
+def show_progress(current, total):
+    progress_percentage = (current / total) * 100
+    print(f"Progress: {progress_percentage:.2f}% ({current}/{total} chunks processed)", end='\r', flush=True)
+
+# Function to replace the 5-pixel border from the original image into the processed image
+def replace_border(original_image, processed_image, border_size=16):
+    # Ensure the dimensions of both images match
+    if original_image.shape != processed_image.shape:
+        raise ValueError("Original image and processed image must have the same dimensions.")
+    
+    # Replace the top border
+    processed_image[:border_size, :] = original_image[:border_size, :]
+    
+    # Replace the bottom border
+    processed_image[-border_size:, :] = original_image[-border_size:, :]
+    
+    # Replace the left border
+    processed_image[:, :border_size] = original_image[:, :border_size]
+    
+    # Replace the right border
+    processed_image[:, -border_size:] = original_image[:, -border_size:]
+    
+    return processed_image
+
+
+
+# Function to stretch an image
+def stretch_image(image):
+    """
+    Perform a linear stretch on the image with unlinked channels for all images.
+    """
+    original_min = np.min(image)
+    stretched_image = image - original_min  # Shift image so that the min is 0
+
+    # Capture the original medians before any further processing
+    original_medians = []
+    for c in range(3):  # Assume 3-channel input
+        channel_median = np.median(stretched_image[..., c])
+        original_medians.append(channel_median)
+
+    # Define the target median for stretching
+    target_median = 0.25
+
+    # Apply the stretch for each channel
+    for c in range(3):
+        channel_median = original_medians[c]
+        if channel_median != 0:
+            stretched_image[..., c] = ((channel_median - 1) * target_median * stretched_image[..., c]) / (
+                channel_median * (target_median + stretched_image[..., c] - 1) - target_median * stretched_image[..., c]
+            )
+
+    # Clip stretched image to [0, 1] range
+    stretched_image = np.clip(stretched_image, 0, 1)
+
+    # Return the stretched image, original min, and original medians
+    return stretched_image, original_min, original_medians
+
+
+
+# Function to unstretch an image
+def unstretch_image(image, original_medians, original_min):
+    """
+    Undo the stretch to return the image to the original linear state with unlinked channels.
+    Handles both single-channel and 3-channel images.
+    """
+    was_single_channel = False  # Local flag to check if the image was originally mono
+
+
+
+    # Check if the image is single-channel
+    if image.ndim == 3 and image.shape[2] == 1:
+        was_single_channel = True  # Mark the image as originally single-channel
+        image = np.repeat(image, 3, axis=2)  # Convert to 3-channel by duplicating
+
+    elif image.ndim == 2:  # True single-channel case
+        was_single_channel = True
+        image = np.stack([image] * 3, axis=-1)  # Convert to 3-channel by duplicating
+
+    # Process each channel independently
+    for c in range(3):  # Assume 3-channel input at this point
+        channel_median = np.median(image[..., c])
+        if channel_median != 0 and original_medians[c] != 0:
+
+            image[..., c] = ((channel_median - 1) * original_medians[c] * image[..., c]) / (
+                channel_median * (original_medians[c] + image[..., c] - 1) - original_medians[c] * image[..., c]
+            )
+
+        else:
+            print(f"Channel {c} - Skipping unstretch due to zero median.")
+
+    image += original_min  # Add back the original minimum
+    image = np.clip(image, 0, 1)  # Clip to [0, 1] range
+
+    # If the image was originally single-channel, convert back to single-channel
+    if was_single_channel:
+        image = np.mean(image, axis=2, keepdims=True)  # Convert back to single-channel
+
+
+
+    return image
+
+
+
+
+
+
+
+
+
+
+# Function to add a border of median value around the image
+def add_border(image, border_size=16):
+    if image.ndim == 2:                                # mono
+        med = np.median(image)
+        return np.pad(image,
+                      ((border_size, border_size), (border_size, border_size)),
+                      mode="constant",
+                      constant_values=med)
+
+    elif image.ndim == 3 and image.shape[2] == 3:       # RGB
+        meds = np.median(image, axis=(0, 1)).astype(image.dtype)  # (3,)
+        padded = [np.pad(image[..., c],
+                         ((border_size, border_size), (border_size, border_size)),
+                         mode="constant",
+                         constant_values=float(meds[c]))
+                  for c in range(3)]
+        return np.stack(padded, axis=-1)
+    else:
+        raise ValueError("add_border expects mono or RGB image.")
+
+# Function to remove the border added around the image
+def remove_border(image, border_size=5):
+    if len(image.shape) == 2:
+        return image[border_size:-border_size, border_size:-border_size]
+    else:
+        return image[border_size:-border_size, border_size:-border_size, :]
+
+# Function to denoise the image
+def denoise_image(image_path: str,
+                  denoise_strength: float,
+                  models: dict,
+                  denoise_mode: str = 'luminance',
+                  separate_channels: bool = False,
+                  color_denoise_strength: float | None = None):
+    """
+    Denoises the input image using the specified model and mode.
+
+    Args:
+        image_path (str): Path to the input image.
+        denoise_strength (float): Strength of the denoising effect.
+        device (torch.device or str): Device for PyTorch or ONNX inference.
+        model: PyTorch model or ONNX session for denoising.
+        denoise_mode (str): 'luminance' or 'full' denoising mode.
+        is_onnx (bool): Whether to use ONNX for inference.
+
+    Returns:
+        tuple: Denoised image, original header, bit depth, file extension, is_mono, and metadata.
+    """
+    device      = models["device"]
+    is_onnx     = models["is_onnx"]
+    mono_model  = models["mono_model"]
+
+
+    # Get file extension
+    file_extension = os.path.splitext(image_path)[1].lower()
+    if file_extension not in ['.png', '.tif', '.tiff', '.fit', '.fits', '.xisf', 'jpg', 'jpeg']:
+        print(f"Ignoring non-image file: {image_path}")
+        return None, None, None, None, None, None, None
+
+    original_header = None
+    image = None
+    bit_depth = "32-bit float"  # Default bit depth
+    is_mono = False
+    file_meta, image_meta = None, None
+    stretch_needed = None
+
+    try:
+        # Load the image based on its extension
+        if file_extension in ['.tif', '.tiff']:
+            image = tiff.imread(image_path)
+            print(f"Loaded TIFF image with dtype: {image.dtype}")
+            if image.dtype == np.uint16:
+                image = image.astype(np.float32) / 65535.0
+                bit_depth = "16-bit"
+
+            elif image.dtype == np.uint8:
+                image = image.astype(np.float32) / 255.0
+                bit_depth = "8-bit"
+            elif image.dtype == np.uint32:
+                image = image.astype(np.float32) / 4294967295.0
+                bit_depth = "32-bit unsigned"
+            else:
+                image.dtype == np.float32
+                bit_depth = "32-bit floating point"  # If dtype is already float
+
+            print(f"Final bit depth set to: {bit_depth}")
+
+            # Handle alpha channels and grayscale TIFFs
+            if image.ndim == 3 and image.shape[-1] == 4:
+                print("Detected alpha channel in TIFF. Removing it.")
+                image = image[:, :, :3]
+            if image.ndim == 2:
+                image = np.stack([image] * 3, axis=-1)  # Convert to RGB
+                is_mono = True
+            else:
+                is_mono = False
+
+        elif file_extension in ['.fits', '.fit']:
+            # Load the FITS image
+            with fits.open(image_path) as hdul:
+                image_data = hdul[0].data
+                original_header = hdul[0].header  # Capture the FITS header
+
+                # Ensure the image data uses the native byte order
+                if image_data.dtype.byteorder not in ('=', '|'):  # Check if byte order is non-native
+                    image_data = image_data.astype(image_data.dtype.newbyteorder('='))  # Convert to native byte order
+
+                # Determine the bit depth based on the data type in the FITS file
+                if image_data.dtype == np.uint16:
+                    bit_depth = "16-bit"
+                    print("Identified 16bit FITS image.")
+                elif image_data.dtype == np.float32:
+                    bit_depth = "32-bit floating point"
+                    print("Identified 32bit floating point FITS image.")
+                elif image_data.dtype == np.uint32:
+                    bit_depth = "32-bit unsigned"
+                    print("Identified 32bit unsigned FITS image.")
+                
+                # Handle 3D FITS data (e.g., RGB or multi-layered data)
+                if image_data.ndim == 3:
+                    print(f"Detected 3D FITS data with shape: {image_data.shape}")
+                    
+                    if image_data.shape[0] == 3:
+                        # Assume RGB in (3, H, W) format
+                        image = np.transpose(image_data, (1, 2, 0))  # Convert to (H, W, 3)
+                        
+                        if bit_depth == "16-bit":
+                            image = image.astype(np.float32) / 65535.0
+                        elif bit_depth == "32-bit unsigned":
+                            bzero = original_header.get('BZERO', 0)
+                            bscale = original_header.get('BSCALE', 1)
+                            image = image.astype(np.float32) * bscale + bzero
+                            image_min = image.min()
+                            image_max = image.max()
+                            image = (image - image_min) / (image_max - image_min)
+                            print(f"Image range after applying BZERO and BSCALE (RGB): min={image_min}, max={image_max}")
+                        # No need to normalize float32
+                        is_mono = False
+
+                    elif image_data.shape[0] == 1:
+                        # Single-channel mono FITS in (1, H, W) format
+                        image = image_data[0]  # Drop channel dim, now (H, W)
+                        image = np.stack([image] * 3, axis=-1)  # Convert to RGB shape (H, W, 3)
+                        is_mono = True
+
+                    elif image_data.shape[2] == 1:
+                        # Shape (H, W, 1) — mono with trailing singleton dim
+                        image = image_data[:, :, 0]
+                        image = np.stack([image] * 3, axis=-1)
+                        is_mono = True
+
+                    else:
+                        raise ValueError(f"Unsupported 3D FITS shape: {image_data.shape}")
+
+                # Handle 2D FITS data (grayscale)
+                elif image_data.ndim == 2:
+                    if bit_depth == "16-bit":
+                        image = image_data.astype(np.float32) / 65535.0  # Normalize to [0, 1] for 16-bit
+                    elif bit_depth == "32-bit unsigned":
+                        # Apply BSCALE and BZERO if present
+                        bzero = original_header.get('BZERO', 0)  # Default to 0 if not present
+                        bscale = original_header.get('BSCALE', 1)  # Default to 1 if not present
+
+                        # Convert to float and apply the scaling and offset
+                        image = image_data.astype(np.float32) * bscale + bzero
+
+                        # Normalize based on the actual data range
+                        image_min = image.min()  # Get the min value after applying BZERO
+                        image_max = image.max()  # Get the max value after applying BZERO
+
+                        # Normalize the image data to the range [0, 1]
+                        image = (image - image_min) / (image_max - image_min)
+                        print(f"Image range after applying BZERO and BSCALE (2D case): min={image_min}, max={image_max}")
+
+
+
+
+                    elif bit_depth == "32-bit floating point":
+                        image = image_data  # No normalization needed for 32-bit float
+                    is_mono = True
+                    image = np.stack([image] * 3, axis=-1)  # Convert to 3-channel for consistency
+                else:
+                    raise ValueError("Unsupported FITS format!")
+
+
+        # Check if file extension is '.xisf'
+        elif file_extension == '.xisf':
+            # Load XISF file
+            xisf = XISF(image_path)
+            image = xisf.read_image(0)  # Assuming the image data is in the first image block
+            image_meta = xisf.get_images_metadata()  # List of metadata blocks for each image
+            file_meta = xisf.get_file_metadata()  # File-level metadata
+
+            # Determine bit depth
+            if image.dtype == np.uint16:
+                image = image.astype(np.float32) / 65535.0
+                bit_depth = "16-bit"
+            elif image.dtype == np.uint8:
+                image = image.astype(np.float32) / 255.0
+                bit_depth = "8-bit"                
+            elif image.dtype == np.uint32:
+                image = image.astype(np.float32) / 4294967295.0
+                bit_depth = "32-bit unsigned"
+            elif image.dtype == np.float32:
+                bit_depth = "32-bit floating point"
+
+            # Check if image is mono (2D array) and stack into 3 channels if mono
+            if len(image.shape) == 2 or (len(image.shape) == 3 and image.shape[2] == 1):
+                image = np.stack([image.squeeze()] * 3, axis=-1)  # Convert mono to RGB by duplicating channels
+                is_mono = True
+            else:
+                is_mono = False
+
+            print(f"Loaded XISF image with bit depth: {bit_depth}, Mono: {is_mono}")
+
+
+        else:  # Assume 8-bit PNG
+            image = np.array(Image.open(image_path).convert('RGB')).astype(np.float32) / 255.0
+            bit_depth = "8-bit"
+            print(f"Loaded {bit_depth} PNG image.")
+
+        if image is None:
+            raise ValueError(f"Failed to load image from: {image_path}")
+
+        # Add a border around the image with the median value
+        stretch_needed = np.median(image - np.min(image)) < 0.05   # decision first
+
+        if stretch_needed:
+            print("normalizing linear data")
+            stretched_core, original_min, original_medians = stretch_image(image)
+        else:
+            stretched_core   = image.astype(np.float32, copy=False)
+            original_min     = np.min(image)
+            original_medians = [np.median(image[..., c]) for c in range(3)] \
+                            if image.ndim == 3 else [np.median(image)]
+
+        # pad AFTER stretch, per-channel median
+        stretched_image = add_border(stretched_core, border_size=16)
+
+        # **Apply TV Denoise on the full image unconditionally**
+        #tv_weight = 0.005  # Fixed weight for TV denoising
+        #print(f"Applying Total Variation Denoising (TVDenoise) on the full image with fixed weight: {tv_weight}")
+        #denoised_full = np.zeros_like(stretched_image)
+        #for c in range(stretched_image.shape[2]):
+        #    denoised_full[:, :, c] = denoise_tv_chambolle(stretched_image[:, :, c], weight=tv_weight, channel_axis=-1)
+        #stretched_image = denoised_full
+        #print("TV Denoising applied.")
+
+        # Process mono or color images
+        if is_mono:
+            # use mono network on channel 0
+            mono_net = models["mono_model"]
+            den_r = denoise_channel(stretched_image[...,0], models["device"], mono_net, True, models["is_onnx"])
+            den = denoise_strength * den_r + (1-denoise_strength)*stretched_image[...,0]
+            denoised_image = den[...,None]
+
+        else:
+            # 1) separate‐channels override
+            if separate_channels:
+                out_ch = []
+                for c in range(3):
+                    dch = denoise_channel(
+                        stretched_image[...,c],
+                        models["device"],
+                        models["mono_model"],
+                        True, models["is_onnx"]
+                    )
+                    blended = blend_images(stretched_image[...,c], dch, denoise_strength)
+                    out_ch.append(blended)
+                denoised_image = np.stack(out_ch, axis=-1)
+
+            # 2) luminance only
+            elif denoise_mode == 'luminance':
+                y, cb, cr = extract_luminance(stretched_image)
+                den_y = denoise_channel(
+                    y,
+                    models["device"],
+                    models["mono_model"],
+                    True, models["is_onnx"]
+                )
+                y2 = blend_images(y, den_y, denoise_strength)
+                denoised_image = merge_luminance(y2, cb, cr)
+
+
+            else:
+                # --- FULL mode: L via NN (AI3_5), chroma via traditional denoise, then recombine ---
+                # Split to YCbCr
+                y, cb, cr = extract_luminance(stretched_image)
+
+                # Denoise L with mono NN
+                # Denoise L with mono NN
+                den_y = denoise_channel(
+                    y,
+                    models["device"],
+                    models["mono_model"],
+                    True, models["is_onnx"]
+                )
+                y2 = blend_images(y, den_y, denoise_strength)
+
+                colorstrength = color_denoise_strength if color_denoise_strength is not None else denoise_strength
+                cb2, cr2 = denoise_chroma(cb, cr,
+                                        strength=colorstrength,
+                                        method="guided",
+                                        guide_y=y)
+
+                # Recombine to RGB
+                denoised_image = merge_luminance(y2, cb2, cr2)
+
+
+        # Unstretch the image
+        if stretch_needed:
+            print("de-normalizing linear data")
+            denoised_image = unstretch_image(denoised_image, original_medians, original_min)
+
+
+        # Remove the border added around the image
+        
+        denoised_image = remove_border(denoised_image, border_size=16)
+        
+        return denoised_image, original_header, bit_depth, file_extension, is_mono, file_meta, image_meta
+
+    except Exception as e:
+        print(f"Error reading image {image_path}: {e}")
+        return None, None, None, None, None, None, None
+    
+
+
+# Function to denoise a single channel
+def denoise_channel(channel, device, model, is_mono=False, is_onnx=False):
+    """
+    Denoises a single image channel using the specified model.
+
+    Args:
+        channel (numpy.ndarray): The input image channel to denoise.
+        device (torch.device or str): Device for PyTorch or ONNX inference.
+        model: PyTorch model or ONNX inference session.
+        is_mono (bool): Whether the input is monochrome.
+        is_onnx (bool): Whether to use ONNX for inference.
+
+    Returns:
+        numpy.ndarray: The denoised image channel.
+    """
+    chunk_size = 256
+    overlap = 64
+    chunks = split_image_into_chunks_with_overlap(channel, chunk_size=chunk_size, overlap=overlap)
+
+    # Debug: validate chunk structure
+    for idx, entry in enumerate(chunks[:3]):
+        print(f"Chunk {idx} type: {type(entry)}, length: {len(entry) if isinstance(entry, tuple) else 'N/A'}")
+
+    denoised_chunks = []
+
+    for idx, entry in enumerate(chunks):
+        try:
+            chunk, i, j = entry
+        except Exception as e:
+            print(f"❌ Failed to unpack chunk {idx}: {entry} ({e})")
+            continue
+
+        original_chunk_shape = chunk.shape
+        if is_onnx:
+            # Prepare ONNX input: Expand to 3 channels
+            chunk_input = chunk[np.newaxis, np.newaxis, :, :].astype(np.float32)  # Shape: (1, 1, H, W)
+            chunk_input = np.tile(chunk_input, (1, 3, 1, 1))  # Shape: (1, 3, H, W)
+
+            # Pad the chunk to 256x256 if necessary
+            if chunk_input.shape[2] != chunk_size or chunk_input.shape[3] != chunk_size:
+                padded_chunk = np.zeros((1, 3, chunk_size, chunk_size), dtype=np.float32)
+                padded_chunk[:, :, :chunk_input.shape[2], :chunk_input.shape[3]] = chunk_input
+                chunk_input = padded_chunk
+
+            # ONNX inference
+            input_name = model.get_inputs()[0].name
+            try:
+                denoised_chunk = model.run(None, {input_name: chunk_input})[0]
+                denoised_chunk = denoised_chunk[0, 0, :original_chunk_shape[0], :original_chunk_shape[1]]
+            except Exception as e:
+                print(f"ONNX inference error for chunk at ({i}, {j}): {e}")
+                continue
+
+        else:
+            # Prepare PyTorch input
+            chunk_tensor = torch.tensor(chunk, dtype=torch.float32).unsqueeze(0).unsqueeze(0).to(device)
+            if not is_mono:
+                chunk_tensor = chunk_tensor.repeat(1, 3, 1, 1)
+            else:
+                chunk_tensor = chunk_tensor.expand(1, 3, chunk_tensor.shape[2], chunk_tensor.shape[3])
+
+            # PyTorch inference
+            try:
+                with torch.no_grad():
+                    if not DISABLE_MIXED_PRECISION and device.type == "cuda":
+                        with torch.cuda.amp.autocast():
+                            denoised_output = model(chunk_tensor).squeeze().cpu().numpy()
+                    else:
+                        denoised_output = model(chunk_tensor).squeeze().cpu().numpy()
+
+                    if denoised_output.ndim == 3:
+                        denoised_chunk = denoised_output[0]  # Take first channel
+                    else:
+                        denoised_chunk = denoised_output
+
+            except Exception as e:
+                print(f"PyTorch inference error for chunk at ({i}, {j}): {e}")
+                continue
+
+        if denoised_chunk is not None:
+            denoised_chunks.append((denoised_chunk, i, j))
+        else:
+            print(f"⚠️ Warning: Denoised chunk at ({i}, {j}) is None — skipping")
+
+        # Show progress update
+        show_progress(idx + 1, len(chunks))
+
+    # Debug: final validation before stitching
+    if not all(isinstance(entry, tuple) and len(entry) == 3 for entry in denoised_chunks):
+        print("❗ Error: One or more denoised chunks are malformed before stitching!")
+
+    # Stitch the chunks back together
+    denoised_channel = stitch_chunks_ignore_border(denoised_chunks, channel.shape, chunk_size=chunk_size, overlap=overlap)
+
+    return denoised_channel
+
+
+
+
+
+# Main process for denoising images
+# Main process for denoising images
+def process_images(input_dir, output_dir,
+                   denoise_strength=None,
+                   color_denoise_strength=None,
+                   use_gpu=True,
+                   denoise_mode='luminance',
+                   separate_channels=False):
+    print((r"""
+ *#        ___     __      ___       __                              #
+ *#       / __/___/ /__   / _ | ___ / /________                      #
+ *#      _\ \/ -_) _ _   / __ |(_-</ __/ __/ _ \                     #
+ *#     /___/\__/\//_/  /_/ |_/___/\__/__/ \___/                     #
+ *#                                                                  #
+ *#              Cosmic Clarity - Denoise V6.6   AI3.6               # 
+ *#                                                                  #
+ *#                         SetiAstro                                #
+ *#                    Copyright © 2025                              #
+ *#                                                                  #
+        """))
+
+    if denoise_strength is None:
+        # GUI path
+        use_gpu, denoise_strength, denoise_mode, separate_channels, gui_color_strength = get_user_input()
+        # if caller didn't supply a color strength, take GUI's; otherwise keep caller's
+        if color_denoise_strength is None:
+            color_denoise_strength = gui_color_strength
+    else:
+        # headless path: if color strength not provided, match luminance strength
+        if color_denoise_strength is None:
+            color_denoise_strength = denoise_strength
+
+    # Create output directory if it doesn't exist
+    if not os.path.exists(output_dir):
+        os.makedirs(output_dir)
+
+    # Load the denoise model
+    models = load_models(exe_dir, use_gpu)
+
+    # Determine whether we're using ONNX or PyTorch
+    print(f"Using {'ONNX' if models['is_onnx'] else 'PyTorch'} models.")
+
+    # Process each image in the input directory
+    for fname in os.listdir(input_dir):
+        if not fname.lower().endswith(('.png', '.tif', '.tiff', '.fits', '.fit', '.xisf', '.jpg', '.jpeg')):
+            continue
+
+        # Unpack exactly what denoise_image returns
+        denoised_image, original_header, bit_depth, file_extension, is_mono, file_meta, image_meta = denoise_image(
+            os.path.join(input_dir, fname),
+            denoise_strength,
+            models,
+            denoise_mode,
+            separate_channels,
+            color_denoise_strength
+        )
+        # Skip on error
+        if denoised_image is None:
+            continue
+
+        # Build output path
+        output_name = os.path.splitext(fname)[0] + "_denoised"
+        output_image_path = os.path.join(output_dir, output_name + file_extension)
+
+        # Save as FITS file with header information if the original was FITS
+        if file_extension in ['.fits', '.fit']:
+            if original_header is not None:
+                if is_mono:  # Grayscale FITS
+                    # Convert the grayscale image back to its original 2D format
+                    if bit_depth == "16-bit":
+                        denoised_image_fits = (denoised_image[:, :, 0] * 65535).astype(np.uint16)
+                        actual_bit_depth = "16-bit"
+                    elif bit_depth == "32-bit unsigned":
+                        bzero = original_header.get('BZERO', 0)
+                        bscale = original_header.get('BSCALE', 1)
+                        denoised_image_fits = (denoised_image[:, :, 0].astype(np.float32) * bscale + bzero).astype(np.uint32)
+                        original_header['BITPIX'] = 32
+                        actual_bit_depth = "32-bit unsigned"
+                    else:  # 32-bit float
+                        denoised_image_fits = denoised_image[:, :, 0].astype(np.float32)
+                        actual_bit_depth = "32-bit"
+                    
+                    # Update header for a 2D (grayscale) image
+                    original_header['NAXIS'] = 2
+                    original_header['NAXIS1'] = denoised_image.shape[1]  # Width
+                    original_header['NAXIS2'] = denoised_image.shape[0]  # Height
+                    if 'NAXIS3' in original_header:
+                        del original_header['NAXIS3']  # Remove if present
+
+                    hdu = fits.PrimaryHDU(denoised_image_fits, header=original_header)
+                
+                else:  # RGB FITS
+                    # Transpose RGB image to FITS-compatible format (channels, height, width)
+                    denoised_image_transposed = np.transpose(denoised_image, (2, 0, 1))
+                    actual_bit_depth=None
+                    if bit_depth == "16-bit":
+                        denoised_image_fits = (denoised_image_transposed * 65535).astype(np.uint16)
+                        actual_bit_depth = "16-bit"
+                    elif bit_depth == "32-bit unsigned":
+                        denoised_image_fits = denoised_image_transposed.astype(np.float32)
+                        original_header['BITPIX'] = -32
+                        actual_bit_depth = "32-bit unsigned"
+                    else:
+                        denoised_image_fits = denoised_image_transposed.astype(np.float32)
+                        actual_bit_depth = "32-bit float"
+
+                    # Update header for a 3D (RGB) image
+                    original_header['NAXIS'] = 3
+                    original_header['NAXIS1'] = denoised_image_transposed.shape[2]  # Width
+                    original_header['NAXIS2'] = denoised_image_transposed.shape[1]  # Height
+                    original_header['NAXIS3'] = denoised_image_transposed.shape[0]  # Channels
+                    
+                    hdu = fits.PrimaryHDU(denoised_image_fits, header=original_header)
+
+                # Write the FITS file
+                hdu.writeto(output_image_path, overwrite=True)
+                print(f"Saved {actual_bit_depth} denoised image to: {output_image_path}")
+
+
+        # Save as TIFF based on the original bit depth if the original was TIFF
+        elif file_extension in ['.tif', '.tiff']:
+            if bit_depth == "16-bit":
+                actual_bit_depth = "16-bit"
+                if is_mono is True:  # Grayscale
+                    tiff.imwrite(output_image_path, (denoised_image[:, :, 0] * 65535).astype(np.uint16))
+                else:  # RGB
+                    tiff.imwrite(output_image_path, (denoised_image * 65535).astype(np.uint16))
+            elif bit_depth == "8-bit":
+                actual_bit_depth = "8-bit"
+                if is_mono:
+                    tiff.imwrite(output_image_path, (denoised_image[:, :, 0] * 255.0).astype(np.uint8))
+                else:
+                    tiff.imwrite(output_image_path, (denoised_image * 255.0).astype(np.uint8))                               
+            elif bit_depth == "32-bit unsigned":
+                actual_bit_depth = "32-bit unsigned"
+                if is_mono is True:  # Grayscale
+                    tiff.imwrite(output_image_path, (denoised_image[:, :, 0] * 4294967295).astype(np.uint32))
+                else:  # RGB
+                    tiff.imwrite(output_image_path, (denoised_image * 4294967295).astype(np.uint32))           
+            else:
+                actual_bit_depth = "32-bit float"
+                if is_mono is True:  # Grayscale
+                    tiff.imwrite(output_image_path, denoised_image[:, :, 0].astype(np.float32))
+                else:  # RGB
+                    tiff.imwrite(output_image_path, denoised_image.astype(np.float32))
+
+            print(f"Saved {actual_bit_depth} denoised image to: {output_image_path}")
+
+        elif file_extension == '.xisf':
+            try:
+                # Debug: Print original image details
+                print(f"Original image shape: {denoised_image.shape}, dtype: {denoised_image.dtype}")
+                print(f"Bit depth: {bit_depth}")
+
+                # Adjust bit depth
+                if bit_depth == "16-bit":
+                    processed_image = (denoised_image * 65535).astype(np.uint16)
+                elif bit_depth == "8-bit":
+                    processed_image = (denoised_image * 255.0).astype(np.uint8)                        
+                elif bit_depth == "32-bit unsigned":
+                    processed_image = (denoised_image * 4294967295).astype(np.uint32)
+                else:  # Default to 32-bit float
+                    processed_image = denoised_image.astype(np.float32)
+
+                # Adjust for mono images
+                if is_mono:
+                    print("Preparing mono image...")
+                    processed_image = processed_image[:, :, 0]  # Take the first channel
+                    processed_image = processed_image[:, :, np.newaxis]  # Add back channel dimension
+                    image_meta[0]['geometry'] = (processed_image.shape[1], processed_image.shape[0], 1)
+                    image_meta[0]['colorSpace'] = 'Gray'  # Update metadata for mono
+
+                # Debug: Print processed image details
+                print(f"Processed image shape: {processed_image.shape}, dtype: {processed_image.dtype}")
+
+                # Save the image
+                XISF.write(
+                    output_image_path,                  # Correct output path
+                    processed_image,                   # Final processed image
+                    creator_app="Seti Astro Cosmic Clarity",
+                    image_metadata=image_meta[0],      # First block of image metadata
+                    xisf_metadata=file_meta,           # File-level metadata
+
+                    shuffle=True
+                )
+                print(f"Saved {bit_depth} XISF denoised image to: {output_image_path}")
+
+            except Exception as e:
+                print(f"Error saving XISF file: {e}")
+
+
+
+
+        # Save as 8-bit PNG if the original was PNG
+        else:
+            output_image_path = os.path.join(output_dir, output_name + ".png")
+            denoised_image_8bit = (denoised_image * 255).astype(np.uint8)
+            denoised_image_pil = Image.fromarray(denoised_image_8bit)
+            actual_bit_depth = "8-bit"
+            denoised_image_pil.save(output_image_path)
+            print(f"Saved {actual_bit_depth} denoised image to: {output_image_path}")
+
+
+
+
+# Define input and output directories
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Cosmic Clarity Denoise Tool")
+    parser.add_argument(
+        '--denoise_strength',
+        type=float,
+        help="Denoise strength (0–1), overrides the GUI slider if provided"
+    )
+    parser.add_argument(
+        '--disable_gpu',
+        action='store_true',
+        help="Disable GPU acceleration and force CPU usage"
+    )
+    parser.add_argument(
+        '--denoise_mode',
+        choices=['luminance','full','separate'],
+        default='luminance',
+        help="Denoise mode: 'luminance', 'full' color model, or 'separate' per-channel"
+    )
+    parser.add_argument(
+        '--separate_channels',
+        action='store_true',
+        help="Alias for --denoise_mode separate"
+    )
+
+    parser.add_argument(
+        '--color_denoise_strength',
+        type=float,
+        help="Color (chroma) denoise strength (0–1). Defaults to --denoise_strength if omitted."
+    )
+
+    args = parser.parse_args()
+
+    # Ensure input/output folders exist
+    input_dir  = os.path.join(exe_dir, 'input')
+    output_dir = os.path.join(exe_dir, 'output')
+    os.makedirs(input_dir, exist_ok=True)
+    os.makedirs(output_dir, exist_ok=True)
+
+    # If the user requested separate_channels, override denoise_mode
+    mode = args.denoise_mode
+    if args.separate_channels:
+        mode = 'separate'
+
+    # Launch processing
+    process_images(
+        input_dir,
+        output_dir,
+        denoise_strength=args.denoise_strength,
+        color_denoise_strength=args.color_denoise_strength,
+        use_gpu=not args.disable_gpu,
+        denoise_mode=mode,
+        separate_channels=args.separate_channels
+    )