画像のハーフトーン化
ハーフトーン化（halftoning）は、グレースケール画像を白黒2値の点の画像に変換する画像処理です。 新聞の写真やモノクロプリンタの印刷など、白か黒かの2値しか表現できない媒体で濃淡を表現するために使われています。 細かい白黒の点は、離れて見ると目の中で平均され、中間の明るさに見えます。
人の目に見える誤差を小さくする2値化の方法として、DBS（Direct Binary Search）がよく知られています。 DBS が小さくする誤差は、そのまま QUBO のエネルギーとして書けます。 このページでは、$256\times 256$ 画素の猫の写真を、DBS と QUBO ソルバー ABS3 で2値化して比べます。
動画のプログラムと章ごとの説明は動画：画像のハーフトーン化のページにあります。
目で見た誤差
$n_r \times n_c$ 画素のグレースケール画像 $I$（画素値 $0 \le I_{ij} \le 255$）が与えられたとします。 2値画像の各画素にバイナリ変数 $x_{ij} \in \lbrace0,1\rbrace$ を割り当てます（$1$ が白、$0$ が黒）。
人の目のぼけは、$(2M{+}1)\times(2M{+}1)$ のフィルタ $G$ との畳み込みで表します （このページでは $M=3$、係数の総和が $255$ の $7\times 7$ 整数ガウシアンを使います）。 2値画像 $x$ をぼかした値は
\[(G * x)_{ij} = \sum_{k=0}^{2M} \sum_{l=0}^{2M} G_{kl}\, x_{i+k-M,\, j+l-M}\]です（画像の外側の画素は $0$ とみなします）。$G$ の係数の総和が $255$ なので、画像の内部では $(G * x)_{ij}$ は $0$〜$255$ の値をとり、目で見た2値画像の明るさを表します。
目標値には元の画像そのものを使います:
\[T_{ij} = \mathrm{round}\!\left(\frac{I_{ij}\, m_{ij}}{255}\right), \qquad m_{ij} = (G * \mathbf{1})_{ij}\]ここで $m_{ij}$ は画素 $(i,j)$ に届くフィルタ係数の合計（$\mathbf{1}$ は全画素が $1$ の画像）です。 画像の内部では $m_{ij} = 255$、すなわち $T_{ij} = I_{ij}$ です。 画像の端ではフィルタの窓が画像からはみ出し、ぼかした値の最大が $m_{ij}$ に下がるので、目標も同じ比率で縮めています。 元の画像もぼかして目標にする方法も考えられますが、ぼかさずにそのまま目標にする方が、輪郭のはっきりした2値画像になります。
目で見た誤差は、ぼかした2値画像と目標値の差の2乗の和です:
\[E(x) = \sum_{i,j} \bigl( (G * x)_{ij} - T_{ij} \bigr)^2\]DBS
DBS は、$E(x)$ を小さくする2値画像を次の手順で求めます。
- 各画素の値を、ランダムなしきい値（$0$〜$254$）より明るければ白、そうでなければ黒にした画像から始めます。
- 画素を左上から右下へ順に（ラスタ順に）見て、その画素の反転と、色の違う隣の 8 画素との入れ替えを試し、$E(x)$ が一番下がるものを選んで実行します。どれでも下がらなければ何もしません。
- 画像全体の走査を、どの画素でも $E(x)$ が下がらなくなるまで繰り返します。
猫の写真では、走査 13 回で止まりました（このページの C++ のプログラムで約 0.2 秒）。 入れ替えを使わず反転だけで同じことをすると、$E(x)$ は $7.0$〜$7.3 \times 10^6$ で止まり、入れ替えありの $4.2 \times 10^6$ よりずっと悪くなります。 白い点を隣へずらすには反転が2回要り、1回目の反転では誤差が増えてしまうからです。
QUBO として表す
$(G * x)_{ij}$ は変数の1次式なので、$E(x)$ は変数の2次式、すなわち制約のない QUBO です。 QUBO のエネルギーが、目で見た誤差の2乗の和そのものになります。 DBS は、この QUBO の解を、1変数（反転）または2変数（入れ替え）ずつ変えて、エネルギーが下がらなくなるまで探す局所探索です。 同じ QUBO を QUBO ソルバーに渡せば、DBS と同じ基準で、より良い2値画像を探せます。
変数と項
$256\times 256$ の画像では、変数は画素の数と同じ $65{,}536$ 個です。 2乗を展開して現れる2つの変数の積 $x_{ij}\, x_{kl}$ は、画素 $(i,j)$ と $(k,l)$ が同じフィルタの窓に入る組だけです。 内部の画素の相手は、上下左右に $2M = 6$ 画素以内のうち $G$ の窓が重なる $128$ 個で、QUBO の項は $4{,}171{,}656$ 個になります。 項の数は画素の数に比例する、疎な QUBO です。
QUBO++ プログラム
入力の画像 halftone_cat.pgm（PGM 形式、$256 \times 256$ 画素）を、実行するディレクトリに置いて実行します。 このプログラムは、DBS の2値画像、ABS3 だけで探した2値画像、DBS の2値画像を hint（ヒント解）として ABS3 に与えて探した2値画像の3つを求めます:
#include <qbpp/qbpp.hpp>
#include <qbpp/abs3_solver.hpp>
#include <algorithm>
#include <fstream>
#include <random>
#include <string>
#include <vector>
constexpr int M = 3; // フィルタ半径（7x7 フィルタ）
// 整数ガウシアン（総和 255）
constexpr int G[2 * M + 1][2 * M + 1] = {{0, 0, 1, 1, 1, 0, 0},
{0, 1, 5, 7, 5, 1, 0},
{1, 5, 14, 21, 14, 5, 1},
{1, 7, 21, 31, 21, 7, 1},
{1, 5, 14, 21, 14, 5, 1},
{0, 1, 5, 7, 5, 1, 0},
{0, 0, 1, 1, 1, 0, 0}};
using Image = std::vector<std::vector<int>>;
int ROWS, COLS;
Image img, target;
bool inside(int i, int j) { return 0 <= i && i < ROWS && 0 <= j && j < COLS; }
// G との畳み込み（画像の外側の画素は 0 とみなす）
Image blur(const Image& a) {
Image out(ROWS, std::vector<int>(COLS, 0));
for (int i = 0; i < ROWS; ++i)
for (int j = 0; j < COLS; ++j)
for (int k = 0; k <= 2 * M; ++k)
for (int l = 0; l <= 2 * M; ++l)
if (inside(i + k - M, j + l - M))
out[i][j] += G[k][l] * a[i + k - M][j + l - M];
return out;
}
// 画素 (i, j) と (i2, j2) のぼけが重なる量: 画像の中の q についての G(q - p) G(q - p2) の和
int overlap(int i, int j, int i2, int j2) {
int s = 0;
for (int k = 0; k <= 2 * M; ++k)
for (int l = 0; l <= 2 * M; ++l) {
int k2 = i + k - i2, l2 = j + l - j2;
if (inside(i + k - M, j + l - M) && 0 <= k2 && k2 <= 2 * M && 0 <= l2 &&
l2 <= 2 * M)
s += G[k][l] * G[k2][l2];
}
return s;
}
// ぼけが画像からはみ出さない画素どうしの overlap は位置の差だけで決まるので、先に求めておく
int INNER[4 * M + 1][4 * M + 1];
int ov(int i, int j, int i2, int j2) {
if (M <= std::min(i, i2) && std::max(i, i2) < ROWS - M &&
M <= std::min(j, j2) && std::max(j, j2) < COLS - M)
return INNER[i2 - i + 2 * M][j2 - j + 2 * M];
return overlap(i, j, i2, j2);
}
// DBS: 反転と隣の 8 画素との入れ替えのうち、誤差の 2 乗の和が一番下がるものを選ぶ。
// どの画素でも下がらなくなるまで、ラスタ順の走査を繰り返す
Image dbs(unsigned seed = 1) {
std::mt19937 rng(seed);
std::uniform_int_distribution<int> threshold(0, 254);
Image b(ROWS, std::vector<int>(COLS));
for (int i = 0; i < ROWS; ++i) // ランダムなしきい値で始める
for (int j = 0; j < COLS; ++j) b[i][j] = img[i][j] > threshold(rng);
Image err = blur(b);
for (int i = 0; i < ROWS; ++i)
for (int j = 0; j < COLS; ++j) err[i][j] -= target[i][j];
Image c = blur(err); // c[i][j]: 画素 (i, j) のぼけと誤差の内積
auto delta = [&](int i, int j) { // 画素 (i, j) を反転したときの E の変化
return 2 * (1 - 2 * b[i][j]) * c[i][j] + ov(i, j, i, j);
};
auto flip = [&](int i, int j) {
int d = 1 - 2 * b[i][j];
b[i][j] ^= 1;
for (int i2 = std::max(0, i - 2 * M); i2 < std::min(ROWS, i + 2 * M + 1); ++i2)
for (int j2 = std::max(0, j - 2 * M); j2 < std::min(COLS, j + 2 * M + 1); ++j2)
c[i2][j2] += d * ov(i, j, i2, j2);
};
for (;;) {
int changed = 0;
for (int i = 0; i < ROWS; ++i)
for (int j = 0; j < COLS; ++j) {
int d0 = delta(i, j);
int best = d0, pi = -1, pj = -1; // 反転
for (int i2 = i - 1; i2 <= i + 1; ++i2) // 色の違う隣の画素との入れ替え
for (int j2 = j - 1; j2 <= j + 1; ++j2)
if (inside(i2, j2) && b[i2][j2] != b[i][j]) {
int de = d0 + delta(i2, j2) - 2 * ov(i, j, i2, j2);
if (de < best) best = de, pi = i2, pj = j2;
}
if (best < 0) {
flip(i, j);
if (pi >= 0) flip(pi, pj);
++changed;
}
}
if (changed == 0) return b;
}
}
// 2 値画像を PGM（P2 形式）で保存
template <typename Bits>
void save(const std::string& name, const Bits& bits) {
std::ofstream ofs(name);
ofs << "P2\n" << COLS << " " << ROWS << "\n255\n";
for (int i = 0; i < ROWS; ++i) {
for (int j = 0; j < COLS; ++j) ofs << 255 * int(bits[i][j]) << " ";
ofs << "\n";
}
}
int main() {
// 入力画像: PGM（P2 形式）のグレースケール画像
std::ifstream in("halftone_cat.pgm");
std::string magic;
int maxval;
in >> magic >> COLS >> ROWS >> maxval;
img.assign(ROWS, std::vector<int>(COLS));
for (auto& row : img)
for (auto& v : row) in >> v;
// 目標値 T: 元画像そのもの（境界は届くフィルタ質量に合わせて縮小）
Image mass = blur(Image(ROWS, std::vector<int>(COLS, 1)));
target.assign(ROWS, std::vector<int>(COLS));
for (int i = 0; i < ROWS; ++i)
for (int j = 0; j < COLS; ++j)
target[i][j] = (img[i][j] * mass[i][j] + 127) / 255;
// E = sum(((G * x) - T)^2)
auto x = qbpp::var("x", ROWS, COLS);
auto zc = qbpp::expr(ROWS, M);
auto zr = qbpp::expr(M, COLS + 2 * M);
auto xw = qbpp::concat(qbpp::concat(zc, x, 1), zc, 1); // 左右にゼロを連結
auto xp = qbpp::concat(qbpp::concat(zr, xw, 0), zr, 0); // 上下にゼロを連結
auto conv = qbpp::expr(ROWS, COLS); // (G * x)
for (int k = 0; k <= 2 * M; ++k)
for (int l = 0; l <= 2 * M; ++l)
if (G[k][l])
conv +=
G[k][l] * xp(qbpp::slice(k, k + ROWS), qbpp::slice(l, l + COLS));
auto f = qbpp::toExpr(0);
for (int i = 0; i < ROWS; ++i)
f += qbpp::sum((conv(i) - qbpp::array(target[i])).sqr());
f.simplify_as_binary();
// DBS の 2 値画像を求め、hint にする
for (int di = -2 * M; di <= 2 * M; ++di)
for (int dj = -2 * M; dj <= 2 * M; ++dj)
INNER[di + 2 * M][dj + 2 * M] = overlap(2 * M, 2 * M, 2 * M + di, 2 * M + dj);
Image dbs_img = dbs();
qbpp::Sol hint(f);
for (int i = 0; i < ROWS; ++i)
for (int j = 0; j < COLS; ++j) hint.set(x[i][j], dbs_img[i][j]);
std::cout << "DBS energy = " << f(hint) << std::endl;
qbpp::ABS3Solver solver(f);
auto sol = solver.search({{"time_limit", 60.0}});
std::cout << "variables = " << sol.info().get("var_count")
<< " terms = " << sol.info().get("term_count") << std::endl;
std::cout << "ABS3 energy = " << sol.energy() << std::endl;
qbpp::Params params({{"time_limit", 60.0}});
params.hint(hint);
auto sol_hint = solver.search(params);
std::cout << "ABS3 + hint energy = " << sol_hint.energy() << std::endl;
save("dbs.pgm", dbs_img);
save("abs3.pgm", sol(x));
save("abs3_hint.pgm", sol_hint(x));
}
QUBO の組み立て
f は配列演算で組み立てています。
qbpp::concatでxの上下左右に幅 $M$ のゼロの配列（qbpp::expr(...)はゼロで初期化されたExprの配列）を連結し、$(n_r{+}2M) \times (n_c{+}2M)$ の配列xpを作ります。xpから縦横に $k, l$ だけずらした $n_r \times n_c$ の部分配列xp(qbpp::slice(k, k + ROWS), qbpp::slice(l, l + COLS))を $G_{kl}$ 倍してconvに足します。二重ループが終わると、convは1次式の配列 $(G * x)$ になります。- 行ごとに目標値
target[i]を引いて要素ごとにsqr()で2乗し、qbpp::sumで合計してfに足すと $E(x)$ になります。simplify_as_binary()が同類項をまとめ、$x^2 = x$ の規則で式を整理します。
DBS
DBS の中では $E(x)$ を毎回計算し直さず、画素を反転したときの変化だけを求めます。 画素 $p$ のぼけ（$p$ だけが白の画像をぼかしたもの）を $g_p$、誤差の画像を $e = (G * x) - T$ とします。 画素 $p$ を $d = \pm 1$ だけ変えると誤差は $d\, g_p$ だけ変わるので、
\[\Delta E = 2d\, \langle g_p, e\rangle + \langle g_p, g_p\rangle\]です（$\langle\cdot,\cdot\rangle$ は画像の中の画素についての積の和）。 c[i][j] が $\langle g_p, e\rangle$、overlap(p, p2) が $\langle g_p, g_{p_2}\rangle$ です。 入れ替えは色の違う2画素の反転なので、2つの反転の変化の和から重なりの分 $2\langle g_p, g_{p_2}\rangle$ を引いたものが変化になります。 変更を実行したら、flip() が上下左右 $2M$ 画素以内の c を直します。
hint
qbpp::Sol hint(f) に hint.set(x[i][j], ...) で DBS の2値画像を変数の値として入れ、解 hint を作ります。 f(hint) は DBS の2値画像の誤差の2乗の和です。 params.hint(hint) を渡した search(params) は、hint を出発点にして探索します。
出力結果
DBS energy = 4232582
variables = 65536 terms = 4171656
ABS3 energy = 3913250
ABS3 + hint energy = 3917986
ABS3 は乱択の探索なので、energy の値は実行ごとに変わります（NVIDIA RTX 3090 を 2 枚積んだ計算機での例）。 DBS は約 0.2 秒で終わります。
3 つの比較
時間と誤差の2乗の和（エネルギー）の変化を、DBS、ABS3 だけ、DBS の2値画像を hint にした ABS3 で比べました（上と同じ計算機。それぞれ何回か実行したうちの1回です）。グラフの DBS は、変化を毎回フィルタの窓の中で計算し直す C++ の実装で測っていて、このページのプログラムより遅く 0.78 秒かかっています。 hint を与えた ABS3 は、DBS が終わった後から探索を始めます:
- DBS は 0.78 秒で止まります。それ以上は誤差を下げられません。
- ABS3 だけの探索は、3 秒あまりで DBS の誤差を下回り、その後も下がり続けます。
- hint を与えると、DBS の解から探し始められるので、最初の数秒は ABS3 だけより良い解が得られます。60 秒後には、どちらもほぼ同じ値になりました。
画像の比較
元の画像と、3 つの2値画像です:
左から、元の画像、DBS（誤差の2乗の和 4,207,803）、ABS3 だけ（3,904,736）、DBS の画像を hint にした ABS3（3,903,575）です（画像とこれらの値は別の実行で求めたもので、上の出力結果とは少し違います）。 目の部分（$40\times 40$ 画素）を拡大すると:
誤差の2乗の和は ABS3 の方が 7% ほど小さいのですが、見た目では3つの差はほとんど分かりません。 実用には、速い DBS で十分です。
まとめ
ハーフトーン化の QUBO は、変数が画素の数だけあり、項の数が画素の数に比例する疎な QUBO です。 画像を大きくすれば規模を自由に変えられ、DBS という速くて良い解を出す方法とも比べられるので、QUBO ソルバーのベンチマークとして興味深い問題です。
点をまとめて打つ クラスタドット・ハーフトーン では、同じ誤差に制約を加えて解きます。