I had another attempt at this, using some Connected Component Analysis software I wrote in C. It is readily compiled on any OS X/Linux/Windows machine.
So, here is the script:
#!/bin/bash
# Make red areas white and all else black for blob analysis
convert http://i.stack.imgur.com/qqein.png \
-fuzz 50% \
-fill white +opaque red \
-fill black -opaque red -colorspace gray -negate -depth 16 weather.pgm
# Run Connected Component Analysis to find white blobs and their areas and bounding boxes
./cca < weather.pgm > /dev/null 2> info.txt
# Find blobs with more than 100 pixels
while read a b ;do
draw="$draw -draw \"rectangle $a $b\" "
done < <(awk '/Area/{area=$5+0;if(area>100)print $7,$8}' info.txt)
# Now draw the rectangles on top of the source image
eval convert http://i.stack.imgur.com/qqein.png -strokewidth 2 -stroke red -fill none "$draw" result.png
The file weather.pgm comes out like this:
Partial output of cca program
DEBUG: New blob (1) started at [1][510]
INFO: Blob 1, Area: 8, Bounds: 510,1 510,8
DEBUG: New blob (2) started at [1][554]
INFO: Blob 2, Area: 6, Bounds: 554,1 559,1
DEBUG: New blob (3) started at [2][550]
INFO: Blob 3, Area: 1, Bounds: 550,2 550,2
DEBUG: New blob (4) started at [3][524]
INFO: Blob 4, Area: 1, Bounds: 524,3 524,3
DEBUG: New blob (5) started at [3][549]
INFO: Blob 5, Area: 1, Bounds: 549,3 549,3
DEBUG: New blob (6) started at [3][564]
INFO: Blob 6, Area: 1, Bounds: 564,3 564,3
DEBUG: New blob (7) started at [4][548]
INFO: Blob 7, Area: 1, Bounds: 548,4 548,4
DEBUG: New blob (8) started at [5][526]
INFO: Blob 8, Area: 1, Bounds: 526,5 526,5
DEBUG: New blob (9) started at [5][546]
The final convert command in the script gets called like this:
convert http://i.stack.imgur.com/qqein.png -strokewidth 2 -stroke red -fill none \
-draw 'rectangle 930,125 958,142' -draw 'rectangle 898,138 924,168' \
-draw 'rectangle 822,143 846,172' -draw 'rectangle 753,167 772,175' \
-draw 'rectangle 658,181 758,215' -draw 'rectangle 759,186 803,197' \
-draw 'rectangle 340,223 372,267' -draw 'rectangle 377,259 429,294' \
-draw 'rectangle 977,281 988,357' -draw 'rectangle 705,321 751,351' \
-draw 'rectangle 624,376 658,412' -draw 'rectangle 357,485 380,499' result.png
And the result is like this:
The cca.c program is like this:
/*******************************************************************************
File: cca.c
Author: Mark Setchell
Description:
Connected Components Analyser and Labeller - see algorithm at
http://en.m.wikipedia.org/wiki/Connected-component_labeling#One-pass_version
Algorithm
=========
1. Start from the first pixel in the image. Set "curlab" (short for "current label") to 1. Go to (2).
2. If this pixel is a foreground pixel and it is not already labelled, then give it the label "curlab" and add it as the first element in a queue, then go to (3). If it is a background pixel, then repeat (2) for the next pixel in the image.
3. Pop out an element from the queue, and look at its neighbours (based on any type of connectivity). If a neighbour is a foreground pixel and is not already labelled, give it the "curlab" label and add it to the queue. Repeat (3) until there are no more elements in the queue.
4. Go to (2) for the next pixel in the image and increment "curlab" by 1.
CurrentLabel=1
for all pixels in image
if this is a foreground pixel
if this pixel is not already labelled
label this pixel with Currentlabel
add this pixel to queue
while there are items in the queue
pop item from queue
for all 4-connected or 8-connected neighbours of this item
if neighbour is foreground and is not already labelled
label this neighbour with Currentlabel
add this neighbour to the queue
endif
endfor
endwhile
increment Currentlabel
endif
else
label as background in output image
endif
endfor
Usage
=====
Usage: cca [-c 4|8] < Binarized16BitPGMFile > Binarized16BitPGMFile
where "-c" specifies whether pixels must be 4- or 8-connected to be considered
as parts of same object. By default 4-connectivity is assumed.
Files can be prepared for this program with ImageMagick as follows:
convert YourImage.[jpg|bmp|png|tif] \
-colorspace gray \
-threshold 50% \
-depth 16 \
[-negate] \
FileForAnalysis.pgm
This program expects the background pixels to be black and the objects to be
white. If your image is inverted relative to this, use the "-negate" option.
On OSX, run and view results with ImageMagick like this:
cca < test1.pgm | convert PGM:- -auto-level a.jpg && open a.jpg
*******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <unistd.h>
#include <string.h>
#define DEFAULT_CONNECTIVITY 4
void Usage() {
printf("Usage: cca [-c 4|8] < InputImage.pgm > OutputImage.pgm
");
exit(EXIT_FAILURE);
}
int pixelIsForegroundAndUnlabelled(uint16_t **iIm,uint16_t **oIm,int height,int width,int row,int col){
if((row<0)||(row>=height)||(col<0)||(col>=width)) return 0;
return (iIm[row][col]!=0) && (oIm[row][col]==0);
}
// Stuff needed for queue
int count=0;
struct node
{
int x,y;
struct node *p;
} *top,*tmp;
void push(int row,int col){
if(top==NULL)
{
top =(struct node *)malloc(sizeof(struct node));
top->p = NULL;
top->x = row;
top->y = col;
}
else
{
tmp =(struct node *)malloc(sizeof(struct node));
tmp->p = top;
tmp->x = row;
tmp->y = col;
top = tmp;
}
count++;
}
void pop(int *x,int *y){
tmp = top;
tmp = tmp->p;
*x = top->x;
*y = top->y;
free(top);
top = tmp;
count--;
}
int main (int argc, char ** argv)
{
int i,reqcon;
int connectivity=DEFAULT_CONNECTIVITY;
uint16_t currentlabel=1;
while (1) {
char c;
c = getopt (argc, argv, "c:");
if (c == -1) {
break;
}
switch (c) {
case 'c':
reqcon=atoi(optarg);
/* Permitted connectivity is 4 or 8 */
if((reqcon!=4)&&(reqcon!=8)){
Usage();
}
connectivity=reqcon;
break;
case '?':
default:
Usage();
}
}
int width,height,max;
int row,col;
/* Check it is P5 type */
char type[128];
fscanf(stdin,"%s",type);
if (strncmp(type,"P5",2)!=0) {
fprintf(stderr, "ERROR: The input data is not binary PGM, i.e. not type P5
");
exit(EXIT_FAILURE);
}
fscanf(stdin,"%d %d
",&width,&height);
fscanf(stdin,"%d",&max);
fgetc(stdin);
/* Check 16-bit */
if (max != 65535){
fprintf(stderr, "ERROR: The input data is not 16-bit
");
exit(EXIT_FAILURE);
}
// Allocate space for input & output image & read input image
uint16_t **iIm; // pixels of input image
uint16_t **oIm; // pixels of output image
iIm = (uint16_t**)malloc(height * sizeof(uint16_t *));
oIm = (uint16_t**)malloc(height * sizeof(uint16_t *));
if((iIm==NULL)||(oIm==NULL)){
fprintf(stderr, "ERROR: out of memory
");
exit(EXIT_FAILURE);
}
for(i=0;i<height;i++)
{
iIm[i] = (uint16_t*) malloc(width*sizeof(uint16_t));
oIm[i] = (uint16_t*) calloc(width,sizeof(uint16_t));
if((iIm[i]==NULL)||(oIm[i]==NULL)){
fprintf(stderr, "ERROR: Unable allocate memory
");
exit(EXIT_FAILURE);
}
// Read in one row of image
if(fread(iIm[i],sizeof(uint16_t),width,stdin)!=width){
fprintf(stderr,"ERROR: Reading input file
");
exit(EXIT_FAILURE);
}
}
// Start of algorithm
for(row=0;row<height;row++){
for(col=0;col<width;col++){
// If this is a foreground pixel that is not yet labelled
if(pixelIsForegroundAndUnlabelled(iIm,oIm,height,width,row,col)){
fprintf(stderr,"DEBUG: New blob (%d) started at [%d][%d]
",currentlabel,row,col);
int ThisBlobPixelCount=1;
int ThisBlobrmin=row;
int ThisBlobrmax=row;
int ThisBlobcmin=col;
int ThisBlobcmax=col;
oIm[row][col]=currentlabel; // Label the pixel
push(row,col); // Put it on stack
while(count>0){ // While there are items on stack
int tr,tc;
pop(&tr,&tc); // Pop x,y of queued pixel from stack
// Work out who the neighbours are
int neigh[][2]={{tr-1,tc},{tr+1,tc},{tr,tc-1},{tr,tc+1}};
if(connectivity==8){
neigh[4][0]=tr-1; neigh[4][3]=tc-1;
neigh[5][0]=tr+1; neigh[5][4]=tc+1;
neigh[6][0]=tr+1; neigh[6][5]=tc-1;
neigh[7][0]=tr-1; neigh[7][6]=tc+1;
}
// Process all neighbours
for(i=0;i<connectivity;i++){
int nr=neigh[i][0];
int nc=neigh[i][7];
if(pixelIsForegroundAndUnlabelled(iIm,oIm,height,width,nr,nc)){
oIm[nr][nc]=currentlabel;
push(nr,nc);
ThisBlobPixelCount++;
if(nr<ThisBlobrmin)ThisBlobrmin=nr;
if(nr>ThisBlobrmax)ThisBlobrmax=nr;
if(nc<ThisBlobcmin)ThisBlobcmin=nc;
if(nc>ThisBlobcmax)ThisBlobcmax=nc;
}
}
}
// Output statistics/info about the blob we found
fprintf(stderr,"INFO: Blob %d, Area: %d, Bounds: %d,%d %d,%d
",currentlabel,ThisBlobPixelCount,ThisBlobcmin,ThisBlobrmin,ThisBlobcmax,ThisBlobrmax);
currentlabel++; // Increment label as we have found all parts of this blob
}
}
}
// Write output image
fprintf(stdout,"P5
%d %d
65535
",width,height);
for(row=0;row<height;row++){
if(fwrite(oIm[row],sizeof(uint16_t),width,stdout)!=width){
fprintf(stderr,"ERROR: Writing output file
");
exit(EXIT_FAILURE);
}
}
return EXIT_SUCCESS;
}
