d460d15579b9c8db4a65a76534b4b77f91f4001a
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>
25 namespace org
26 {
28 namespace siox
29 {
33 //########################################################################
34 //# C L A B
35 //########################################################################
37 static std::map<unsigned long, CLAB> clabLookupTable;
39 /**
40 * Convert integer A, R, G, B values into an pixel value.
41 */
42 static unsigned long getRGB(int a, int r, int g, int b)
43 {
44 if (a<0) a=0;
45 else if (a>255) a=255;
47 if (r<0) r=0;
48 else if (r>255) r=255;
50 if (g<0) g=0;
51 else if (g>255) g=255;
53 if (b<0) b=0;
54 else if (b>255) b=255;
56 return (a<<24)|(r<<16)|(g<<8)|b;
57 }
61 /**
62 * Convert float A, R, G, B values (0.0-1.0) into an pixel value.
63 */
64 static unsigned long getRGB(float a, float r, float g, float b)
65 {
66 return getRGB((int)(a * 256.0),
67 (int)(r * 256.0),
68 (int)(g * 256.0),
69 (int)(b * 256.0));
70 }
74 /**
75 * Construct this CLAB from a packed-pixel ARGB value
76 */
77 CLAB::CLAB(unsigned long rgb)
78 {
79 //First try looking up in the cache
80 std::map<unsigned long, CLAB>::iterator iter;
81 iter = clabLookupTable.find(rgb);
82 if (iter != clabLookupTable.end())
83 {
84 CLAB res = iter->second;
85 C = res.C;
86 L = res.L;
87 A = res.A;
88 B = res.B;
89 }
92 int ir = (rgb>>16) & 0xff;
93 int ig = (rgb>> 8) & 0xff;
94 int ib = (rgb ) & 0xff;
96 float fr = ((float)ir) / 255.0;
97 float fg = ((float)ig) / 255.0;
98 float fb = ((float)ib) / 255.0;
100 if (fr > 0.04045)
101 fr = (float) pow((fr + 0.055) / 1.055, 2.4);
102 else
103 fr = fr / 12.92;
105 if (fg > 0.04045)
106 fg = (float) pow((fg + 0.055) / 1.055, 2.4);
107 else
108 fg = fg / 12.92;
110 if (fb > 0.04045)
111 fb = (float) pow((fb + 0.055) / 1.055, 2.4);
112 else
113 fb = fb / 12.92;
115 fr = fr * 100.0;
116 fg = fg * 100.0;
117 fb = fb * 100.0;
119 // Use white = D65
120 float x = fr * 0.4124 + fg * 0.3576 + fb * 0.1805;
121 float y = fr * 0.2126 + fg * 0.7152 + fb * 0.0722;
122 float z = fr * 0.0193 + fg * 0.1192 + fb * 0.9505;
124 float vx = x / 95.047;
125 float vy = y / 100.000;
126 float vz = z / 108.883;
128 if (vx > 0.008856)
129 vx = (float) pow(vx, 0.3333);
130 else
131 vx = (7.787 * vx) + (16.0 / 116.0);
133 if (vy > 0.008856)
134 vy = (float) pow(vy, 0.3333);
135 else
136 vy = (7.787 * vy) + (16.0 / 116.0);
138 if (vz > 0.008856)
139 vz = (float) pow(vz, 0.3333);
140 else
141 vz = (7.787 * vz) + (16.0 / 116.0);
143 C = 0;
144 L = 116.0 * vy - 16.0;
145 A = 500.0 * (vx - vy);
146 B = 200.0 * (vy - vz);
148 // Cache for next time
149 clabLookupTable[rgb] = *this;
151 }
155 /**
156 * Return this CLAB's value a a packed-pixel ARGB value
157 */
158 unsigned long CLAB::toRGB()
159 {
160 float vy = (L + 16.0) / 116.0;
161 float vx = A / 500.0 + vy;
162 float vz = vy - B / 200.0;
164 float vx3 = vx * vx * vx;
165 float vy3 = vy * vy * vy;
166 float vz3 = vz * vz * vz;
168 if (vy3 > 0.008856)
169 vy = vy3;
170 else
171 vy = (vy - 16.0 / 116.0) / 7.787;
173 if (vx3 > 0.008856)
174 vx = vx3;
175 else
176 vx = (vx - 16.0 / 116.0) / 7.787;
178 if (vz3 > 0.008856)
179 vz = vz3;
180 else
181 vz = (vz - 16.0 / 116.0) / 7.787;
183 float x = 95.047 * vx; //use white = D65
184 float y = 100.000 * vy;
185 float z = 108.883 * vz;
187 vx = x / 100.0;
188 vy = y / 100.0;
189 vz = z / 100.0;
191 float vr =(float)(vx * 3.2406 + vy * -1.5372 + vz * -0.4986);
192 float vg =(float)(vx * -0.9689 + vy * 1.8758 + vz * 0.0415);
193 float vb =(float)(vx * 0.0557 + vy * -0.2040 + vz * 1.0570);
195 if (vr > 0.0031308)
196 vr = (float)(1.055 * pow(vr, (1.0 / 2.4)) - 0.055);
197 else
198 vr = 12.92 * vr;
200 if (vg > 0.0031308)
201 vg = (float)(1.055 * pow(vg, (1.0 / 2.4)) - 0.055);
202 else
203 vg = 12.92 * vg;
205 if (vb > 0.0031308)
206 vb = (float)(1.055 * pow(vb, (1.0 / 2.4)) - 0.055);
207 else
208 vb = 12.92 * vb;
210 return getRGB(0.0, vr, vg, vb);
211 }
214 /**
215 * Computes squared euclidian distance in CLAB space for two colors
216 * given as RGB values.
217 */
218 float CLAB::diffSq(unsigned int rgb1, unsigned int rgb2)
219 {
220 CLAB c1(rgb1);
221 CLAB c2(rgb2);
222 float euclid=0.0f;
223 euclid += (c1.L - c2.L) * (c1.L - c2.L);
224 euclid += (c1.A - c2.A) * (c1.A - c2.A);
225 euclid += (c1.B - c2.B) * (c1.B - c2.B);
226 return euclid;
227 }
230 /**
231 * Computes squared euclidian distance in CLAB space for two colors
232 * given as RGB values.
233 */
234 float CLAB::diff(unsigned int rgb0, unsigned int rgb1)
235 {
236 return (float) sqrt(diffSq(rgb0, rgb1));
237 }
240 //########################################################################
241 //# T U P E L
242 //########################################################################
244 /**
245 * Helper class for storing the minimum distances to a cluster centroid
246 * in background and foreground and the index to the centroids in each
247 * signature for a given color.
248 */
249 class Tupel {
250 public:
252 Tupel()
253 {
254 minBgDist = 0.0f;
255 indexMinBg = 0;
256 minFgDist = 0.0f;
257 indexMinFg = 0;
258 }
259 Tupel(float minBgDistArg, long indexMinBgArg,
260 float minFgDistArg, long indexMinFgArg)
261 {
262 minBgDist = minBgDistArg;
263 indexMinBg = indexMinBgArg;
264 minFgDist = minFgDistArg;
265 indexMinFg = indexMinFgArg;
266 }
267 Tupel(const Tupel &other)
268 {
269 minBgDist = other.minBgDist;
270 indexMinBg = other.indexMinBg;
271 minFgDist = other.minFgDist;
272 indexMinFg = other.indexMinFg;
273 }
274 Tupel &operator=(const Tupel &other)
275 {
276 minBgDist = other.minBgDist;
277 indexMinBg = other.indexMinBg;
278 minFgDist = other.minFgDist;
279 indexMinFg = other.indexMinFg;
280 return *this;
281 }
282 virtual ~Tupel()
283 {}
285 float minBgDist;
286 long indexMinBg;
287 float minFgDist;
288 long indexMinFg;
290 };
294 //########################################################################
295 //# S I O X I M A G E
296 //########################################################################
298 /**
299 * Create an image with the given width and height
300 */
301 SioxImage::SioxImage(unsigned int widthArg, unsigned int heightArg)
302 {
303 init(width, height);
304 }
306 /**
307 * Copy constructor
308 */
309 SioxImage::SioxImage(const SioxImage &other)
310 {
311 pixdata = NULL;
312 cmdata = NULL;
313 assign(other);
314 }
316 /**
317 * Assignment
318 */
319 SioxImage &SioxImage::operator=(const SioxImage &other)
320 {
321 assign(other);
322 return *this;
323 }
326 /**
327 * Clean up after use.
328 */
329 SioxImage::~SioxImage()
330 {
331 if (pixdata) delete[] pixdata;
332 if (cmdata) delete[] cmdata;
333 }
335 /**
336 * Returns true if the previous operation on this image
337 * was successful, else false.
338 */
339 bool SioxImage::isValid()
340 {
341 return valid;
342 }
344 /**
345 * Sets whether an operation was successful, and whether
346 * this image should be considered a valid one.
347 * was successful, else false.
348 */
349 void SioxImage::setValid(bool val)
350 {
351 valid = val;
352 }
355 /**
356 * Set a pixel at the x,y coordinates to the given value.
357 * If the coordinates are out of range, do nothing.
358 */
359 void SioxImage::setPixel(unsigned int x,
360 unsigned int y,
361 unsigned int pixval)
362 {
363 if (x > width || y > height)
364 return;
365 unsigned long offset = width * y + x;
366 pixdata[offset] = pixval;
367 }
369 /**
370 * Set a pixel at the x,y coordinates to the given r, g, b values.
371 * If the coordinates are out of range, do nothing.
372 */
373 void SioxImage::setPixel(unsigned int x, unsigned int y,
374 unsigned int a,
375 unsigned int r,
376 unsigned int g,
377 unsigned int b)
378 {
379 if (x > width || y > height)
380 return;
381 unsigned long offset = width * y + x;
382 unsigned int pixval = ((a << 24) & 0xff000000) |
383 ((r << 16) & 0x00ff0000) |
384 ((g << 8) & 0x0000ff00) |
385 ((b ) & 0x000000ff);
386 pixdata[offset] = pixval;
387 }
391 /**
392 * Get a pixel at the x,y coordinates given. If
393 * the coordinates are out of range, return 0;
394 */
395 unsigned int SioxImage::getPixel(unsigned int x, unsigned int y)
396 {
397 if (x > width || y > height)
398 return 0L;
399 unsigned long offset = width * y + x;
400 return pixdata[offset];
401 }
403 /**
404 * Return the image data buffer
405 */
406 unsigned int *SioxImage::getImageData()
407 {
408 return pixdata;
409 }
411 /**
412 * Set a confidence value at the x,y coordinates to the given value.
413 * If the coordinates are out of range, do nothing.
414 */
415 void SioxImage::setConfidence(unsigned int x,
416 unsigned int y,
417 float confval)
418 {
419 if (x > width || y > height)
420 return;
421 unsigned long offset = width * y + x;
422 cmdata[offset] = confval;
423 }
425 /**
426 * Get a confidence valueat the x,y coordinates given. If
427 * the coordinates are out of range, return 0;
428 */
429 float SioxImage::getConfidence(unsigned int x, unsigned int y)
430 {
431 if (x > width || y > height)
432 return 0.0;
433 unsigned long offset = width * y + x;
434 return cmdata[offset];
435 }
437 /**
438 * Return the confidence data buffer
439 */
440 float *SioxImage::getConfidenceData()
441 {
442 return cmdata;
443 }
446 /**
447 * Return the width of this image
448 */
449 int SioxImage::getWidth()
450 {
451 return width;
452 }
454 /**
455 * Return the height of this image
456 */
457 int SioxImage::getHeight()
458 {
459 return height;
460 }
462 /**
463 * Initialize values. Used by constructors
464 */
465 void SioxImage::init(unsigned int widthArg, unsigned int heightArg)
466 {
467 valid = true;
468 width = widthArg;
469 height = heightArg;
470 imageSize = width * height;
471 pixdata = new unsigned int[imageSize];
472 cmdata = new float[imageSize];
473 for (unsigned long i=0 ; i<imageSize ; i++)
474 {
475 pixdata[i] = 0;
476 cmdata[i] = 0.0;
477 }
478 }
480 /**
481 * Assign values to that of another
482 */
483 void SioxImage::assign(const SioxImage &other)
484 {
485 if (pixdata) delete[] pixdata;
486 if (cmdata) delete[] cmdata;
487 valid = other.valid;
488 width = other.width;
489 height = other.height;
490 imageSize = width * height;
491 pixdata = new unsigned int[imageSize];
492 cmdata = new float[imageSize];
493 for (unsigned long i=0 ; i<imageSize ; i++)
494 {
495 pixdata[i] = other.pixdata[i];
496 cmdata[i] = other.cmdata[i];
497 }
498 }
501 /**
502 * Write the image to a PPM file
503 */
504 bool SioxImage::writePPM(const std::string fileName)
505 {
507 FILE *f = fopen(fileName.c_str(), "wb");
508 if (!f)
509 return false;
511 fprintf(f, "P6 %d %d 255\n", width, height);
513 for (unsigned int y=0 ; y<height; y++)
514 {
515 for (unsigned int x=0 ; x<width ; x++)
516 {
517 unsigned int rgb = getPixel(x, y);
518 //unsigned int alpha = (rgb>>24) & 0xff;
519 unsigned int r = ((rgb>>16) & 0xff);
520 unsigned int g = ((rgb>> 8) & 0xff);
521 unsigned int b = ((rgb ) & 0xff);
522 fputc((unsigned char) r, f);
523 fputc((unsigned char) g, f);
524 fputc((unsigned char) b, f);
525 }
526 }
527 fclose(f);
528 return true;
529 }
532 #ifdef HAVE_GLIB
534 /**
535 * Create an image from a GdkPixbuf
536 */
537 SioxImage::SioxImage(GdkPixbuf *buf)
538 {
539 if (!buf)
540 return;
542 unsigned int width = gdk_pixbuf_get_width(buf);
543 unsigned int height = gdk_pixbuf_get_height(buf);
544 init(width, height); //DO THIS NOW!!
547 guchar *pixldata = gdk_pixbuf_get_pixels(buf);
548 int rowstride = gdk_pixbuf_get_rowstride(buf);
549 int n_channels = gdk_pixbuf_get_n_channels(buf);
551 //### Fill in the cells with RGB values
552 int row = 0;
553 for (unsigned int y=0 ; y<height ; y++)
554 {
555 guchar *p = pixldata + row;
556 for (unsigned int x=0 ; x<width ; x++)
557 {
558 int r = (int)p[0];
559 int g = (int)p[1];
560 int b = (int)p[2];
561 int alpha = (int)p[3];
563 setPixel(x, y, alpha, r, g, b);
564 p += n_channels;
565 }
566 row += rowstride;
567 }
569 }
572 /**
573 * Create a GdkPixbuf from this image
574 */
575 GdkPixbuf *SioxImage::getGdkPixbuf()
576 {
577 guchar *pixdata = (guchar *)
578 malloc(sizeof(guchar) * width * height * 3);
579 if (!pixdata)
580 return NULL;
582 int n_channels = 3;
583 int rowstride = width * 3;
585 GdkPixbuf *buf = gdk_pixbuf_new_from_data(pixdata, GDK_COLORSPACE_RGB,
586 0, 8, width, height,
587 rowstride, NULL, NULL);
589 //### Fill in the cells with RGB values
590 int row = 0;
591 for (unsigned int y=0 ; y < height ; y++)
592 {
593 guchar *p = pixdata + row;
594 for (unsigned x=0 ; x < width ; x++)
595 {
596 unsigned int rgb = getPixel(x, y);
597 p[0] = (rgb >> 16) & 0xff;
598 p[1] = (rgb >> 8) & 0xff;
599 p[2] = (rgb ) & 0xff;
600 p += n_channels;
601 }
602 row += rowstride;
603 }
605 return buf;
606 }
608 #endif /* GLIB */
613 //########################################################################
614 //# S I O X
615 //########################################################################
617 //##############
618 //## PUBLIC
619 //##############
621 /**
622 * Confidence corresponding to a certain foreground region (equals one).
623 */
624 const float Siox::CERTAIN_FOREGROUND_CONFIDENCE=1.0f;
626 /**
627 * Confidence for a region likely being foreground.
628 */
629 const float Siox::FOREGROUND_CONFIDENCE=0.8f;
631 /**
632 * Confidence for foreground or background type being equally likely.
633 */
634 const float Siox::UNKNOWN_REGION_CONFIDENCE=0.5f;
636 /**
637 * Confidence for a region likely being background.
638 */
639 const float Siox::BACKGROUND_CONFIDENCE=0.1f;
641 /**
642 * Confidence corresponding to a certain background reagion (equals zero).
643 */
644 const float Siox::CERTAIN_BACKGROUND_CONFIDENCE=0.0f;
646 /**
647 * Construct a Siox engine
648 */
649 Siox::Siox()
650 {
651 init();
652 }
655 /**
656 *
657 */
658 Siox::~Siox()
659 {
660 cleanup();
661 }
664 /**
665 * Error logging
666 */
667 void Siox::error(char *fmt, ...)
668 {
669 char msgbuf[256];
670 va_list args;
671 va_start(args, fmt);
672 vsnprintf(msgbuf, 255, fmt, args);
673 va_end(args) ;
674 #ifdef HAVE_GLIB
675 g_warning("Siox error: %s\n", msgbuf);
676 #else
677 fprintf(stderr, "Siox error: %s\n", msgbuf);
678 #endif
679 }
681 /**
682 * Trace logging
683 */
684 void Siox::trace(char *fmt, ...)
685 {
686 char msgbuf[256];
687 va_list args;
688 va_start(args, fmt);
689 vsnprintf(msgbuf, 255, fmt, args);
690 va_end(args) ;
691 #ifdef HAVE_GLIB
692 g_message("Siox: %s\n", msgbuf);
693 #else
694 fprintf(stdout, "Siox: %s\n", msgbuf);
695 #endif
696 }
701 /**
702 * Extract the foreground of the original image, according
703 * to the values in the confidence matrix.
704 */
705 SioxImage Siox::extractForeground(const SioxImage &originalImage,
706 unsigned int backgroundFillColor)
707 {
708 init();
710 SioxImage workImage = originalImage;
712 //# fetch some info from the image
713 width = workImage.getWidth();
714 height = workImage.getHeight();
715 pixelCount = width * height;
716 image = workImage.getImageData();
717 cm = workImage.getConfidenceData();
718 labelField = new int[pixelCount];
720 trace("### Creating signatures");
722 //#### create color signatures
723 std::vector<CLAB> knownBg;
724 std::vector<CLAB> knownFg;
725 for (int x = 0 ; x < workImage.getWidth() ; x++)
726 for (int y = 0 ; y < workImage.getHeight() ; y++)
727 {
728 float cm = workImage.getConfidence(x, y);
729 unsigned int pix = workImage.getPixel(x, y);
730 if (cm <= BACKGROUND_CONFIDENCE)
731 knownBg.push_back(pix); //note: uses CLAB(rgb)
732 else if (cm >= FOREGROUND_CONFIDENCE)
733 knownFg.push_back(pix);
734 }
736 trace("knownBg:%d knownFg:%d", knownBg.size(), knownFg.size());
739 std::vector<CLAB> bgSignature ;
740 if (!colorSignature(knownBg, bgSignature, 3))
741 {
742 error("Could not create background signature");
743 workImage.setValid(false);
744 return workImage;
745 }
746 std::vector<CLAB> fgSignature ;
747 if (!colorSignature(knownFg, fgSignature, 3))
748 {
749 error("Could not create foreground signature");
750 delete[] labelField;
751 workImage.setValid(false);
752 return workImage;
753 }
755 //trace("### bgSignature:%d", bgSignature.size());
757 if (bgSignature.size() < 1)
758 {
759 // segmentation impossible
760 error("Signature size is < 1. Segmentation is impossible");
761 delete[] labelField;
762 workImage.setValid(false);
763 return workImage;
764 }
766 // classify using color signatures,
767 // classification cached in hashmap for drb and speedup purposes
768 std::map<unsigned int, Tupel> hs;
770 for (unsigned int i=0; i<pixelCount; i++)
771 {
772 if (cm[i] >= FOREGROUND_CONFIDENCE)
773 {
774 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
775 }
776 else if (cm[i] <= BACKGROUND_CONFIDENCE)
777 {
778 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
779 }
780 else // somewhere in between
781 {
782 bool isBackground = true;
783 std::map<unsigned int, Tupel>::iterator iter = hs.find(i);
784 if (iter != hs.end()) //found
785 {
786 Tupel tupel = iter->second;
787 isBackground = tupel.minBgDist <= tupel.minFgDist;
788 }
789 else
790 {
791 CLAB lab(image[i]);
792 float minBg = sqrEuclidianDist(lab, bgSignature[0]);
793 int minIndex=0;
794 for (unsigned int j=1; j<bgSignature.size() ; j++)
795 {
796 float d = sqrEuclidianDist(lab, bgSignature[j]);
797 if (d<minBg)
798 {
799 minBg = d;
800 minIndex = j;
801 }
802 }
803 Tupel tupel(0.0f, 0, 0.0f, 0);
804 tupel.minBgDist = minBg;
805 tupel.indexMinBg = minIndex;
806 float minFg = 1.0e6f;
807 minIndex = -1;
808 for (unsigned int j = 0 ; j < fgSignature.size() ; j++)
809 {
810 float d = sqrEuclidianDist(lab, fgSignature[j]);
811 if (d < minFg)
812 {
813 minFg = d;
814 minIndex = j;
815 }
816 }
817 tupel.minFgDist = minFg;
818 tupel.indexMinFg = minIndex;
819 if (fgSignature.size() == 0)
820 {
821 isBackground=(minBg <= clusterSize);
822 // remove next line to force behaviour of old algorithm
823 //error("foreground signature does not exist");
824 //delete[] labelField;
825 //workImage.setValid(false);
826 //return workImage;
827 }
828 else
829 {
830 isBackground = minBg < minFg;
831 }
832 hs[image[i]] = tupel;
833 }
835 if (isBackground)
836 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
837 else
838 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
839 }
840 }
842 trace("### postProcessing");
844 //## postprocessing
845 smooth(cm, width, height, 0.33f, 0.33f, 0.33f); // average
846 normalizeMatrix(cm, pixelCount);
847 erode(cm, width, height);
848 keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.0/*sizeFactorToKeep*/);
850 for (int i=0; i < 2/*smoothness*/; i++)
851 smooth(cm, width, height, 0.33f, 0.33f, 0.33f); // average
853 normalizeMatrix(cm, pixelCount);
855 for (unsigned int i=0; i < pixelCount; i++)
856 {
857 if (cm[i] >= UNKNOWN_REGION_CONFIDENCE)
858 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
859 else
860 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
861 }
863 keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.5/*sizeFactorToKeep*/);
864 fillColorRegions();
865 dilate(cm, width, height);
867 delete[] labelField;
869 //#### Yaay. We are done. Now clear everything but the background
870 for (unsigned int y = 0 ; y < height ; y++)
871 for (unsigned int x = 0 ; x < width ; x++)
872 {
873 float conf = workImage.getConfidence(x, y);
874 if (conf < FOREGROUND_CONFIDENCE)
875 {
876 workImage.setPixel(x, y, backgroundFillColor);
877 }
878 }
880 trace("### Done");
881 return workImage;
882 }
886 //##############
887 //## PRIVATE
888 //##############
890 /**
891 * Initialize the Siox engine to its 'pristine' state.
892 * Performed at the beginning of extractForeground().
893 */
894 void Siox::init()
895 {
896 limits[0] = 0.64f;
897 limits[1] = 1.28f;
898 limits[2] = 2.56f;
900 float negLimits[3];
901 negLimits[0] = -limits[0];
902 negLimits[1] = -limits[1];
903 negLimits[2] = -limits[2];
905 clusterSize = sqrEuclidianDist(limits, 3, negLimits);
906 }
909 /**
910 * Clean up any debris from processing.
911 */
912 void Siox::cleanup()
913 {
914 }
919 /**
920 * Stage 1 of the color signature work. 'dims' will be either
921 * 2 for grays, or 3 for colors
922 */
923 void Siox::colorSignatureStage1(CLAB *points,
924 unsigned int leftBase,
925 unsigned int rightBase,
926 unsigned int recursionDepth,
927 unsigned int *clusterCount,
928 const unsigned int dims)
929 {
931 unsigned int currentDim = recursionDepth % dims;
932 CLAB point = points[leftBase];
933 float min = point(currentDim);
934 float max = min;
936 for (unsigned int i = leftBase + 1; i < rightBase ; i++)
937 {
938 point = points[i];
939 float curval = point(currentDim);
940 if (curval < min) min = curval;
941 if (curval > max) max = curval;
942 }
944 //Do the Rubner-rule split (sounds like a dance)
945 if (max - min > limits[currentDim])
946 {
947 float pivotPoint = (min + max) / 2.0; //average
948 unsigned int left = leftBase;
949 unsigned int right = rightBase - 1;
951 //# partition points according to the dimension
952 while (true)
953 {
954 while ( true )
955 {
956 point = points[left];
957 if (point(currentDim) > pivotPoint)
958 break;
959 left++;
960 }
961 while ( true )
962 {
963 point = points[right];
964 if (point(currentDim) <= pivotPoint)
965 break;
966 right--;
967 }
969 if (left > right)
970 break;
972 point = points[left];
973 points[left] = points[right];
974 points[right] = point;
976 left++;
977 right--;
978 }
980 //# Recurse and create sub-trees
981 colorSignatureStage1(points, leftBase, left,
982 recursionDepth + 1, clusterCount, dims);
983 colorSignatureStage1(points, left, rightBase,
984 recursionDepth + 1, clusterCount, dims);
985 }
986 else
987 {
988 //create a leaf
989 CLAB newpoint;
991 newpoint.C = rightBase - leftBase;
993 for (; leftBase < rightBase ; leftBase++)
994 {
995 newpoint.add(points[leftBase]);
996 }
998 //printf("clusters:%d\n", *clusters);
1000 if (newpoint.C != 0)
1001 newpoint.mul(1.0 / (float)newpoint.C);
1002 points[*clusterCount] = newpoint;
1003 (*clusterCount)++;
1004 }
1005 }
1009 /**
1010 * Stage 2 of the color signature work
1011 */
1012 void Siox::colorSignatureStage2(CLAB *points,
1013 unsigned int leftBase,
1014 unsigned int rightBase,
1015 unsigned int recursionDepth,
1016 unsigned int *clusterCount,
1017 const float threshold,
1018 const unsigned int dims)
1019 {
1022 unsigned int currentDim = recursionDepth % dims;
1023 CLAB point = points[leftBase];
1024 float min = point(currentDim);
1025 float max = min;
1027 for (unsigned int i = leftBase+ 1; i < rightBase; i++)
1028 {
1029 point = points[i];
1030 float curval = point(currentDim);
1031 if (curval < min) min = curval;
1032 if (curval > max) max = curval;
1033 }
1035 //Do the Rubner-rule split (sounds like a dance)
1036 if (max - min > limits[currentDim])
1037 {
1038 float pivotPoint = (min + max) / 2.0; //average
1039 unsigned int left = leftBase;
1040 unsigned int right = rightBase - 1;
1042 //# partition points according to the dimension
1043 while (true)
1044 {
1045 while ( true )
1046 {
1047 point = points[left];
1048 if (point(currentDim) > pivotPoint)
1049 break;
1050 left++;
1051 }
1052 while ( true )
1053 {
1054 point = points[right];
1055 if (point(currentDim) <= pivotPoint)
1056 break;
1057 right--;
1058 }
1060 if (left > right)
1061 break;
1063 point = points[left];
1064 points[left] = points[right];
1065 points[right] = point;
1067 left++;
1068 right--;
1069 }
1071 //# Recurse and create sub-trees
1072 colorSignatureStage2(points, leftBase, left,
1073 recursionDepth + 1, clusterCount, threshold, dims);
1074 colorSignatureStage2(points, left, rightBase,
1075 recursionDepth + 1, clusterCount, threshold, dims);
1076 }
1077 else
1078 {
1079 //### Create a leaf
1080 unsigned int sum = 0;
1081 for (unsigned int i = leftBase; i < rightBase; i++)
1082 sum += points[i].C;
1084 if ((float)sum >= threshold)
1085 {
1086 float scale = (float)(rightBase - leftBase);
1087 CLAB newpoint;
1089 for (; leftBase < rightBase; leftBase++)
1090 newpoint.add(points[leftBase]);
1092 if (scale != 0.0)
1093 newpoint.mul(1.0 / scale);
1094 points[*clusterCount] = newpoint;
1095 (*clusterCount)++;
1096 }
1097 }
1098 }
1102 /**
1103 * Main color signature method
1104 */
1105 bool Siox::colorSignature(const std::vector<CLAB> &inputVec,
1106 std::vector<CLAB> &result,
1107 const unsigned int dims)
1108 {
1110 unsigned int length = inputVec.size();
1112 if (length < 1) // no error. just don't do anything
1113 return true;
1115 CLAB *input = (CLAB *) malloc(length * sizeof(CLAB));
1117 if (!input)
1118 {
1119 error("Could not allocate buffer for signature");
1120 return false;
1121 }
1122 for (unsigned int i=0 ; i < length ; i++)
1123 input[i] = inputVec[i];
1125 unsigned int stage1length = 0;
1126 colorSignatureStage1(input, 0, length, 0,
1127 &stage1length, dims);
1129 unsigned int stage2length = 0;
1130 colorSignatureStage2(input, 0, stage1length, 0,
1131 &stage2length, length * 0.001, dims);
1133 result.clear();
1134 for (unsigned int i=0 ; i < stage2length ; i++)
1135 result.push_back(input[i]);
1137 free(input);
1139 return true;
1140 }
1144 /**
1145 *
1146 */
1147 void Siox::keepOnlyLargeComponents(float threshold,
1148 double sizeFactorToKeep)
1149 {
1150 for (unsigned long idx = 0 ; idx<pixelCount ; idx++)
1151 labelField[idx] = -1;
1153 int curlabel = 0;
1154 int maxregion= 0;
1155 int maxblob = 0;
1157 // slow but easy to understand:
1158 std::vector<int> labelSizes;
1159 for (unsigned long i=0 ; i<pixelCount ; i++)
1160 {
1161 int regionCount = 0;
1162 if (labelField[i] == -1 && cm[i] >= threshold)
1163 {
1164 regionCount = depthFirstSearch(i, threshold, curlabel++);
1165 labelSizes.push_back(regionCount);
1166 }
1168 if (regionCount>maxregion)
1169 {
1170 maxregion = regionCount;
1171 maxblob = curlabel-1;
1172 }
1173 }
1175 for (unsigned int i=0 ; i<pixelCount ; i++)
1176 {
1177 if (labelField[i] != -1)
1178 {
1179 // remove if the component is to small
1180 if (labelSizes[labelField[i]] * sizeFactorToKeep < maxregion)
1181 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
1183 // add maxblob always to foreground
1184 if (labelField[i] == maxblob)
1185 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
1186 }
1187 }
1189 }
1193 int Siox::depthFirstSearch(int startPos,
1194 float threshold, int curLabel)
1195 {
1196 // stores positions of labeled pixels, where the neighbours
1197 // should still be checked for processing:
1199 //trace("startPos:%d threshold:%f curLabel:%d",
1200 // startPos, threshold, curLabel);
1202 std::vector<int> pixelsToVisit;
1203 int componentSize = 0;
1205 if (labelField[startPos]==-1 && cm[startPos]>=threshold)
1206 {
1207 labelField[startPos] = curLabel;
1208 componentSize++;
1209 pixelsToVisit.push_back(startPos);
1210 }
1213 while (pixelsToVisit.size() > 0)
1214 {
1215 int pos = pixelsToVisit[pixelsToVisit.size() - 1];
1216 pixelsToVisit.erase(pixelsToVisit.end() - 1);
1217 unsigned int x = pos % width;
1218 unsigned int y = pos / width;
1220 // check all four neighbours
1221 int left = pos - 1;
1222 if (x-1>=0 && labelField[left]==-1 && cm[left]>=threshold)
1223 {
1224 labelField[left]=curLabel;
1225 componentSize++;
1226 pixelsToVisit.push_back(left);
1227 }
1229 int right = pos + 1;
1230 if (x+1 < width && labelField[right]==-1 && cm[right]>=threshold)
1231 {
1232 labelField[right]=curLabel;
1233 componentSize++;
1234 pixelsToVisit.push_back(right);
1235 }
1237 int top = pos - width;
1238 if (y-1>=0 && labelField[top]==-1 && cm[top]>=threshold)
1239 {
1240 labelField[top]=curLabel;
1241 componentSize++;
1242 pixelsToVisit.push_back(top);
1243 }
1245 int bottom = pos + width;
1246 if (y+1 < height && labelField[bottom]==-1
1247 && cm[bottom]>=threshold)
1248 {
1249 labelField[bottom]=curLabel;
1250 componentSize++;
1251 pixelsToVisit.push_back(bottom);
1252 }
1254 }
1255 return componentSize;
1256 }
1260 /**
1261 *
1262 */
1263 void Siox::fillColorRegions()
1264 {
1265 for (unsigned long idx = 0 ; idx<pixelCount ; idx++)
1266 labelField[idx] = -1;
1268 //int maxRegion=0; // unused now
1269 std::vector<int> pixelsToVisit;
1270 for (unsigned long i=0; i<pixelCount; i++)
1271 { // for all pixels
1272 if (labelField[i]!=-1 || cm[i]<UNKNOWN_REGION_CONFIDENCE)
1273 {
1274 continue; // already visited or bg
1275 }
1277 unsigned int origColor = image[i];
1278 unsigned long curLabel = i+1;
1279 labelField[i] = curLabel;
1280 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
1282 // int componentSize = 1;
1283 pixelsToVisit.push_back(i);
1284 // depth first search to fill region
1285 while (pixelsToVisit.size() > 0)
1286 {
1287 int pos = pixelsToVisit[pixelsToVisit.size() - 1];
1288 pixelsToVisit.erase(pixelsToVisit.end() - 1);
1289 unsigned int x=pos % width;
1290 unsigned int y=pos / width;
1291 // check all four neighbours
1292 int left = pos-1;
1293 if (x-1 >= 0 && labelField[left] == -1
1294 && CLAB::diff(image[left], origColor)<1.0)
1295 {
1296 labelField[left]=curLabel;
1297 cm[left]=CERTAIN_FOREGROUND_CONFIDENCE;
1298 // ++componentSize;
1299 pixelsToVisit.push_back(left);
1300 }
1301 int right = pos+1;
1302 if (x+1 < width && labelField[right]==-1
1303 && CLAB::diff(image[right], origColor)<1.0)
1304 {
1305 labelField[right]=curLabel;
1306 cm[right]=CERTAIN_FOREGROUND_CONFIDENCE;
1307 // ++componentSize;
1308 pixelsToVisit.push_back(right);
1309 }
1310 int top = pos - width;
1311 if (y-1>=0 && labelField[top]==-1
1312 && CLAB::diff(image[top], origColor)<1.0)
1313 {
1314 labelField[top]=curLabel;
1315 cm[top]=CERTAIN_FOREGROUND_CONFIDENCE;
1316 // ++componentSize;
1317 pixelsToVisit.push_back(top);
1318 }
1319 int bottom = pos + width;
1320 if (y+1 < height && labelField[bottom]==-1
1321 && CLAB::diff(image[bottom], origColor)<1.0)
1322 {
1323 labelField[bottom]=curLabel;
1324 cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE;
1325 // ++componentSize;
1326 pixelsToVisit.push_back(bottom);
1327 }
1328 }
1329 //if (componentSize>maxRegion) {
1330 // maxRegion=componentSize;
1331 //}
1332 }
1333 }
1338 /**
1339 * Applies the morphological dilate operator.
1340 *
1341 * Can be used to close small holes in the given confidence matrix.
1342 */
1343 void Siox::dilate(float *cm, int xres, int yres)
1344 {
1346 for (int y=0; y<yres; y++)
1347 {
1348 for (int x=0; x<xres-1; x++)
1349 {
1350 int idx=(y*xres)+x;
1351 if (cm[idx+1]>cm[idx])
1352 cm[idx]=cm[idx+1];
1353 }
1354 }
1356 for (int y=0; y<yres; y++)
1357 {
1358 for (int x=xres-1; x>=1; x--)
1359 {
1360 int idx=(y*xres)+x;
1361 if (cm[idx-1]>cm[idx])
1362 cm[idx]=cm[idx-1];
1363 }
1364 }
1366 for (int y=0; y<yres-1; y++)
1367 {
1368 for (int x=0; x<xres; x++)
1369 {
1370 int idx=(y*xres)+x;
1371 if (cm[((y+1)*xres)+x] > cm[idx])
1372 cm[idx]=cm[((y+1)*xres)+x];
1373 }
1374 }
1376 for (int y=yres-1; y>=1; y--)
1377 {
1378 for (int x=0; x<xres; x++)
1379 {
1380 int idx=(y*xres)+x;
1381 if (cm[((y-1)*xres)+x] > cm[idx])
1382 cm[idx]=cm[((y-1)*xres)+x];
1383 }
1384 }
1385 }
1387 /**
1388 * Applies the morphological erode operator.
1389 */
1390 void Siox::erode(float *cm, int xres, int yres)
1391 {
1392 for (int y=0; y<yres; y++)
1393 {
1394 for (int x=0; x<xres-1; x++)
1395 {
1396 int idx=(y*xres)+x;
1397 if (cm[idx+1] < cm[idx])
1398 cm[idx]=cm[idx+1];
1399 }
1400 }
1401 for (int y=0; y<yres; y++)
1402 {
1403 for (int x=xres-1; x>=1; x--)
1404 {
1405 int idx=(y*xres)+x;
1406 if (cm[idx-1] < cm[idx])
1407 cm[idx]=cm[idx-1];
1408 }
1409 }
1410 for (int y=0; y<yres-1; y++)
1411 {
1412 for (int x=0; x<xres; x++)
1413 {
1414 int idx=(y*xres)+x;
1415 if (cm[((y+1)*xres)+x] < cm[idx])
1416 cm[idx]=cm[((y+1)*xres)+x];
1417 }
1418 }
1419 for (int y=yres-1; y>=1; y--)
1420 {
1421 for (int x=0; x<xres; x++)
1422 {
1423 int idx=(y*xres)+x;
1424 if (cm[((y-1)*xres)+x] < cm[idx])
1425 cm[idx]=cm[((y-1)*xres)+x];
1426 }
1427 }
1428 }
1433 /**
1434 * Normalizes the matrix to values to [0..1].
1435 */
1436 void Siox::normalizeMatrix(float *cm, int cmSize)
1437 {
1438 float max= -1000000.0f;
1439 for (int i=0; i<cmSize; i++)
1440 if (max<cm[i] > max)
1441 max=cm[i];
1443 if (max<=0.0 || max==1.0)
1444 return;
1446 float alpha=1.00f/max;
1447 premultiplyMatrix(alpha, cm, cmSize);
1448 }
1450 /**
1451 * Multiplies matrix with the given scalar.
1452 */
1453 void Siox::premultiplyMatrix(float alpha, float *cm, int cmSize)
1454 {
1455 for (int i=0; i<cmSize; i++)
1456 cm[i]=alpha*cm[i];
1457 }
1459 /**
1460 * Blurs confidence matrix with a given symmetrically weighted kernel.
1461 *
1462 * In the standard case confidence matrix entries are between 0...1 and
1463 * the weight factors sum up to 1.
1464 */
1465 void Siox::smooth(float *cm, int xres, int yres,
1466 float f1, float f2, float f3)
1467 {
1468 for (int y=0; y<yres; y++)
1469 {
1470 for (int x=0; x<xres-2; x++)
1471 {
1472 int idx=(y*xres)+x;
1473 cm[idx]=f1*cm[idx]+f2*cm[idx+1]+f3*cm[idx+2];
1474 }
1475 }
1476 for (int y=0; y<yres; y++)
1477 {
1478 for (int x=xres-1; x>=2; x--)
1479 {
1480 int idx=(y*xres)+x;
1481 cm[idx]=f3*cm[idx-2]+f2*cm[idx-1]+f1*cm[idx];
1482 }
1483 }
1484 for (int y=0; y<yres-2; y++)
1485 {
1486 for (int x=0; x<xres; x++)
1487 {
1488 int idx=(y*xres)+x;
1489 cm[idx]=f1*cm[idx]+f2*cm[((y+1)*xres)+x]+f3*cm[((y+2)*xres)+x];
1490 }
1491 }
1492 for (int y=yres-1; y>=2; y--)
1493 {
1494 for (int x=0; x<xres; x++)
1495 {
1496 int idx=(y*xres)+x;
1497 cm[idx]=f3*cm[((y-2)*xres)+x]+f2*cm[((y-1)*xres)+x]+f1*cm[idx];
1498 }
1499 }
1500 }
1502 /**
1503 * Squared Euclidian distance of p and q.
1504 */
1505 float Siox::sqrEuclidianDist(float *p, int pSize, float *q)
1506 {
1507 float sum=0.0;
1508 for (int i=0; i<pSize; i++)
1509 {
1510 float v = p[i] - q[i];
1511 sum += v*v;
1512 }
1513 return sum;
1514 }
1516 /**
1517 * Squared Euclidian distance of p and q.
1518 */
1519 float Siox::sqrEuclidianDist(const CLAB &p, const CLAB &q)
1520 {
1521 float sum=0;
1522 sum += (p.L - q.L) * (p.L - q.L);
1523 sum += (p.A - q.A) * (p.A - q.A);
1524 sum += (p.B - q.B) * (p.B - q.B);
1525 return sum;
1526 }
1535 } // namespace siox
1536 } // namespace org
1538 //########################################################################
1539 //# E N D O F F I L E
1540 //########################################################################