1//
2//  PREqualize.m
3//  PRICE
4//  Image level equalization
5//
6//  Created by Riccardo Mottola on Fri Dec 05 2003.
7//  Copyright (c) 2003-2014 Carduus. All rights reserved.
8//
9// This application is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
10// This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
11
12#include <math.h>
13#include <limits.h>
14
15#import "PREqualize.h"
16
17#if defined (__SVR4) && defined (__sun)
18#define rintf(x) rint(x)
19#endif
20
21
22@implementation PREqualize
23
24- (PRImage *)filterImage:(PRImage *)image with:(NSArray *)parameters progressPanel:(PRCProgress *)progressPanel
25{
26    int space;
27
28    /* interpret the parameters */
29    space = [[parameters objectAtIndex:0] intValue];
30
31    return [self equalizeImage:image :space];
32}
33
34- (NSString *)actionName
35{
36    return @"Equalize";
37}
38
39- (PRImage *)equalizeImage:(PRImage *)srcImage :(int)space
40{
41  NSBitmapImageRep *srcImageRep;
42  PRImage          *destImage;
43  NSBitmapImageRep *destImageRep;
44  NSInteger        w, h;
45  NSInteger        x, y;
46  NSInteger        i;
47  unsigned char    *srcData;
48  unsigned char    *destData;
49  NSInteger srcSamplesPerPixel;
50  NSInteger destSamplesPerPixel;
51  register NSInteger srcBytesPerPixel;
52  register NSInteger destBytesPerPixel;
53  register NSInteger srcBytesPerRow;
54  register NSInteger destBytesPerRow;
55  int              pixNum;
56  BOOL             hasAlpha;
57
58  /* some trace */
59  NSLog(@"levels: %d", UCHAR_MAX);
60  NSLog(@"space: %d", space);
61
62    /* get source image representation and associated information */
63    srcImageRep = [srcImage bitmapRep];
64
65    w = [srcImageRep pixelsWide];
66    h = [srcImageRep pixelsHigh];
67    pixNum = h * w;
68    NSLog(@"pixels: %d", pixNum);
69    srcBytesPerRow = [srcImageRep bytesPerRow];
70    srcSamplesPerPixel = [srcImageRep samplesPerPixel];
71    srcBytesPerPixel = [srcImageRep bitsPerPixel] / 8;
72    destSamplesPerPixel = srcSamplesPerPixel;
73
74    hasAlpha = [srcImageRep hasAlpha];
75
76    /* allocate destination image and its representation */
77    destImage = [[PRImage alloc] initWithSize:NSMakeSize(w, h)];
78    destImageRep = [[NSBitmapImageRep alloc]
79                initWithBitmapDataPlanes:NULL
80                              pixelsWide:w
81                              pixelsHigh:h
82                           bitsPerSample:[srcImageRep bitsPerSample]
83                         samplesPerPixel:destSamplesPerPixel
84                                hasAlpha:[srcImage hasAlpha]
85                                isPlanar:NO
86                          colorSpaceName:[srcImageRep colorSpaceName]
87                             bytesPerRow:w*destSamplesPerPixel
88                            bitsPerPixel:0];
89    srcData = [srcImageRep bitmapData];
90    destData = [destImageRep bitmapData];
91    destBytesPerRow = [destImageRep bytesPerRow];
92    destBytesPerPixel = [destImageRep bitsPerPixel] / 8;
93
94    if ([srcImage hasColor])
95    {
96        if (space == COLOR_SPACE_RGB)
97        {
98            unsigned long int histogramDenormR[UCHAR_MAX+1]; /* not normalized pixel count for each level */
99            unsigned long int histogramDenormG[UCHAR_MAX+1]; /* not normalized pixel count for each level */
100            unsigned long int histogramDenormB[UCHAR_MAX+1]; /* not normalized pixel count for each level */
101            float histogramR[UCHAR_MAX+1];                   /* normalized histogram */
102            float histogramG[UCHAR_MAX+1];                   /* normalized histogram */
103            float histogramB[UCHAR_MAX+1];                   /* normalized histogram */
104            float cumulativeHistogramR[UCHAR_MAX+1];         /* cumulative histogram */
105            float cumulativeHistogramG[UCHAR_MAX+1];         /* cumulative histogram */
106            float cumulativeHistogramB[UCHAR_MAX+1];         /* cumulative histogram */
107
108            /* calculate the histogram */
109            for (i = 0; i <= UCHAR_MAX; i++)
110                histogramDenormR[i] = histogramDenormG[i] = histogramDenormB[i] =  0;
111            for (y = 0; y < h; y++)
112                for (x = 0; x < w*3; x += 3)
113                {
114                    histogramDenormR[srcData[y*srcBytesPerRow + srcBytesPerPixel*x]]++;
115                    histogramDenormG[srcData[y*srcBytesPerRow + srcBytesPerPixel*x + 1]]++;
116                    histogramDenormB[srcData[y*srcBytesPerRow + srcBytesPerPixel*x + 2]]++;
117                }
118
119            /* normalize histogram */
120            for (i = 0; i <= UCHAR_MAX; i++)
121            {
122                histogramR[i] = (float)histogramDenormR[i] / (float)pixNum;
123                histogramG[i] = (float)histogramDenormG[i] / (float)pixNum;
124                histogramB[i] = (float)histogramDenormB[i] / (float)pixNum;
125            }
126
127            /* cumulative histogram */
128            cumulativeHistogramR[0] = histogramR[0];
129            cumulativeHistogramG[0] = histogramG[0];
130            cumulativeHistogramB[0] = histogramB[0];
131            for (i = 1; i <= UCHAR_MAX; i++)
132            {
133                cumulativeHistogramR[i] = cumulativeHistogramR[i-1] + histogramR[i];
134                cumulativeHistogramG[i] = cumulativeHistogramG[i-1] + histogramG[i];
135                cumulativeHistogramB[i] = cumulativeHistogramB[i-1] + histogramB[i];
136            }
137
138            /* equalize */
139            for (y = 0; y < h; y++)
140                for (x = 0; x < w; x++)
141                {
142                    destData[y*destBytesPerRow + destBytesPerPixel*x]     = floor((UCHAR_MAX+0.9)*cumulativeHistogramR[srcData[y*srcBytesPerRow + srcBytesPerPixel*x]]);
143                    destData[y*destBytesPerRow + destBytesPerPixel*x + 1] = floor((UCHAR_MAX+0.9)*cumulativeHistogramG[srcData[y*srcBytesPerRow + srcBytesPerPixel*x + 1]]);
144                    destData[y*destBytesPerRow + destBytesPerPixel*x + 2] = floor((UCHAR_MAX+0.9)*cumulativeHistogramB[srcData[y*srcBytesPerRow + srcBytesPerPixel*x + 2]]);
145                    if (hasAlpha)
146                      destData[y*destBytesPerRow + destBytesPerPixel*x + 3] = srcData[y*srcBytesPerRow + srcBytesPerPixel*x + 3];
147                }
148        } else if (space == COLOR_SPACE_YUV)
149        {
150            unsigned long int histogramDenormY[UCHAR_MAX+1]; /* not normalized pixel count for each level */
151            float histogramY[UCHAR_MAX+1];                   /* normalized histogram */
152            float cumulativeHistogramY[UCHAR_MAX+1];         /* cumulative histogram */
153            register int r, g, b;
154            unsigned char yy, cb, cr;
155
156/*
157 *            JPEG-YCbCr (601) from "digital 8-bit R'G'B'  "
158 *            ========================================================================
159 *            Y' =       + 0.299    * R'd + 0.587    * G'd + 0.114    * B'd
160 *            Cb = 128   - 0.168736 * R'd - 0.331264 * G'd + 0.5      * B'd
161 *            Cr = 128   + 0.5      * R'd - 0.418688 * G'd - 0.081312 * B'd
162 *            ........................................................................
163 *            R'd, G'd, B'd   in {0, 1, 2, ..., 255}
164 *            Y', Cb, Cr      in {0, 1, 2, ..., 255}
165 * R = Y                    + 1.402   (Cr-128)
166 * G = Y - 0.34414 (Cb-128) - 0.71414 (Cr-128)
167 * B = Y + 1.772   (Cb-128)
168 */
169            /* first we convert the whole image to YCC-d */
170            for (y = 0; y < h; y++)
171                for (x = 0; x < w; x++)
172                {
173                    r = srcData[y*srcBytesPerRow + srcBytesPerPixel*x];
174                    g = srcData[y*srcBytesPerRow + srcBytesPerPixel*x + 1];
175                    b = srcData[y*srcBytesPerRow + srcBytesPerPixel*x + 2];
176                    yy = rintf(       + 0.2990f*r + 0.5870f*g + 0.1140f*b);
177                    cb = rintf(128.0f - 0.1687f*r - 0.3313f*g + 0.5000f*b);
178                    cr = rintf(128.0f + 0.5000f*r - 0.4187f*g - 0.0813f*b);
179
180                    destData[y*destBytesPerRow + destBytesPerPixel*x]     = yy;
181                    destData[y*destBytesPerRow + destBytesPerPixel*x + 1] = cb;
182                    destData[y*destBytesPerRow + destBytesPerPixel*x + 2] = cr;
183
184                    if (hasAlpha)
185                      destData[y*destBytesPerRow + destBytesPerPixel*x + 3] = srcData[y*srcBytesPerRow + srcBytesPerPixel*x + 3];
186                }
187
188            /* calculate the histogram */
189            for (i = 0; i <= UCHAR_MAX; i++)
190                histogramDenormY[i] =  0;
191            for (y = 0; y < h; y++)
192                for (x = 0; x < w; x++)
193                    histogramDenormY[destData[y*destBytesPerRow + destBytesPerPixel*x]]++;
194
195            /* normalize histogram */
196            for (i = 0; i <= UCHAR_MAX; i++)
197                histogramY[i] = (float)histogramDenormY[i] / (float)pixNum;
198
199            /* cumulative histogram */
200            cumulativeHistogramY[0] = histogramY[0];
201            for (i = 1; i <= UCHAR_MAX; i++)
202                cumulativeHistogramY[i] = cumulativeHistogramY[i-1] + histogramY[i];
203
204
205            /* equalize */
206            for (y = 0; y < h; y++)
207                for (x = 0; x < w; x++)
208                {
209                    destData[y*destBytesPerRow + destBytesPerPixel*x] = floor((UCHAR_MAX+0.9)*cumulativeHistogramY[destData[y*destBytesPerRow + destBytesPerPixel*x]]);
210                }
211
212            /* now we convert back to RGB */
213            for (y = 0; y < h; y++)
214                for (x = 0; x < w; x++)
215                {
216                    yy = destData[y*destBytesPerRow + destBytesPerPixel*x];
217                    cb = destData[y*destBytesPerRow + destBytesPerPixel*x + 1];
218                    cr = destData[y*destBytesPerRow + destBytesPerPixel*x + 2];
219                    r = yy                     + (int)rintf(1.40200f*(cr-128));
220                    g = yy - (int)rintf(0.34414f*(cb-128) + 0.71414f*(cr-128));
221                    b = yy + (int)rintf(1.77200f*(cb-128));
222                    r = r > UCHAR_MAX ? UCHAR_MAX : r;
223                    g = g > UCHAR_MAX ? UCHAR_MAX : g;
224                    b = b > UCHAR_MAX ? UCHAR_MAX : b;
225                    r = r < 0 ? 0 : r;
226                    g = g < 0 ? 0 : g;
227                    b = b < 0 ? 0 : b;
228                    destData[y*destBytesPerRow + destBytesPerPixel*x]     = r;
229                    destData[y*destBytesPerRow + destBytesPerPixel*x + 1] = g;
230                    destData[y*destBytesPerRow + destBytesPerPixel*x + 2] = b;
231
232                    /* no need to convert back alpha channel */
233                }
234        }
235    } else
236    {
237        unsigned long int histogramDenorm[UCHAR_MAX+1]; /* not normalized pixel count for each level */
238        float histogram[UCHAR_MAX+1];                   /* normalized histogram */
239        float cumulativeHistogram[UCHAR_MAX+1];         /* cumulative histogram */
240        /* calculate the histogram */
241        for (i = 0; i <= UCHAR_MAX; i++)
242            histogramDenorm[i] = 0;
243        for (y = 0; y < h; y++)
244            for (x = 0; x < w; x++)
245                histogramDenorm[srcData[y*srcBytesPerRow + x]]++;
246
247        /* normalize histogram */
248        for (i = 0; i <= UCHAR_MAX; i++)
249            histogram[i] = (float)histogramDenorm[i] / (float)pixNum;
250
251        /* cumulative histogram */
252        cumulativeHistogram[0] = histogram[0];
253        for (i = 1; i <= UCHAR_MAX; i++)
254            cumulativeHistogram[i] = cumulativeHistogram[i-1] + histogram[i];
255
256        /* equalize */
257        for (y = 0; y < h; y++)
258            for (x = 0; x < w; x++)
259              {
260                destData[y*destBytesPerRow + destBytesPerPixel*x] = floor((UCHAR_MAX+0.9)*cumulativeHistogram[srcData[y*srcBytesPerRow + srcBytesPerPixel*x]]);
261                if (hasAlpha)
262                  destData[y*destBytesPerRow + destBytesPerPixel*x + 1] = srcData[y*srcBytesPerRow + srcBytesPerPixel*x + 1];
263              }
264    }
265
266    [destImage setBitmapRep:destImageRep];
267    [destImageRep release];
268    [destImage autorelease];
269    return destImage;
270}
271
272
273@end
274