-
Notifications
You must be signed in to change notification settings - Fork 7
/
CommentedExample.m
242 lines (211 loc) · 7.82 KB
/
CommentedExample.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
%
% Hello,
%
% This file is a script that you can use to try out the S-CIELAB
% calculations we have developed. The sequence of steps in this file
% will (1) illustrate how to calculate the S-CIELAB representation for an
% image, and (2) illustrate how to calculate the S-CIELAB errors
% (delta Es) for a pair of images.
%
% This file is written as if you are going to evaluate the error
% between a pair of images that are presently coded as RGB data. The
% initial calculation converts the RGB data to CIE XYZ data.
% The calculations are similar if your images are coded as CMYK,
% but you will need to convert from CMYK to XYZ instead of converting
% from RGB to XYZ.
%
%
%%END INTRODUCTION
%%%%EXAMPLE BEGINS
% 1. Reading in your image data
%
% First, load in an RGB image. The example image below consists of
% three matrices, rHats,gHats and bHats, each of size 128 x 192.
load images/hats
load images/hatsCompressed
% If you would like to view the image, we suggest you use a tiff
% viewer, like xv, and type
%
% unix('xv -perfect images/hats.tiff &')
%
% If you have the Matlab image toolbox, you can read in the tiff files
% directly (which is what we usually do) via
%
% [rHats gHats bHats] = tiffread('images/hats.tiff');
% [rHatsc gHatsc bHatsc] = tiffread('images/hatsCompressed.tiff');
%
% To see the different color planes, you could try looking at
%
% colormap( gray(128)*diag([1,0,0])), imagesc(rHats)
% colormap( gray(128)*diag([0,1,0])), imagesc(gHats)
% colormap( gray(128)*diag([0,0,1])), imagesc(bHats)
% 2. Defining the calibration information
%
% 2.1: Spatial calibration information.
%
% Now, we need to specify the viewing conditions. Specifically, we
% need to know the relevant spatial parameters for when the digital
% image is viewed. This parameter tells us how many digital samples
% there are per degree of viewing angle. We call it sampPerDeg. For
% a typical monitor (72 dots per inch) viewed at 18 inches, one inch
% sweeps out 3.1798 = (180/pi) * atan(1/18) deg. Hence, there are
% 72/3.1798 samples in a single degree of visual angle for this monitor.
%
% round(72 / ((180/pi)*atan(1/18)))
sampPerDeg = 23;
% 2.2: Color calibration information
%
% Next, we need to specify two aspects of the display characteristcs
% so that we can specify how the displayed image affects the cone
% photoreceptors. To make this estimate, we need to know something
% about
% 1) the effect each display primary has on your cones, and
% 2) the relationship between the frame-buffer values and the
% intensity of the display primaries
%
% For this example, we assumed that your display is pretty standard
% (and like some of ours). These assumptions are not going to be
% preicsely correct for your CRTs; they are not even close for
% printers or LCD displays.
%
% A broad description of the issues concerning color calibration may be
% found in the following references:
%
% Brainard -- CR and A paper
% ???
% Wandell-- Appendix B
%
% 2.2.1 -- Effects of the display primaries on the human cones
%
% To compute the effect of the display primaries on the cones, we
% need to know
% 1) The spectral power distribution (SPD) of the display, and
% 2) The relative absorptions of the human cones.
%
% We include a file containing the display SPDs of one of our
% monitors in the file displaySPD.mat. The file contains the
% variables:
%
% wavelength: a vector describing the wavelengths of the cone
% sensitivities (370:73))
% displaySPD: a 361x3 matrix of the spectral power distributions
% of the red, green, and blue primaries
%
load displaySPD
% We usually compute with respect to the human cone estimates from
% Smith and Pokorny (19XX). These are represented in the file named
% SmithPokornyCones.mat
%
% The file contains the following variables:
% wavelength: a vector describing the wavelengths of the cone
% sensitivities (370:73))
% cones: a 361 x 3 matrix matrix describing the cone
% sensitivities at the corresponding wavelengths
%
load SmithPokornyCones
% If you would like to see the cone sensitivities or the display
% spectral power distribution we have assumed, type:
%
% plot(wavelength,cones)
% plot(wavelength, displaySPD)
%
% REFERENCE:
%
% Now, we can compute the 3 x 3 transformation that maps the linear
% intensity of the display r,g,b signals into the cone absorptions
% (l,m,s).
rgb2lms = cones'* displaySPD;
% 2.2.2 -- Gamma correction
%
% There is a nonlinear relationship between the frame-buffer values
% in the image data and the intensity of the light emitted by each of
% the primaries on your display. This relationship is captured by a
% function commonly called the ``display gamma curve.''
%
% The file displayGamma.mat contains the variables
% gamma: is a look-up table that maps
% framebuffer values -> relative linear display intensity
% where the maximum display intensity is 1.0.
%
% invGamma: is the inverse of gamma and maps
% relative lin. disp. intensity -> frame-buffer value
% This table has been interpolated to 1024 intensity levels.
%
load displayGamma
% 2.2.3 --- The CIELAB computation requires specifiying the
% white-point of the data set. In our view, each image should have
% its own white point. In a generally abominable practice, which we
% follow here, the white point is commonly set to be the white of the
% monitor. We don't like this much because the white point is
% supposed to depend on the image data, not the display device. But,
% nevermind.
%
rgbWhite = [1 1 1];
whitepoint = rgbWhite * rgb2lms';
% We have helped you out here by writing the little routine
% ``cmatrix'' that returns some of the most frequently used matrices,
% such as opp2lms, opp2xyz, and so forth. The entries returned by
% that routine were computed using the same method shown above for
% computing rgb2lms.
%
% Setting up the image data for S-CIELAB
% 3.1 -- Convert the RGB data to LMS (or XYZ if you like).
%
% Now, we convert the rgb data into the lms format. This takes place
% in two steps. First, we correct the r,g,b frame-buffer values for
% the gamma curve. Second, we convert the linearized r,g,b values
% into cone absorptions using the matrix rgb2lms
%
img = [ rHats gHats bHats];
imgRGB = dac2rgb(img,gammaTable);
img1LMS = changeColorSpace(imgRGB,rgb2lms);
img = [ rHatsc gHatsc bHatsc];
imgRGB = dac2rgb(img,gammaTable);
img2LMS = changeColorSpace(imgRGB,rgb2lms);
% 3.2 -- Identify the color space
%
% The scielab function accepts input in a few different formats. We
% need to tell it whether the input data are in lms, or xyz format.
%
imageformat = 'lms';
% 4. -- Run the scielab function.
% See the paper to understand what it does more fully.
%
% REFERENCE:
% SID paper
%
errorImage = scielab(sampPerDeg, img1LMS, img2LMS, whitepoint, imageformat);
% 5. -- Examining and interpreting the results.
%
max(img1LMS(:)), min(img1LMS(:))
max(img2LMS(:)), min(img2LMS(:))
% For this example, most of the errors are less than 10, as this
% histogram shows.
%
hist(errorImage(:),[1:2:14])
sum(errorImage(:) > 20) % We think this is 173
% Let's study the spatial distribution of the errors and just mark the
% ones that are 10 or larger using green
%
errorTruncated = min(128*(errorImage/10),128*ones(size(errorImage)));
figure(1)
colormap([gray(127); [0 1 0]])
colormap([gray(127); [1 1 1]])
image(errorTruncated)
% If you have the image processing toolbox, you can find out where the
% edges are and overlay the edges with the locations of the scielab
% errors
%
%% For Matlab version 5, use the following command:
edgeImage = 129 * double(edge(rHats,'prewitt'));
%% For Matlab version 4 or 3, use the following command:
edgeImage = 129 * edge(rHats,'prewitt');
comparison = max(edgeImage,errorTruncated);
mp = [gray(127); [0 1 0]; [1 0 0] ];
colormap(mp)
image(comparison)
figure(2)
colormap(gray(128));
imagesc([rHats + gHats + bHats] / 3)
%%%%%EXAMPLE ENDS
close all;