1 /**********
2 *
3 * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
4 * University of Manchester, UK. All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without modification,
7 * are permitted provided that the following conditions are met:
8 *
9 * . Redistributions of source code must retain the above copyright notice,
10 * this list of conditions and the following disclaimer.
11 *
12 * . Redistributions in binary form must reproduce the above copyright notice,
13 * this list of conditions and the following disclaimer in the documentation
14 * and/or other materials provided with the distribution.
15 *
16 * . Neither the name of the University of Manchester nor the names of its
17 * contributors may be used to endorse or promote products derived from this
18 * software without specific prior written permission.
19 *
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 *
33 **********
34 *
35 * Program : TINA
36 * File :
37 * Date :
38 * Version :
39 * CVS Id :
40 *
41 * Author : paul.bromiley@manchester.ac.uk
42 *
43 * Notes :
44 *
45 * Colour image conversion routines from RGB to YIQ, IHS, and IJK (my own
46 * invention: see Tina Memo no. 2001-015.
47 *
48 * Each function takes eight arguments: RGB fields of the image to be converted,
49 * three pointers (that should be NULL) for the outpur fields, an input error matrix,
50 * and a pointer for an output error matrix (whih shouls also be NULL).
51 *
52 * The covar matrix should be a (void *) pointer to a 2D (void *) array, for which
53 * each element is a 3 element fvector, containing the standard deviation of the noise
54 * estimate for the R, G and B values in the input images. If this argment is NULL, no
55 * error propagation will be performed, allowing these functions to be called for
56 * conversions of the stack.
57 *
58 * At the moment, the error propagation is only performed for the colour fields
59 * i.e. not for the intensity fields. Therefore, the output error matrices from the
60 * conversion functions contain only two elements on the fvector.
61 *
62 *********
63 */
64
65 #include "imgColour_convert.h"
66
67 #if HAVE_CONFIG_H
68 # include <config.h>
69 #endif
70
71 #include<math.h>
72 #include<tina/math/mathPro.h>
73 #include<tina/math/mathDef.h>
74 #include<tina/image/imgDef.h>
75 #include<tina/image/imgPro.h>
76
77 static double pi = 3.1415926535;
78
79
80 static void hsi_scaler(Imrect *im_red, Imrect *im_green, Imrect *im_blue, void ***covar_matrix)
81 {
82 int i, j, k;
83 Imregion *roi=NULL;
84 float *covar_row=NULL, scale=1.0;
85 float pix_r, pix_g, pix_b;
86 float min_r, max_r, min_g, max_g, min_b, max_b;
87
88 roi = im_red->region;
89 imf_minmax(im_red, &min_r, &max_r);
90 imf_minmax(im_green, &min_g, &max_g);
91 imf_minmax(im_blue, &min_b, &max_b);
92
93 scale = ((float)fmax3(((double)(max_r)), ((double)(max_g)), ((double)(max_b))));
94
95 for(j = roi->ly; j<roi->uy; j++)
96 {
97 for(k=roi->lx; k<roi->ux; k++)
98 {
99 pix_r = im_get_pixf(im_red, j, k);
100 pix_g = im_get_pixf(im_green, j, k);
101 pix_b = im_get_pixf(im_blue, j, k);
102
103 im_put_pixf((pix_r/scale), im_red, j, k);
104 im_put_pixf((pix_g/scale), im_green, j, k);
105 im_put_pixf((pix_b/scale), im_blue, j, k);
106 }
107 }
108 if(covar_matrix!=NULL)
109 {
110 for(j = roi->ly; j<roi->uy; j++)
111 {
112 for(k=roi->lx; k<roi->ux; k++)
113 {
114 covar_row = covar_matrix[j][k];
115 for(i=0; i<3; i++)
116 {
117 covar_row[i] = covar_row[i]/scale;
118 }
119 }
120 }
121 }
122 }
123
124
125 void **hsi_convert(Imrect *im_r, Imrect *im_g, Imrect *im_b, void ***covar_matrix)
126 {
127 int j, k;
128 float red_pix, green_pix, blue_pix, int_pix, hue_pix, sat_pix, top_line, bottom_line;
129 float red_noise, blue_noise, green_noise;
130 float sigma_s, sigma_h, alpha, inter_1, inter_2, inter_3, inter_4;
131 float *covar_row = NULL, *covar_row_out=NULL;
132 Imregion *roi=NULL;
133 Imrect *im_int=NULL, *im_sat=NULL, *im_hue=NULL;
134 Imrect *im_red=NULL, *im_green=NULL, *im_blue=NULL;
135 void **image_vector=NULL, ***covar_matrix_out=NULL, ***covar_matrix_int=NULL;
136
137 /* Convert to hsi */
138
139 roi = im_r->region;
140 if(covar_matrix!=NULL)
141 {
142 covar_matrix_int = parray_alloc(roi->ly, roi->lx, roi->uy, roi->ux);
143 covar_matrix_out = parray_alloc(roi->ly, roi->lx, roi->uy, roi->ux);
144 for(j = roi->ly; j<roi->uy; j++)
145 {
146 for(k=roi->lx; k<roi->ux; k++)
147 {
148 covar_row = covar_matrix[j][k];
149 covar_matrix_int[j][k] = fvector_copy((void *)covar_row, 0, 3);
150 }
151 }
152 }
153 im_red = im_cast(im_r, float_v);
154 im_green = im_cast(im_g, float_v);
155 im_blue = im_cast(im_b, float_v);
156 hsi_scaler(im_red, im_green, im_blue, covar_matrix_int);
157
158 im_hue = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), NULL, float_v);
159 im_sat = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), NULL, float_v);
160 im_int = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), NULL, float_v);
161 image_vector=pvector_alloc(0, 4);
162 image_vector[0] = im_hue;
163 image_vector[1] = im_sat;
164 image_vector[2] = im_int;
165 image_vector[3] = covar_matrix_out;
166
167 for(j = roi->ly; j<roi->uy; j++)
168 {
169 for(k=roi->lx; k<roi->ux; k++)
170 {
171 red_pix = im_get_pixf(im_red, j, k);
172 green_pix = im_get_pixf(im_green, j, k);
173 blue_pix = im_get_pixf(im_blue, j, k);
174
175 int_pix = (red_pix + blue_pix + green_pix) / 3;
176 sat_pix = 1-((3/(red_pix + blue_pix + green_pix))
177 *fmin3(red_pix, green_pix, blue_pix));
178
179 if(red_pix==green_pix&&red_pix==blue_pix)
180 {
181 sat_pix = 0;
182 }
183 if(sat_pix<0.0) sat_pix=0.0;
184 if(sat_pix>1.0) sat_pix=1.0;
185
186 if(sat_pix==0)
187 {
188 hue_pix = 0;
189 }
190 else
191 {
192 top_line = (0.5*((red_pix-green_pix)+(red_pix-blue_pix)));
193 bottom_line = sqrt(( sqr((red_pix-green_pix))+((red_pix-blue_pix)*
194 (green_pix-blue_pix)) ));
195 hue_pix = acos((top_line/bottom_line));
196 }
197 if(blue_pix>green_pix)
198 {
199 hue_pix = 2.0*pi - hue_pix;
200 }
201 im_put_pixf(int_pix, im_int, j, k);
202 im_put_pixf(sat_pix, im_sat, j, k);
203 im_put_pixf(hue_pix, im_hue, j, k);
204
205 if(covar_matrix!=NULL)
206 {
207 /* The next few lines contain the error propagation results
208 for the RGB to IHS conversion equations */
209
210 covar_row = covar_matrix_int[j][k];
211 red_noise = covar_row[0];
212 green_noise = covar_row[1];
213 blue_noise = covar_row[2];
214
215 if(fmin3(red_pix, blue_pix, green_pix )==red_pix)
216 {
217 sigma_s = (3.0/(sqr(red_pix+green_pix+blue_pix)))*
218 sqrt(sqr(green_pix+blue_pix)*sqr(red_noise) +
219 sqr(red_pix)*(sqr(green_noise)+sqr(blue_noise)));
220 }
221 else if(fmin3(red_pix, blue_pix, green_pix)==green_pix)
222 {
223 sigma_s = (3.0/(sqr(red_pix+green_pix+blue_pix)))*
224 sqrt(sqr(blue_pix+red_pix)*sqr(green_noise)
225 + sqr(green_pix)*(sqr(blue_noise)+sqr(red_noise)));
226 }
227 else
228 {
229 sigma_s = (3.0/(sqr(red_pix+green_pix+blue_pix)))*
230 sqrt(sqr(green_pix+red_pix)*sqr(blue_noise)
231 + sqr(blue_pix)*(sqr(green_noise)+sqr(red_noise)));
232 }
233
234 if(sat_pix!=0)
235 {
236 alpha = (sqr(blue_pix - green_pix) + sqr(blue_pix-red_pix)
237 + sqr(green_pix-red_pix))/2.0;
238
239 inter_1 = 1.0/ (1.0-( (sqr(blue_pix-red_pix)
240 *sqr(green_pix-red_pix))/(4.0*alpha)));
241
242 inter_2 = (green_pix-red_pix) * ( (1.0/sqrt(alpha))
243 -( (0.5*(blue_pix-red_pix)*(2.0*blue_pix-green_pix-red_pix))
244 / (pow(alpha, 1.5))) );
245
246 inter_3 = (blue_pix-red_pix) * ( (1.0/sqrt(alpha))
247 -( (0.5*(green_pix-red_pix)*(2.0*green_pix-blue_pix-red_pix))
248 / (pow(alpha, 1.5))) );
249
250 inter_4 = ((red_pix-blue_pix)/sqrt(alpha)) +
251 ((red_pix-green_pix)/sqrt(alpha)) +
252 ((0.5*(green_pix-red_pix)*(blue_pix-red_pix)*(blue_pix +
253 green_pix - 2.0*red_pix))/pow(alpha, 1.5));
254
255 sigma_h = 0.5*sqrt((inter_1*( sqr(inter_2)*sqr(blue_noise)
256 + sqr(inter_2)*sqr(green_noise)
257 + sqr(inter_3)*sqr(red_noise) ) ));
258 }
259 else
260 {
261 sigma_h = 0.0;
262 }
263
264 /* End of error propagation */
265
266 covar_row_out = fvector_alloc(0, 2);
267 covar_row_out[0] = sigma_s;
268 covar_row_out[1] = sigma_h;
269 covar_matrix_out[j][k] = covar_row_out;
270 }
271 }
272 }
273 im_free(im_red);
274 im_free(im_green);
275 im_free(im_blue);
276
277 if(covar_matrix!=NULL)
278 {
279 for(j=roi->ly; j<roi->uy; j++)
280 {
281 for(k=roi->lx; k<roi->ux; k++)
282 {
283 pvector_free(covar_matrix_int[j][k], 0);
284 }
285 }
286 parray_free(covar_matrix_int, roi->ly, roi->lx, roi->uy, roi->ux);
287 }
288 return(image_vector);
289 }
290
291 void **yiq_convert(Imrect *im_red, Imrect *im_green, Imrect *im_blue, void ***covar_matrix)
292 {
293 int j, k;
294 float red_pix, green_pix, blue_pix, y_pix, i_pix, q_pix;
295 float *covar_row=NULL, *covar_row_out=NULL;
296 float sigma_i, sigma_q;
297 Imregion *roi=NULL;
298 float red_noise, blue_noise, green_noise;
299 Imrect *im_y=NULL, *im_i=NULL, *im_q=NULL;
300 void **image_vector=NULL, ***covar_matrix_out=NULL;
301
302 /* Convert to yiq */
303
304 roi = im_red->region;
305 if(covar_matrix!=NULL)
306 {
307 covar_matrix_out = parray_alloc(roi->ly, roi->lx, roi->uy, roi->ux);
308 }
309 im_y = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), NULL, float_v);
310 im_i = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), NULL, float_v);
311 im_q = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), NULL, float_v);
312 image_vector = pvector_alloc(0, 4);
313 image_vector[0] = im_y;
314 image_vector[1] = im_i;
315 image_vector[2] = im_q;
316 image_vector[3] = covar_matrix_out;
317
318 for(j = roi->ly; j<roi->uy; j++)
319 {
320 for(k=roi->lx; k<roi->ux; k++)
321 {
322 red_pix = im_get_pixf(im_red, j, k);
323 green_pix = im_get_pixf(im_green, j, k);
324 blue_pix = im_get_pixf(im_blue, j, k);
325
326 y_pix = 0.299*red_pix + 0.587*green_pix + 0.114*blue_pix;
327 i_pix = 0.569*red_pix - 0.274*green_pix - 0.322*blue_pix;
328 q_pix = 0.211*red_pix - 0.523*green_pix + 0.312*blue_pix;
329
330 im_put_pixf(y_pix, im_y, j, k);
331 im_put_pixf(i_pix, im_i, j, k);
332 im_put_pixf(q_pix, im_q, j, k);
333
334 if(covar_matrix!=NULL)
335 {
336 /* The next few lines contain the error propagation results
337 for the RGB to YIQ conversion equations */
338
339 covar_row = covar_matrix[j][k];
340 red_noise = covar_row[0];
341 green_noise = covar_row[1];
342 blue_noise = covar_row[2];
343
344 sigma_i = sqrt(sqr(0.596)*sqr(red_noise) +
345 sqr(0.274)*sqr(green_noise) + sqr(0.322)*sqr(blue_noise));
346 sigma_q = sqrt(sqr(0.211)*sqr(red_noise) +
347 sqr(0.523)*sqr(green_noise) +sqr(0.312)*sqr(blue_noise));
348
349 covar_row_out = fvector_alloc(0, 2);
350 covar_row_out[0] = sigma_i;
351 covar_row_out[1] = sigma_q;
352 covar_matrix_out[j][k] = covar_row_out;
353 }
354 }
355 }
356 return(image_vector);
357 }
358
359
360 void **ijk_convert(Imrect *im_red, Imrect *im_green, Imrect *im_blue, void ***covar_matrix)
361 {
362 int j, k;
363 float red_pix, green_pix, blue_pix, i2_pix, j_pix, k_pix;
364 float sigma_j, sigma_k;
365 float *covar_row=NULL, *covar_row_out=NULL;
366 Imregion *roi=NULL;
367 float red_noise, green_noise, blue_noise;
368 Imrect *im_i2=NULL, *im_j=NULL, *im_k=NULL;
369 void **image_vector=NULL, ***covar_matrix_out=NULL;
370
371 roi = im_red->region;
372 if(covar_matrix!=NULL)
373 {
374 covar_matrix_out = parray_alloc(roi->ly, roi->lx, roi->uy, roi->ux);
375 }
376 im_i2 = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), NULL, float_v);
377 im_j = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), NULL, float_v);
378 im_k = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), NULL, float_v);
379 image_vector = pvector_alloc(0, 4);
380 image_vector[0] = im_i2;
381 image_vector[1] = im_j;
382 image_vector[2] = im_k;
383 image_vector[3] = covar_matrix_out;
384
385 for(j = roi->ly; j<roi->uy; j++)
386 {
387 for(k=roi->lx; k<roi->ux; k++)
388 {
389 red_pix = im_get_pixf(im_red, j, k);
390 green_pix = im_get_pixf(im_green, j, k);
391 blue_pix = im_get_pixf(im_blue, j, k);
392
393 i2_pix = (red_pix + green_pix + blue_pix)/sqrt(3.0);
394 j_pix = (2.0*red_pix - green_pix - blue_pix)/sqrt(6.0);
395 k_pix = (blue_pix-green_pix)/sqrt(2.0);
396
397 im_put_pixf(i2_pix, im_i2, j, k);
398 im_put_pixf(j_pix, im_j, j, k);
399 im_put_pixf(k_pix, im_k, j, k);
400
401 if(covar_matrix!=NULL)
402 {
403 /* The next few lines contain the error propagation results
404 for the RGB to IJK conversion equations */
405
406 covar_row = covar_matrix[j][k];
407 red_noise = covar_row[0];
408 green_noise = covar_row[1];
409 blue_noise = covar_row[2];
410
411 sigma_j = (4.0*sqr(red_noise) + sqr(blue_noise) + sqr(green_noise))/6.0;
412 sigma_j = sqrt(sigma_j);
413 sigma_k = (sqr(green_noise) + sqr(blue_noise))/2.0;
414 sigma_k = sqrt(sigma_k);
415
416 covar_row_out = fvector_alloc(0, 2);
417 covar_row_out[0] = sigma_j;
418 covar_row_out[1] = sigma_k;
419 covar_matrix_out[j][k] = covar_row_out;
420 }
421 }
422 }
423 return(image_vector);
424 }
425
426
427 void **rgb_normalise(Imrect *im_red, Imrect *im_green, Imrect *im_blue, void ***covar_matrix)
428 {
429 Imregion *roi=NULL;
430 Imrect *im_r2=NULL, *im_g2=NULL, *im_b2=NULL;
431 int i, j;
432 float red_pix, green_pix, blue_pix, tot_pix;
433 float red_noise, green_noise, blue_noise;
434 float sigma_r, sigma_g, sigma_b;
435 float *covar_row = NULL, *covar_row_out=NULL;
436 void ***covar_matrix_out=NULL, **image_vector=NULL;
437
438 /* Subroutine to normalise RGB images */
439
440 roi=im_red->region;
441 if(covar_matrix!=NULL)
442 {
443 covar_matrix_out = parray_alloc(roi->ly, roi->lx, roi->uy, roi->ux);
444 }
445 im_r2 = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), NULL, float_v);
446 im_g2 = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), NULL, float_v);
447 im_b2 = im_alloc((roi->uy-roi->ly), (roi->ux-roi->lx), NULL, float_v);
448 image_vector = pvector_alloc(0, 4);
449 image_vector[0] = im_r2;
450 image_vector[1] = im_g2;
451 image_vector[2] = im_b2;
452 image_vector[3] = covar_matrix_out;
453
454 for(i=roi->ly; i<roi->uy; i++)
455 {
456 for(j = roi->lx; j<roi->ux; j++)
457 {
458 red_pix = (float)im_get_pixf(im_red, i, j);
459 green_pix = (float)im_get_pixf(im_green, i, j);
460 blue_pix = (float)im_get_pixf(im_blue, i, j);
461 tot_pix = red_pix+green_pix+blue_pix;
462
463 im_put_pixf((red_pix/tot_pix), im_r2, i, j);
464 im_put_pixf((green_pix/tot_pix), im_g2, i, j);
465 im_put_pixf((blue_pix/tot_pix), im_b2, i, j);
466
467 if(covar_matrix!=NULL)
468 {
469 /* Do error propagation on the normalisation */
470
471 covar_row = covar_matrix[i][j];
472 red_noise = covar_row[0];
473 green_noise = covar_row[1];
474 blue_noise = covar_row[2];
475
476 sigma_r = (sqr(red_noise)*sqr(blue_pix + green_pix)
477 + sqr(blue_noise)*sqr(red_pix) + sqr(green_noise)*sqr(red_pix))
478 / pow(tot_pix, 4.0);
479 sigma_g = (sqr(green_noise)*sqr(red_pix + blue_pix)
480 + sqr(blue_noise)*sqr(green_pix) + sqr(red_noise)*sqr(green_pix))
481 / pow(tot_pix, 4.0);
482 sigma_b = (sqr(blue_noise)*sqr(red_pix + green_pix)
483 + sqr(red_noise)*sqr(blue_pix) + sqr(green_noise)*sqr(blue_pix))
484 / pow(tot_pix, 4.0);
485 sigma_r = sqrt(sigma_r);
486 sigma_g = sqrt(sigma_g);
487 sigma_b = sqrt(sigma_b);
488
489 covar_row_out = fvector_alloc(0, 3);
490 covar_row_out[0] = sigma_r;
491 covar_row_out[1] = sigma_g;
492 covar_row_out[2] = sigma_b;
493 covar_matrix_out[i][j] = covar_row_out;
494 }
495 }
496 }
497 return(image_vector);
498 }
499
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.