source: libpipi/trunk/pipi/dither/dbs.c @ 2902

Last change on this file since 2902 was 2902, checked in by Sam Hocevar, 11 years ago

Support C99 types on Win32 through the same hacks as in libcaca.

File size: 7.0 KB
Line 
1/*
2 *  libpipi       Pathetic image processing interface library
3 *  Copyright (c) 2004-2008 Sam Hocevar <sam@zoy.org>
4 *                All Rights Reserved
5 *
6 *  $Id$
7 *
8 *  This library is free software. It comes without any warranty, to
9 *  the extent permitted by applicable law. You can redistribute it
10 *  and/or modify it under the terms of the Do What The Fuck You Want
11 *  To Public License, Version 2, as published by Sam Hocevar. See
12 *  http://sam.zoy.org/wtfpl/COPYING for more details.
13 */
14
15/*
16 * dbs.c: Direct Binary Search dithering functions
17 */
18
19#include "config.h"
20
21#include <string.h>
22#include <stdlib.h>
23
24#include <math.h>
25
26#include "pipi.h"
27#include "pipi_internals.h"
28
29#define CELL 16
30
31#define N 7
32#define NN ((N * 2 + 1))
33
34/* FIXME: though the algorithm is supposed to stop, we do not have a real,
35 * guaranteed stop condition here. */
36
37pipi_image_t *pipi_dither_dbs(pipi_image_t *img)
38{
39    double kernel[NN * NN];
40    double t = 0.;
41    pipi_image_t *src, *dst, *tmp1, *tmp2;
42    pipi_pixels_t *dstp, *tmp1p, *tmp2p;
43    int *changelist;
44    float *dstdata, *tmp1data, *tmp2data;
45    int i, j, x, y, w, h, cw, ch;
46
47    /* Build our human visual system kernel. */
48    for(j = 0; j < NN; j++)
49        for(i = 0; i < NN; i++)
50        {
51            double a = (double)(i - N);
52            double b = (double)(j - N);
53            kernel[j * NN + i] =
54                    exp(-(a * a + b * b) / (2. * 1.6 * 1.6))
55                  + exp(-(a * a + b * b) / (2. * 0.6 * 0.6));
56            t += kernel[j * NN + i];
57        }
58
59    for(j = 0; j < NN; j++)
60        for(i = 0; i < NN; i++)
61            kernel[j * NN + i] /= t;
62
63    w = img->w;
64    h = img->h;
65
66    cw = (w + CELL - 1) / CELL;
67    ch = (h + CELL - 1) / CELL;
68    changelist = malloc(cw * ch * sizeof(int));
69    memset(changelist, 0, cw * ch * sizeof(int));
70
71    src = pipi_copy(img);
72    pipi_getpixels(src, PIPI_PIXELS_Y_F);
73
74    tmp1 = pipi_convolution(src, NN, NN, kernel);
75    tmp1p = pipi_getpixels(tmp1, PIPI_PIXELS_Y_F);
76    tmp1data = (float *)tmp1p->pixels;
77
78    dst = pipi_dither_random(src);
79    dstp = pipi_getpixels(dst, PIPI_PIXELS_Y_F);
80    dstdata = (float *)dstp->pixels;
81
82    pipi_free(src);
83
84    tmp2 = pipi_convolution(dst, NN, NN, kernel);
85    tmp2p = pipi_getpixels(tmp2, PIPI_PIXELS_Y_F);
86    tmp2data = (float *)tmp2p->pixels;
87
88    for(;;)
89    {
90        int cx, cy, n;
91        int allchanges = 0;
92
93        for(cy = 0; cy < ch; cy++)
94        for(cx = 0; cx < cw; cx++)
95        {
96            int changes = 0;
97
98            if(changelist[cy * cw + cx] >= 2)
99                continue;
100
101            for(y = cy * CELL; y < (cy + 1) * CELL; y++)
102            for(x = cx * CELL; x < (cx + 1) * CELL; x++)
103            {
104                double d, d2, e, best = 0.;
105                int opx = -1, opy = -1;
106
107                d = dstdata[y * w + x];
108
109                /* Compute the effect of a toggle */
110                e = 0.;
111                for(j = -N; j < N + 1; j++)
112                {
113                    if(y + j < 0 || y + j >= h)
114                        continue;
115
116                    for(i = -N; i < N + 1; i++)
117                    {
118                        double m, p, q1, q2;
119
120                        if(x + i < 0 || x + i >= w)
121                            continue;
122
123                        m = kernel[(j + N) * NN + i + N];
124                        p = tmp1data[(y + j) * w + x + i];
125                        q1 = tmp2data[(y + j) * w + x + i];
126                        q2 = q1 - m * d + m * (1. - d);
127                        e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
128                    }
129                }
130                if(e > best)
131                {
132                    best = e;
133                    opx = opy = 0;
134                }
135
136                /* Compute the effect of swaps */
137                for(n = 0; n < 8; n++)
138                {
139                    static int const step[] =
140                      { 0, 1, 0, -1, -1, 0, 1, 0, -1, -1, -1, 1, 1, -1, 1, 1 };
141                    int idx = step[n * 2], idy = step[n * 2 + 1];
142                    if(y + idy < 0 || y + idy >= h
143                         || x + idx < 0 || x + idx >= w)
144                        continue;
145                    d2 = dstdata[(y + idy) * w + x + idx];
146                    if(d2 == d)
147                        continue;
148                    e = 0.;
149                    for(j = -N; j < N + 1; j++)
150                    {
151                        if(y + j < 0 || y + j >= h)
152                            continue;
153                        if(j - idy + N < 0 || j - idy + N >= NN)
154                            continue;
155                        for(i = -N; i < N + 1; i++)
156                        {
157                            double ma, mb, p, q1, q2;
158                            if(x + i < 0 || x + i >= w)
159                                continue;
160                            if(i - idx + N < 0 || i - idx + N >= NN)
161                                continue;
162                            ma = kernel[(j + N) * NN + i + N];
163                            mb = kernel[(j - idy + N) * NN + i - idx + N];
164                            p = tmp1data[(y + j) * w + x + i];
165                            q1 = tmp2data[(y + j) * w + x + i];
166                            q2 = q1 - ma * d + ma * d2 - mb * d2 + mb * d;
167                            e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
168                        }
169                    }
170                    if(e > best)
171                    {
172                        best = e;
173                        opx = idx;
174                        opy = idy;
175                    }
176                }
177
178                /* Apply the change if interesting */
179                if(best <= 0.)
180                    continue;
181                if(opx || opy)
182                {
183                    d2 = dstdata[(y + opy) * w + x + opx];
184                    dstdata[(y + opy) * w + x + opx] = d;
185                }
186                else
187                    d2 = 1. - d;
188                dstdata[y * w + x] = d2;
189                for(j = -N; j < N + 1; j++)
190                for(i = -N; i < N + 1; i++)
191                {
192                    double m = kernel[(j + N) * NN + i + N];
193                    if(y + j >= 0 && y + j < h
194                        && x + i >= 0 && x + i < w)
195                    {
196                        t = tmp2data[(y + j) * w + x + i];
197                        tmp2data[(y + j) * w + x + i] = t + m * (d2 - d);
198                    }
199                    if((opx || opy) && y + opy + j >= 0 && y + opy + j < h
200                                    && x + opx + i >= 0 && x + opx + i < w)
201                    {
202                        t = tmp2data[(y + opy + j) * w + x + opx + i];
203                        tmp2data[(y + opy + j) * w + x + opx + i]
204                                = t + m * (d - d2);
205                    }
206                }
207
208                changes++;
209            }
210
211            if(changes == 0)
212                changelist[cy * cw + cx]++;
213
214            allchanges += changes;
215        }
216
217        if(allchanges == 0)
218            break;
219    }
220
221    free(changelist);
222
223    pipi_free(tmp1);
224    pipi_free(tmp2);
225
226    return dst;
227}
228
Note: See TracBrowser for help on using the repository browser.