31d99ca038aa357dc2f08b42ab01c079dbb2b738
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-segmentator.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<long, CLAB> RGB_TO_LAB;
50 //forward decls
51 static void premultiplyMatrix(float alpha, float *cm, int cmSize);
52 //static float colordiffsq(long rgb0, long rgb1);
53 //static int getAlpha(long argb);
54 static int getRed(long rgb);
55 static int getGreen(long rgb);
56 static int getBlue(long rgb);
57 //static long packPixel(int a, int r, int g, int b);
58 static CLAB rgbToClab(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(long rgb1, 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(long rgb0, 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(long rgb)
378 {
379 std::map<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, 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, 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(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(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(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(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-segmentator.h
932 SioxSegmentator::SioxSegmentator(int w, int h, float *limitsArg, int limitsSize)
933 {
935 imgWidth = w;
936 imgHeight = h;
937 labelField = new int[imgWidth*imgHeight];
939 if (!limitsArg) {
940 limits = new float[3];
941 limits[0] = 0.64f;
942 limits[1] = 1.28f;
943 limits[2] = 2.56f;
944 } else {
945 limits = new float[limitsSize];
946 for (int i=0 ; i<limitsSize ; i++)
947 limits[i] = limitsArg[i];
948 }
950 float negLimits[3];
951 negLimits[0] = -limits[0];
952 negLimits[1] = -limits[1];
953 negLimits[2] = -limits[2];
954 clusterSize = sqrEuclidianDist(limits, 3, negLimits);
956 segmentated=false;
958 origImage = NULL;
959 }
961 SioxSegmentator::~SioxSegmentator()
962 {
963 delete labelField;
964 delete limits;
965 delete origImage;
966 }
969 /**
970 * Error logging
971 */
972 void SioxSegmentator::error(char *fmt, ...)
973 {
974 va_list args;
975 fprintf(stderr, "SioxSegmentator error:");
976 va_start(args, fmt);
977 vfprintf(stderr, fmt, args);
978 va_end(args) ;
979 fprintf(stderr, "\n");
980 }
982 /**
983 * Trace logging
984 */
985 void SioxSegmentator::trace(char *fmt, ...)
986 {
987 va_list args;
988 fprintf(stderr, "SioxSegmentator:");
989 va_start(args, fmt);
990 vfprintf(stderr, fmt, args);
991 va_end(args) ;
992 fprintf(stderr, "\n");
993 }
997 bool SioxSegmentator::segmentate(long *image, int imageSize,
998 float *cm, int cmSize,
999 int smoothness, double sizeFactorToKeep)
1000 {
1001 segmentated=false;
1003 hs.clear();
1005 // save image for drb
1006 origImage=new long[imageSize];
1007 for (int i=0 ; i<imageSize ; i++)
1008 origImage[i] = image[i];
1010 // create color signatures
1011 for (int i=0; i<cmSize; i++) {
1012 if (cm[i]<=BACKGROUND_CONFIDENCE)
1013 knownBg.push_back(rgbToClab(image[i]));
1014 else if (cm[i]>=FOREGROUND_CONFIDENCE)
1015 knownFg.push_back(rgbToClab(image[i]));
1016 }
1018 bgSignature = createSignature(knownBg, limits, BACKGROUND_CONFIDENCE);
1019 fgSignature = createSignature(knownFg, limits, BACKGROUND_CONFIDENCE);
1021 if (bgSignature.size() < 1) {
1022 // segmentation impossible
1023 return false;
1024 }
1026 // classify using color signatures,
1027 // classification cached in hashmap for drb and speedup purposes
1028 for (int i=0; i<cmSize; i++) {
1029 if (cm[i]>=FOREGROUND_CONFIDENCE) {
1030 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1031 continue;
1032 }
1033 if (cm[i]>BACKGROUND_CONFIDENCE) {
1034 bool isBackground=true;
1035 std::map<long, Tupel>::iterator iter = hs.find(i);
1036 Tupel tupel(0.0f, 0, 0.0f, 0);
1037 if (iter == hs.end()) {
1038 CLAB lab = rgbToClab(image[i]);
1039 float minBg = sqrEuclidianDist(lab, bgSignature[0]);
1040 int minIndex=0;
1041 for (unsigned int j=1; j<bgSignature.size(); j++) {
1042 float d = sqrEuclidianDist(lab, bgSignature[j]);
1043 if (d<minBg) {
1044 minBg=d;
1045 minIndex=j;
1046 }
1047 }
1048 tupel.minBgDist=minBg;
1049 tupel.indexMinBg=minIndex;
1050 float minFg = 1.0e6f;
1051 minIndex=-1;
1052 for (unsigned int j=0 ; j<fgSignature.size() ; j++) {
1053 float d = sqrEuclidianDist(lab, fgSignature[j]);
1054 if (d<minFg) {
1055 minFg=d;
1056 minIndex=j;
1057 }
1058 }
1059 tupel.minFgDist=minFg;
1060 tupel.indexMinFg=minIndex;
1061 if (fgSignature.size()==0) {
1062 isBackground=(minBg<=clusterSize);
1063 // remove next line to force behaviour of old algorithm
1064 error("foreground signature does not exist");
1065 return false;
1066 } else {
1067 isBackground=minBg<minFg;
1068 }
1069 hs[image[i]] = tupel;
1070 } else {
1071 isBackground=tupel.minBgDist<=tupel.minFgDist;
1072 }
1073 if (isBackground) {
1074 cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1075 } else {
1076 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1077 }
1078 } else {
1079 cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1080 }
1081 }
1083 // postprocessing
1084 smoothcm(cm, imgWidth, imgHeight, 0.33f, 0.33f, 0.33f); // average
1085 normalizeMatrix(cm, cmSize);
1086 erode(cm, imgWidth, imgHeight);
1087 keepOnlyLargeComponents(cm, cmSize, UNKNOWN_REGION_CONFIDENCE, sizeFactorToKeep);
1089 for (int i=0; i<smoothness; i++) {
1090 smoothcm(cm, imgWidth, imgHeight, 0.33f, 0.33f, 0.33f); // average
1091 }
1093 normalizeMatrix(cm, cmSize);
1095 for (int i=0; i<cmSize; i++) {
1096 if (cm[i]>=UNKNOWN_REGION_CONFIDENCE) {
1097 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1098 } else {
1099 cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1100 }
1101 }
1103 keepOnlyLargeComponents(cm, cmSize, UNKNOWN_REGION_CONFIDENCE, sizeFactorToKeep);
1104 fillColorRegions(cm, cmSize, image);
1105 dilate(cm, imgWidth, imgHeight);
1107 segmentated=true;
1108 return true;
1109 }
1113 void SioxSegmentator::keepOnlyLargeComponents(float *cm, int cmSize,
1114 float threshold,
1115 double sizeFactorToKeep)
1116 {
1117 int idx = 0;
1118 for (int i=0 ; i<imgHeight ; i++)
1119 for (int j=0 ; j<imgWidth ; j++)
1120 labelField[idx++] = -1;
1122 int curlabel = 0;
1123 int maxregion= 0;
1124 int maxblob = 0;
1126 // slow but easy to understand:
1127 std::vector<int> labelSizes;
1128 for (int i=0 ; i<cmSize ; i++) {
1129 regionCount=0;
1130 if (labelField[i]==-1 && cm[i]>=threshold) {
1131 regionCount=depthFirstSearch(cm, i, threshold, curlabel++);
1132 labelSizes.push_back(regionCount);
1133 }
1135 if (regionCount>maxregion) {
1136 maxregion=regionCount;
1137 maxblob=curlabel-1;
1138 }
1139 }
1141 for (int i=0 ; i<cmSize ; i++) {
1142 if (labelField[i]!=-1) {
1143 // remove if the component is to small
1144 if (labelSizes[labelField[i]]*sizeFactorToKeep < maxregion)
1145 cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1147 // add maxblob always to foreground
1148 if (labelField[i]==maxblob)
1149 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1150 }
1151 }
1152 }
1156 int SioxSegmentator::depthFirstSearch(float *cm, int i, float threshold, int curLabel)
1157 {
1158 // stores positions of labeled pixels, where the neighbours
1159 // should still be checked for processing:
1160 std::vector<int> pixelsToVisit;
1161 int componentSize=0;
1162 if (labelField[i]==-1 && cm[i]>=threshold) { // label #i
1163 labelField[i] = curLabel;
1164 ++componentSize;
1165 pixelsToVisit.push_back(i);
1166 }
1167 while (pixelsToVisit.size() > 0) {
1168 int pos=pixelsToVisit[pixelsToVisit.size()-1];
1169 pixelsToVisit.erase(pixelsToVisit.end()-1);
1170 int x=pos%imgWidth;
1171 int y=pos/imgWidth;
1172 // check all four neighbours
1173 int left = pos-1;
1174 if (x-1>=0 && labelField[left]==-1 && cm[left]>=threshold) {
1175 labelField[left]=curLabel;
1176 ++componentSize;
1177 pixelsToVisit.push_back(left);
1178 }
1179 int right = pos+1;
1180 if (x+1<imgWidth && labelField[right]==-1 && cm[right]>=threshold) {
1181 labelField[right]=curLabel;
1182 ++componentSize;
1183 pixelsToVisit.push_back(right);
1184 }
1185 int top = pos-imgWidth;
1186 if (y-1>=0 && labelField[top]==-1 && cm[top]>=threshold) {
1187 labelField[top]=curLabel;
1188 ++componentSize;
1189 pixelsToVisit.push_back(top);
1190 }
1191 int bottom = pos+imgWidth;
1192 if (y+1<imgHeight && labelField[bottom]==-1
1193 && cm[bottom]>=threshold) {
1194 labelField[bottom]=curLabel;
1195 ++componentSize;
1196 pixelsToVisit.push_back(bottom);
1197 }
1198 }
1199 return componentSize;
1200 }
1202 void SioxSegmentator::subpixelRefine(int x, int y, int brushmode,
1203 float threshold, float *cf, int brushsize)
1204 {
1205 subpixelRefine(x-brushsize, y-brushsize,
1206 2*brushsize, 2*brushsize,
1207 brushmode, threshold, cf);
1208 }
1211 bool SioxSegmentator::subpixelRefine(int xa, int ya, int dx, int dy,
1212 int brushmode,
1213 float threshold, float *cf)
1214 {
1215 if (!segmentated) {
1216 error("no segmentation yet");
1217 return false;
1218 }
1220 int x0 = (xa > 0) ? xa : 0;
1221 int y0 = (ya > 0) ? ya : 0;
1223 int xTo = (imgWidth - 1 < xa+dx ) ? imgWidth-1 : xa+dx;
1224 int yTo = (imgHeight - 1 < ya+dy ) ? imgHeight-1 : ya+dy;
1226 for (int ey=y0; ey<yTo; ++ey) {
1227 for (int ex=x0; ex<xTo; ++ex) {
1228 /* we are using a rect, not necessary
1229 if (!area.contains(ex, ey)) {
1230 continue;
1231 }
1232 */
1233 long val=origImage[ey*imgWidth+ex];
1234 long orig=val;
1235 float minDistBg = 0.0f;
1236 float minDistFg = 0.0f;
1237 std::map<long, Tupel>::iterator iter = hs.find(val);
1238 if (iter != hs.end()) {
1239 minDistBg=(float) sqrt((float)iter->second.minBgDist);
1240 minDistFg=(float) sqrt((float)iter->second.minFgDist);
1241 } else {
1242 continue;
1243 }
1244 if (ADD_EDGE == brushmode) { // handle adder
1245 if (cf[ey*imgWidth+ex]<FOREGROUND_CONFIDENCE) { // postprocessing wins
1246 float alpha;
1247 if (minDistFg==0) {
1248 alpha=CERTAIN_FOREGROUND_CONFIDENCE;
1249 } else {
1250 alpha = (minDistBg/minDistFg < CERTAIN_FOREGROUND_CONFIDENCE) ?
1251 minDistBg/minDistFg : CERTAIN_FOREGROUND_CONFIDENCE;
1252 }
1253 if (alpha<threshold) { // background with certain confidence decided by user.
1254 alpha=CERTAIN_BACKGROUND_CONFIDENCE;
1255 }
1256 val = setAlpha(alpha, orig);
1257 cf[ey*imgWidth+ex]=alpha;
1258 }
1259 } else if (SUB_EDGE == brushmode) { // handle substractor
1260 if (cf[ey*imgWidth+ex]>FOREGROUND_CONFIDENCE) {
1261 // foreground, we want to take something away
1262 float alpha;
1263 if (minDistBg==0) {
1264 alpha=CERTAIN_BACKGROUND_CONFIDENCE;
1265 } else {
1266 alpha=CERTAIN_FOREGROUND_CONFIDENCE-
1267 (minDistFg/minDistBg < CERTAIN_FOREGROUND_CONFIDENCE) ? // more background -> >1
1268 minDistFg/minDistBg : CERTAIN_FOREGROUND_CONFIDENCE;
1269 // bg = gf -> 1
1270 // more fg -> <1
1271 }
1272 if (alpha<threshold) { // background with certain confidence decided by user
1273 alpha=CERTAIN_BACKGROUND_CONFIDENCE;
1274 }
1275 val = setAlpha(alpha, orig);
1276 cf[ey*imgWidth+ex]=alpha;
1277 }
1278 } else {
1279 error("unknown brush mode: %d", brushmode);
1280 return false;
1281 }
1282 }
1283 }
1284 return true;
1285 }
1289 void SioxSegmentator::fillColorRegions(float *cm, int cmSize, long *image)
1290 {
1291 int idx = 0;
1292 for (int i=0 ; i<imgHeight ; i++)
1293 for (int j=0 ; i<imgWidth ; j++)
1294 labelField[idx++] = -1;
1296 //int maxRegion=0; // unused now
1297 std::vector<int> pixelsToVisit;
1298 for (int i=0; i<cmSize; i++) { // for all pixels
1299 if (labelField[i]!=-1 || cm[i]<UNKNOWN_REGION_CONFIDENCE) {
1300 continue; // already visited or bg
1301 }
1302 int origColor=image[i];
1303 int curLabel=i+1;
1304 labelField[i]=curLabel;
1305 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1306 // int componentSize = 1;
1307 pixelsToVisit.push_back(i);
1308 // depth first search to fill region
1309 while (pixelsToVisit.size() > 0) {
1310 int pos=pixelsToVisit[pixelsToVisit.size()-1];
1311 pixelsToVisit.erase(pixelsToVisit.end()-1);
1312 int x=pos%imgWidth;
1313 int y=pos/imgWidth;
1314 // check all four neighbours
1315 int left = pos-1;
1316 if (x-1>=0 && labelField[left]==-1
1317 && labcolordiff(image[left], origColor)<1.0) {
1318 labelField[left]=curLabel;
1319 cm[left]=CERTAIN_FOREGROUND_CONFIDENCE;
1320 // ++componentSize;
1321 pixelsToVisit.push_back(left);
1322 }
1323 int right = pos+1;
1324 if (x+1<imgWidth && labelField[right]==-1
1325 && labcolordiff(image[right], origColor)<1.0) {
1326 labelField[right]=curLabel;
1327 cm[right]=CERTAIN_FOREGROUND_CONFIDENCE;
1328 // ++componentSize;
1329 pixelsToVisit.push_back(right);
1330 }
1331 int top = pos-imgWidth;
1332 if (y-1>=0 && labelField[top]==-1
1333 && labcolordiff(image[top], origColor)<1.0) {
1334 labelField[top]=curLabel;
1335 cm[top]=CERTAIN_FOREGROUND_CONFIDENCE;
1336 // ++componentSize;
1337 pixelsToVisit.push_back(top);
1338 }
1339 int bottom = pos+imgWidth;
1340 if (y+1<imgHeight && labelField[bottom]==-1
1341 && labcolordiff(image[bottom], origColor)<1.0) {
1342 labelField[bottom]=curLabel;
1343 cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE;
1344 // ++componentSize;
1345 pixelsToVisit.push_back(bottom);
1346 }
1347 }
1348 //if (componentSize>maxRegion) {
1349 // maxRegion=componentSize;
1350 //}
1351 }
1352 }
1367 } //namespace siox
1368 } //namespace org