ColorMoments — Color Moments¶

def transforms():
    """(name, x-axis label, [(strength, PIL image -> (image, (file extension, quality)))]).

    Each edit is applied alone to the original. A variant is saved as PPM, which is
    lossless (and quick to write for a corpus of thousands of variants), so the edit is
    the only difference; JPEG quality is the one edit that is the saving itself.
    """
    def jpeg(q):
        return lambda im: (im, ("jpg", q))

    def scale(s):
        return lambda im: (im.resize((max(1, round(im.width * s)), max(1, round(im.height * s))),
                                     Image.LANCZOS), ("ppm", None))

    def rotate(deg):
        # The corners the turn uncovers are filled with the image's mean color, so they
        # add no edge of their own.
        def f(im):
            fill = tuple(int(v) for v in np.asarray(im).reshape(-1, 3).mean(axis=0))
            return im.rotate(deg, resample=Image.BICUBIC, fillcolor=fill), ("ppm", None)
        return f

    def brightness(k):
        return lambda im: (ImageEnhance.Brightness(im).enhance(k), ("ppm", None))

    def contrast(k):
        return lambda im: (ImageEnhance.Contrast(im).enhance(k), ("ppm", None))

    def gamma(g):
        def f(im):
            a = (np.asarray(im).astype(np.float64) / 255.0) ** g
            return Image.fromarray((a * 255.0).round().astype(np.uint8)), ("ppm", None)
        return f

    def blur(r):
        return lambda im: (im.filter(ImageFilter.GaussianBlur(r)), ("ppm", None))

    def noise(sigma):
        # A fixed seed: the same noise on every build, so the chart does not move.
        def f(im):
            rng = np.random.default_rng(1)
            a = np.asarray(im).astype(np.float64) + rng.normal(0, sigma, (im.height, im.width, 3))
            return Image.fromarray(np.clip(a, 0, 255).round().astype(np.uint8)), ("ppm", None)
        return f

    def crop(p):
        def f(im):
            dx, dy = round(im.width * p / 2), round(im.height * p / 2)
            return im.crop((dx, dy, im.width - dx, im.height - dy)), ("ppm", None)
        return f

    return [
        ("JPEG quality", "quality", [(q, jpeg(q)) for q in (95, 75, 50, 30, 15, 5)]),
        ("Downscale", "scale", [(s, scale(s)) for s in (0.75, 0.5, 0.25, 0.125)]),
        ("Rotation", "degrees", [(d, rotate(d)) for d in (1, 2, 5, 10, 20, 45)]),
        ("Brightness", "factor", [(k, brightness(k)) for k in (0.5, 0.7, 0.85, 1.15, 1.3, 1.5)]),
        ("Contrast", "factor", [(k, contrast(k)) for k in (0.5, 0.7, 0.85, 1.15, 1.3, 1.5)]),
        ("Gamma", "exponent", [(g, gamma(g)) for g in (0.5, 0.7, 0.85, 1.2, 1.5, 2.0)]),
        ("Gaussian blur", "radius, px", [(r, blur(r)) for r in (0.5, 1, 2, 4, 8)]),
        ("Noise", "σ, gray levels", [(s, noise(s)) for s in (2, 5, 10, 20, 40)]),
        ("Crop", "border removed, %", [(p, crop(p / 100)) for p in (5, 10, 20, 30)]),
    ]

def measure_variants(tool, base, ref):
    """{algorithm: {transform: [(strength, value)]}}: every edit of `base` (a Pillow
    image), of both groups of transforms.py, compared with `ref`, the same pixels saved
    losslessly.

    Every variant is written to a file, and `site_stages measure` hashes them with the
    library and compares each variant with the reference, by each algorithm's own
    metric; `value` is None where the comparison does not apply.
    """
    with tempfile.TemporaryDirectory() as tmp:
        files, keys = [], []
        for name, _, steps in transforms() + content_edits():
            for strength, op in steps:
                files.append(write_edit(op, base, os.path.join(tmp, str(len(files)))))
                keys.append((name, strength))
        rows = run_lines(tool, "measure", ref, *files)
    result = {}
    for (name, strength), row in zip(keys, rows):
        for algo, value in row.items():
            if algo != "file":
                result.setdefault(algo, {}).setdefault(name, []).append((strength, value))
    return result


def measure_robustness(tool, image):
    """measure_variants() for one image file."""
    base = Image.open(image).convert("RGB")
    with tempfile.TemporaryDirectory() as tmp:
        ref = os.path.join(tmp, "reference.ppm")
        base.save(ref)
        return measure_variants(tool, base, ref)

def download(row, unreachable):
    """The photograph's bytes as the manifest names them, from the copy of the corpus or
    else from Commons; (None, why) when neither has them. A source that cannot be reached
    is added to `unreachable` and not tried again."""
    why = "no source"
    for source in (MIRROR + row["sha256"] + extension(row), row["url"]):
        host = urllib.parse.urlsplit(source).netloc
        if host in unreachable:
            continue
        try:
            data = get(source)
        except urllib.error.HTTPError as e:
            why = f"{host}: {e}"
            continue
        except OSError as e:
            unreachable[host] = str(e)
            why = f"cannot reach {host} ({e})"
            continue
        finally:
            if source == row["url"]:
                time.sleep(0.2)  # Commons asks clients to pace their requests
        if hashlib.sha256(data).hexdigest() == row["sha256"]:
            return data, None
        why = f"{host}: the download differs from the manifest"
    return None, why


def unpack_archive(rows, d):
    """Writes the rows' photographs that the copy's archive holds, verified, into d.
    An archive that cannot be downloaded or read writes nothing; the files are then
    downloaded one by one."""
    wanted = {row["sha256"] + extension(row): row for row in rows}
    try:
        req = urllib.request.Request(MIRROR_ARCHIVE, headers={"User-Agent": USER_AGENT})
        with urllib.request.urlopen(req, timeout=300) as r, \
                tarfile.open(fileobj=r, mode="r|gz") as tar:
            for member in tar:
                row = wanted.get(os.path.basename(member.name))
                if row is None or not member.isfile():
                    continue
                data = tar.extractfile(member).read()
                if hashlib.sha256(data).hexdigest() != row["sha256"]:
                    continue
                path = os.path.join(d, local_name(row))
                with open(path + ".part", "wb") as f:
                    f.write(data)
                os.replace(path + ".part", path)
    except (OSError, tarfile.TarError) as e:
        print(f"warning: photo corpus: the archive of the copy: {e}", file=sys.stderr)


_fetched = {}


def fetch(cache=None):
    """The verified local paths of the manifest's photographs, downloading what is missing.
    Returns (paths, problems); a photograph that is unavailable is a problem, not an error.
    One process downloads once: a second call returns the first one's answer."""
    d = cache_dir(cache)
    if d in _fetched:
        return _fetched[d]
    os.makedirs(d, exist_ok=True)

    def cached(row):
        path = os.path.join(d, local_name(row))
        if not os.path.exists(path):
            return False
        with open(path, "rb") as f:
            return hashlib.sha256(f.read()).hexdigest() == row["sha256"]

    rows = read_manifest()
    missing = [row for row in rows if not cached(row)]
    if len(missing) > 10:
        unpack_archive(missing, d)
    paths, problems, unreachable = [], [], {}
    for row in rows:
        path = os.path.join(d, local_name(row))
        if row not in missing or cached(row):
            paths.append(path)
            continue
        data, why = download(row, unreachable)
        if data is None:
            problems.append(f"{row['file']}: {why}")
            continue
        with open(path + ".part", "wb") as f:
            f.write(data)
        os.replace(path + ".part", path)
        paths.append(path)
    _fetched[d] = (paths, problems)
    return paths, problems

def measure_corpus(tool, corpus):
    """{"n", "label", "robust": {algo: {transform: [[strength, [value per image]]]}},
    "edits": {the same, for the edits that change the picture}, "different": {algo:
    [value per pair]}} for one corpus, from the cache when nothing that decides it has
    changed. Values are each algorithm's own metric, as `site_stages measure` and
    `site_stages pairs` print it; None where it does not apply."""
    paths = images(tool, corpus)

    def measure():
        with tempfile.TemporaryDirectory() as tmp:
            # The originals are decoded once, by Pillow, and written losslessly, so that
            # the reference an edit is compared with and the image that is edited are the
            # same pixels.
            refs = [os.path.join(tmp, f"{i:04d}.ppm") for i in range(len(paths))]
            with concurrent.futures.ProcessPoolExecutor() as pool:
                per_image = list(pool.map(_measure_one, [(tool, p, r) for p, r in
                                                         zip(paths, refs)]))
            rows = run_lines(tool, "pairs", *refs) if len(refs) >= 2 else []

        robust, edits = {}, {}
        content = {name for name, _, _ in content_edits()}
        for result in per_image:
            for algo, by_transform in result.items():
                for name, points in by_transform.items():
                    group = edits if name in content else robust
                    slot = group.setdefault(algo, {}).setdefault(
                        name, [[s, []] for s, _ in points])
                    for (strength, value), cell in zip(points, slot):
                        cell[1].append(value)
        different = {}
        for row in rows:
            for algo, value in row.items():
                if algo not in ("a", "b"):
                    different.setdefault(algo, []).append(value)
        return {"n": len(paths), "label": label(corpus, len(paths)), "robust": robust,
                "edits": edits, "different": different}

    return cached(corpus, cache_key(paths), measure, f"{corpus} corpus, {len(paths)} images")

# A "copy" is an original after one moderate edit: one strength of each transform, the
# kind of change a picture goes through between two places it is stored.
COPY_STRENGTHS = {
    "JPEG quality": 50,
    "Downscale": 0.5,
    "Rotation": 2,
    "Brightness": 1.15,
    "Contrast": 1.15,
    "Gamma": 1.2,
    "Gaussian blur": 1,
    "Noise": 5,
    "Crop": 5,
}

def separability(copies, different, lower_is_closer, recall=RECALL):
    """d', the threshold that keeps `recall` of the copies (RECALL unless a page asks for
    another share), and the share of different pairs that threshold also accepts.

    d' is the gap between the two means in units of their pooled spread, as
    tests/src/test_hash_properties.c defines it (separability()), signed so that a
    positive d' means copies are closer than different images.
    """
    c, d = np.asarray(copies, float), np.asarray(different, float)
    sign = 1.0 if lower_is_closer else -1.0
    pooled = np.sqrt((c.std() ** 2 + d.std() ** 2) / 2.0)
    dprime = sign * (d.mean() - c.mean()) / pooled if pooled > 0 else float("inf")
    # The threshold: the closest value that still accepts `recall` of the copies.
    t = sign * np.quantile(sign * c, recall, method="inverted_cdf")
    accepted = (d <= t) if lower_is_closer else (d >= t)
    return dprime, t, accepted.mean()

def content_edits():
    """The edits that change the picture, in the shape of transforms(); a strength is a
    number, or a word where the edit has no scale.

    A patch replaces a square at the center, covering the given share of the frame, with
    the square of the same size at the top-left corner of the same image: a local edit
    with content the image itself supplies. A hue rotation turns every pixel's hue by the
    given angle and keeps its saturation and value, a recoloring. The turns and the mirror
    change only the orientation.
    """
    def patch(area):
        def f(im):
            side = area ** 0.5
            w, h = round(im.width * side), round(im.height * side)
            x0, y0 = (im.width - w) // 2, (im.height - h) // 2
            out = im.copy()
            out.paste(im.crop((0, 0, w, h)), (x0, y0))
            return out, ("ppm", None)
        return f

    def hue(degrees):
        def f(im):
            hsv = np.asarray(im.convert("HSV")).copy()
            hsv[..., 0] = (hsv[..., 0].astype(int) + round(degrees * 256 / 360)) % 256
            return Image.fromarray(hsv, "HSV").convert("RGB"), ("ppm", None)
        return f

    def turn(method):
        return lambda im: (im.transpose(method), ("ppm", None))

    return [
        ("Local patch", "share of the frame, %", [(a, patch(a / 100)) for a in (1, 4, 9, 16)]),
        ("Hue rotation", "degrees", [(d, hue(d)) for d in (30, 90, 180)]),
        ("Turn and mirror", "", [("90°", turn(Image.Transpose.ROTATE_90)),
                                 ("180°", turn(Image.Transpose.ROTATE_180)),
                                 ("mirror", turn(Image.Transpose.FLIP_LEFT_RIGHT))]),
    ]

def moments(hexstr):
    """A ColorMoments digest as its nine numbers, [channel][moment]: each pair of bytes a
    signed 16-bit big-endian number, divided by 128 (PH_VECTOR16_SCALE)."""
    return (np.frombuffer(bytes.fromhex(hexstr), ">i2").astype(float) / 128.0).reshape(3, 3)


def l2(a, b):
    """ph_l2_distance() on two decoded digests."""
    return float(np.sqrt(((np.asarray(a) - np.asarray(b)) ** 2).sum()))

def standardized(v):
    """Per channel, the third central moment over σ³ (the standardized skewness, which
    does not depend on the channel's spread), from the decoded skewness s = cbrt(m3);
    `v` is [...][channel][moment]."""
    with np.errstate(divide="ignore", invalid="ignore"):
        return (v[..., 2] / v[..., 1]) ** 3


def skew_measurement(tool):
    """Per corpus: (originals [image][channel][moment], copies [image][edit][channel]
    [moment]), from the digests the library computes for every original and its nine
    copies (corpus.measure_settings)."""
    measured = corpus.measure_settings(tool, "color_moments-digests")
    out = {}
    for key, d in measured.items():
        if d["n"]:
            v = np.array([[moments(row["digest"]) for row in img] for img in d["images"]])
            out[key] = (v[:, 0], v[:, 1:])
    return out

def write_copies(base, tmp, stem):
    """The copies of the separability measurement: one moderate strength of each of the
    nine edits (COPY_STRENGTHS), saved as measure/robustness.py saves them; their paths."""
    return [write_edit(dict(steps)[COPY_STRENGTHS[name]], base,
                       os.path.join(tmp, f"{stem}-{k}"))
            for k, (name, _, steps) in enumerate(transforms())]


def _variants_one(args):
    """One original and its copies, through `site_stages <mode>`: one row per file."""
    tool, mode, path, tmp, stem = args
    base = Image.open(path).convert("RGB")
    ref = os.path.join(tmp, f"{stem}.ppm")
    base.save(ref)
    files = [ref] + write_copies(base, tmp, stem)
    rows = run_lines(tool, mode, *files)
    for f in files:
        os.remove(f)
    return [{k: v for k, v in row.items() if k != "file"} for row in rows]


def measure_settings(tool, mode):
    """{corpus: {"n", "label", "images": [[original, copy…] rows per image]}}: every
    original of both corpora and its nine copies through `site_stages <mode>`, which
    prints one JSON row of digests per file, one digest per setting it tries. Cached like
    the corpora themselves, as <mode>-<corpus>."""
    result = {}
    for corpus in CORPORA:
        paths = images(tool, corpus)
        rows = []
        if paths:
            def measure():
                with tempfile.TemporaryDirectory() as tmp, \
                        concurrent.futures.ProcessPoolExecutor() as pool:
                    return list(pool.map(_variants_one, [(tool, mode, p, tmp, f"{i:04d}")
                                                         for i, p in enumerate(paths)]))
            rows = cached(f"{mode}-{corpus}", cache_key(paths), measure,
                          f"{mode}, {len(paths)} images and their copies")
        result[corpus] = {"n": len(paths), "label": label(corpus, len(paths)), "images": rows}
    return result