How to label a pixel of color image in region growing ?

I have two image one of which is RGB color image and other is binary image containing seed pixels. I want to label(color) neighboring pixel by the seed pixel label(color) , if the euclidean distance between seed pixel and neighboring pixel less than by Threshold value (T) else by white color. Temp is used to store label(color) value of seed pixel....

   clear all
  img=imread('ldy.jpg');
  figure, imshow(img), title('Input Image');
  imggray=rgb2gray(img);
  bin=imread('seedbinary.jpg');
  figure, imshow(bin);
  T = graythresh(imggray);
  [rows, columns, numberOfColorBands] = size(img);
  temp=zeros(rows, columns, numberOfColorBands);
  for i=1:rows;
      for j=1:columns;
          if bin(i, j)==1;
              temp(i, j, 1)=img(i, j, 1);
              temp(i, j, 2)=img(i, j, 2);
              temp(i, j, 3)=img(i, j, 3);
          end
      end
  end
  figure, imshow(temp), title('Seed Pixel Label');
  for i=2:rows-1;
      for j=2:columns-1;
          if bin(i, j)==1;
              d1=sqrt(double((img(i, j, 1)-img(i, j+1, 1))^2+(img(i, j, 2)-img(i, j+1, 2))^2+(img(i, j, 3)-img(i, j+1, 3))^2));
              if d1<=T
              img(i, j+1, 1)=temp(i, j, 1);
              img(i, j+1, 2)=temp(i, j, 2);
              img(i, j+1, 3)=temp(i, j, 3);
              else
              img(i, j+1, 1)=255;
              img(i, j+1, 2)=255;
              img(i, j+1, 3)=255;    
              end
          end
      end
  end
  figure, imshow(img), title('Output Image');
There is no error generated but not providing desired result.... Your effort will be appreciated ... Tell me where I am doing mistake....

Answer 1

[David Young]

I suspect that the value of T is one problem. T is based on the image gray-levels, but it is used to threshold squared gray-level differences. T is probably not a good threshold for differences.

I suggest removing the call to graythresh, and simply setting T to a constant instead. Experiment with the value of the constant until you get a result that looks reasonable. After that worry about how to set T automatically if you need to.

If that doesn't help, please explain what is actually wrong with the results - in what way are they not what you expect?

[Rajesh Gothwal]

This is not an issue. I want to subdivide an color image into different regions separated by white line, But The above code provide output image with some white pixels. Region color should be the color of seed pixel.

[David Young]

It's quite difficult to relate your code to your description of the desired output. Please could you explain why you only take differences in the row direction (between positions i,j and i,j+1) and not in the column direction? And are you aiming to have each region, between the white lines, all the same colour?

[Rajesh Gothwal]

This is the code actually i am using. There are two problem with output image 1) Regions are not growing after one iteration and 2)Regions are not labeled by seed pixel color. No, I do not want same color. I want to color each region by the color of that region's seed pixel.

clear all
img=imread('ldy.jpg');
figure, imshow(img), title('Input Image');
imggray=rgb2gray(img);
bin=imread('seedbinary.jpg');
figure, imshow(bin), title('Seed Pixel in Binary');
T = graythresh(imggray);
[rows, columns, numberOfColorBands] = size(img);
temp=zeros(rows, columns, numberOfColorBands);

for i=1:rows;
    for j=1:columns;
        if bin(i, j)==1

              temp(i, j, 1)=img(i, j, 1);
              temp(i, j, 2)=img(i, j, 2);
              temp(i, j, 3)=img(i, j, 3);

          end
      end
  end
  figure, imshow(temp), title('Seed Pixel Label');
  for k=1:1000;
  for i=2:rows-1;
      for j=2:columns-1;
          if bin(i, j)==1;
              d1=sqrt(double((img(i, j, 1)-img(i, j+1, 1))^2+(img(i, j, 2)-img(i, j+1, 2))^2+(img(i, j, 3)-img(i, j+1, 3))^2));
              if d1<=T
              img(i, j+1, 1)=temp(i, j, 1);
              img(i, j+1, 2)=temp(i, j, 2);
              img(i, j+1, 3)=temp(i, j, 3);
              else
              img(i, j+1, 1)=255;
              img(i, j+1, 2)=255;
              img(i, j+1, 3)=255;    
              end

              d2=sqrt(double((img(i, j, 1)-img(i-1, j+1, 1))^2+(img(i, j, 2)-img(i-1, j+1, 2))^2+(img(i, j, 3)-img(i-1, j+1, 3))^2));
              if d2<=T
              img(i-1, j+1, 1)=temp(i, j, 1);
              img(i-1, j+1, 2)=temp(i, j, 2);
              img(i-1, j+1, 3)=temp(i, j, 3);
              else
              img(i-1, j+1, 1)=255;
              img(i-1, j+1, 2)=255;
              img(i-1, j+1, 3)=255;
              end

              d3=sqrt(double((img(i, j, 1)-img(i-1, j, 1))^2+(img(i, j, 2)-img(i-1, j, 2))^2+(img(i, j, 3)-img(i-1, j, 3))^2));
              if d3<=T
              img(i-1, j, 1)=temp(i, j, 1);
              img(i-1, j, 2)=temp(i, j, 2);
              img(i-1, j, 3)=temp(i, j, 3);
              else
              img(i-1, j, 1)=255;
              img(i-1, j, 2)=255;
              img(i-1, j, 3)=255;
              end

              d4=sqrt(double((img(i, j, 1)-img(i-1, j-1, 1))^2+(img(i, j, 2)-img(i-1, j-1, 2))^2+(img(i, j, 3)-img(i-1, j-1, 3))^2));
              if d4<=T
              img(i-1, j-1, 1)=temp(i, j, 1);
              img(i-1, j-1, 2)=temp(i, j, 2);
              img(i-1, j-1, 3)=temp(i, j, 3);
              else
              img(i-1, j-1, 1)=255;
              img(i-1, j-1, 2)=255;
              img(i-1, j-1, 3)=255;
              end

              d5=sqrt(double((img(i, j, 1)-img(i, j-1, 1))^2+(img(i, j, 2)-img(i, j-1, 2))^2+(img(i, j, 3)-img(i, j-1, 3))^2));
              if d5<=T
              img(i, j-1, 1)=temp(i, j, 1);
              img(i, j-1, 2)=temp(i, j, 2);
              img(i, j-1, 3)=temp(i, j, 3);
              else
              img(i, j-1, 1)=255;
              img(i, j-1, 2)=255;
              img(i, j-1, 3)=255;
              end

              d6=sqrt(double((img(i, j, 1)-img(i+1, j-1, 1))^2+(img(i, j, 2)-img(i+1, j-1, 2))^2+(img(i, j, 3)-img(i+1, j-1, 3))^2));
              if d6<=T
              img(i+1, j-1, 1)=temp(i, j, 1);
              img(i+1, j-1, 2)=temp(i, j, 2);
              img(i+1, j-1, 3)=temp(i, j, 3);
              else
              img(i+1, j-1, 1)=255;
              img(i+1, j-1, 2)=255;
              img(i+1, j-1, 3)=255;
              end

              d7=sqrt(double((img(i, j, 1)-img(i+1, j, 1))^2+(img(i, j, 2)-img(i+1, j, 2))^2+(img(i, j, 3)-img(i+1, j, 3))^2));
              if d7<=T
              img(i+1, j, 1)=temp(i, j, 1);
              img(i+1, j, 2)=temp(i, j, 2);
              img(i+1, j, 3)=temp(i, j, 3);
              else
              img(i+1, j, 1)=255;
              img(i+1, j, 2)=255;
              img(i+1, j, 3)=255;
              end

              d8=sqrt(double((img(i, j, 1)-img(i+1, j+1, 1))^2+(img(i, j, 2)-img(i+1, j+1, 2))+(img(i, j, 3)-img(i+1, j+1, 3))^2));
              if d8<=T
              img(i+1, j+1, 1)=temp(i, j, 1);
              img(i+1, j+1, 2)=temp(i, j, 2);
              img(i+1, j+1, 3)=temp(i, j, 3);
              else
              img(i+1, j+1, 1)=255;
              img(i+1, j+1, 2)=255;
              img(i+1, j+1, 3)=255;
              end
          end
      end
  end
  end
  figure, imshow(img), title('Output Image');
  % imwrite(img,'output.jpg','jpg')

[Image Analyst]

I didn't go over all of it but one thing that struck me immediately was clipping performed by this:

double((img(i, j, 1)-img(i, j+1, 1)

You can't subtract and then make it double. After subtraction the clipping has already happened so casting to double does not help at all. You need to cast each to double first, before the subtraction.

double((img(i, j, 1)) - double(img(i, j+1, 1))

Answer 2

[Image Analyst]

Please post your images. This does not seem like a particularly intelligent or robust algorithm but it might be okay for certain images. We can advise you better once we see your image(s).

You might also take a look at algorithms like http://www.mathworks.com/matlabcentral/fileexchange/37197-dem-diffused-expectation-maximisation-for-image-segmentation

[David Young]

I was thinking he might be heading towards the watershed transform, but there's too little information so far.

[Rajesh Gothwal]

Yes, The desired output image will be like a watershed transform but watershed transform suffers with over-segmentation so I am using this code and region color will be the original image color not the color generated by colormap. I have posted images and complete code. Please look into the matter.

[David Young]

As far as I can see, your code only ever modifies the image at the seed locations - the array bin is never modified, so each iteration makes changes at the same place. So there's no way that a value can propagate across a region.

If your problem is that the watershed oversegments, there's a demo that addresses it exactly. Look at this: http://www.mathworks.co.uk/products/demos/image/watershed/ipexwatershed.html - looks like just the thing, assuming you have the Image Processing Toolbox.

[Image Analyst]

That's one of my favorite demos. Rajesh's code reminded me a little bit of anisotropic diffusion the way he calculates the difference squared in the 8 directions. I've attached a demo of that, it's really nice.

But what I don't understand is why the code does an Otsu threshold on the grayscale version and then seems to only process within that region. And then when it does you get some hysteresis or something because when you're processing the (i,j) pixel, you make changes to a neighbor of that which is then checked as the 3x3 window moves away. I guess that is how it's propagating the region. Some comments would help.

Another problem is T. T is number between 0 and 1 while d1 through d8 are in the range 0-255 so the condition d1<=T will only ever be met for pixels that are exactly 0. I'm not really sure what the intent there is. So why are you setting the pixel to 0 (which is what temp is at those locations) when the pixel is already zero there? Again adding comments may help you talk this out with yourself to discover errors. Also, displaying intermediate images may help reveal what is going wrong.