00001 #include "sadie.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 void
00013 HISTEQ (IMAGE * in,
00014 short nlev,
00015
00016 short inc,
00017 IMAGE ** out
00018
00019
00020 )
00021 {
00022 register short i, j, k;
00023 char msg[SLEN];
00024 long *hist = 0;
00025 double *dist = 0, factor, width;
00026
00027 if (TIMES)
00028 TIMING (T_HISTEQ);
00029 if (NAMES)
00030 {
00031 MESSAGE ('I', "");
00032 MESSAGE ('I', "HISTEQ");
00033 MESSAGE ('I', "");
00034 sprintf (msg, " Input image: %s", in->text);
00035 MESSAGE ('I', msg);
00036 sprintf (msg, " Histogram resolution: %d", nlev);
00037 MESSAGE ('I', msg);
00038 sprintf (msg, " Line and pixel subsample increment: %d", inc);
00039 MESSAGE ('I', msg);
00040 MESSAGE ('I', " ...............");
00041 }
00042
00043
00044 if (!CHECKIMG (in))
00045 {
00046 MESSAGE ('E', " Can't identify image.");
00047 goto the_end;
00048 }
00049 else if (inc <= 0)
00050 {
00051 MESSAGE ('E', " Subsample increment must be greater than zero.");
00052 goto the_end;
00053 }
00054 else if (nlev <= 0)
00055 {
00056 MESSAGE ('E', " Number of histogram bins must be greater than zero.");
00057 goto the_end;
00058 }
00059
00060
00061 if (!CHECKIMG (*out))
00062 GETMEM (in->nbnd, in->nlin, in->npix, out);
00063 if (!*out)
00064 goto the_end;
00065
00066
00067 if (!(hist = (long *) malloc (in->nbnd * nlev * sizeof (long))))
00068 {
00069 MESSAGE ('E', " Memory request failed.");
00070 goto the_end;
00071 }
00072 if (!(dist = (double *) malloc (in->nbnd * nlev * sizeof (double))))
00073 {
00074 MESSAGE ('E', " Memory request failed.");
00075 goto the_end;
00076 }
00077
00078
00079 RANGE (in);
00080 HISTOGRM (in, inc, in->gmin, in->gmax, nlev, hist);
00081 for (i = 0; i < in->nbnd; i++)
00082 {
00083 for (j = 1; j < nlev; j++)
00084 {
00085 hist[i * nlev + j] += hist[i * nlev + j - 1];
00086 }
00087 }
00088 factor = (double) (in->gmax - in->gmin) / (double) hist[nlev - 1];
00089 for (i = 0; i < in->nbnd; i++)
00090 {
00091 for (j = 0; j < nlev; j++)
00092 {
00093 dist[i * nlev + j] = hist[i * nlev + j] * factor;
00094 }
00095 }
00096 width = (double) (in->gmax - in->gmin) / (double) (nlev - 1);
00097
00098
00099 for (i = 0; i < in->nbnd; i++)
00100 {
00101 for (j = 0; j < in->nlin; j++)
00102 {
00103 for (k = 0; k < in->npix; k++)
00104 {
00105 (*out)->data[i][j][k] =
00106 (PIXEL) (dist
00107 [i * nlev +
00108 (short) ((double) (in->data[i][j][k] - in->gmin) /
00109 width)] + (double) in->gmin);
00110 }
00111 }
00112 }
00113 (*out)->gmin = (PIXEL) ((double) in->gmin + dist[0]);
00114 (*out)->gmax = (PIXEL) ((double) in->gmin + dist[nlev - 1]);
00115
00116 the_end:
00117 if (hist)
00118 free (hist);
00119 if (dist)
00120 free (dist);
00121 if (TIMES)
00122 TIMING (T_EXIT);
00123 }