f06f943e054c22b4f34e76d04fe59a7bf3fb2a6b
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 <math.h>
21 #include <stdarg.h>
22 #include <map>
23 #include <algorithm>
26 namespace org
27 {
29 namespace siox
30 {
34 //########################################################################
35 //# C L A B
36 //########################################################################
38 /**
39 * Convert integer A, R, G, B values into an pixel value.
40 */
41 static unsigned long getRGB(int a, int r, int g, int b)
42 {
43 if (a<0) a=0;
44 else if (a>255) a=255;
46 if (r<0) r=0;
47 else if (r>255) r=255;
49 if (g<0) g=0;
50 else if (g>255) g=255;
52 if (b<0) b=0;
53 else if (b>255) b=255;
55 return (a<<24)|(r<<16)|(g<<8)|b;
56 }
60 /**
61 * Convert float A, R, G, B values (0.0-1.0) into an pixel value.
62 */
63 static unsigned long getRGB(float a, float r, float g, float b)
64 {
65 return getRGB((int)(a * 256.0),
66 (int)(r * 256.0),
67 (int)(g * 256.0),
68 (int)(b * 256.0));
69 }
73 //#########################################
74 //# Root approximations for large speedup.
75 //# By njh!
76 //#########################################
77 static const int ROOT_TAB_SIZE = 16;
78 static float cbrt_table[ROOT_TAB_SIZE +1];
80 double CieLab::cbrt(double x)
81 {
82 double y = cbrt_table[int(x*ROOT_TAB_SIZE )]; // assuming x \in [0, 1]
83 y = (2.0 * y + x/(y*y))/3.0;
84 y = (2.0 * y + x/(y*y))/3.0; // polish twice
85 return y;
86 }
88 static float qn_table[ROOT_TAB_SIZE +1];
90 double CieLab::qnrt(double x)
91 {
92 double y = qn_table[int(x*ROOT_TAB_SIZE )]; // assuming x \in [0, 1]
93 double Y = y*y;
94 y = (4.0*y + x/(Y*Y))/5.0;
95 Y = y*y;
96 y = (4.0*y + x/(Y*Y))/5.0; // polish twice
97 return y;
98 }
100 double CieLab::pow24(double x)
101 {
102 double onetwo = x*qnrt(x);
103 return onetwo*onetwo;
104 }
107 static bool _clab_inited_ = false;
108 void CieLab::init()
109 {
110 if (!_clab_inited_)
111 {
112 cbrt_table[0] = pow(float(1)/float(ROOT_TAB_SIZE*2), 0.3333);
113 qn_table[0] = pow(float(1)/float(ROOT_TAB_SIZE*2), 0.2);
114 for(int i = 1; i < ROOT_TAB_SIZE +1; i++)
115 {
116 cbrt_table[i] = pow(float(i)/float(ROOT_TAB_SIZE), 0.3333);
117 qn_table[i] = pow(float(i)/float(ROOT_TAB_SIZE), 0.2);
118 }
119 _clab_inited_ = true;
120 }
121 }
125 /**
126 * Construct this CieLab from a packed-pixel ARGB value
127 */
128 CieLab::CieLab(unsigned long rgb)
129 {
130 init();
132 int ir = (rgb>>16) & 0xff;
133 int ig = (rgb>> 8) & 0xff;
134 int ib = (rgb ) & 0xff;
136 float fr = ((float)ir) / 255.0;
137 float fg = ((float)ig) / 255.0;
138 float fb = ((float)ib) / 255.0;
140 //printf("fr:%f fg:%f fb:%f\n", fr, fg, fb);
141 if (fr > 0.04045)
142 //fr = (float) pow((fr + 0.055) / 1.055, 2.4);
143 fr = (float) pow24((fr + 0.055) / 1.055);
144 else
145 fr = fr / 12.92;
147 if (fg > 0.04045)
148 //fg = (float) pow((fg + 0.055) / 1.055, 2.4);
149 fg = (float) pow24((fg + 0.055) / 1.055);
150 else
151 fg = fg / 12.92;
153 if (fb > 0.04045)
154 //fb = (float) pow((fb + 0.055) / 1.055, 2.4);
155 fb = (float) pow24((fb + 0.055) / 1.055);
156 else
157 fb = fb / 12.92;
159 // Use white = D65
160 const float x = fr * 0.4124 + fg * 0.3576 + fb * 0.1805;
161 const float y = fr * 0.2126 + fg * 0.7152 + fb * 0.0722;
162 const float z = fr * 0.0193 + fg * 0.1192 + fb * 0.9505;
164 float vx = x / 0.95047;
165 float vy = y;
166 float vz = z / 1.08883;
168 //printf("vx:%f vy:%f vz:%f\n", vx, vy, vz);
169 if (vx > 0.008856)
170 //vx = (float) pow(vx, 0.3333);
171 vx = (float) cbrt(vx);
172 else
173 vx = (7.787 * vx) + (16.0 / 116.0);
175 if (vy > 0.008856)
176 //vy = (float) pow(vy, 0.3333);
177 vy = (float) cbrt(vy);
178 else
179 vy = (7.787 * vy) + (16.0 / 116.0);
181 if (vz > 0.008856)
182 //vz = (float) pow(vz, 0.3333);
183 vz = (float) cbrt(vz);
184 else
185 vz = (7.787 * vz) + (16.0 / 116.0);
187 C = 0;
188 L = 116.0 * vy - 16.0;
189 A = 500.0 * (vx - vy);
190 B = 200.0 * (vy - vz);
191 }
195 /**
196 * Return this CieLab's value converted to a packed-pixel ARGB value
197 */
198 unsigned long CieLab::toRGB()
199 {
200 float vy = (L + 16.0) / 116.0;
201 float vx = A / 500.0 + vy;
202 float vz = vy - B / 200.0;
204 float vx3 = vx * vx * vx;
205 float vy3 = vy * vy * vy;
206 float vz3 = vz * vz * vz;
208 if (vy3 > 0.008856)
209 vy = vy3;
210 else
211 vy = (vy - 16.0 / 116.0) / 7.787;
213 if (vx3 > 0.008856)
214 vx = vx3;
215 else
216 vx = (vx - 16.0 / 116.0) / 7.787;
218 if (vz3 > 0.008856)
219 vz = vz3;
220 else
221 vz = (vz - 16.0 / 116.0) / 7.787;
223 vx *= 0.95047; //use white = D65
224 vz *= 1.08883;
226 float vr =(float)(vx * 3.2406 + vy * -1.5372 + vz * -0.4986);
227 float vg =(float)(vx * -0.9689 + vy * 1.8758 + vz * 0.0415);
228 float vb =(float)(vx * 0.0557 + vy * -0.2040 + vz * 1.0570);
230 if (vr > 0.0031308)
231 vr = (float)(1.055 * pow(vr, (1.0 / 2.4)) - 0.055);
232 else
233 vr = 12.92 * vr;
235 if (vg > 0.0031308)
236 vg = (float)(1.055 * pow(vg, (1.0 / 2.4)) - 0.055);
237 else
238 vg = 12.92 * vg;
240 if (vb > 0.0031308)
241 vb = (float)(1.055 * pow(vb, (1.0 / 2.4)) - 0.055);
242 else
243 vb = 12.92 * vb;
245 return getRGB(0.0, vr, vg, vb);
246 }
249 /**
250 * Squared Euclidian distance between this and another color
251 */
252 float CieLab::diffSq(const CieLab &other)
253 {
254 float sum=0.0;
255 sum += (L - other.L) * (L - other.L);
256 sum += (A - other.A) * (A - other.A);
257 sum += (B - other.B) * (B - other.B);
258 return sum;
259 }
261 /**
262 * Computes squared euclidian distance in CieLab space for two colors
263 * given as RGB values.
264 */
265 float CieLab::diffSq(unsigned int rgb1, unsigned int rgb2)
266 {
267 CieLab c1(rgb1);
268 CieLab c2(rgb2);
269 float euclid = c1.diffSq(c2);
270 return euclid;
271 }
274 /**
275 * Computes squared euclidian distance in CieLab space for two colors
276 * given as RGB values.
277 */
278 float CieLab::diff(unsigned int rgb0, unsigned int rgb1)
279 {
280 return (float) sqrt(diffSq(rgb0, rgb1));
281 }
285 //########################################################################
286 //# T U P E L
287 //########################################################################
289 /**
290 * Helper class for storing the minimum distances to a cluster centroid
291 * in background and foreground and the index to the centroids in each
292 * signature for a given color.
293 */
294 class Tupel {
295 public:
297 Tupel()
298 {
299 minBgDist = 0.0f;
300 indexMinBg = 0;
301 minFgDist = 0.0f;
302 indexMinFg = 0;
303 }
304 Tupel(float minBgDistArg, long indexMinBgArg,
305 float minFgDistArg, long indexMinFgArg)
306 {
307 minBgDist = minBgDistArg;
308 indexMinBg = indexMinBgArg;
309 minFgDist = minFgDistArg;
310 indexMinFg = indexMinFgArg;
311 }
312 Tupel(const Tupel &other)
313 {
314 minBgDist = other.minBgDist;
315 indexMinBg = other.indexMinBg;
316 minFgDist = other.minFgDist;
317 indexMinFg = other.indexMinFg;
318 }
319 Tupel &operator=(const Tupel &other)
320 {
321 minBgDist = other.minBgDist;
322 indexMinBg = other.indexMinBg;
323 minFgDist = other.minFgDist;
324 indexMinFg = other.indexMinFg;
325 return *this;
326 }
327 virtual ~Tupel()
328 {}
330 float minBgDist;
331 long indexMinBg;
332 float minFgDist;
333 long indexMinFg;
335 };
339 //########################################################################
340 //# S I O X I M A G E
341 //########################################################################
343 /**
344 * Create an image with the given width and height
345 */
346 SioxImage::SioxImage(unsigned int widthArg, unsigned int heightArg)
347 {
348 init(width, height);
349 }
351 /**
352 * Copy constructor
353 */
354 SioxImage::SioxImage(const SioxImage &other)
355 {
356 pixdata = NULL;
357 cmdata = NULL;
358 assign(other);
359 }
361 /**
362 * Assignment
363 */
364 SioxImage &SioxImage::operator=(const SioxImage &other)
365 {
366 assign(other);
367 return *this;
368 }
371 /**
372 * Clean up after use.
373 */
374 SioxImage::~SioxImage()
375 {
376 if (pixdata) delete[] pixdata;
377 if (cmdata) delete[] cmdata;
378 }
380 /**
381 * Error logging
382 */
383 void SioxImage::error(char *fmt, ...)
384 {
385 char msgbuf[256];
386 va_list args;
387 va_start(args, fmt);
388 vsnprintf(msgbuf, 255, fmt, args);
389 va_end(args) ;
390 #ifdef HAVE_GLIB
391 g_warning("SioxImage error: %s\n", msgbuf);
392 #else
393 fprintf(stderr, "SioxImage error: %s\n", msgbuf);
394 #endif
395 }
398 /**
399 * Returns true if the previous operation on this image
400 * was successful, else false.
401 */
402 bool SioxImage::isValid()
403 {
404 return valid;
405 }
407 /**
408 * Sets whether an operation was successful, and whether
409 * this image should be considered a valid one.
410 * was successful, else false.
411 */
412 void SioxImage::setValid(bool val)
413 {
414 valid = val;
415 }
418 /**
419 * Set a pixel at the x,y coordinates to the given value.
420 * If the coordinates are out of range, do nothing.
421 */
422 void SioxImage::setPixel(unsigned int x,
423 unsigned int y,
424 unsigned int pixval)
425 {
426 if (x >= width || y >= height)
427 {
428 error("setPixel: out of bounds (%d,%d)/(%d,%d)",
429 x, y, width, height);
430 return;
431 }
432 unsigned long offset = width * y + x;
433 pixdata[offset] = pixval;
434 }
436 /**
437 * Set a pixel at the x,y coordinates to the given r, g, b values.
438 * If the coordinates are out of range, do nothing.
439 */
440 void SioxImage::setPixel(unsigned int x, unsigned int y,
441 unsigned int a,
442 unsigned int r,
443 unsigned int g,
444 unsigned int b)
445 {
446 if (x >= width || y >= height)
447 {
448 error("setPixel: out of bounds (%d,%d)/(%d,%d)",
449 x, y, width, height);
450 return;
451 }
452 unsigned long offset = width * y + x;
453 unsigned int pixval = ((a << 24) & 0xff000000) |
454 ((r << 16) & 0x00ff0000) |
455 ((g << 8) & 0x0000ff00) |
456 ((b ) & 0x000000ff);
457 pixdata[offset] = pixval;
458 }
462 /**
463 * Get a pixel at the x,y coordinates given. If
464 * the coordinates are out of range, return 0;
465 */
466 unsigned int SioxImage::getPixel(unsigned int x, unsigned int y)
467 {
468 if (x >= width || y >= height)
469 {
470 error("getPixel: out of bounds (%d,%d)/(%d,%d)",
471 x, y, width, height);
472 return 0L;
473 }
474 unsigned long offset = width * y + x;
475 return pixdata[offset];
476 }
478 /**
479 * Return the image data buffer
480 */
481 unsigned int *SioxImage::getImageData()
482 {
483 return pixdata;
484 }
486 /**
487 * Set a confidence value at the x,y coordinates to the given value.
488 * If the coordinates are out of range, do nothing.
489 */
490 void SioxImage::setConfidence(unsigned int x,
491 unsigned int y,
492 float confval)
493 {
494 if (x >= width || y >= height)
495 {
496 error("setConfidence: out of bounds (%d,%d)/(%d,%d)",
497 x, y, width, height);
498 return;
499 }
500 unsigned long offset = width * y + x;
501 cmdata[offset] = confval;
502 }
504 /**
505 * Get a confidence valueat the x,y coordinates given. If
506 * the coordinates are out of range, return 0;
507 */
508 float SioxImage::getConfidence(unsigned int x, unsigned int y)
509 {
510 if (x >= width || y >= height)
511 {
512 g_warning("getConfidence: out of bounds (%d,%d)/(%d,%d)",
513 x, y, width, height);
514 return 0.0;
515 }
516 unsigned long offset = width * y + x;
517 return cmdata[offset];
518 }
520 /**
521 * Return the confidence data buffer
522 */
523 float *SioxImage::getConfidenceData()
524 {
525 return cmdata;
526 }
529 /**
530 * Return the width of this image
531 */
532 int SioxImage::getWidth()
533 {
534 return width;
535 }
537 /**
538 * Return the height of this image
539 */
540 int SioxImage::getHeight()
541 {
542 return height;
543 }
545 /**
546 * Initialize values. Used by constructors
547 */
548 void SioxImage::init(unsigned int widthArg, unsigned int heightArg)
549 {
550 valid = true;
551 width = widthArg;
552 height = heightArg;
553 imageSize = width * height;
554 pixdata = new unsigned int[imageSize];
555 cmdata = new float[imageSize];
556 for (unsigned long i=0 ; i<imageSize ; i++)
557 {
558 pixdata[i] = 0;
559 cmdata[i] = 0.0;
560 }
561 }
563 /**
564 * Assign values to that of another
565 */
566 void SioxImage::assign(const SioxImage &other)
567 {
568 if (pixdata) delete[] pixdata;
569 if (cmdata) delete[] cmdata;
570 valid = other.valid;
571 width = other.width;
572 height = other.height;
573 imageSize = width * height;
574 pixdata = new unsigned int[imageSize];
575 cmdata = new float[imageSize];
576 for (unsigned long i=0 ; i<imageSize ; i++)
577 {
578 pixdata[i] = other.pixdata[i];
579 cmdata[i] = other.cmdata[i];
580 }
581 }
584 /**
585 * Write the image to a PPM file
586 */
587 bool SioxImage::writePPM(const std::string fileName)
588 {
590 FILE *f = fopen(fileName.c_str(), "wb");
591 if (!f)
592 return false;
594 fprintf(f, "P6 %d %d 255\n", width, height);
596 for (unsigned int y=0 ; y<height; y++)
597 {
598 for (unsigned int x=0 ; x<width ; x++)
599 {
600 unsigned int rgb = getPixel(x, y);
601 //unsigned int alpha = (rgb>>24) & 0xff;
602 unsigned int r = ((rgb>>16) & 0xff);
603 unsigned int g = ((rgb>> 8) & 0xff);
604 unsigned int b = ((rgb ) & 0xff);
605 fputc((unsigned char) r, f);
606 fputc((unsigned char) g, f);
607 fputc((unsigned char) b, f);
608 }
609 }
610 fclose(f);
611 return true;
612 }
615 #ifdef HAVE_GLIB
617 /**
618 * Create an image from a GdkPixbuf
619 */
620 SioxImage::SioxImage(GdkPixbuf *buf)
621 {
622 if (!buf)
623 return;
625 unsigned int width = gdk_pixbuf_get_width(buf);
626 unsigned int height = gdk_pixbuf_get_height(buf);
627 init(width, height); //DO THIS NOW!!
630 guchar *pixldata = gdk_pixbuf_get_pixels(buf);
631 int rowstride = gdk_pixbuf_get_rowstride(buf);
632 int n_channels = gdk_pixbuf_get_n_channels(buf);
634 //### Fill in the cells with RGB values
635 int row = 0;
636 for (unsigned int y=0 ; y<height ; y++)
637 {
638 guchar *p = pixldata + row;
639 for (unsigned int x=0 ; x<width ; x++)
640 {
641 int r = (int)p[0];
642 int g = (int)p[1];
643 int b = (int)p[2];
644 int alpha = (int)p[3];
646 setPixel(x, y, alpha, r, g, b);
647 p += n_channels;
648 }
649 row += rowstride;
650 }
652 }
655 /**
656 * Create a GdkPixbuf from this image
657 */
658 GdkPixbuf *SioxImage::getGdkPixbuf()
659 {
660 bool has_alpha = true;
661 int n_channels = has_alpha ? 4 : 3;
663 guchar *pixdata = (guchar *)
664 malloc(sizeof(guchar) * width * height * n_channels);
665 if (!pixdata)
666 return NULL;
668 int rowstride = width * n_channels;
670 GdkPixbuf *buf = gdk_pixbuf_new_from_data(pixdata,
671 GDK_COLORSPACE_RGB,
672 has_alpha, 8, width, height,
673 rowstride, NULL, NULL);
675 //### Fill in the cells with RGB values
676 int row = 0;
677 for (unsigned int y=0 ; y < height ; y++)
678 {
679 guchar *p = pixdata + row;
680 for (unsigned x=0 ; x < width ; x++)
681 {
682 unsigned int rgb = getPixel(x, y);
683 p[0] = (rgb >> 16) & 0xff;//r
684 p[1] = (rgb >> 8) & 0xff;//g
685 p[2] = (rgb ) & 0xff;//b
686 if ( n_channels > 3 )
687 {
688 p[3] = (rgb >> 24) & 0xff;//a
689 }
690 p += n_channels;
691 }
692 row += rowstride;
693 }
695 return buf;
696 }
698 #endif /* GLIB */
703 //########################################################################
704 //# S I O X
705 //########################################################################
707 //##############
708 //## PUBLIC
709 //##############
711 /**
712 * Confidence corresponding to a certain foreground region (equals one).
713 */
714 const float Siox::CERTAIN_FOREGROUND_CONFIDENCE=1.0f;
716 /**
717 * Confidence for a region likely being foreground.
718 */
719 const float Siox::FOREGROUND_CONFIDENCE=0.8f;
721 /**
722 * Confidence for foreground or background type being equally likely.
723 */
724 const float Siox::UNKNOWN_REGION_CONFIDENCE=0.5f;
726 /**
727 * Confidence for a region likely being background.
728 */
729 const float Siox::BACKGROUND_CONFIDENCE=0.1f;
731 /**
732 * Confidence corresponding to a certain background reagion (equals zero).
733 */
734 const float Siox::CERTAIN_BACKGROUND_CONFIDENCE=0.0f;
736 /**
737 * Construct a Siox engine
738 */
739 Siox::Siox()
740 {
741 sioxObserver = NULL;
742 init();
743 }
745 /**
746 * Construct a Siox engine
747 */
748 Siox::Siox(SioxObserver *observer)
749 {
750 init();
751 sioxObserver = observer;
752 }
755 /**
756 *
757 */
758 Siox::~Siox()
759 {
760 cleanup();
761 }
764 /**
765 * Error logging
766 */
767 void Siox::error(char *fmt, ...)
768 {
769 char msgbuf[256];
770 va_list args;
771 va_start(args, fmt);
772 vsnprintf(msgbuf, 255, fmt, args);
773 va_end(args) ;
774 #ifdef HAVE_GLIB
775 g_warning("Siox error: %s\n", msgbuf);
776 #else
777 fprintf(stderr, "Siox error: %s\n", msgbuf);
778 #endif
779 }
781 /**
782 * Trace logging
783 */
784 void Siox::trace(char *fmt, ...)
785 {
786 char msgbuf[256];
787 va_list args;
788 va_start(args, fmt);
789 vsnprintf(msgbuf, 255, fmt, args);
790 va_end(args) ;
791 #ifdef HAVE_GLIB
792 g_message("Siox: %s\n", msgbuf);
793 #else
794 fprintf(stdout, "Siox: %s\n", msgbuf);
795 #endif
796 }
800 /**
801 * Progress reporting
802 */
803 bool Siox::progressReport(float percentCompleted)
804 {
805 if (!sioxObserver)
806 return true;
808 bool ret = sioxObserver->progress(percentCompleted);
810 if (!ret)
811 {
812 trace("User selected abort");
813 keepGoing = false;
814 }
816 return ret;
817 }
822 /**
823 * Extract the foreground of the original image, according
824 * to the values in the confidence matrix.
825 */
826 SioxImage Siox::extractForeground(const SioxImage &originalImage,
827 unsigned int backgroundFillColor)
828 {
829 trace("### Start");
831 init();
832 keepGoing = true;
834 SioxImage workImage = originalImage;
836 //# fetch some info from the image
837 width = workImage.getWidth();
838 height = workImage.getHeight();
839 pixelCount = width * height;
840 image = workImage.getImageData();
841 cm = workImage.getConfidenceData();
842 labelField = new int[pixelCount];
844 trace("### Creating signatures");
846 //#### create color signatures
847 std::vector<CieLab> knownBg;
848 std::vector<CieLab> knownFg;
849 CieLab *imageClab = new CieLab[pixelCount];
850 for (unsigned long i=0 ; i<pixelCount ; i++)
851 {
852 float conf = cm[i];
853 unsigned int pix = image[i];
854 CieLab lab(pix);
855 imageClab[i] = lab;
856 if (conf <= BACKGROUND_CONFIDENCE)
857 knownBg.push_back(lab);
858 else if (conf >= FOREGROUND_CONFIDENCE)
859 knownFg.push_back(lab);
860 }
862 /*
863 std::vector<CieLab> imageClab;
864 for (int y = 0 ; y < workImage.getHeight() ; y++)
865 for (int x = 0 ; x < workImage.getWidth() ; x++)
866 {
867 float cm = workImage.getConfidence(x, y);
868 unsigned int pix = workImage.getPixel(x, y);
869 CieLab lab(pix);
870 imageClab.push_back(lab);
871 if (cm <= BACKGROUND_CONFIDENCE)
872 knownBg.push_back(lab); //note: uses CieLab(rgb)
873 else if (cm >= FOREGROUND_CONFIDENCE)
874 knownFg.push_back(lab);
875 }
876 */
878 if (!progressReport(10.0))
879 {
880 error("User aborted");
881 workImage.setValid(false);
882 delete[] imageClab;
883 delete[] labelField;
884 return workImage;
885 }
887 trace("knownBg:%d knownFg:%d", knownBg.size(), knownFg.size());
890 std::vector<CieLab> bgSignature ;
891 if (!colorSignature(knownBg, bgSignature, 3))
892 {
893 error("Could not create background signature");
894 workImage.setValid(false);
895 delete[] imageClab;
896 delete[] labelField;
897 return workImage;
898 }
900 if (!progressReport(30.0))
901 {
902 error("User aborted");
903 workImage.setValid(false);
904 delete[] imageClab;
905 delete[] labelField;
906 return workImage;
907 }
910 std::vector<CieLab> fgSignature ;
911 if (!colorSignature(knownFg, fgSignature, 3))
912 {
913 error("Could not create foreground signature");
914 workImage.setValid(false);
915 delete[] imageClab;
916 delete[] labelField;
917 return workImage;
918 }
920 //trace("### bgSignature:%d", bgSignature.size());
922 if (bgSignature.size() < 1)
923 {
924 // segmentation impossible
925 error("Signature size is < 1. Segmentation is impossible");
926 workImage.setValid(false);
927 delete[] imageClab;
928 delete[] labelField;
929 return workImage;
930 }
932 if (!progressReport(30.0))
933 {
934 error("User aborted");
935 workImage.setValid(false);
936 delete[] imageClab;
937 delete[] labelField;
938 return workImage;
939 }
942 // classify using color signatures,
943 // classification cached in hashmap for drb and speedup purposes
944 trace("### Analyzing image");
946 std::map<unsigned int, Tupel> hs;
948 unsigned int progressResolution = pixelCount / 10;
950 for (unsigned int i=0; i<pixelCount; i++)
951 {
952 if (i % progressResolution == 0)
953 {
954 float progress =
955 30.0 + 60.0 * (float)i / (float)pixelCount;
956 //trace("### progress:%f", progress);
957 if (!progressReport(progress))
958 {
959 error("User aborted");
960 delete[] imageClab;
961 delete[] labelField;
962 workImage.setValid(false);
963 return workImage;
964 }
965 }
967 if (cm[i] >= FOREGROUND_CONFIDENCE)
968 {
969 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
970 }
971 else if (cm[i] <= BACKGROUND_CONFIDENCE)
972 {
973 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
974 }
975 else // somewhere in between
976 {
977 bool isBackground = true;
978 std::map<unsigned int, Tupel>::iterator iter = hs.find(i);
979 if (iter != hs.end()) //found
980 {
981 Tupel tupel = iter->second;
982 isBackground = tupel.minBgDist <= tupel.minFgDist;
983 }
984 else
985 {
986 CieLab lab = imageClab[i];
987 float minBg = lab.diffSq(bgSignature[0]);
988 int minIndex = 0;
989 for (unsigned int j=1; j<bgSignature.size() ; j++)
990 {
991 float d = lab.diffSq(bgSignature[j]);
992 if (d<minBg)
993 {
994 minBg = d;
995 minIndex = j;
996 }
997 }
998 Tupel tupel(0.0f, 0, 0.0f, 0);
999 tupel.minBgDist = minBg;
1000 tupel.indexMinBg = minIndex;
1001 float minFg = 1.0e6f;
1002 minIndex = -1;
1003 for (unsigned int j = 0 ; j < fgSignature.size() ; j++)
1004 {
1005 float d = lab.diffSq(fgSignature[j]);
1006 if (d < minFg)
1007 {
1008 minFg = d;
1009 minIndex = j;
1010 }
1011 }
1012 tupel.minFgDist = minFg;
1013 tupel.indexMinFg = minIndex;
1014 if (fgSignature.size() == 0)
1015 {
1016 isBackground = (minBg <= clusterSize);
1017 // remove next line to force behaviour of old algorithm
1018 //error("foreground signature does not exist");
1019 //delete[] labelField;
1020 //workImage.setValid(false);
1021 //return workImage;
1022 }
1023 else
1024 {
1025 isBackground = minBg < minFg;
1026 }
1027 hs[image[i]] = tupel;
1028 }
1030 if (isBackground)
1031 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
1032 else
1033 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
1034 }
1035 }
1038 delete[] imageClab;
1040 trace("### postProcessing");
1043 //#### postprocessing
1044 smooth(cm, width, height, 0.333f, 0.333f, 0.333f); // average
1045 normalizeMatrix(cm, pixelCount);
1046 erode(cm, width, height);
1047 keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.0/*sizeFactorToKeep*/);
1049 //for (int i=0; i < 2/*smoothness*/; i++)
1050 // smooth(cm, width, height, 0.333f, 0.333f, 0.333f); // average
1052 normalizeMatrix(cm, pixelCount);
1054 for (unsigned int i=0; i < pixelCount; i++)
1055 {
1056 if (cm[i] >= UNKNOWN_REGION_CONFIDENCE)
1057 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
1058 else
1059 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
1060 }
1062 keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.5/*sizeFactorToKeep*/);
1063 fillColorRegions();
1064 dilate(cm, width, height);
1066 if (!progressReport(100.0))
1067 {
1068 error("User aborted");
1069 delete[] labelField;
1070 workImage.setValid(false);
1071 return workImage;
1072 }
1075 //#### We are done. Now clear everything but the background
1076 for (unsigned long i = 0; i<pixelCount ; i++)
1077 {
1078 float conf = cm[i];
1079 if (conf < FOREGROUND_CONFIDENCE)
1080 image[i] = backgroundFillColor;
1081 }
1083 delete[] labelField;
1085 trace("### Done");
1086 keepGoing = false;
1087 return workImage;
1088 }
1092 //##############
1093 //## PRIVATE
1094 //##############
1096 /**
1097 * Initialize the Siox engine to its 'pristine' state.
1098 * Performed at the beginning of extractForeground().
1099 */
1100 void Siox::init()
1101 {
1102 limits[0] = 0.64f;
1103 limits[1] = 1.28f;
1104 limits[2] = 2.56f;
1106 float negLimits[3];
1107 negLimits[0] = -limits[0];
1108 negLimits[1] = -limits[1];
1109 negLimits[2] = -limits[2];
1111 clusterSize = sqrEuclidianDist(limits, 3, negLimits);
1112 }
1115 /**
1116 * Clean up any debris from processing.
1117 */
1118 void Siox::cleanup()
1119 {
1120 }
1125 /**
1126 * Stage 1 of the color signature work. 'dims' will be either
1127 * 2 for grays, or 3 for colors
1128 */
1129 void Siox::colorSignatureStage1(CieLab *points,
1130 unsigned int leftBase,
1131 unsigned int rightBase,
1132 unsigned int recursionDepth,
1133 unsigned int *clusterCount,
1134 const unsigned int dims)
1135 {
1137 unsigned int currentDim = recursionDepth % dims;
1138 CieLab point = points[leftBase];
1139 float min = point(currentDim);
1140 float max = min;
1142 for (unsigned int i = leftBase + 1; i < rightBase ; i++)
1143 {
1144 point = points[i];
1145 float curval = point(currentDim);
1146 if (curval < min) min = curval;
1147 if (curval > max) max = curval;
1148 }
1150 //Do the Rubner-rule split (sounds like a dance)
1151 if (max - min > limits[currentDim])
1152 {
1153 float pivotPoint = (min + max) / 2.0; //average
1154 unsigned int left = leftBase;
1155 unsigned int right = rightBase - 1;
1157 //# partition points according to the dimension
1158 while (true)
1159 {
1160 while ( true )
1161 {
1162 point = points[left];
1163 if (point(currentDim) > pivotPoint)
1164 break;
1165 left++;
1166 }
1167 while ( true )
1168 {
1169 point = points[right];
1170 if (point(currentDim) <= pivotPoint)
1171 break;
1172 right--;
1173 }
1175 if (left > right)
1176 break;
1178 point = points[left];
1179 points[left] = points[right];
1180 points[right] = point;
1182 left++;
1183 right--;
1184 }
1186 //# Recurse and create sub-trees
1187 colorSignatureStage1(points, leftBase, left,
1188 recursionDepth + 1, clusterCount, dims);
1189 colorSignatureStage1(points, left, rightBase,
1190 recursionDepth + 1, clusterCount, dims);
1191 }
1192 else
1193 {
1194 //create a leaf
1195 CieLab newpoint;
1197 newpoint.C = rightBase - leftBase;
1199 for (; leftBase < rightBase ; leftBase++)
1200 {
1201 newpoint.add(points[leftBase]);
1202 }
1204 //printf("clusters:%d\n", *clusters);
1206 if (newpoint.C != 0)
1207 newpoint.mul(1.0 / (float)newpoint.C);
1208 points[*clusterCount] = newpoint;
1209 (*clusterCount)++;
1210 }
1211 }
1215 /**
1216 * Stage 2 of the color signature work
1217 */
1218 void Siox::colorSignatureStage2(CieLab *points,
1219 unsigned int leftBase,
1220 unsigned int rightBase,
1221 unsigned int recursionDepth,
1222 unsigned int *clusterCount,
1223 const float threshold,
1224 const unsigned int dims)
1225 {
1228 unsigned int currentDim = recursionDepth % dims;
1229 CieLab point = points[leftBase];
1230 float min = point(currentDim);
1231 float max = min;
1233 for (unsigned int i = leftBase+ 1; i < rightBase; i++)
1234 {
1235 point = points[i];
1236 float curval = point(currentDim);
1237 if (curval < min) min = curval;
1238 if (curval > max) max = curval;
1239 }
1241 //Do the Rubner-rule split (sounds like a dance)
1242 if (max - min > limits[currentDim])
1243 {
1244 float pivotPoint = (min + max) / 2.0; //average
1245 unsigned int left = leftBase;
1246 unsigned int right = rightBase - 1;
1248 //# partition points according to the dimension
1249 while (true)
1250 {
1251 while ( true )
1252 {
1253 point = points[left];
1254 if (point(currentDim) > pivotPoint)
1255 break;
1256 left++;
1257 }
1258 while ( true )
1259 {
1260 point = points[right];
1261 if (point(currentDim) <= pivotPoint)
1262 break;
1263 right--;
1264 }
1266 if (left > right)
1267 break;
1269 point = points[left];
1270 points[left] = points[right];
1271 points[right] = point;
1273 left++;
1274 right--;
1275 }
1277 //# Recurse and create sub-trees
1278 colorSignatureStage2(points, leftBase, left,
1279 recursionDepth + 1, clusterCount, threshold, dims);
1280 colorSignatureStage2(points, left, rightBase,
1281 recursionDepth + 1, clusterCount, threshold, dims);
1282 }
1283 else
1284 {
1285 //### Create a leaf
1286 unsigned int sum = 0;
1287 for (unsigned int i = leftBase; i < rightBase; i++)
1288 sum += points[i].C;
1290 if ((float)sum >= threshold)
1291 {
1292 float scale = (float)(rightBase - leftBase);
1293 CieLab newpoint;
1295 for (; leftBase < rightBase; leftBase++)
1296 newpoint.add(points[leftBase]);
1298 if (scale != 0.0)
1299 newpoint.mul(1.0 / scale);
1300 points[*clusterCount] = newpoint;
1301 (*clusterCount)++;
1302 }
1303 }
1304 }
1308 /**
1309 * Main color signature method
1310 */
1311 bool Siox::colorSignature(const std::vector<CieLab> &inputVec,
1312 std::vector<CieLab> &result,
1313 const unsigned int dims)
1314 {
1316 unsigned int length = inputVec.size();
1318 if (length < 1) // no error. just don't do anything
1319 return true;
1321 CieLab *input = (CieLab *) malloc(length * sizeof(CieLab));
1323 if (!input)
1324 {
1325 error("Could not allocate buffer for signature");
1326 return false;
1327 }
1328 for (unsigned int i=0 ; i < length ; i++)
1329 input[i] = inputVec[i];
1331 unsigned int stage1length = 0;
1332 colorSignatureStage1(input, 0, length, 0,
1333 &stage1length, dims);
1335 unsigned int stage2length = 0;
1336 colorSignatureStage2(input, 0, stage1length, 0,
1337 &stage2length, length * 0.001, dims);
1339 result.clear();
1340 for (unsigned int i=0 ; i < stage2length ; i++)
1341 result.push_back(input[i]);
1343 free(input);
1345 return true;
1346 }
1350 /**
1351 *
1352 */
1353 void Siox::keepOnlyLargeComponents(float threshold,
1354 double sizeFactorToKeep)
1355 {
1356 for (unsigned long idx = 0 ; idx<pixelCount ; idx++)
1357 labelField[idx] = -1;
1359 int curlabel = 0;
1360 int maxregion= 0;
1361 int maxblob = 0;
1363 // slow but easy to understand:
1364 std::vector<int> labelSizes;
1365 for (unsigned long i=0 ; i<pixelCount ; i++)
1366 {
1367 int regionCount = 0;
1368 if (labelField[i] == -1 && cm[i] >= threshold)
1369 {
1370 regionCount = depthFirstSearch(i, threshold, curlabel++);
1371 labelSizes.push_back(regionCount);
1372 }
1374 if (regionCount>maxregion)
1375 {
1376 maxregion = regionCount;
1377 maxblob = curlabel-1;
1378 }
1379 }
1381 for (unsigned int i=0 ; i<pixelCount ; i++)
1382 {
1383 if (labelField[i] != -1)
1384 {
1385 // remove if the component is to small
1386 if (labelSizes[labelField[i]] * sizeFactorToKeep < maxregion)
1387 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
1389 // add maxblob always to foreground
1390 if (labelField[i] == maxblob)
1391 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
1392 }
1393 }
1395 }
1399 int Siox::depthFirstSearch(int startPos,
1400 float threshold, int curLabel)
1401 {
1402 // stores positions of labeled pixels, where the neighbours
1403 // should still be checked for processing:
1405 //trace("startPos:%d threshold:%f curLabel:%d",
1406 // startPos, threshold, curLabel);
1408 std::vector<int> pixelsToVisit;
1409 int componentSize = 0;
1411 if (labelField[startPos]==-1 && cm[startPos]>=threshold)
1412 {
1413 labelField[startPos] = curLabel;
1414 componentSize++;
1415 pixelsToVisit.push_back(startPos);
1416 }
1419 while (pixelsToVisit.size() > 0)
1420 {
1421 int pos = pixelsToVisit[pixelsToVisit.size() - 1];
1422 pixelsToVisit.erase(pixelsToVisit.end() - 1);
1423 unsigned int x = pos % width;
1424 unsigned int y = pos / width;
1426 // check all four neighbours
1427 int left = pos - 1;
1428 if (((int)x)-1>=0 && labelField[left]==-1 && cm[left]>=threshold)
1429 {
1430 labelField[left]=curLabel;
1431 componentSize++;
1432 pixelsToVisit.push_back(left);
1433 }
1435 int right = pos + 1;
1436 if (x+1 < width && labelField[right]==-1 && cm[right]>=threshold)
1437 {
1438 labelField[right]=curLabel;
1439 componentSize++;
1440 pixelsToVisit.push_back(right);
1441 }
1443 int top = pos - width;
1444 if (((int)y)-1>=0 && labelField[top]==-1 && cm[top]>=threshold)
1445 {
1446 labelField[top]=curLabel;
1447 componentSize++;
1448 pixelsToVisit.push_back(top);
1449 }
1451 int bottom = pos + width;
1452 if (y+1 < height && labelField[bottom]==-1
1453 && cm[bottom]>=threshold)
1454 {
1455 labelField[bottom]=curLabel;
1456 componentSize++;
1457 pixelsToVisit.push_back(bottom);
1458 }
1460 }
1461 return componentSize;
1462 }
1466 /**
1467 *
1468 */
1469 void Siox::fillColorRegions()
1470 {
1471 for (unsigned long idx = 0 ; idx<pixelCount ; idx++)
1472 labelField[idx] = -1;
1474 //int maxRegion=0; // unused now
1475 std::vector<int> pixelsToVisit;
1476 for (unsigned long i=0; i<pixelCount; i++)
1477 { // for all pixels
1478 if (labelField[i]!=-1 || cm[i]<UNKNOWN_REGION_CONFIDENCE)
1479 {
1480 continue; // already visited or bg
1481 }
1483 unsigned int origColor = image[i];
1484 unsigned long curLabel = i+1;
1485 labelField[i] = curLabel;
1486 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
1488 // int componentSize = 1;
1489 pixelsToVisit.push_back(i);
1490 // depth first search to fill region
1491 while (pixelsToVisit.size() > 0)
1492 {
1493 int pos = pixelsToVisit[pixelsToVisit.size() - 1];
1494 pixelsToVisit.erase(pixelsToVisit.end() - 1);
1495 unsigned int x=pos % width;
1496 unsigned int y=pos / width;
1497 // check all four neighbours
1498 int left = pos-1;
1499 if (((int)x)-1 >= 0 && labelField[left] == -1
1500 && CieLab::diff(image[left], origColor)<1.0)
1501 {
1502 labelField[left]=curLabel;
1503 cm[left]=CERTAIN_FOREGROUND_CONFIDENCE;
1504 // ++componentSize;
1505 pixelsToVisit.push_back(left);
1506 }
1507 int right = pos+1;
1508 if (x+1 < width && labelField[right]==-1
1509 && CieLab::diff(image[right], origColor)<1.0)
1510 {
1511 labelField[right]=curLabel;
1512 cm[right]=CERTAIN_FOREGROUND_CONFIDENCE;
1513 // ++componentSize;
1514 pixelsToVisit.push_back(right);
1515 }
1516 int top = pos - width;
1517 if (((int)y)-1>=0 && labelField[top]==-1
1518 && CieLab::diff(image[top], origColor)<1.0)
1519 {
1520 labelField[top]=curLabel;
1521 cm[top]=CERTAIN_FOREGROUND_CONFIDENCE;
1522 // ++componentSize;
1523 pixelsToVisit.push_back(top);
1524 }
1525 int bottom = pos + width;
1526 if (y+1 < height && labelField[bottom]==-1
1527 && CieLab::diff(image[bottom], origColor)<1.0)
1528 {
1529 labelField[bottom]=curLabel;
1530 cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE;
1531 // ++componentSize;
1532 pixelsToVisit.push_back(bottom);
1533 }
1534 }
1535 //if (componentSize>maxRegion) {
1536 // maxRegion=componentSize;
1537 //}
1538 }
1539 }
1544 /**
1545 * Applies the morphological dilate operator.
1546 *
1547 * Can be used to close small holes in the given confidence matrix.
1548 */
1549 void Siox::dilate(float *cm, int xres, int yres)
1550 {
1552 for (int y=0; y<yres; y++)
1553 {
1554 for (int x=0; x<xres-1; x++)
1555 {
1556 int idx=(y*xres)+x;
1557 if (cm[idx+1]>cm[idx])
1558 cm[idx]=cm[idx+1];
1559 }
1560 }
1562 for (int y=0; y<yres; y++)
1563 {
1564 for (int x=xres-1; x>=1; x--)
1565 {
1566 int idx=(y*xres)+x;
1567 if (cm[idx-1]>cm[idx])
1568 cm[idx]=cm[idx-1];
1569 }
1570 }
1572 for (int y=0; y<yres-1; y++)
1573 {
1574 for (int x=0; x<xres; x++)
1575 {
1576 int idx=(y*xres)+x;
1577 if (cm[((y+1)*xres)+x] > cm[idx])
1578 cm[idx]=cm[((y+1)*xres)+x];
1579 }
1580 }
1582 for (int y=yres-1; y>=1; y--)
1583 {
1584 for (int x=0; x<xres; x++)
1585 {
1586 int idx=(y*xres)+x;
1587 if (cm[((y-1)*xres)+x] > cm[idx])
1588 cm[idx]=cm[((y-1)*xres)+x];
1589 }
1590 }
1591 }
1593 /**
1594 * Applies the morphological erode operator.
1595 */
1596 void Siox::erode(float *cm, int xres, int yres)
1597 {
1598 for (int y=0; y<yres; y++)
1599 {
1600 for (int x=0; x<xres-1; x++)
1601 {
1602 int idx=(y*xres)+x;
1603 if (cm[idx+1] < cm[idx])
1604 cm[idx]=cm[idx+1];
1605 }
1606 }
1607 for (int y=0; y<yres; y++)
1608 {
1609 for (int x=xres-1; x>=1; x--)
1610 {
1611 int idx=(y*xres)+x;
1612 if (cm[idx-1] < cm[idx])
1613 cm[idx]=cm[idx-1];
1614 }
1615 }
1616 for (int y=0; y<yres-1; y++)
1617 {
1618 for (int x=0; x<xres; x++)
1619 {
1620 int idx=(y*xres)+x;
1621 if (cm[((y+1)*xres)+x] < cm[idx])
1622 cm[idx]=cm[((y+1)*xres)+x];
1623 }
1624 }
1625 for (int y=yres-1; y>=1; y--)
1626 {
1627 for (int x=0; x<xres; x++)
1628 {
1629 int idx=(y*xres)+x;
1630 if (cm[((y-1)*xres)+x] < cm[idx])
1631 cm[idx]=cm[((y-1)*xres)+x];
1632 }
1633 }
1634 }
1638 /**
1639 * Normalizes the matrix to values to [0..1].
1640 */
1641 void Siox::normalizeMatrix(float *cm, int cmSize)
1642 {
1643 float max= -1000000.0f;
1644 for (int i=0; i<cmSize; i++)
1645 if (cm[i] > max) max=cm[i];
1647 //good to use STL, but max() is not iterative
1648 //float max = *std::max(cm, cm + cmSize);
1650 if (max<=0.0 || max==1.0)
1651 return;
1653 float alpha=1.00f/max;
1654 premultiplyMatrix(alpha, cm, cmSize);
1655 }
1657 /**
1658 * Multiplies matrix with the given scalar.
1659 */
1660 void Siox::premultiplyMatrix(float alpha, float *cm, int cmSize)
1661 {
1662 for (int i=0; i<cmSize; i++)
1663 cm[i]=alpha*cm[i];
1664 }
1666 /**
1667 * Blurs confidence matrix with a given symmetrically weighted kernel.
1668 *
1669 * In the standard case confidence matrix entries are between 0...1 and
1670 * the weight factors sum up to 1.
1671 */
1672 void Siox::smooth(float *cm, int xres, int yres,
1673 float f1, float f2, float f3)
1674 {
1675 for (int y=0; y<yres; y++)
1676 {
1677 for (int x=0; x<xres-2; x++)
1678 {
1679 int idx=(y*xres)+x;
1680 cm[idx]=f1*cm[idx]+f2*cm[idx+1]+f3*cm[idx+2];
1681 }
1682 }
1683 for (int y=0; y<yres; y++)
1684 {
1685 for (int x=xres-1; x>=2; x--)
1686 {
1687 int idx=(y*xres)+x;
1688 cm[idx]=f3*cm[idx-2]+f2*cm[idx-1]+f1*cm[idx];
1689 }
1690 }
1691 for (int y=0; y<yres-2; y++)
1692 {
1693 for (int x=0; x<xres; x++)
1694 {
1695 int idx=(y*xres)+x;
1696 cm[idx]=f1*cm[idx]+f2*cm[((y+1)*xres)+x]+f3*cm[((y+2)*xres)+x];
1697 }
1698 }
1699 for (int y=yres-1; y>=2; y--)
1700 {
1701 for (int x=0; x<xres; x++)
1702 {
1703 int idx=(y*xres)+x;
1704 cm[idx]=f3*cm[((y-2)*xres)+x]+f2*cm[((y-1)*xres)+x]+f1*cm[idx];
1705 }
1706 }
1707 }
1709 /**
1710 * Squared Euclidian distance of p and q.
1711 */
1712 float Siox::sqrEuclidianDist(float *p, int pSize, float *q)
1713 {
1714 float sum=0.0;
1715 for (int i=0; i<pSize; i++)
1716 {
1717 float v = p[i] - q[i];
1718 sum += v*v;
1719 }
1720 return sum;
1721 }
1728 } // namespace siox
1729 } // namespace org
1731 //########################################################################
1732 //# E N D O F F I L E
1733 //########################################################################