00001
00002
00003 #include <queue>
00004 #include <stdlib.h>
00005
00006 #include "Watershed.h"
00007
00008 static const int MAXGREYVAL = 256;
00009
00010 static void sortPixels(ByteImage &im, int *sorted_pixels, int *hist)
00011 {
00012 int i, s, *index[MAXGREYVAL];
00013 byte *tmp;
00014
00015 for (i=0, tmp=im[0]; i<im.getSize(); i++, tmp++)
00016 hist[*tmp]++;
00017
00018 s=0;
00019 for (i=0;i<MAXGREYVAL; i++)
00020 {
00021 index[i] = sorted_pixels+s;
00022 s+=hist[i];
00023 }
00024
00025 for (i=0, tmp=im[0]; i<im.getSize(); i++, tmp++)
00026 *(index[*tmp]++) = i;
00027 }
00028
00029 #define MASK -2
00030 #define WSHED 0
00031 #define INIT -1
00032 #define RANGE(i, j, h, w) ( (i>=0) && (j>=0) && (i<h) && (j<w) )
00033
00034 void Watershed(ByteImage &im, CONN_TYPE connect, IntImage &out)
00035 {
00036 int h, x, y, gval=0, pixel, cur_dist, cur_label=0, imsize = im.getSize(), *dist_im=0, *SortPixels=0, *hist=0;
00037 int *current, *out_im, width=im.getWidth(), height=im.getHeight();
00038 queue<int> q;
00039
00040 try
00041 {
00042 dist_im = (int *)calloc(imsize, sizeof(int));
00043 Assert<Fatal_error>(dist_im);
00044 SortPixels = (int *)malloc(imsize * sizeof(int));
00045 Assert<Fatal_error>(SortPixels);
00046 hist = (int *)calloc(MAXGREYVAL, sizeof(int));
00047 Assert<Fatal_error>(hist);
00048 }
00049 catch(Fatal_error e) { e.print("Unable to allocate memory"); }
00050
00051 out.makeImage(width, height);
00052 out.clearImage(INIT);
00053 out_im = out[0];
00054
00055 sortPixels(im, SortPixels, hist);
00056
00057 for(h=0;h<MAXGREYVAL;h++) {
00058 for(current=SortPixels+gval; current<SortPixels+gval+hist[h]; current++) {
00059 pixel = *current;
00060 out_im[pixel]=MASK;
00061 x=pixel % width;
00062 y=pixel / width;
00063 for(int k=-1;k<=1;k++)
00064 for(int m=-1;m<=1;m++) {
00065 if(RANGE(y+k, x+m, height, width))
00066 if(((connect==FOUR) && (((k==0) && (m!=0)) || ((k!=0) && (m==0)))) ||
00067 ((connect==EIGHT) && ((k!=0) || (m!=0)) )) {
00068 if((out[y+k][x+m]>0) || (out[y+k][x+m]==WSHED)) {
00069 q.push(pixel);
00070 dist_im[pixel]=1;
00071 k=1; break;
00072 }
00073 }
00074 }
00075 }
00076
00077 cur_dist=1;
00078 q.push(-1);
00079
00080 #define POP { pixel=q.front(); q.pop(); }
00081
00082 for(;;) {
00083 POP
00084 if(pixel==-1) {
00085 if(q.empty()) break;
00086 else {
00087 q.push(-1);
00088 cur_dist++;
00089 POP
00090 }
00091 }
00092 x=pixel % width;
00093 y=pixel / width;
00094 for(int k=-1;k<=1;k++)
00095 for(int m=-1;m<=1;m++) {
00096 if(RANGE(y+k, x+m, height, width))
00097 if(((connect==FOUR) && (((k==0) && (m!=0)) || ((k!=0) && (m==0)))) ||
00098 ((connect==EIGHT) && ((k!=0) || (m!=0)) )) {
00099 int neigh = (y+k)*width + (x+m);
00100 if((dist_im[neigh]<cur_dist) && ((out_im[neigh]>0) || (out_im[neigh]==WSHED)) ) {
00101 if(out_im[neigh]>0) {
00102 if((out_im[pixel]==MASK) || (out_im[pixel]==WSHED))
00103 out_im[pixel]=out_im[neigh];
00104 else if(out_im[pixel]!=out_im[neigh])
00105 out_im[pixel]=WSHED;
00106 }
00107 else if(out_im[pixel]==MASK) out_im[pixel]=WSHED;
00108 }
00109 else if((out_im[neigh]==MASK) && (dist_im[neigh]==0)) {
00110 dist_im[neigh]=cur_dist+1;
00111 q.push(neigh);
00112 }
00113 }
00114 }
00115 }
00116
00117 for(current=SortPixels+gval; current<SortPixels+gval+hist[h]; current++) {
00118 pixel = *current;
00119 dist_im[pixel]=0;
00120 if(out_im[pixel]==MASK) {
00121 cur_label++;
00122 q.push(pixel);
00123 out_im[pixel]=cur_label;
00124 while(!q.empty()) {
00125 POP
00126 x=pixel % width;
00127 y=pixel / width;
00128 for(int k=-1;k<=1;k++)
00129 for(int m=-1;m<=1;m++) {
00130 if(RANGE(y+k, x+m, height, width))
00131 if(((connect==FOUR) && (((k==0) && (m!=0)) || ((k!=0) && (m==0)))) ||
00132 ((connect==EIGHT) && ((k!=0) || (m!=0)) )) {
00133 int neigh = (y+k)*width + (x+m);
00134 if(out_im[neigh]==MASK) {
00135 q.push(neigh);
00136 out_im[neigh]=cur_label;
00137 }
00138 }
00139 }
00140 }
00141 }
00142 }
00143 gval+=hist[h];
00144 }
00145
00146 for(y=0;y<height;y++)
00147 for(x=0;x<width;x++) {
00148 if(out[y][x]>0) {
00149 for(int k=-1;k<=1;k++)
00150 for(int m=-1;m<=1;m++) {
00151 if(RANGE(y+k, x+m, height, width))
00152 if(((connect==FOUR) && (((k==0) && (m!=0)) || ((k!=0) && (m==0)))) ||
00153 ((connect==EIGHT) && ((k!=0) || (m!=0)) )) {
00154 if((out[y][x]>out[y+k][x+m]) && (out[y+k][x+m]>0)) {
00155 out[y][x]=WSHED;
00156 k=1; break;
00157 }
00158 }
00159 }
00160 }
00161 }
00162
00163 free(dist_im);
00164 free(SortPixels);
00165 free(hist);
00166 }
00167
00168
00169
00170
00171