7975d940066516578d56c68f833d073ce2be5c42
1 /*
2 Copyright 2005, 2006 by Gerald Friedland, Kristian Jantz and Lars Knipping
4 Conversion to C++ for Inkscape by Bob Jamison
6 Licensed under the Apache License, Version 2.0 (the "License");
7 you may not use this file except in compliance with the License.
8 You may obtain a copy of the License at
10 http://www.apache.org/licenses/LICENSE-2.0
12 Unless required by applicable law or agreed to in writing, software
13 distributed under the License is distributed on an "AS IS" BASIS,
14 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 See the License for the specific language governing permissions and
16 limitations under the License.
17 */
18 #include "siox.h"
20 #include <string>
22 #include <stdarg.h> //for error() and trace()
23 #include <math.h> //sqrt(), pow(), round(), etc
26 namespace org
27 {
28 namespace siox
29 {
31 //########################################################################
32 //## U T I L S (originally Utils.java)
33 //########################################################################
35 /**
36 * Collection of auxiliary image processing methods used by the
37 * SioxSegmentator mainly for postprocessing.
38 *
39 * @author G. Friedland, K. Jantz, L. Knipping
40 * @version 1.05
41 *
42 * Conversion to C++ by Bob Jamison
43 *
44 */
47 /** Caches color conversion values to speed up RGB->CIELAB conversion.*/
48 static std::map<unsigned long, CLAB> RGB_TO_LAB;
50 //forward decls
51 static void premultiplyMatrix(float alpha, float *cm, int cmSize);
52 //static float colordiffsq(unsigned long rgb0, unsigned long rgb1);
53 //static int getAlpha(unsigned long argb);
54 static int getRed(unsigned long rgb);
55 static int getGreen(unsigned long rgb);
56 static int getBlue(unsigned long rgb);
57 //static unsigned long packPixel(int a, int r, int g, int b);
58 static CLAB rgbToClab(unsigned long rgb);
60 /**
61 * Applies the morphological dilate operator.
62 *
63 * Can be used to close small holes in the given confidence matrix.
64 *
65 * @param cm Confidence matrix to be processed.
66 * @param xres Horizontal resolution of the matrix.
67 * @param yres Vertical resolution of the matrix.
68 */
69 static void dilate(float *cm, int xres, int yres)
70 {
71 for (int y=0; y<yres; y++) {
72 for (int x=0; x<xres-1; x++) {
73 int idx=(y*xres)+x;
74 if (cm[idx+1]>cm[idx])
75 cm[idx]=cm[idx+1];
76 }
77 }
78 for (int y=0; y<yres; y++) {
79 for (int x=xres-1; x>=1; x--) {
80 int idx=(y*xres)+x;
81 if (cm[idx-1]>cm[idx])
82 cm[idx]=cm[idx-1];
83 }
84 }
85 for (int y=0; y<yres-1; y++) {
86 for (int x=0; x<xres; x++) {
87 int idx=(y*xres)+x;
88 if (cm[((y+1)*xres)+x] > cm[idx])
89 cm[idx]=cm[((y+1)*xres)+x];
90 }
91 }
92 for (int y=yres-1; y>=1; y--) {
93 for (int x=0; x<xres; x++) {
94 int idx=(y*xres)+x;
95 if (cm[((y-1)*xres)+x] > cm[idx])
96 cm[idx]=cm[((y-1)*xres)+x];
97 }
98 }
99 }
101 /**
102 * Applies the morphological erode operator.
103 *
104 * @param cm Confidence matrix to be processed.
105 * @param xres Horizontal resolution of the matrix.
106 * @param yres Vertical resolution of the matrix.
107 */
108 static void erode(float *cm, int xres, int yres)
109 {
110 for (int y=0; y<yres; y++) {
111 for (int x=0; x<xres-1; x++) {
112 int idx=(y*xres)+x;
113 if (cm[idx+1] < cm[idx])
114 cm[idx]=cm[idx+1];
115 }
116 }
117 for (int y=0; y<yres; y++) {
118 for (int x=xres-1; x>=1; x--) {
119 int idx=(y*xres)+x;
120 if (cm[idx-1] < cm[idx])
121 cm[idx]=cm[idx-1];
122 }
123 }
124 for (int y=0; y<yres-1; y++) {
125 for (int x=0; x<xres; x++) {
126 int idx=(y*xres)+x;
127 if (cm[((y+1)*xres)+x] < cm[idx])
128 cm[idx]=cm[((y+1)*xres)+x];
129 }
130 }
131 for (int y=yres-1; y>=1; y--) {
132 for (int x=0; x<xres; x++) {
133 int idx=(y*xres)+x;
134 if (cm[((y-1)*xres)+x] < cm[idx])
135 cm[idx]=cm[((y-1)*xres)+x];
136 }
137 }
138 }
140 /**
141 * Normalizes the matrix to values to [0..1].
142 *
143 * @param cm The matrix to be normalized.
144 */
145 static void normalizeMatrix(float *cm, int cmSize)
146 {
147 float max=0.0f;
148 for (int i=0; i<cmSize; i++) {
149 if (max<cm[i])
150 max=cm[i];
151 }
152 if (max<=0.0)
153 return;
154 else if (max==1.00)
155 return;
157 float alpha=1.00f/max;
158 premultiplyMatrix(alpha, cm, cmSize);
159 }
161 /**
162 * Multiplies matrix with the given scalar.
163 *
164 * @param alpha The scalar value.
165 * @param cm The matrix of values be multiplied with alpha.
166 * @param cmSize The matrix length.
167 */
168 static void premultiplyMatrix(float alpha, float *cm, int cmSize)
169 {
170 for (int i=0; i<cmSize; i++)
171 cm[i]=alpha*cm[i];
172 }
174 /**
175 * Blurs confidence matrix with a given symmetrically weighted kernel.
176 * <P>
177 * In the standard case confidence matrix entries are between 0...1 and
178 * the weight factors sum up to 1.
179 *
180 * @param cm The matrix to be smoothed.
181 * @param xres Horizontal resolution of the matrix.
182 * @param yres Vertical resolution of the matrix.
183 * @param f1 Weight factor for the first pixel.
184 * @param f2 Weight factor for the mid-pixel.
185 * @param f3 Weight factor for the last pixel.
186 */
187 static void smoothcm(float *cm, int xres, int yres,
188 float f1, float f2, float f3)
189 {
190 for (int y=0; y<yres; y++) {
191 for (int x=0; x<xres-2; x++) {
192 int idx=(y*xres)+x;
193 cm[idx]=f1*cm[idx]+f2*cm[idx+1]+f3*cm[idx+2];
194 }
195 }
196 for (int y=0; y<yres; y++) {
197 for (int x=xres-1; x>=2; x--) {
198 int idx=(y*xres)+x;
199 cm[idx]=f3*cm[idx-2]+f2*cm[idx-1]+f1*cm[idx];
200 }
201 }
202 for (int y=0; y<yres-2; y++) {
203 for (int x=0; x<xres; x++) {
204 int idx=(y*xres)+x;
205 cm[idx]=f1*cm[idx]+f2*cm[((y+1)*xres)+x]+f3*cm[((y+2)*xres)+x];
206 }
207 }
208 for (int y=yres-1; y>=2; y--) {
209 for (int x=0; x<xres; x++) {
210 int idx=(y*xres)+x;
211 cm[idx]=f3*cm[((y-2)*xres)+x]+f2*cm[((y-1)*xres)+x]+f1*cm[idx];
212 }
213 }
214 }
216 /**
217 * Squared Euclidian distance of p and q.
218 * <P>
219 * Usage hint: When only comparisons between Euclidian distances without
220 * actual values are needed, the squared distance can be preferred
221 * for being faster to compute.
222 *
223 * @param p First euclidian point coordinates.
224 * @param pSize Length of coordinate array.
225 * @param q Second euclidian point coordinates.
226 * Dimension must not be smaller than that of p.
227 * Any extra dimensions will be ignored.
228 * @return Squared euclidian distance of p and q.
229 * @see #euclid
230 */
231 static float sqrEuclidianDist(float *p, int pSize, float *q)
232 {
233 float sum=0;
234 for (int i=0; i<pSize; i++)
235 sum+=(p[i]-q[i])*(p[i]-q[i]);
236 return sum;
237 }
239 /**
240 * Squared Euclidian distance of p and q.
241 * <P>
242 * Usage hint: When only comparisons between Euclidian distances without
243 * actual values are needed, the squared distance can be preferred
244 * for being faster to compute.
245 *
246 * @param p CLAB value
247 * @param q second CLAB value
248 * @return Squared euclidian distance of p and q.
249 * @see #euclid
250 */
251 static float sqrEuclidianDist(const CLAB &p, const CLAB &q)
252 {
253 float sum=0;
254 sum += (p.C - q.C) * (p.C - q.C);
255 sum += (p.L - q.L) * (p.L - q.L);
256 sum += (p.A - q.A) * (p.A - q.A);
257 sum += (p.B - q.B) * (p.B - q.B);
258 return sum;
259 }
261 /**
262 * Euclidian distance of p and q.
263 *
264 * @param p First euclidian point coordinates.
265 * @param pSize Length of coordinate array.
266 * @param q Second euclidian point coordinates.
267 * Dimension must not be smaller than that of p.
268 * Any extra dimensions will be ignored.
269 * @return Squared euclidian distance of p and q.
270 * @see #sqrEuclidianDist
271 */
272 /*
273 static float euclid(float *p, int pSize, float *q)
274 {
275 return (float)sqrt(sqrEuclidianDist(p, pSize, q));
276 }
277 */
279 /**
280 * Computes Euclidian distance of two RGB color values.
281 *
282 * @param rgb0 First color value.
283 * @param rgb1 Second color value.
284 * @return Euclidian distance between the two color values.
285 */
286 /*
287 static float colordiff(long rgb0, long rgb1)
288 {
289 return (float)sqrt(colordiffsq(rgb0, rgb1));
290 }
291 */
293 /**
294 * Computes squared euclidian distance of two RGB color values.
295 * <P>
296 * Note: Faster to compute than colordiff
297 *
298 * @param rgb0 First color value.
299 * @param rgb1 Second color value.
300 * @return Squared Euclidian distance between the two color values.
301 */
302 /*
303 static float colordiffsq(long rgb0, long rgb1)
304 {
305 int rDist=getRed(rgb0) - getRed(rgb1);
306 int gDist=getGreen(rgb0) - getGreen(rgb1);
307 int bDist=getBlue(rgb0) - getBlue(rgb1);
309 return (float)(rDist*rDist+gDist*gDist+bDist*bDist);
310 }
311 */
313 /**
314 * Averages two ARGB colors.
315 *
316 * @param argb0 First color value.
317 * @param argb1 Second color value.
318 * @return The averaged ARGB color.
319 */
320 /*
321 static long average(long argb0, long argb1)
322 {
323 long ret = packPixel(
324 (getAlpha(argb0) + getAlpha(argb1))/2,
325 (getRed(argb0) + getRed(argb1) )/2,
326 (getGreen(argb0) + getGreen(argb1))/2,
327 (getBlue(argb0) + getBlue(argb1) )/2);
328 return ret;
329 }
330 */
332 /**
333 * Computes squared euclidian distance in CLAB space for two colors
334 * given as RGB values.
335 *
336 * @param rgb0 First color value.
337 * @param rgb1 Second color value.
338 * @return Squared Euclidian distance in CLAB space.
339 */
340 static float labcolordiffsq(unsigned long rgb1, unsigned long rgb2)
341 {
342 CLAB c1 = rgbToClab(rgb1);
343 CLAB c2 = rgbToClab(rgb2);
344 float euclid=0.0f;
345 euclid += (c1.L - c2.L) * (c1.L - c2.L);
346 euclid += (c1.A - c2.A) * (c1.A - c2.A);
347 euclid += (c1.B - c2.B) * (c1.B - c2.B);
348 return euclid;
349 }
352 /**
353 * Computes squared euclidian distance in CLAB space for two colors
354 * given as RGB values.
355 *
356 * @param rgb0 First color value.
357 * @param rgb1 Second color value.
358 * @return Euclidian distance in CLAB space.
359 */
360 static float labcolordiff(unsigned long rgb0, unsigned long rgb1)
361 {
362 return (float)sqrt(labcolordiffsq(rgb0, rgb1));
363 }
366 /**
367 * Converts 24-bit RGB values to {l,a,b} float values.
368 * <P>
369 * The conversion used is decribed at
370 * <a href="http://www.easyrgb.com/math.php?MATH=M7#text7">CLAB Conversion</a>
371 * for reference white D65. Note that that the conversion is computational
372 * expensive. Result are cached to speed up further conversion calls.
373 *
374 * @param rgb RGB color value,
375 * @return CLAB color value tripel.
376 */
377 static CLAB rgbToClab(unsigned long rgb)
378 {
379 std::map<unsigned long, CLAB>::iterator iter = RGB_TO_LAB.find(rgb);
380 if (iter != RGB_TO_LAB.end())
381 {
382 CLAB res = iter->second;
383 return res;
384 }
386 int R=getRed(rgb);
387 int G=getGreen(rgb);
388 int B=getBlue(rgb);
390 float var_R=(R/255.0f); //R = From 0 to 255
391 float var_G=(G/255.0f); //G = From 0 to 255
392 float var_B=(B/255.0f); //B = From 0 to 255
394 if (var_R>0.04045)
395 var_R=(float) pow((var_R+0.055f)/1.055f, 2.4);
396 else
397 var_R=var_R/12.92f;
399 if (var_G>0.04045)
400 var_G=(float) pow((var_G+0.055f)/1.055f, 2.4);
401 else
402 var_G=var_G/12.92f;
404 if (var_B>0.04045)
405 var_B=(float) pow((var_B+0.055f)/1.055f, 2.4);
406 else
407 var_B=var_B/12.92f;
409 var_R=var_R*100.0f;
410 var_G=var_G*100.0f;
411 var_B=var_B*100.0f;
413 // Observer. = 2�, Illuminant = D65
414 float X=var_R*0.4124f + var_G*0.3576f + var_B*0.1805f;
415 float Y=var_R*0.2126f + var_G*0.7152f + var_B*0.0722f;
416 float Z=var_R*0.0193f + var_G*0.1192f + var_B*0.9505f;
418 float var_X=X/95.047f;
419 float var_Y=Y/100.0f;
420 float var_Z=Z/108.883f;
422 if (var_X>0.008856f)
423 var_X=(float) pow(var_X, 0.3333f);
424 else
425 var_X=(7.787f*var_X)+(16.0f/116.0f);
427 if (var_Y>0.008856f)
428 var_Y=(float) pow(var_Y, 0.3333f);
429 else
430 var_Y=(7.787f*var_Y)+(16.0f/116.0f);
432 if (var_Z>0.008856f)
433 var_Z=(float) pow(var_Z, 0.3333f);
434 else
435 var_Z=(7.787f*var_Z)+(16.0f/116.0f);
437 CLAB lab((116.0f*var_Y)-16.0f , 500.0f*(var_X-var_Y), 200.0f*(var_Y-var_Z));
439 RGB_TO_LAB[rgb] = lab;
441 return lab;
442 }
444 /**
445 * Converts an CLAB value to a RGB color value.
446 * <P>
447 * This is the reverse operation to rgbToClab.
448 * @param clab CLAB value.
449 * @return RGB value.
450 * @see #rgbToClab
451 */
452 /*
453 static long clabToRGB(const CLAB &clab)
454 {
455 float L=clab.L;
456 float a=clab.A;
457 float b=clab.B;
459 float var_Y=(L+16.0f)/116.0f;
460 float var_X=a/500.0f+var_Y;
461 float var_Z=var_Y-b/200.0f;
463 float var_yPow3=(float)pow(var_Y, 3.0);
464 float var_xPow3=(float)pow(var_X, 3.0);
465 float var_zPow3=(float)pow(var_Z, 3.0);
467 if (var_yPow3>0.008856f)
468 var_Y=var_yPow3;
469 else
470 var_Y=(var_Y-16.0f/116.0f)/7.787f;
472 if (var_xPow3>0.008856f)
473 var_X=var_xPow3;
474 else
475 var_X=(var_X-16.0f/116.0f)/7.787f;
477 if (var_zPow3>0.008856f)
478 var_Z=var_zPow3;
479 else
480 var_Z=(var_Z-16.0f/116.0f)/7.787f;
482 float X=95.047f * var_X; //ref_X= 95.047 Observer=2�, Illuminant=D65
483 float Y=100.0f * var_Y; //ref_Y=100.000
484 float Z=108.883f * var_Z; //ref_Z=108.883
486 var_X=X/100.0f; //X = From 0 to ref_X
487 var_Y=Y/100.0f; //Y = From 0 to ref_Y
488 var_Z=Z/100.0f; //Z = From 0 to ref_Y
490 float var_R=(float)(var_X * 3.2406f + var_Y * -1.5372f + var_Z * -0.4986f);
491 float var_G=(float)(var_X * -0.9689f + var_Y * 1.8758f + var_Z * 0.0415f);
492 float var_B=(float)(var_X * 0.0557f + var_Y * -0.2040f + var_Z * 1.0570f);
494 if (var_R>0.0031308f)
495 var_R=(float)(1.055f*pow(var_R, (1.0f/2.4f))-0.055f);
496 else
497 var_R=12.92f*var_R;
499 if (var_G>0.0031308f)
500 var_G=(float)(1.055f*pow(var_G, (1.0f/2.4f))-0.055f);
501 else
502 var_G=12.92f*var_G;
504 if (var_B>0.0031308f)
505 var_B=(float)(1.055f*pow(var_B, (1.0f/2.4f))-0.055f);
506 else
507 var_B=12.92f*var_B;
509 int R = (int)lround(var_R*255.0f);
510 int G = (int)lround(var_G*255.0f);
511 int B = (int)lround(var_B*255.0f);
513 return packPixel(0xFF, R, G, B);
514 }
515 */
517 /**
518 * Sets the alpha byte of a pixel.
519 *
520 * Constructs alpha to values from 0 to 255.
521 * @param alpha Alpha value from 0 (transparent) to 255 (opaque).
522 * @param rgb The 24bit rgb color to be combined with the alpga value.
523 * @return An ARBG calor value.
524 */
525 static long setAlpha(int alpha, unsigned long rgb)
526 {
527 if (alpha>255)
528 alpha=0;
529 else if (alpha<0)
530 alpha=0;
531 return (alpha<<24)|(rgb&0xFFFFFF);
532 }
534 /**
535 * Sets the alpha byte of a pixel.
536 *
537 * Constricts alpha to values from 0 to 255.
538 * @param alpha Alpha value from 0.0f (transparent) to 1.0f (opaque).
539 * @param rgb The 24bit rgb color to be combined with the alpga value.
540 * @return An ARBG calor value.
541 */
542 static long setAlpha(float alpha, unsigned long rgb)
543 {
544 return setAlpha((int)(255.0f*alpha), rgb);
545 }
547 /**
548 * Limits the values of a,r,g,b to values from 0 to 255 and puts them
549 * together into an 32 bit integer.
550 *
551 * @param a Alpha part, the first byte.
552 * @param r Red part, the second byte.
553 * @param g Green part, the third byte.
554 * @param b Blue part, the fourth byte.
555 * @return A ARBG value packed to an int.
556 */
557 /*
558 static long packPixel(int a, int r, int g, int b)
559 {
560 if (a<0)
561 a=0;
562 else if (a>255)
563 a=255;
565 if (r<0)
566 r=0;
567 else if (r>255)
568 r=255;
570 if (g<0)
571 g=0;
572 else if (g>255)
573 g=255;
575 if (b<0)
576 b=0;
577 else if (b>255)
578 b=255;
580 return (a<<24)|(r<<16)|(g<<8)|b;
581 }
582 */
584 /**
585 * Returns the alpha component of an ARGB value.
586 *
587 * @param argb An ARGB color value.
588 * @return The alpha component, ranging from 0 to 255.
589 */
590 /*
591 static int getAlpha(unsigned long argb)
592 {
593 return (argb>>24)&0xFF;
594 }
595 */
597 /**
598 * Returns the red component of an (A)RGB value.
599 *
600 * @param rgb An (A)RGB color value.
601 * @return The red component, ranging from 0 to 255.
602 */
603 static int getRed(unsigned long rgb)
604 {
605 return (rgb>>16)&0xFF;
606 }
609 /**
610 * Returns the green component of an (A)RGB value.
611 *
612 * @param rgb An (A)RGB color value.
613 * @return The green component, ranging from 0 to 255.
614 */
615 static int getGreen(unsigned long rgb)
616 {
617 return (rgb>>8)&0xFF;
618 }
620 /**
621 * Returns the blue component of an (A)RGB value.
622 *
623 * @param rgb An (A)RGB color value.
624 * @return The blue component, ranging from 0 to 255.
625 */
626 static int getBlue(unsigned long rgb)
627 {
628 return (rgb)&0xFF;
629 }
631 /**
632 * Returns a string representation of a CLAB value.
633 *
634 * @param clab The CLAB value.
635 * @param clabSize Size of the CLAB value.
636 * @return A string representation of the CLAB value.
637 */
638 /*
639 static std::string clabToString(const CLAB &clab)
640 {
641 std::string buff;
642 char nbuf[60];
643 snprintf(nbuf, 59, "%5.3f, %5.3f, %5.3f", clab.L, clab.A, clab.B);
644 buff = nbuf;
645 return buff;
646 }
647 */
649 //########################################################################
650 //## C O L O R S I G N A T U R E (originally ColorSignature.java)
651 //########################################################################
653 /**
654 * Representation of a color signature.
655 * <br><br>
656 * This class implements a clustering algorithm based on a modified kd-tree.
657 * The splitting rule is to simply divide the given interval into two equally
658 * sized subintervals.
659 * In the <code>stageone()</code>, approximate clusters are found by building
660 * up such a tree and stopping when an interval at a node has become smaller
661 * than the allowed cluster diameter, which is given by <code>limits</code>.
662 * At this point, clusters may be split in several nodes.<br>
663 * Therefore, in <code>stagetwo()</code>, nodes that belong to several clusters
664 * are recombined by another k-d tree clustering on the prior cluster
665 * centroids. To guarantee a proper level of abstraction, clusters that contain
666 * less than 0.01% of the pixels of the entire sample are removed. Please
667 * refer to the file NOTICE to get links to further documentation.
668 *
669 * @author Gerald Friedland, Lars Knipping
670 * @version 1.02
671 *
672 * Conversion to C++ by Bob Jamison
673 *
674 */
676 /**
677 * Stage one of clustering.
678 * @param points float[][] the input points in LAB space
679 * @param depth int used for recursion, start with 0
680 * @param clusters ArrayList an Arraylist to store the clusters
681 * @param limits float[] the cluster diameters
682 */
683 static void stageone(std::vector<CLAB> &points,
684 int depth,
685 std::vector< std::vector<CLAB> > &clusters,
686 float *limits)
687 {
688 if (points.size() < 1)
689 return;
691 int dims=3;
692 int curdim=depth%dims;
693 float min = 0.0f;
694 float max = 0.0f;
695 if (curdim == 0)
696 {
697 min=points[0].C;
698 max=points[0].C;
699 // find maximum and minimum
700 for (unsigned int i=1; i<points.size(); i++) {
701 if (min>points[i].C)
702 min=points[i].C;
703 if (max<points[i].C)
704 max=points[i].C;
705 }
706 }
707 else if (curdim == 1)
708 {
709 min=points[0].L;
710 max=points[0].L;
711 // find maximum and minimum
712 for (unsigned int i=1; i<points.size(); i++) {
713 if (min>points[i].L)
714 min=points[i].L;
715 if (max<points[i].L)
716 max=points[i].L;
717 }
718 }
719 else if (curdim == 2)
720 {
721 min=points[0].A;
722 max=points[0].A;
723 // find maximum and minimum
724 for (unsigned int i=1; i<points.size(); i++) {
725 if (min>points[i].A)
726 min=points[i].A;
727 if (max<points[i].A)
728 max=points[i].A;
729 }
730 }
731 else if (curdim == 3)
732 {
733 min=points[0].B;
734 max=points[0].B;
735 // find maximum and minimum
736 for (unsigned int i=1; i<points.size(); i++) {
737 if (min>points[i].B)
738 min=points[i].B;
739 if (max<points[i].B)
740 max=points[i].B;
741 }
742 }
744 if (max-min>limits[curdim]) { // Split according to Rubner-Rule
745 // split
746 float pivotvalue=((max-min)/2.0f)+min;
748 std::vector<CLAB> smallerpoints; // allocate mem
749 std::vector<CLAB> biggerpoints;
751 for (unsigned int i=0; i<points.size(); i++) { // do actual split
752 float v = 0.0f;
753 if (curdim==0)
754 v = points[i].C;
755 else if (curdim==1)
756 v = points[i].L;
757 else if (curdim==2)
758 v = points[i].A;
759 else if (curdim==3)
760 v = points[i].B;
761 if (v <= pivotvalue) {
762 smallerpoints.push_back(points[i]);
763 } else {
764 biggerpoints.push_back(points[i]);
765 }
766 } // create subtrees
767 stageone(smallerpoints, depth+1, clusters, limits);
768 stageone(biggerpoints, depth+1, clusters, limits);
769 } else { // create leave
770 clusters.push_back(points);
771 }
772 }
774 /**
775 * Stage two of clustering.
776 * @param points float[][] the input points in LAB space
777 * @param depth int used for recursion, start with 0
778 * @param clusters ArrayList an Arraylist to store the clusters
779 * @param limits float[] the cluster diameters
780 * @param total int the total number of points as given to stageone
781 * @param threshold should be 0.01 - abstraction threshold
782 */
783 static void stagetwo(std::vector<CLAB> &points,
784 int depth,
785 std::vector< std::vector<CLAB> > &clusters,
786 float *limits, int total, float threshold)
787 {
788 if (points.size() < 1)
789 return;
791 int curdim=depth%3; // without cardinality
792 float min = 0.0f;
793 float max = 0.0f;
794 if (curdim == 0)
795 {
796 min=points[0].L;
797 max=points[0].L;
798 // find maximum and minimum
799 for (unsigned int i=1; i<points.size(); i++) {
800 if (min>points[i].L)
801 min=points[i].L;
802 if (max<points[i].L)
803 max=points[i].L;
804 }
805 }
806 else if (curdim == 1)
807 {
808 min=points[0].A;
809 max=points[0].A;
810 // find maximum and minimum
811 for (unsigned int i=1; i<points.size(); i++) {
812 if (min>points[i].A)
813 min=points[i].A;
814 if (max<points[i].A)
815 max=points[i].A;
816 }
817 }
818 else if (curdim == 2)
819 {
820 min=points[0].B;
821 max=points[0].B;
822 // find maximum and minimum
823 for (unsigned int i=1; i<points.size(); i++) {
824 if (min>points[i].B)
825 min=points[i].B;
826 if (max<points[i].B)
827 max=points[i].B;
828 }
829 }
832 if (max-min>limits[curdim]) { // Split according to Rubner-Rule
833 // split
834 float pivotvalue=((max-min)/2.0f)+min;
836 std::vector<CLAB> smallerpoints; // allocate mem
837 std::vector<CLAB> biggerpoints;
839 for (unsigned int i=0; i<points.size(); i++) { // do actual split
840 float v = 0.0f;
841 if (curdim==0)
842 v = points[i].L;
843 else if (curdim==1)
844 v = points[i].A;
845 else if (curdim==2)
846 v = points[i].B;
848 if (v <= pivotvalue) {
849 smallerpoints.push_back(points[i]);
850 } else {
851 biggerpoints.push_back(points[i]);
852 }
853 } // create subtrees
854 stagetwo(smallerpoints, depth+1, clusters, limits, total, threshold);
855 stagetwo(biggerpoints, depth+1, clusters, limits, total, threshold);
856 } else { // create leave
857 float sum=0.0;
858 for (unsigned int i=0; i<points.size(); i++) {
859 sum+=points[i].B;
860 }
861 if (((sum*100.0)/total)>=threshold) {
862 CLAB point;
863 for (unsigned int i=0; i<points.size(); i++) {
864 point.C += points[i].C;
865 point.L += points[i].L;
866 point.A += points[i].A;
867 point.B += points[i].B;
868 }
869 point.C /= points.size();
870 point.L /= points.size();
871 point.A /= points.size();
872 point.B /= points.size();
874 std::vector<CLAB> newCluster;
875 newCluster.push_back(point);
876 clusters.push_back(newCluster);
877 }
878 }
879 }
881 /**
882 * Create a color signature for the given set of pixels.
883 * @param input float[][] a set of pixels in LAB space
884 * @param length int the number of pixels that should be processed from the input
885 * @param limits float[] the cluster diameters for each dimension
886 * @param threshold float the abstraction threshold
887 * @return float[][] a color siganture containing cluster centroids in LAB space
888 */
889 static std::vector<CLAB> createSignature(std::vector<CLAB> &input,
890 float *limits, float threshold)
891 {
892 std::vector< std::vector<CLAB> > clusters1;
893 std::vector< std::vector<CLAB> > clusters2;
895 stageone(input, 0, clusters1, limits);
897 std::vector<CLAB> centroids;
898 for (unsigned int i=0; i<clusters1.size(); i++) {
899 std::vector<CLAB> cluster = clusters1[i];
900 CLAB centroid; // +1 for the cardinality
901 for (unsigned int k=0; k<cluster.size(); k++) {
902 centroid.L += cluster[k].L;
903 centroid.A += cluster[k].A;
904 centroid.B += cluster[k].B;
905 }
906 centroid.C = cluster.size();
907 centroid.L /= cluster.size();
908 centroid.A /= cluster.size();
909 centroid.B /= cluster.size();
911 centroids.push_back(centroid);
912 }
913 stagetwo(centroids, 0, clusters2, limits, input.size(), threshold); // 0.1 -> see paper by tomasi
915 std::vector<CLAB> res;
916 for (unsigned int i=0 ; i<clusters2.size() ; i++)
917 for (unsigned int j=0 ; j<clusters2[i].size() ; j++)
918 res.push_back(clusters2[i][j]);
920 return res;
921 }
925 //########################################################################
926 //## S I O X S E G M E N T A T O R (originally SioxSegmentator.java)
927 //########################################################################
929 //### NOTE: Doxygen comments are in siox.h
931 /** Confidence corresponding to a certain foreground region (equals one). */
932 const float SioxSegmentator::CERTAIN_FOREGROUND_CONFIDENCE=1.0f;
934 /** Confidence for a region likely being foreground.*/
935 const float SioxSegmentator::FOREGROUND_CONFIDENCE=0.8f;
937 /** Confidence for foreground or background type being equally likely.*/
938 const float SioxSegmentator::UNKNOWN_REGION_CONFIDENCE=0.5f;
940 /** Confidence for a region likely being background.*/
941 const float SioxSegmentator::BACKGROUND_CONFIDENCE=0.1f;
943 /** Confidence corresponding to a certain background reagion (equals zero). */
944 const float SioxSegmentator::CERTAIN_BACKGROUND_CONFIDENCE=0.0f;
947 SioxSegmentator::SioxSegmentator(int w, int h, float *limitsArg, int limitsSize)
948 {
950 imgWidth = w;
951 imgHeight = h;
952 pixelCount = imgWidth * imgHeight;
954 labelField = new int[pixelCount];
956 if (!limitsArg) {
957 limits = new float[3];
958 limits[0] = 0.64f;
959 limits[1] = 1.28f;
960 limits[2] = 2.56f;
961 } else {
962 limits = new float[limitsSize];
963 for (int i=0 ; i<limitsSize ; i++)
964 limits[i] = limitsArg[i];
965 }
967 float negLimits[3];
968 negLimits[0] = -limits[0];
969 negLimits[1] = -limits[1];
970 negLimits[2] = -limits[2];
971 clusterSize = sqrEuclidianDist(limits, 3, negLimits);
973 segmentated=false;
975 origImage = NULL;
976 }
978 SioxSegmentator::~SioxSegmentator()
979 {
980 delete labelField;
981 delete limits;
982 delete origImage;
983 }
986 /**
987 * Error logging
988 */
989 void SioxSegmentator::error(char *fmt, ...)
990 {
991 va_list args;
992 fprintf(stderr, "SioxSegmentator error:");
993 va_start(args, fmt);
994 vfprintf(stderr, fmt, args);
995 va_end(args) ;
996 fprintf(stderr, "\n");
997 }
999 /**
1000 * Trace logging
1001 */
1002 void SioxSegmentator::trace(char *fmt, ...)
1003 {
1004 va_list args;
1005 fprintf(stderr, "SioxSegmentator:");
1006 va_start(args, fmt);
1007 vfprintf(stderr, fmt, args);
1008 va_end(args) ;
1009 fprintf(stderr, "\n");
1010 }
1014 bool SioxSegmentator::segmentate(unsigned long *image,
1015 float *cm,
1016 int smoothness, double sizeFactorToKeep)
1017 {
1018 segmentated=false;
1020 hs.clear();
1022 // save image for drb
1023 origImage=new long[pixelCount];
1024 for (int i=0 ; i<pixelCount ; i++)
1025 origImage[i] = image[i];
1027 // create color signatures
1028 for (int i=0; i<pixelCount; i++) {
1029 if (cm[i]<=BACKGROUND_CONFIDENCE)
1030 knownBg.push_back(rgbToClab(image[i]));
1031 else if (cm[i]>=FOREGROUND_CONFIDENCE)
1032 knownFg.push_back(rgbToClab(image[i]));
1033 }
1035 bgSignature = createSignature(knownBg, limits, BACKGROUND_CONFIDENCE);
1036 fgSignature = createSignature(knownFg, limits, BACKGROUND_CONFIDENCE);
1038 if (bgSignature.size() < 1) {
1039 // segmentation impossible
1040 return false;
1041 }
1043 // classify using color signatures,
1044 // classification cached in hashmap for drb and speedup purposes
1045 for (int i=0; i<pixelCount; i++) {
1046 if (cm[i]>=FOREGROUND_CONFIDENCE) {
1047 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1048 continue;
1049 }
1050 if (cm[i]>BACKGROUND_CONFIDENCE) {
1051 bool isBackground=true;
1052 std::map<unsigned long, Tupel>::iterator iter = hs.find(i);
1053 Tupel tupel(0.0f, 0, 0.0f, 0);
1054 if (iter == hs.end()) {
1055 CLAB lab = rgbToClab(image[i]);
1056 float minBg = sqrEuclidianDist(lab, bgSignature[0]);
1057 int minIndex=0;
1058 for (unsigned int j=1; j<bgSignature.size(); j++) {
1059 float d = sqrEuclidianDist(lab, bgSignature[j]);
1060 if (d<minBg) {
1061 minBg=d;
1062 minIndex=j;
1063 }
1064 }
1065 tupel.minBgDist=minBg;
1066 tupel.indexMinBg=minIndex;
1067 float minFg = 1.0e6f;
1068 minIndex=-1;
1069 for (unsigned int j=0 ; j<fgSignature.size() ; j++) {
1070 float d = sqrEuclidianDist(lab, fgSignature[j]);
1071 if (d<minFg) {
1072 minFg=d;
1073 minIndex=j;
1074 }
1075 }
1076 tupel.minFgDist=minFg;
1077 tupel.indexMinFg=minIndex;
1078 if (fgSignature.size()==0) {
1079 isBackground=(minBg<=clusterSize);
1080 // remove next line to force behaviour of old algorithm
1081 error("foreground signature does not exist");
1082 return false;
1083 } else {
1084 isBackground=minBg<minFg;
1085 }
1086 hs[image[i]] = tupel;
1087 } else {
1088 isBackground=tupel.minBgDist<=tupel.minFgDist;
1089 }
1090 if (isBackground) {
1091 cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1092 } else {
1093 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1094 }
1095 } else {
1096 cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1097 }
1098 }
1100 // postprocessing
1101 smoothcm(cm, imgWidth, imgHeight, 0.33f, 0.33f, 0.33f); // average
1102 normalizeMatrix(cm, pixelCount);
1103 erode(cm, imgWidth, imgHeight);
1104 keepOnlyLargeComponents(cm, UNKNOWN_REGION_CONFIDENCE, sizeFactorToKeep);
1106 for (int i=0; i<smoothness; i++) {
1107 smoothcm(cm, imgWidth, imgHeight, 0.33f, 0.33f, 0.33f); // average
1108 }
1110 normalizeMatrix(cm, pixelCount);
1112 for (int i=0; i<pixelCount; i++) {
1113 if (cm[i]>=UNKNOWN_REGION_CONFIDENCE) {
1114 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1115 } else {
1116 cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1117 }
1118 }
1120 keepOnlyLargeComponents(cm, UNKNOWN_REGION_CONFIDENCE, sizeFactorToKeep);
1121 fillColorRegions(cm, image);
1122 dilate(cm, imgWidth, imgHeight);
1124 segmentated=true;
1125 return true;
1126 }
1130 void SioxSegmentator::keepOnlyLargeComponents(float *cm,
1131 float threshold,
1132 double sizeFactorToKeep)
1133 {
1134 int idx = 0;
1135 for (int i=0 ; i<imgHeight ; i++)
1136 for (int j=0 ; j<imgWidth ; j++)
1137 labelField[idx++] = -1;
1139 int curlabel = 0;
1140 int maxregion= 0;
1141 int maxblob = 0;
1143 // slow but easy to understand:
1144 std::vector<int> labelSizes;
1145 for (int i=0 ; i<pixelCount ; i++) {
1146 regionCount=0;
1147 if (labelField[i]==-1 && cm[i]>=threshold) {
1148 regionCount=depthFirstSearch(cm, i, threshold, curlabel++);
1149 labelSizes.push_back(regionCount);
1150 }
1152 if (regionCount>maxregion) {
1153 maxregion=regionCount;
1154 maxblob=curlabel-1;
1155 }
1156 }
1158 for (int i=0 ; i<pixelCount ; i++) {
1159 if (labelField[i]!=-1) {
1160 // remove if the component is to small
1161 if (labelSizes[labelField[i]]*sizeFactorToKeep < maxregion)
1162 cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1164 // add maxblob always to foreground
1165 if (labelField[i]==maxblob)
1166 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1167 }
1168 }
1169 }
1173 int SioxSegmentator::depthFirstSearch(float *cm, int i, float threshold, int curLabel)
1174 {
1175 // stores positions of labeled pixels, where the neighbours
1176 // should still be checked for processing:
1177 std::vector<int> pixelsToVisit;
1178 int componentSize=0;
1179 if (labelField[i]==-1 && cm[i]>=threshold) { // label #i
1180 labelField[i] = curLabel;
1181 ++componentSize;
1182 pixelsToVisit.push_back(i);
1183 }
1184 while (pixelsToVisit.size() > 0) {
1185 int pos=pixelsToVisit[pixelsToVisit.size()-1];
1186 pixelsToVisit.erase(pixelsToVisit.end()-1);
1187 int x=pos%imgWidth;
1188 int y=pos/imgWidth;
1189 // check all four neighbours
1190 int left = pos-1;
1191 if (x-1>=0 && labelField[left]==-1 && cm[left]>=threshold) {
1192 labelField[left]=curLabel;
1193 ++componentSize;
1194 pixelsToVisit.push_back(left);
1195 }
1196 int right = pos+1;
1197 if (x+1<imgWidth && labelField[right]==-1 && cm[right]>=threshold) {
1198 labelField[right]=curLabel;
1199 ++componentSize;
1200 pixelsToVisit.push_back(right);
1201 }
1202 int top = pos-imgWidth;
1203 if (y-1>=0 && labelField[top]==-1 && cm[top]>=threshold) {
1204 labelField[top]=curLabel;
1205 ++componentSize;
1206 pixelsToVisit.push_back(top);
1207 }
1208 int bottom = pos+imgWidth;
1209 if (y+1<imgHeight && labelField[bottom]==-1
1210 && cm[bottom]>=threshold) {
1211 labelField[bottom]=curLabel;
1212 ++componentSize;
1213 pixelsToVisit.push_back(bottom);
1214 }
1215 }
1216 return componentSize;
1217 }
1219 void SioxSegmentator::subpixelRefine(int x, int y, int brushmode,
1220 float threshold, float *cf, int brushsize)
1221 {
1222 subpixelRefine(x-brushsize, y-brushsize,
1223 2*brushsize, 2*brushsize,
1224 brushmode, threshold, cf);
1225 }
1228 bool SioxSegmentator::subpixelRefine(int xa, int ya, int dx, int dy,
1229 int brushmode,
1230 float threshold, float *cf)
1231 {
1232 if (!segmentated) {
1233 error("no segmentation yet");
1234 return false;
1235 }
1237 int x0 = (xa > 0) ? xa : 0;
1238 int y0 = (ya > 0) ? ya : 0;
1240 int xTo = (imgWidth - 1 < xa+dx ) ? imgWidth-1 : xa+dx;
1241 int yTo = (imgHeight - 1 < ya+dy ) ? imgHeight-1 : ya+dy;
1243 for (int ey=y0; ey<yTo; ++ey) {
1244 for (int ex=x0; ex<xTo; ++ex) {
1245 /* we are using a rect, not necessary
1246 if (!area.contains(ex, ey)) {
1247 continue;
1248 }
1249 */
1250 unsigned long val=origImage[ey*imgWidth+ex];
1251 unsigned long orig=val;
1252 float minDistBg = 0.0f;
1253 float minDistFg = 0.0f;
1254 std::map<unsigned long, Tupel>::iterator iter = hs.find(val);
1255 if (iter != hs.end()) {
1256 minDistBg=(float) sqrt((float)iter->second.minBgDist);
1257 minDistFg=(float) sqrt((float)iter->second.minFgDist);
1258 } else {
1259 continue;
1260 }
1261 if (ADD_EDGE == brushmode) { // handle adder
1262 if (cf[ey*imgWidth+ex]<FOREGROUND_CONFIDENCE) { // postprocessing wins
1263 float alpha;
1264 if (minDistFg==0) {
1265 alpha=CERTAIN_FOREGROUND_CONFIDENCE;
1266 } else {
1267 alpha = (minDistBg/minDistFg < CERTAIN_FOREGROUND_CONFIDENCE) ?
1268 minDistBg/minDistFg : CERTAIN_FOREGROUND_CONFIDENCE;
1269 }
1270 if (alpha<threshold) { // background with certain confidence decided by user.
1271 alpha=CERTAIN_BACKGROUND_CONFIDENCE;
1272 }
1273 val = setAlpha(alpha, orig);
1274 cf[ey*imgWidth+ex]=alpha;
1275 }
1276 } else if (SUB_EDGE == brushmode) { // handle substractor
1277 if (cf[ey*imgWidth+ex]>FOREGROUND_CONFIDENCE) {
1278 // foreground, we want to take something away
1279 float alpha;
1280 if (minDistBg==0) {
1281 alpha=CERTAIN_BACKGROUND_CONFIDENCE;
1282 } else {
1283 alpha=CERTAIN_FOREGROUND_CONFIDENCE-
1284 (minDistFg/minDistBg < CERTAIN_FOREGROUND_CONFIDENCE) ? // more background -> >1
1285 minDistFg/minDistBg : CERTAIN_FOREGROUND_CONFIDENCE;
1286 // bg = gf -> 1
1287 // more fg -> <1
1288 }
1289 if (alpha<threshold) { // background with certain confidence decided by user
1290 alpha=CERTAIN_BACKGROUND_CONFIDENCE;
1291 }
1292 val = setAlpha(alpha, orig);
1293 cf[ey*imgWidth+ex]=alpha;
1294 }
1295 } else {
1296 error("unknown brush mode: %d", brushmode);
1297 return false;
1298 }
1299 }
1300 }
1301 return true;
1302 }
1306 void SioxSegmentator::fillColorRegions(float *cm, unsigned long *image)
1307 {
1308 int idx = 0;
1309 for (int i=0 ; i<imgHeight ; i++)
1310 for (int j=0 ; i<imgWidth ; j++)
1311 labelField[idx++] = -1;
1313 //int maxRegion=0; // unused now
1314 std::vector<int> pixelsToVisit;
1315 for (int i=0; i<pixelCount; i++) { // for all pixels
1316 if (labelField[i]!=-1 || cm[i]<UNKNOWN_REGION_CONFIDENCE) {
1317 continue; // already visited or bg
1318 }
1319 int origColor=image[i];
1320 int curLabel=i+1;
1321 labelField[i]=curLabel;
1322 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1323 // int componentSize = 1;
1324 pixelsToVisit.push_back(i);
1325 // depth first search to fill region
1326 while (pixelsToVisit.size() > 0) {
1327 int pos=pixelsToVisit[pixelsToVisit.size()-1];
1328 pixelsToVisit.erase(pixelsToVisit.end()-1);
1329 int x=pos%imgWidth;
1330 int y=pos/imgWidth;
1331 // check all four neighbours
1332 int left = pos-1;
1333 if (x-1>=0 && labelField[left]==-1
1334 && labcolordiff(image[left], origColor)<1.0) {
1335 labelField[left]=curLabel;
1336 cm[left]=CERTAIN_FOREGROUND_CONFIDENCE;
1337 // ++componentSize;
1338 pixelsToVisit.push_back(left);
1339 }
1340 int right = pos+1;
1341 if (x+1<imgWidth && labelField[right]==-1
1342 && labcolordiff(image[right], origColor)<1.0) {
1343 labelField[right]=curLabel;
1344 cm[right]=CERTAIN_FOREGROUND_CONFIDENCE;
1345 // ++componentSize;
1346 pixelsToVisit.push_back(right);
1347 }
1348 int top = pos-imgWidth;
1349 if (y-1>=0 && labelField[top]==-1
1350 && labcolordiff(image[top], origColor)<1.0) {
1351 labelField[top]=curLabel;
1352 cm[top]=CERTAIN_FOREGROUND_CONFIDENCE;
1353 // ++componentSize;
1354 pixelsToVisit.push_back(top);
1355 }
1356 int bottom = pos+imgWidth;
1357 if (y+1<imgHeight && labelField[bottom]==-1
1358 && labcolordiff(image[bottom], origColor)<1.0) {
1359 labelField[bottom]=curLabel;
1360 cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE;
1361 // ++componentSize;
1362 pixelsToVisit.push_back(bottom);
1363 }
1364 }
1365 //if (componentSize>maxRegion) {
1366 // maxRegion=componentSize;
1367 //}
1368 }
1369 }
1384 } //namespace siox
1385 } //namespace org