Metadata-Version: 2.4
Name: with-simple-color
Version: 0.1.0
Summary: Palette-constrained color halftoning optimized for Gaussian reconstruction.
License-Expression: MIT
License-File: LICENSE
Author: GGN_2015
Author-email: neko@jlulug.org
Requires-Python: >=3.10
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Requires-Dist: matplotlib (>=3.10.9)
Requires-Dist: numba (>=0.67.0)
Requires-Dist: numpy (>=2.2.6)
Requires-Dist: pillow (>=10.0.0)
Requires-Dist: scipy (>=1.15.3)
Description-Content-Type: text/markdown

# with_simple_color

Create RGB halftones using only a specified palette. Rather than minimizing
the difference between individual output and input pixels, this project
minimizes the difference between the **Gaussian-blurred output** and the input.
It is intended for images viewed through a spatial averaging process.

**The returned and saved result is always the unblurred palette image.**
Gaussian blur is used internally to evaluate the reconstruction error, not as
an output effect. The CLI does not save or display a blurred reconstruction.

## Installation

Python 3.10 or newer is required. From the project directory:

```sh
uv sync
uv run with-simple-color --help
```

Alternatively, install with `python -m pip install .` and use
`python -m with_simple_color`. NumPy and SciPy provide numerical operations,
Numba accelerates the sequential optimization and initialization, Pillow handles
image files, and Matplotlib provides optional comparison plots. The first run
includes Numba compilation overhead; compiled kernels are cached.

## Quick start

Run the built-in 64-by-64 demo without opening a window:

```sh
uv run with-simple-color --sigma 2 --output halftone.png
```

Reproduce the original larger synthetic example with a fixed seed:

```sh
uv run with-simple-color --size 256 --sigma 10 --max-iter 15 --seed 42 --output halftone.png
```

With the default eight-color CLI palette, this example reduces total squared
reconstruction error from about `62.3408` after initialization to `1.36979`
after refinement, compared with `12555.8` for nearest-color quantization alone.

Process an image using the eight RGB cube corners:

```sh
uv run with-simple-color --input photo.png --palette cube --sigma 2 --max-iter 30 --output halftone.png
```

Available built-in CLI palettes are `cube` (the default, eight colors), `rgbk`
(red, green, blue, black), `cmyk` (white, black, yellow, magenta, cyan at
0.75 intensity), and `bw` (black and white). A JSON file can specify an
arbitrary nonempty list of normalized RGB triples, including fractional colors:

```json
[
  [0.10, 0.15, 0.20],
  [0.85, 0.30, 0.20],
  [0.95, 0.90, 0.80]
]
```

```sh
uv run with-simple-color --input photo.png --palette palette.json --output halftone.npy
```

Use `--show` for an interactive comparison, or `--comparison comparison.png`
to save a two-panel comparison of the original and unblurred palette image.
Neither option displays a blurred image. `--seed 0` is the CLI default.
`--boundary`, `--init`, `--switch-fraction`, and `--tolerance` expose optimization
settings; `--quiet` suppresses progress. `--size` changes the synthetic demo size.
All output paths must be distinct and must not overwrite an input image or palette.

### Example outputs

Unblurred palette halftones of the same image, rendered with the built-in
palettes:

| Original | `cube` | `rgbk` | `cmyk` | `bw` |
| --- | --- | --- | --- | --- |
| ![](https://github.com/GGN-2015/with_simple_color/blob/main/img/raw.webp?raw=true) | ![](https://github.com/GGN-2015/with_simple_color/blob/main/img/cube.png?raw=true) | ![](https://github.com/GGN-2015/with_simple_color/blob/main/img/rgbk.png?raw=true) | ![](https://github.com/GGN-2015/with_simple_color/blob/main/img/cmyk.png?raw=true) | ![](https://github.com/GGN-2015/with_simple_color/blob/main/img/bw.png?raw=true) |

### Image and color conventions

- The Python API and `.npy` inputs require finite arrays of shape
  `(height, width, 3)` with values in `[0, 1]`. Integer 0/1 arrays are valid;
  divide 8-bit image arrays by 255 before calling the API.
- Ordinary image files are loaded as 8-bit RGB, with EXIF orientation applied.
  Grayscale is expanded to RGB; transparency is composited over white.
- API results and `.npy` outputs are `float64`, preserving the supplied palette
  values after conversion to `float64`. Use `.npy` for exact fractional colors.
- PNG output rounds to 8-bit RGB. Fractional palette entries not representable
  as multiples of `1/255` cannot be preserved exactly in PNG. Reported loss
  refers to the in-memory result, before PNG rounding.
- The objective operates in the supplied RGB coordinates. Image loading does
  not perform ICC color management or linear-light conversion. For a physical
  light-averaging model, convert both the image and palette to linear-light RGB
  before using the API, and apply the appropriate display conversion afterward.

## Python API

```python
import numpy as np

from with_simple_color import (
    compute_loss,
    quantize,
    validate_B_general,
)

image = np.full((32, 32, 3), 0.4)
palette = np.array([[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]])
quantized, losses = quantize(
    image,
    sigma=2.0,
    palette=palette,
    max_iter=30,
    seed=42,
    verbose=False,
    return_history=True,
)
np.save("halftone.npy", quantized)
assert validate_B_general(quantized, palette)[0]
assert all(after <= before for before, after in zip(losses, losses[1:]))
print(compute_loss(quantized, image, sigma=2.0))
```

Without `return_history=True`, `quantize` returns only the unblurred image. The API's
default palette remains red, green, blue, and black, unlike the CLI demo's
eight-color palette. Existing helper names remain available from
`with_simple_color.main` and the package root. Use `image=` and `quantized=`
for named image arguments; the old positional argument order is unchanged.

Important options:

| Argument | Meaning |
| --- | --- |
| `sigma` | Finite positive Gaussian standard deviation in pixels. |
| `init="auto"` | Select the lower actual reconstruction loss of nearest-color and serpentine Floyd-Steinberg initializations. |
| `init="nearest"` | Start with the closest palette color at each pixel. |
| `init="error_diffusion"` | Use serpentine Floyd-Steinberg diffusion as a starting point, not as the optimizer. |
| `init="probabilistic"` | Sample distance-weighted palette colors. This is only a heuristic, not an unbiased color mixture. |
| `max_iter=30` | Maximum sequential sweeps; zero returns the initialization unchanged. |
| `switch_fraction=1.0` | Cap accepted changes per sweep at `ceil(switch_fraction * height * width)`, in `(0, 1]`. Smaller values may need more sweeps. |
| `boundary_mode="reflect"` | `reflect`, `mirror`, `nearest`, `wrap`, or `constant`. Constant padding is zero; SciPy's grid aliases are also accepted by the API. |
| `seed=None` | Seed for random coordinate order and probabilistic initialization; set an integer for repeatability. |
| `tolerance=1e-12` | Minimum absolute decrease in total squared error needed to accept one recoloring. |
| `return_history=False` | Return `(quantized, losses)` when true. History starts with the initialization loss. |

`refine_general(quantized, image, sigma, palette, ...)` improves an existing
palette image. Inputs are not mutated. `plot_comparison(..., show=False)`
returns a two-panel Matplotlib figure containing only the original and unblurred
palette image, suitable for saving without an interactive window. Plotting is
imported lazily; the internal reconstruction is not clipped before computing
the objective.

## Objective and algorithm

Let `image` be the target and `quantized[p]` be one palette color at every pixel.
Let `H` be the spatial Gaussian blur, applied independently to each RGB channel.
The objective is:

```text
minimize    L(quantized) = ||H(quantized) - image||_F^2
subject to  quantized[p] belongs to palette for every pixel p
```

The target is **not** blurred. `H` uses SciPy's discrete, normalized Gaussian
with `truncate=4.0` and radius `int(4 * sigma + 0.5)`, including the selected
boundary extension. The channel axis is never blurred.

For a recoloring at pixel `p`, define `delta = new_color - old_color`, the current
residual `residual = H(quantized) - image`, and `column_p = H(unit_impulse_at_p)`.
Its exact change in loss is:

```text
delta_loss = 2 * dot(delta, sum_over_pixels(column_p * residual))
             + dot(delta, delta) * sum_over_pixels(column_p ** 2)
```

The optimizer evaluates all palette colors for one pixel, accepts the best
strictly improving choice, and **immediately updates the residual** before
considering another pixel. This is discrete coordinate descent with exact
single-pixel moves. It avoids treating overlapping pixel changes as independent.

Boundary-aware, sparse one-dimensional operator columns are constructed and
combined separably. Their squared norms depend on pixel position near edges;
they are not replaced by an infinite-domain kernel norm. This also avoids
incorrectly assuming that every boundary extension makes `H` self-adjoint.
No dense two-dimensional pixel-to-pixel blur matrix is constructed.

Each completed sweep recomputes the residual using SciPy. A sweep that fails
to reduce the measured objective is rolled back to protect against numerical
drift. Returned loss history therefore never increases. A full sweep with no
accepted changes certifies only that no single-pixel move improves the objective
by more than `tolerance`.

### Why the previous implementation failed

- Simultaneous switches were scored separately, omitting cross terms between
  overlapping Gaussian footprints. The combined update could increase error.
  For an 8-by-8 RGB target of `0.25`, a black initialization, a black/white
  palette, and `sigma=1`, one old sweep increased loss from `12` to `108`.
- A single interior kernel norm was used at every pixel, despite boundary
  extensions changing the effective impulse response. Its radius also differed
  from SciPy's rounding rule for some sigma values.
- Applying the forward blur as an adjoint was not valid for every supported
  boundary mode.
- Casting colors to `int8` destroyed fractional palettes, and floating-point
  diffusion offsets could not be used as array indices by Numba.
- Distance-based random sampling was not a mean-preserving palette mixture.

### Limitations and cost

This is a discrete, nonconvex optimization problem. The implementation does
**not** promise the global minimum. It uses recolorings, not multi-pixel swaps;
escaping a one-pixel local optimum can require coordinated changes. Reaching
`max_iter` does not certify convergence. Different seeds and initializations
can yield different local solutions.

For normalized boundary modes, every reconstructed color is a convex combination
of palette colors. A target outside that convex hull cannot be reproduced
exactly. In particular, red/green/blue/black cannot represent white by spatial
mixing. The eight RGB cube corners span the full unit RGB cube, but spatial
coupling and discrete pixel counts still prevent arbitrary exact reconstruction.
Zero padding in `constant` mode additionally mixes black into border pixels.

With `pixel_count = height * width`, `color_count` palette colors, and Gaussian
radius `radius`, a sweep costs approximately
`O(pixel_count * ((2 * radius + 1)**2 + color_count))`. Large blur radii and
high-resolution images can be expensive; begin with a small image and a modest
sweep limit. Sparse operator construction uses
`O((height + width) * (2 * radius + 1))` temporary storage, in addition to image
buffers and the palette. There is no `height * width * color_count` cost volume.

## Development checks

```sh
uv run python -m compileall -q with_simple_color
uv run python -m with_simple_color --size 16 --max-iter 3 --output smoke.npy --comparison smoke.png
uv build
```

The numerical implementation can be independently checked by constructing each
Gaussian impulse with SciPy and comparing the predicted `delta_loss` with two
full objective evaluations. For tiny images, enumerate every recoloring after
convergence to verify one-pixel optimality. Check all five boundary modes,
singleton image dimensions, kernels larger than the image, fractional palettes,
and nonincreasing loss histories rather than judging only a displayed picture.

## License

MIT. See `LICENSE`.

