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

Last change on this file since 2672 was 2672, checked in by Sam Hocevar, 12 years ago
  • dbs.c: generate the initial halftone using random dithering instead of Floyd-Steinberg in order to avoid the energy level getting locally stuck around error diffusion structure artifacts.
File size: 7.0 KB
Line 
1/*
2 *  libpipi       Proper image processing implementation 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#include "common.h"
21
22#include <string.h>
23#include <stdlib.h>
24
25#include <math.h>
26
27#include "pipi.h"
28#include "pipi_internals.h"
29
30#define CELL 16
31
32#define N 7
33#define NN ((N * 2 + 1))
34
35/* FIXME: though the algorithm is supposed to stop, we do not have a real,
36 * guaranteed stop condition here. */
37
38pipi_image_t *pipi_dither_dbs(pipi_image_t *img)
39{
40    double kernel[NN * NN];
41    double t = 0.;
42    pipi_image_t *src, *dst, *tmp1, *tmp2;
43    pipi_pixels_t *dstp, *tmp1p, *tmp2p;
44    int *changelist;
45    float *dstdata, *tmp1data, *tmp2data;
46    int i, j, x, y, w, h, cw, ch;
47
48    /* Build our human visual system kernel. */
49    for(j = 0; j < NN; j++)
50        for(i = 0; i < NN; i++)
51        {
52            double a = (double)(i - N);
53            double b = (double)(j - N);
54            kernel[j * NN + i] =
55                    exp(-(a * a + b * b) / (2. * 1.6 * 1.6))
56                  + exp(-(a * a + b * b) / (2. * 0.6 * 0.6));
57            t += kernel[j * NN + i];
58        }
59
60    for(j = 0; j < NN; j++)
61        for(i = 0; i < NN; i++)
62            kernel[j * NN + i] /= t;
63
64    w = img->w;
65    h = img->h;
66
67    cw = (w + CELL - 1) / CELL;
68    ch = (h + CELL - 1) / CELL;
69    changelist = malloc(cw * ch * sizeof(int));
70    memset(changelist, 0, cw * ch * sizeof(int));
71
72    src = pipi_copy(img);
73    pipi_getpixels(src, PIPI_PIXELS_Y_F);
74
75    tmp1 = pipi_convolution(src, NN, NN, kernel);
76    tmp1p = pipi_getpixels(tmp1, PIPI_PIXELS_Y_F);
77    tmp1data = (float *)tmp1p->pixels;
78
79    dst = pipi_dither_random(src);
80    dstp = pipi_getpixels(dst, PIPI_PIXELS_Y_F);
81    dstdata = (float *)dstp->pixels;
82
83    pipi_free(src);
84
85    tmp2 = pipi_convolution(dst, NN, NN, kernel);
86    tmp2p = pipi_getpixels(tmp2, PIPI_PIXELS_Y_F);
87    tmp2data = (float *)tmp2p->pixels;
88
89    for(;;)
90    {
91        int cx, cy, n;
92        int allchanges = 0;
93
94        for(cy = 0; cy < ch; cy++)
95        for(cx = 0; cx < cw; cx++)
96        {
97            int changes = 0;
98
99            if(changelist[cy * cw + cx] >= 2)
100                continue;
101
102            for(y = cy * CELL; y < (cy + 1) * CELL; y++)
103            for(x = cx * CELL; x < (cx + 1) * CELL; x++)
104            {
105                double d, d2, e, best = 0.;
106                int opx = -1, opy = -1;
107
108                d = dstdata[y * w + x];
109
110                /* Compute the effect of a toggle */
111                e = 0.;
112                for(j = -N; j < N + 1; j++)
113                {
114                    if(y + j < 0 || y + j >= h)
115                        continue;
116
117                    for(i = -N; i < N + 1; i++)
118                    {
119                        double m, p, q1, q2;
120
121                        if(x + i < 0 || x + i >= w)
122                            continue;
123
124                        m = kernel[(j + N) * NN + i + N];
125                        p = tmp1data[(y + j) * w + x + i];
126                        q1 = tmp2data[(y + j) * w + x + i];
127                        q2 = q1 - m * d + m * (1. - d);
128                        e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
129                    }
130                }
131                if(e > best)
132                {
133                    best = e;
134                    opx = opy = 0;
135                }
136
137                /* Compute the effect of swaps */
138                for(n = 0; n < 8; n++)
139                {
140                    static int const step[] =
141                      { 0, 1, 0, -1, -1, 0, 1, 0, -1, -1, -1, 1, 1, -1, 1, 1 };
142                    int idx = step[n * 2], idy = step[n * 2 + 1];
143                    if(y + idy < 0 || y + idy >= h
144                         || x + idx < 0 || x + idx >= w)
145                        continue;
146                    d2 = dstdata[(y + idy) * w + x + idx];
147                    if(d2 == d)
148                        continue;
149                    e = 0.;
150                    for(j = -N; j < N + 1; j++)
151                    {
152                        if(y + j < 0 || y + j >= h)
153                            continue;
154                        if(j - idy + N < 0 || j - idy + N >= NN)
155                            continue;
156                        for(i = -N; i < N + 1; i++)
157                        {
158                            double ma, mb, p, q1, q2;
159                            if(x + i < 0 || x + i >= w)
160                                continue;
161                            if(i - idx + N < 0 || i - idx + N >= NN)
162                                continue;
163                            ma = kernel[(j + N) * NN + i + N];
164                            mb = kernel[(j - idy + N) * NN + i - idx + N];
165                            p = tmp1data[(y + j) * w + x + i];
166                            q1 = tmp2data[(y + j) * w + x + i];
167                            q2 = q1 - ma * d + ma * d2 - mb * d2 + mb * d;
168                            e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
169                        }
170                    }
171                    if(e > best)
172                    {
173                        best = e;
174                        opx = idx;
175                        opy = idy;
176                    }
177                }
178
179                /* Apply the change if interesting */
180                if(best <= 0.)
181                    continue;
182                if(opx || opy)
183                {
184                    d2 = dstdata[(y + opy) * w + x + opx];
185                    dstdata[(y + opy) * w + x + opx] = d;
186                }
187                else
188                    d2 = 1. - d;
189                dstdata[y * w + x] = d2;
190                for(j = -N; j < N + 1; j++)
191                for(i = -N; i < N + 1; i++)
192                {
193                    double m = kernel[(j + N) * NN + i + N];
194                    if(y + j >= 0 && y + j < h
195                        && x + i >= 0 && x + i < w)
196                    {
197                        t = tmp2data[(y + j) * w + x + i];
198                        tmp2data[(y + j) * w + x + i] = t + m * (d2 - d);
199                    }
200                    if((opx || opy) && y + opy + j >= 0 && y + opy + j < h
201                                    && x + opx + i >= 0 && x + opx + i < w)
202                    {
203                        t = tmp2data[(y + opy + j) * w + x + opx + i];
204                        tmp2data[(y + opy + j) * w + x + opx + i]
205                                = t + m * (d - d2);
206                    }
207                }
208
209                changes++;
210            }
211
212            if(changes == 0)
213                changelist[cy * cw + cx]++;
214
215            allchanges += changes;
216        }
217
218        if(allchanges == 0)
219            break;
220    }
221
222    free(changelist);
223
224    pipi_free(tmp1);
225    pipi_free(tmp2);
226
227    return dst;
228}
229
Note: See TracBrowser for help on using the repository browser.