summaryrefslogtreecommitdiff
path: root/src/process/im_statistics.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/process/im_statistics.cpp')
-rw-r--r--src/process/im_statistics.cpp341
1 files changed, 341 insertions, 0 deletions
diff --git a/src/process/im_statistics.cpp b/src/process/im_statistics.cpp
new file mode 100644
index 0000000..b9f086d
--- /dev/null
+++ b/src/process/im_statistics.cpp
@@ -0,0 +1,341 @@
+/** \file
+ * \brief Image Statistics Calculations
+ *
+ * See Copyright Notice in im_lib.h
+ * $Id: im_statistics.cpp,v 1.1 2008/10/17 06:16:33 scuri Exp $
+ */
+
+
+#include <im.h>
+#include <im_util.h>
+#include <im_math_op.h>
+
+#include "im_process_ana.h"
+
+#include <stdlib.h>
+#include <memory.h>
+#include <math.h>
+
+static unsigned long count_map(const imImage* image)
+{
+ unsigned long histo[256];
+ int size = image->width * image->height;
+ imCalcHistogram((imbyte*)image->data[0], size, histo, 0);
+ unsigned long numcolor = 0;
+
+ for (int i = 0; i < 256; i++)
+ {
+ if(histo[i] != 0)
+ numcolor++;
+ }
+
+ return numcolor;
+}
+
+// will count also all the 3 components color spaces
+static unsigned long count_rgb(const imImage* image)
+{
+ imbyte *count = (imbyte*)calloc(sizeof(imbyte), 1 << 21 ); /* (2^24)/8=2^21 ~ 2Mb */
+ if (!count)
+ return (unsigned long)-1;
+
+ int size = image->width * image->height;
+ imbyte *red = (imbyte*)image->data[0];
+ imbyte *green = (imbyte*)image->data[1];
+ imbyte *blue = (imbyte*)image->data[2];
+
+ int index;
+ unsigned long numcolor = 0;
+
+ for(int i = 0; i < size; i++)
+ {
+ index = red[i] << 16 | green[i] << 8 | blue[i];
+
+ if(imDataBitGet(count, index) == 0)
+ numcolor++;
+
+ imDataBitSet(count, index, 1);
+ }
+
+ free(count);
+
+ return numcolor;
+}
+
+unsigned long imCalcCountColors(const imImage* image)
+{
+ if (imColorModeDepth(image->color_space) > 1)
+ return count_rgb(image);
+ else
+ return count_map(image);
+}
+
+void imCalcHistogram(const imbyte* map, int size, unsigned long* histo, int cumulative)
+{
+ int i;
+
+ memset(histo, 0, 256 * sizeof(unsigned long));
+
+ for (i = 0; i < size; i++)
+ histo[*map++]++;
+
+ if (cumulative)
+ {
+ /* make cumulative histogram */
+ for (i = 1; i < 256; i++)
+ histo[i] += histo[i-1];
+ }
+}
+
+void imCalcUShortHistogram(const imushort* map, int size, unsigned long* histo, int cumulative)
+{
+ int i;
+
+ memset(histo, 0, 65535 * sizeof(unsigned long));
+
+ for (i = 0; i < size; i++)
+ histo[*map++]++;
+
+ if (cumulative)
+ {
+ /* make cumulative histogram */
+ for (i = 1; i < 65535; i++)
+ histo[i] += histo[i-1];
+ }
+}
+
+void imCalcGrayHistogram(const imImage* image, unsigned long* histo, int cumulative)
+{
+ int i;
+
+ memset(histo, 0, 256 * sizeof(unsigned long));
+
+ if (image->color_space == IM_GRAY)
+ {
+ imbyte* map = (imbyte*)image->data[0];
+ for (i = 0; i < image->count; i++)
+ histo[*map++]++;
+ }
+ else if (image->color_space == IM_MAP || image->color_space == IM_BINARY)
+ {
+ imbyte* map = (imbyte*)image->data[0];
+ imbyte gray_map[256], r, g, b;
+
+ for (i = 0; i < image->palette_count; i++)
+ {
+ imColorDecode(&r, &g, &b, image->palette[i]);
+ gray_map[i] = (imbyte)((299*r + 587*g + 114*b) / 1000);
+ }
+
+ for (i = 0; i < image->count; i++)
+ {
+ int index = *map++;
+ histo[gray_map[index]]++;
+ }
+ }
+ else
+ {
+ imbyte gray;
+ imbyte* r = (imbyte*)image->data[0];
+ imbyte* g = (imbyte*)image->data[1];
+ imbyte* b = (imbyte*)image->data[2];
+ for (i = 0; i < image->count; i++)
+ {
+ gray = (imbyte)((299*(*r++) + 587*(*g++) + 114*(*b++)) / 1000);
+ histo[gray]++;
+ }
+ }
+
+ if (cumulative)
+ {
+ /* make cumulative histogram */
+ for (i = 1; i < 256; i++)
+ histo[i] += histo[i-1];
+ }
+}
+
+template <class T>
+static void DoStats(T* data, int count, imStats* stats)
+{
+ memset(stats, 0, sizeof(imStats));
+
+ stats->min = (float)data[0];
+ stats->max = (float)data[0];
+
+ for (int i = 0; i < count; i++)
+ {
+ if (data[i] < stats->min)
+ stats->min = (float)data[i];
+
+ if (data[i] > stats->max)
+ stats->max = (float)data[i];
+
+ if (data[i] > 0)
+ stats->positive++;
+
+ if (data[i] < 0)
+ stats->negative++;
+
+ if (data[i] == 0)
+ stats->zeros++;
+
+ stats->mean += (float)data[i];
+ stats->stddev += ((float)data[i])*((float)data[i]);
+ }
+
+ stats->mean /= float(count);
+ stats->stddev = (float)sqrt((stats->stddev - count * stats->mean*stats->mean)/(count-1.0));
+}
+
+void imCalcImageStatistics(const imImage* image, imStats* stats)
+{
+ int count = image->width * image->height;
+
+ for (int i = 0; i < image->depth; i++)
+ {
+ switch(image->data_type)
+ {
+ case IM_BYTE:
+ DoStats((imbyte*)image->data[i], count, &stats[i]);
+ break;
+ case IM_USHORT:
+ DoStats((imushort*)image->data[i], count, &stats[i]);
+ break;
+ case IM_INT:
+ DoStats((int*)image->data[i], count, &stats[i]);
+ break;
+ case IM_FLOAT:
+ DoStats((float*)image->data[i], count, &stats[i]);
+ break;
+ }
+ }
+}
+
+void imCalcHistogramStatistics(const imImage* image, imStats* stats)
+{
+ int image_size = image->width * image->height;
+ unsigned long histo[256];
+
+ for (int d = 0; d < image->depth; d++)
+ {
+ imCalcHistogram((imbyte*)image->data[d], image_size, histo, 0);
+ DoStats((unsigned long*)histo, 256, &stats[d]);
+ }
+}
+
+void imCalcHistoImageStatistics(const imImage* image, int* median, int* mode)
+{
+ unsigned long histo[256];
+
+ for (int d = 0; d < image->depth; d++)
+ {
+ int i;
+ imCalcHistogram((imbyte*)image->data[d], image->count, histo, 0);
+
+ unsigned long half = image->count/2;
+ unsigned long count = histo[0];
+ for (i = 1; i < 256; i++)
+ {
+ if (count > half)
+ {
+ median[d] = i-1;
+ break;
+ }
+
+ count += histo[i];
+ }
+
+ unsigned long max = histo[0];
+ for (i = 1; i < 256; i++)
+ {
+ if (max < histo[i])
+ max = histo[i];
+ }
+
+ int found_mode = 0;
+ for (i = 0; i < 256; i++)
+ {
+ if (histo[i] == max)
+ {
+ if (found_mode)
+ {
+ mode[d] = -1;
+ break;
+ }
+
+ mode[d] = i;
+ found_mode = 1;
+ }
+ }
+ }
+}
+
+float imCalcSNR(const imImage* image, const imImage* noise_image)
+{
+ imStats stats[3];
+ imCalcImageStatistics((imImage*)image, stats);
+
+ imStats noise_stats[3];
+ imCalcImageStatistics((imImage*)noise_image, noise_stats);
+
+ if (image->color_space == IM_RGB)
+ {
+ noise_stats[0].stddev += noise_stats[1].stddev;
+ noise_stats[0].stddev += noise_stats[2].stddev;
+ noise_stats[0].stddev /= 3;
+ stats[0].stddev += stats[1].stddev;
+ stats[0].stddev += stats[2].stddev;
+ stats[0].stddev /= 3;
+ }
+
+ if (noise_stats[0].stddev == 0)
+ return 0;
+
+ return float(20.*log10(stats[0].stddev / noise_stats[0].stddev));
+}
+
+template <class T>
+static float DoRMSOp(T *map1, T *map2, int count)
+{
+ float rmserror = 0.0f;
+ float diff;
+
+ for (int i = 0; i < count; i++)
+ {
+ diff = float(map1[i] - map2[i]);
+ rmserror += diff * diff;
+ }
+
+ return rmserror;
+}
+
+float imCalcRMSError(const imImage* image1, const imImage* image2)
+{
+ float rmserror = 0.0f;
+
+ int count = image1->count*image1->depth;
+
+ switch(image1->data_type)
+ {
+ case IM_BYTE:
+ rmserror = DoRMSOp((imbyte*)image1->data[0], (imbyte*)image2->data[0], count);
+ break;
+ case IM_USHORT:
+ rmserror = DoRMSOp((imushort*)image1->data[0], (imushort*)image2->data[0], count);
+ break;
+ case IM_INT:
+ rmserror = DoRMSOp((int*)image1->data[0], (int*)image2->data[0], count);
+ break;
+ case IM_FLOAT:
+ rmserror = DoRMSOp((float*)image1->data[0], (float*)image2->data[0], count);
+ break;
+ case IM_CFLOAT:
+ rmserror = DoRMSOp((float*)image1->data[0], (float*)image2->data[0], 2*count);
+ break;
+ }
+
+ rmserror = float(sqrt(rmserror / float((count * image1->depth))));
+
+ return rmserror;
+}
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-21 22:30:19 +0000