Main Page   Modules   Class Hierarchy   Compound List   File List   Compound Members   Related Pages   Examples  

Watershed.cpp

00001 /* Copyright (c) 2001 A. Jalba */
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     }  // endfor(;;)                   
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   }  //endfor h                 
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