1// Aseprite Render Library
2// Copyright (c) 2001-2017 David Capello
3//
4// This file is released under the terms of the MIT license.
5// Read LICENSE.txt for more information.
6
7#ifndef RENDER_MEDIAN_CUT_H_INCLUDED
8#define RENDER_MEDIAN_CUT_H_INCLUDED
9#pragma once
10
11#include "doc/color.h"
12
13#include <list>
14#include <queue>
15
16namespace render {
17
18 template<class Histogram>
19 class Box {
20
21 // These classes are used as parameters for some Box's generic
22 // member functions, so we can access to a different axis using
23 // the same generic function (i=Red channel in RAxisGetter, etc.).
24 struct RAxisGetter { static std::size_t at(const Histogram& h, int i, int j, int k, int l) { return h.at(i, j, k, l); } };
25 struct GAxisGetter { static std::size_t at(const Histogram& h, int i, int j, int k, int l) { return h.at(j, i, k, l); } };
26 struct BAxisGetter { static std::size_t at(const Histogram& h, int i, int j, int k, int l) { return h.at(j, k, i, l); } };
27 struct AAxisGetter { static std::size_t at(const Histogram& h, int i, int j, int k, int l) { return h.at(j, k, l, i); } };
28
29 // These classes are used as template parameter to split a Box
30 // along an axis (see splitAlongAxis)
31 struct RAxisSplitter {
32 static Box box1(const Box& box, int r) { return Box(box.r1, box.g1, box.b1, box.a1, r, box.g2, box.b2, box.a2); }
33 static Box box2(const Box& box, int r) { return Box(r, box.g1, box.b1, box.a1, box.r2, box.g2, box.b2, box.a2); }
34 };
35 struct GAxisSplitter {
36 static Box box1(const Box& box, int g) { return Box(box.r1, box.g1, box.b1, box.a1, box.r2, g, box.b2, box.a2); }
37 static Box box2(const Box& box, int g) { return Box(box.r1, g, box.b1, box.a1, box.r2, box.g2, box.b2, box.a2); }
38 };
39 struct BAxisSplitter {
40 static Box box1(const Box& box, int b) { return Box(box.r1, box.g1, box.b1, box.a1, box.r2, box.g2, b, box.a2); }
41 static Box box2(const Box& box, int b) { return Box(box.r1, box.g1, b, box.a1, box.r2, box.g2, box.b2, box.a2); }
42 };
43 struct AAxisSplitter {
44 static Box box1(const Box& box, int a) { return Box(box.r1, box.g1, box.b1, box.a1, box.r2, box.g2, box.b2, a ); }
45 static Box box2(const Box& box, int a) { return Box(box.r1, box.g1, box.b1, a, box.r2, box.g2, box.b2, box.a2); }
46 };
47
48 public:
49 Box(int r1, int g1, int b1, int a1,
50 int r2, int g2, int b2, int a2)
51 : r1(r1), g1(g1), b1(b1), a1(a1)
52 , r2(r2), g2(g2), b2(b2), a2(a2)
53 , points(0)
54 , volume(calculateVolume()) {
55 }
56
57 // Shrinks each plane of the box to a position where there are
58 // points in the histogram.
59 void shrink(const Histogram& histogram) {
60 axisShrink<RAxisGetter>(histogram, r1, r2, g1, g2, b1, b2, a1, a2);
61 axisShrink<GAxisGetter>(histogram, g1, g2, r1, r2, b1, b2, a1, a2);
62 axisShrink<BAxisGetter>(histogram, b1, b2, r1, r2, g1, g2, a1, a2);
63 axisShrink<AAxisGetter>(histogram, a1, a2, r1, r2, g1, g2, b1, b2);
64
65 // Calculate number of points inside the box (this is done by
66 // first time here, because the Box ctor didn't calculate it).
67 points = countPoints(histogram);
68
69 // Recalculate the volume (used in operator<).
70 volume = calculateVolume();
71 }
72
73 bool split(const Histogram& histogram, std::priority_queue<Box>& boxes) const {
74 // Split along the largest dimension of the box.
75 if ((r2-r1) >= (g2-g1) &&
76 (r2-r1) >= (b2-b1) &&
77 (r2-r1) >= (a2-a1)) {
78 return splitAlongAxis<RAxisGetter, RAxisSplitter>(histogram, boxes, r1, r2, g1, g2, b1, b2, a1, a2);
79 }
80
81 if ((g2-g1) >= (r2-r1) &&
82 (g2-g1) >= (b2-b1) &&
83 (g2-g1) >= (a2-a1)) {
84 return splitAlongAxis<GAxisGetter, GAxisSplitter>(histogram, boxes, g1, g2, r1, r2, b1, b2, a1, a2);
85 }
86
87 if ((b2-b1) >= (r2-r1) &&
88 (b2-b1) >= (g2-g1) &&
89 (b2-b1) >= (a2-a1)) {
90 return splitAlongAxis<BAxisGetter, BAxisSplitter>(histogram, boxes, b1, b2, r1, r2, g1, g2, a1, a2);
91 }
92
93 return splitAlongAxis<AAxisGetter, AAxisSplitter>(histogram, boxes, a1, a2, r1, r2, g1, g2, b1, b2);
94 }
95
96 // Returns the color enclosed by the box calculating the mean of
97 // all histogram's points inside the box.
98 uint32_t meanColor(const Histogram& histogram) const {
99 std::size_t r = 0, g = 0, b = 0, a = 0;
100 std::size_t count = 0;
101 int i, j, k, l;
102
103 for (i=r1; i<=r2; ++i)
104 for (j=g1; j<=g2; ++j)
105 for (k=b1; k<=b2; ++k)
106 for (l=a1; l<=a2; ++l) {
107 int c = histogram.at(i, j, k, l);
108 r += c * i;
109 g += c * j;
110 b += c * k;
111 a += c * l;
112 count += c;
113 }
114
115 // No colors in the box? This should not be possible.
116 ASSERT(count > 0 && "Box without histogram points, you must fill the histogram before using this function.");
117 if (count == 0)
118 return doc::rgba(0, 0, 0, 255);
119
120 // Calculate the mean. We have to do this before the *255
121 // multiplication to avoid a 32-bit overflow. E.g. Alpha channel
122 // is the most proper to overflow the 32-bit capacity in case
123 // all pixels are opaque.
124 r /= count;
125 g /= count;
126 b /= count;
127 a /= count;
128
129 return doc::rgba(int(255 * r / (Histogram::RElements-1)),
130 int(255 * g / (Histogram::GElements-1)),
131 int(255 * b / (Histogram::BElements-1)),
132 int(255 * a / (Histogram::AElements-1)));
133 }
134
135 // The boxes will be sort in the priority_queue by volume.
136 bool operator<(const Box& other) const {
137 return volume < other.volume;
138 }
139
140 private:
141
142 // Calculates the volume from the current box's dimensions. The
143 // value returned by this function is cached in the "volume"
144 // variable member of Box class to avoid multiplying several
145 // times.
146 int calculateVolume() const {
147 return (r2-r1+1) * (g2-g1+1) * (b2-b1+1) * (a2-a1+1);
148 }
149
150 // Returns the number of histogram's points inside the box bounds.
151 std::size_t countPoints(const Histogram& histogram) const {
152 std::size_t count = 0;
153 int i, j, k, l;
154
155 for (i=r1; i<=r2; ++i)
156 for (j=g1; j<=g2; ++j)
157 for (k=b1; k<=b2; ++k)
158 for (l=a1; l<=a2; ++l)
159 count += histogram.at(i, j, k, l);
160
161 return count;
162 }
163
164 // Reduces the specified side of the box (i1/i2) along the
165 // specified axis (if AxisGetter is RAxisGetter, then i1=r1,
166 // i2=r2; if AxisGetter is GAxisGetter, then i1=g1, i2=g2).
167 template<class AxisGetter>
168 static void axisShrink(const Histogram& histogram,
169 int& i1, int& i2,
170 const int& j1, const int& j2,
171 const int& k1, const int& k2,
172 const int& l1, const int& l2)
173 {
174 int j, k, l;
175
176 // Shrink i1.
177 for (; i1<i2; ++i1) {
178 for (j=j1; j<=j2; ++j) {
179 for (k=k1; k<=k2; ++k) {
180 for (l=l1; l<=l2; ++l) {
181 if (AxisGetter::at(histogram, i1, j, k, l) > 0)
182 goto doneA;
183 }
184 }
185 }
186 }
187
188 doneA:;
189
190 for (; i2>i1; --i2) {
191 for (j=j1; j<=j2; ++j) {
192 for (k=k1; k<=k2; ++k) {
193 for (l=l1; l<=l2; ++l) {
194 if (AxisGetter::at(histogram, i2, j, k, l) > 0)
195 goto doneB;
196 }
197 }
198 }
199 }
200
201 doneB:;
202 }
203
204 // Splits the box in two sub-boxes (if it's possible) along the
205 // specified axis by AxisGetter template parameter and "i1/i2"
206 // arguments. Returns true if the split was done and the "boxes"
207 // queue contains the new two sub-boxes resulting from the split
208 // operation.
209 template<class AxisGetter, class AxisSplitter>
210 bool splitAlongAxis(const Histogram& histogram,
211 std::priority_queue<Box>& boxes,
212 const int& i1, const int& i2,
213 const int& j1, const int& j2,
214 const int& k1, const int& k2,
215 const int& l1, const int& l2) const {
216 // These two variables will be used to count how many points are
217 // in each side of the box if we split it in "i" position.
218 std::size_t totalPoints1 = 0;
219 std::size_t totalPoints2 = this->points;
220 int i, j, k, l;
221
222 // We will try to split the box along the "i" axis. Imagine a
223 // plane which its normal vector is "i" axis, so we will try to
224 // move this plane from "i1" to "i2" to find the median, where
225 // the number of points in both sides of the plane are
226 // approximated the same.
227 for (i=i1; i<=i2; ++i) {
228 std::size_t planePoints = 0;
229
230 // We count all points in "i" plane.
231 for (j=j1; j<=j2; ++j)
232 for (k=k1; k<=k2; ++k)
233 for (l=l1; l<=l2; ++l)
234 planePoints += AxisGetter::at(histogram, i, j, k, l);
235
236 // As we move the plane to split through "i" axis One side is getting more points,
237 totalPoints1 += planePoints;
238 totalPoints2 -= planePoints;
239
240 if (totalPoints1 > totalPoints2) {
241 if (totalPoints2 > 0) {
242 Box box1(AxisSplitter::box1(*this, i));
243 Box box2(AxisSplitter::box2(*this, i+1));
244 box1.points = totalPoints1;
245 box2.points = totalPoints2;
246 boxes.push(box1);
247 boxes.push(box2);
248 return true;
249 }
250 else if (totalPoints1-planePoints > 0) {
251 Box box1(AxisSplitter::box1(*this, i-1));
252 Box box2(AxisSplitter::box2(*this, i));
253 box1.points = totalPoints1-planePoints;
254 box2.points = totalPoints2+planePoints;
255 boxes.push(box1);
256 boxes.push(box2);
257 return true;
258 }
259 else
260 return false;
261 }
262 }
263 return false;
264 }
265
266 int r1, g1, b1, a1; // Min point (closest to origin)
267 int r2, g2, b2, a2; // Max point
268 std::size_t points; // Number of points in the space which enclose this box
269 int volume;
270 }; // end of class Box
271
272 // Median Cut Algorithm as described in P. Heckbert, "Color image
273 // quantization for frame buffer display,", Computer Graphics,
274 // 16(3), pp. 297-307 (1982)
275 template<class Histogram>
276 void median_cut(const Histogram& histogram, std::size_t maxBoxes, std::vector<uint32_t>& result) {
277 // We need a priority queue to split bigger boxes first (see Box::operator<).
278 std::priority_queue<Box<Histogram> > boxes;
279
280 // First we start with one big box containing all histogram's samples.
281 boxes.push(Box<Histogram>(0, 0, 0, 0,
282 Histogram::RElements-1,
283 Histogram::GElements-1,
284 Histogram::BElements-1,
285 Histogram::AElements-1));
286
287 // Then we split each box until we reach the maximum specified by
288 // the user (maxBoxes) or until there aren't more boxes to split.
289 while (!boxes.empty() && boxes.size() < maxBoxes) {
290 // Get and remove the first (bigger) box to process from "boxes" queue.
291 Box<Histogram> box(boxes.top());
292 boxes.pop();
293
294 // Shrink the box to the minimum, to enclose the same points in
295 // the histogram.
296 box.shrink(histogram);
297
298 // Try to split the box along the largest axis.
299 if (!box.split(histogram, boxes)) {
300 // If we were not able to split the box (maybe because it is
301 // too small or there are not enough points to split it), then
302 // we add the box's color to the "result" vector directly (the
303 // box is not in the queue anymore).
304 if (result.size() < maxBoxes)
305 result.push_back(box.meanColor(histogram));
306 else
307 return;
308 }
309 }
310
311 // When we reach the maximum number of boxes, we convert each box
312 // to a color for the "result" vector.
313 while (!boxes.empty() && result.size() < maxBoxes) {
314 const Box<Histogram>& box(boxes.top());
315 doc::color_t color = box.meanColor(histogram);
316 result.push_back(color);
317 boxes.pop();
318 }
319 }
320
321} // namespace render
322
323#endif
324