---

ImageGrey.cpp

Go to the documentation of this file.

/**************************************************************************
 *
 * Source code for Real Time Image Library (libRTImage)
 *
 * Leeds Vision Group give permission for this code to be copied, modified 
 * and distributed within the University of Leeds subject to the following 
 * conditions:-
 *
 * - The code is not to be used for commercial gain.
 * - The code and use thereof will be attributed to the author where
 *   appropriate (inluding demonstrations which rely on it's use).
 * - All modified, distributions of the source files will retain this header.
 *
 ****************************************************************************/

#include <new>          // bad_alloc exception
#include <string.h>     // memcpy
#include <MicrosoftBMP_RT.h>
#include <MicrosoftRCBMP_RT.h>
#include <MicrosoftMonoBMP_RT.h>
#include <math.h>
#include <limits.h>
#include <stdlib.h>
#include <iostream>

using namespace std;
 
#include "ImageGrey.h"

static int int_compare(const void *i1, const void *i2) ;

ImageGrey::ImageGrey(ImageGrey& img) : ImageBase(img)
{
        unsigned int array_size = img.width * img.height;
        brightness = new PixelGrey[array_size];

        *this = img;
}


ImageGrey::ImageGrey(unsigned int w, unsigned int h) : ImageBase(w,h)
{
        unsigned int array_size = width * height;
        brightness = new PixelGrey[array_size];

//      clear();
}

ImageGrey::ImageGrey(char *filename, FileFormat format)
{

    if(format == FORMAT_BMP_24){
        MicrosoftBMP_RT bmp(filename) ;

        if(!bmp.read()) throw("Invalid Bitmap file specified") ;

        width = bmp.ImageWidth ;
        height = bmp.ImageHeight ;

        brightness = new int[width*height] ;

        bmp.extract_data_grey((int *)brightness) ;
    }
    else if(format == FORMAT_BMP_RC){
        MicrosoftRCBMP_RT bmp(filename) ;

        if(!bmp.read()) throw("Invalid Bitmap file specified") ;

        width = bmp.ImageWidth ;
        height = bmp.ImageHeight ;

        brightness = new int[width*height] ;

        bmp.extract_data_grey((int *)brightness) ;
    }
    else if(format == FORMAT_BMP_MONO){
        MicrosoftMonoBMP_RT bmp(filename) ;

        if(!bmp.read()) throw("Invalid Bitmap file specified") ;

        width = bmp.ImageWidth ;
        height = bmp.ImageHeight ;

        brightness = new int[width*height] ;

        bmp.extract_data_grey((int *)brightness) ;
    }
    else if(format == FORMAT_BMP_ANY){
    // Check all formats
         MicrosoftBMP_RT* bmp ;
         MicrosoftRCBMP_RT *bmp_rc;
         MicrosoftMonoBMP_RT *bmp_mono;

         if((bmp = new  MicrosoftBMP_RT(filename)) != NULL &&
            bmp->read()){
         // 24 bit bmp
            width = bmp->ImageWidth ;
            height = bmp->ImageHeight ;

            brightness = new int[width*height] ;

            bmp->extract_data_grey((int *)brightness) ;
         }
         else if((bmp_rc = new  MicrosoftRCBMP_RT(filename)) != NULL &&
            bmp_rc->read()){
         // Reduced clour (palette) bmp
            width = bmp_rc->ImageWidth ;
            height = bmp_rc->ImageHeight ;

            brightness = new int[width*height] ;

            bmp_rc->extract_data_grey((int *)brightness) ;

         }
         else if((bmp_mono = new  MicrosoftMonoBMP_RT(filename)) != NULL &&
            bmp_mono->read()){
         // Monochrome (b+w) bmp
            width = bmp_mono->ImageWidth ;
            height = bmp_mono->ImageHeight ;

            brightness = new int[width*height] ;

            bmp_mono->extract_data_grey((int *)brightness) ;
         }
    }




}


ImageGrey::ImageGrey(ImageGrey&    img,
                     unsigned int min_x,
                     unsigned int min_y,
                     unsigned int max_x,
                     unsigned int max_y) : 
                                        ImageBase(max_x-min_x+1, max_y-min_y+1)
{
        brightness = new int[width*height];
        copy(img, min_x, min_y, max_x, max_y, 0, 0) ;
}


ImageGrey::~ImageGrey()
{
    delete [] brightness;
}


void ImageGrey::set_pixel(unsigned int x, unsigned int y, PixelGrey level)
{
        // Check that co-ords are valid. 
        if (x > width || y > height)
                throw "Invalid Coordinates in method set_rgb";

        unsigned int offset = y * width + x;
        brightness[offset]  = level;
}


void ImageGrey::get_pixel(unsigned int x, unsigned int y, PixelGrey& p)
{
        p = brightness[y * width + x];
}



ImageGrey& ImageGrey::clear(int k)
{
        int *start = brightness;
        int *end   = start + width * height;
                
        for (int *d = start; d < end; d++) *d = k;

        return *this;
}


ImageGrey& ImageGrey::copy(ImageGrey& img)
{       
        /* check we're not trying to copy ourself */
        if (this == &img)
                throw "Source and Destination are same Image";
        
        /* First check images are the same size and valid */
        if (width  != img.width || height != img.height) 
                throw "Difference in image dimensions"; 
        
        /* Perform copy */
        memcpy(brightness, img.brightness, width * height * sizeof(int));

        return *this;
}

ImageGrey& ImageGrey::copy(ImageGrey&    img,
                           unsigned int min_x,
                           unsigned int min_y,
                           unsigned int max_x,
                           unsigned int max_y,
                           unsigned int start_x,
                           unsigned int start_y)
{
        unsigned int sx ;
        unsigned int sy ;
        unsigned int cntx ;
        unsigned int cnty ;
        int          *s_ptr ;
        int          *d_ptr ;
        unsigned int s_width ;

        s_width = img.get_width() ;

        /* check we're not trying to copy ourself */
        if (this == &img)
                throw "Source and Destination are same Image";

        /* check sanity */
        if(min_x > max_x || min_y > max_y)
                throw "Invalid arguments to ImageRGB::copy()" ;

        /* Check images are big enough */
        sx = max_x - min_x + 1 ;
        sy = max_y - min_y + 1 ;

        if(sx+start_x > width ||
           sy+start_y > height)
                throw "Dest. image is not big enough" ;

         if(max_x > s_width ||
            max_y > img.get_height() )
                throw "Source images is not big enough" ;

         for(cnty=0 ; cnty<sy ; cnty++){

                 s_ptr = (int*)img.brightness + (min_x + (min_y+cnty)*s_width) ;
                 d_ptr = (int*)brightness + (start_x + (start_y+cnty)*width) ;

                 for(cntx=0 ; cntx<sx ; cntx++){
                                     
                         *(d_ptr++) = *(s_ptr++) ;
                 }
         }

         return *this;
}

ImageGrey& ImageGrey::operator += (int fact) 
{
    unsigned int cnt ;
    unsigned int sz ;
    int          *p1 ;

    sz = width * height ;
    p1 = brightness ;

    for(cnt=0 ; cnt<sz ; cnt++, p1++){
        *p1 += fact ;
    }

    return *this;
}

ImageGrey& ImageGrey::operator -= (int fact)
{
    unsigned int cnt ;
    unsigned int sz ;
    int          *p1 ;

    sz = width * height ;
    p1 = brightness ;

    for(cnt=0 ; cnt<sz ; cnt++, p1++){
        *p1 -= fact ;
    }

    return *this;
}

ImageGrey& ImageGrey::operator *= (int fact)
{
    unsigned int cnt ;
    unsigned int sz ;
    int          *p1 ;

    sz = width * height ;
    p1 = brightness ;

    for(cnt=0 ; cnt<sz ; cnt++, p1++){
        *p1 *= fact ;
    }

    return *this;
}

ImageGrey& ImageGrey::operator /= (int fact)
{
    unsigned int cnt ;
    unsigned int sz ;
    int          *p1 ;

    sz = width * height ;
    p1 = brightness ;

    for(cnt=0 ; cnt<sz ; cnt++, p1++){
        *p1 /= fact ;
    }

    return *this;
}


ImageGrey& ImageGrey::operator -= (ImageGrey& img) 
{
    unsigned int cnt ;
    unsigned int sz ;
    int          *p1 ;
    int          *p2 ;

    sz = width * height ;
    p1 = brightness ;
    p2 = img.brightness ;

    for(cnt=0 ; cnt<sz ; cnt++, p1++, p2++){
        *p1 -= *p2 ;
    } 

    return *this;
}

void ImageGrey::absolute()
{
    unsigned int cnt ; 
    unsigned int sz ;
    int          *p1 ;

    p1 = brightness ;
    sz = width * height ;

    for(cnt=0 ; cnt<sz ; cnt++, p1++){
        if(*p1<0) *p1 = -(*p1) ;
    }

    return ;
}

void ImageGrey::threshold(int thresh)
{
    unsigned int cnt ;
    unsigned int sz ;
    int          *p1 ;

    p1 = brightness ;
    sz = width * height ;

    for(cnt=0 ; cnt<sz ; cnt++, p1++){
        if(*p1<thresh) *p1 = 0 ;
        else           *p1 = 255 ;
    }

    return ;
}

bool ImageGrey::save(char *filename, FileFormat format)
{
 
    bool            ret_val = true ;

    if(format == FORMAT_BMP_24){
    // 24 Bit (full colour) bitmap 
        MicrosoftBMP_RT *btmp;
        btmp = new MicrosoftBMP_RT(filename);
        if(btmp==NULL) ret_val = false ;
        else if(!btmp->set_size(width,height)) ret_val = false ;
        else{
            btmp->put_data_grey((int *)brightness) ;
            ret_val = btmp->write();
        }
        delete btmp ;
    }
    if(format == FORMAT_BMP_RC || format == FORMAT_BMP_ANY){
    // Greyscale bitmap (+default)
        MicrosoftRCBMP_RT *btmp;
        btmp = new MicrosoftRCBMP_RT(filename);
        if(btmp==NULL) ret_val = false ;
        else if(!btmp->set_size(width,height)) ret_val = false ;
        else if(!btmp->set_palette_grey256()) ret_val = false ;
        else{
            btmp->put_data_grey((int *)brightness) ;
            ret_val = btmp->write();
        }
        delete btmp ;
    }
    if(format == FORMAT_BMP_MONO){
    // Greyscale bitmap
        MicrosoftMonoBMP_RT *btmp;
        btmp = new MicrosoftMonoBMP_RT(filename);
        if(btmp==NULL) ret_val = false ;
        else if(!btmp->set_size(width,height)) ret_val = false ;
        else{
            btmp->put_data_grey((int *)brightness) ;
            ret_val = btmp->write();
        }
        delete btmp ;
    }
 

 
    return ret_val ;
}

bool ImageGrey::sobel_horizontal(ImageGrey &result)
{
    bool     ret_val = true ;
    int      *ob_ptr ;
    int      *rs_ptr ;
    int      *tl_ptr ;
    int      *tr_ptr ;
    int      *b_ptr ;
    int      *t_ptr ;
    int      *bl_ptr ;
    int      *br_ptr ;
    unsigned int cnt_x ;
    unsigned int cnt_y ;
    register int pix_val ;

    if(result.get_width() != width   ||
       result.get_height() != height ){
        ret_val = false ;
    }
    else{
        ob_ptr = brightness ;
        rs_ptr = result.brightness ;
        tl_ptr = ob_ptr - width - 1 ;
        tr_ptr = ob_ptr - width + 1 ;
        t_ptr  = ob_ptr - width ;
        b_ptr  = ob_ptr + width ;
        bl_ptr = ob_ptr + width - 1 ;
        br_ptr = ob_ptr + width + 1 ;
        for(cnt_y=0 ; cnt_y<height ; cnt_y++){
            for(cnt_x=0 ; cnt_x<width ; cnt_x++, ob_ptr++, rs_ptr++, tl_ptr++,
                tr_ptr++,b_ptr++,t_ptr++,bl_ptr++,br_ptr++){
                if(cnt_y>0 &&
                   cnt_x>0 &&
                   cnt_y<height-1 &&
                   cnt_x<width-1 ){
                    pix_val = -1 * (*tl_ptr) +
                              -2 * (*t_ptr)  +
                              -1 * (*tr_ptr) +
                               1 * (*bl_ptr) +
                               2 * (*b_ptr)  +
                               1 * (*br_ptr)  ;
                    *rs_ptr = pix_val / 8 ;
                }
                else{
                    *rs_ptr = 0 ;
                }
            }
        }
    }

    return ret_val ;
 
}

bool ImageGrey::sobel_vertical(ImageGrey &result)
{
    bool     ret_val = true ;
    int      *ob_ptr ;
    int      *rs_ptr ;
    int      *tl_ptr ;
    int      *tr_ptr ;
    int      *l_ptr ;
    int      *r_ptr ;
    int      *bl_ptr ;
    int      *br_ptr ;
    unsigned int cnt_x ;
    unsigned int cnt_y ;
    register int pix_val ;

    if(result.get_width() != width   ||
       result.get_height() != height ){
        ret_val = false ;
    }
    else{
        ob_ptr = brightness ;
        rs_ptr = result.brightness ;
        l_ptr = ob_ptr - 1 ;
        r_ptr = ob_ptr + 1 ;
        tl_ptr = ob_ptr - width - 1 ;
        tr_ptr = ob_ptr - width + 1 ;
        bl_ptr = ob_ptr + width - 1 ;
        br_ptr = ob_ptr + width + 1 ;
        for(cnt_y=0 ; cnt_y<height ; cnt_y++){
            for(cnt_x=0 ; cnt_x<width ; cnt_x++, ob_ptr++, rs_ptr++, tl_ptr++,
                tr_ptr++,l_ptr++,r_ptr++,bl_ptr++,br_ptr++){
               if(cnt_y>0 &&
                   cnt_x>0 &&
                   cnt_y<height-1 &&
                   cnt_x<width-1 ){

                    pix_val = -1 * *tl_ptr +
                              -2 * *l_ptr  +
                              -1 * *bl_ptr +
                               1 * *tr_ptr +
                               2 * *r_ptr  +
                               1 * *br_ptr ;
                    *rs_ptr = pix_val / 8 ;
                }
                else{
                    *rs_ptr = 0 ;
                }
            }
        }
    }

    return ret_val ;
}

bool ImageGrey::sobel(ImageGrey &result)
{
    bool     ret_val = true ;
    int      *ob_ptr ;
    int      *rs_ptr ;
    int      *tl_ptr ;
    int      *tr_ptr ;
    int      *l_ptr ;
    int      *r_ptr ;
    int      *t_ptr ;
    int      *b_ptr ;
    int      *bl_ptr ;
    int      *br_ptr ;
    unsigned int cnt_x ;
    unsigned int cnt_y ;
    register int pix_val_h ;
    register int pix_val_v ;

    if(result.get_width() != width   ||
       result.get_height() != height ){
        ret_val = false ;
    }
    else{
        ob_ptr = brightness ;
        rs_ptr = result.brightness ;
        l_ptr = ob_ptr - 1 ;
        r_ptr = ob_ptr + 1 ;
        t_ptr  = ob_ptr - width ;
        b_ptr  = ob_ptr + width ;
        tl_ptr = ob_ptr - width - 1 ;
        tr_ptr = ob_ptr - width + 1 ;
        bl_ptr = ob_ptr + width - 1 ;
        br_ptr = ob_ptr + width + 1 ;
        for(cnt_y=0 ; cnt_y<height ; cnt_y++){
            for(cnt_x=0 ; cnt_x<width ; cnt_x++, ob_ptr++, rs_ptr++, tl_ptr++,
                tr_ptr++,l_ptr++,r_ptr++,bl_ptr++,br_ptr++, t_ptr++, b_ptr++){
               if(cnt_y>0 &&
                   cnt_x>0 &&
                   cnt_y<height-1 &&
                   cnt_x<width-1 ){

                    pix_val_h = -1 * (*tl_ptr) +
                                -2 * (*t_ptr)  +
                                -1 * (*tr_ptr) +
                                 1 * (*bl_ptr) +
                                 2 * (*b_ptr)  +
                                 1 * (*br_ptr)  ;

                    pix_val_v = -1 * *tl_ptr +
                                -2 * *l_ptr  +
                                -1 * *bl_ptr +
                                 1 * *tr_ptr +
                                 2 * *r_ptr  +
                                 1 * *br_ptr ;

                    *rs_ptr = (int)sqrt((double)(pix_val_h*pix_val_h + 
                                                 pix_val_v*pix_val_v)) ;
                }
                else{
                    *rs_ptr = 0 ;
                }
            }
        }
    }

    return ret_val ;
}

bool ImageGrey::dilate(ImageGrey &result)
{
    bool         ret_val = true ;
    unsigned int cnt_x ;
    unsigned int cnt_y ;
    int          *ob_ptr ;
    int          *rs_ptr ;
    int          *ne_ptr ;

    ob_ptr = brightness ;
    rs_ptr = result.brightness ;

    if(result.get_width() != width   ||
       result.get_height() != height ){
    /* Inputs are invalid */
        ret_val = false ;
    }
    else{
        for(cnt_y=0; cnt_y<height; cnt_y++){
            for(cnt_x=0; cnt_x<width; cnt_x++, ob_ptr++, rs_ptr++){
                *rs_ptr = *ob_ptr ;
                if( cnt_y>0 &&
                    cnt_x>0 &&
                    cnt_x<width-1 &&
                    cnt_y<height-1){
                /* Pixel is NOT a border pixel */
                    if(!(*ob_ptr)){
                    /* Pixel is black */
                        ne_ptr = ob_ptr - width - 1 ;
                        if(*(ne_ptr++)) *rs_ptr = 255 ;
                        if(*(ne_ptr++)) *rs_ptr = 255 ;
                        if(*ne_ptr)     *rs_ptr = 255 ;
                        ne_ptr += width ;
                        if(*(ne_ptr--)) *rs_ptr = 255 ;
                        ne_ptr-- ;
                        if(*ne_ptr)     *rs_ptr = 255 ;
                        ne_ptr += width ;
                        if(*(ne_ptr++)) *rs_ptr = 255 ;
                        if(*(ne_ptr++)) *rs_ptr = 255 ;
                        if(*ne_ptr)     *rs_ptr = 255 ;
                    }
                }
            }
        }
    }

    return ret_val ;

}

bool ImageGrey::errode(ImageGrey &result)
{
    return erode(result) ;
}

bool ImageGrey::erode(ImageGrey &result)
{
    bool         ret_val = true ;
    unsigned int cnt_x ;
    unsigned int cnt_y ;
    int          *ob_ptr ;
    int          *rs_ptr ;
    int          *ne_ptr ;

    ob_ptr = brightness ;
    rs_ptr = result.brightness ;

    if(result.get_width() != width   ||
       result.get_height() != height ){
    /* Inputs are invalid */
        ret_val = false ;
    }
    else{
        for(cnt_y=0; cnt_y<height; cnt_y++){
            for(cnt_x=0; cnt_x<width; cnt_x++, ob_ptr++, rs_ptr++){
                *rs_ptr = *ob_ptr ;
                if( cnt_y>0 &&
                    cnt_x>0 &&
                    cnt_x<width-1 &&
                    cnt_y<height-1){
                /* Pixel is NOT a border pixel */
                    if(*ob_ptr){
                    /* Pixel is white */
                        ne_ptr = ob_ptr - width - 1 ;
                        if(!*(ne_ptr++)) *rs_ptr = 0 ;
                        if(!*(ne_ptr++)) *rs_ptr = 0 ;
                        if(!*ne_ptr)     *rs_ptr = 0 ;
                        ne_ptr += width ;
                        if(!*(ne_ptr--)) *rs_ptr = 0 ;
                        ne_ptr-- ;
                        if(!*ne_ptr)     *rs_ptr = 0 ;
                        ne_ptr += width ;
                        if(!*(ne_ptr++)) *rs_ptr = 0 ;
                        if(!*(ne_ptr++)) *rs_ptr = 0 ;

                        if(!*ne_ptr)     *rs_ptr = 0 ;
                    }
                }
            }
        }
    }

    return ret_val ;
}

bool ImageGrey::extract_region(int region_no, ImageGrey &result)
{
    int    *ob_ptr ;
    int    *rs_ptr ;
    bool   ret_val = true ;

    if(result.get_width() != width    ||
       result.get_height() != height  ){
    /* Inputs are invalid */
        ret_val = false ;
    }
    else{
        for( ob_ptr = brightness, rs_ptr = result.brightness ;
             ob_ptr < brightness + width*height ;
             ob_ptr++, rs_ptr++){
        /* Loop through all pixels */
            if(*ob_ptr == region_no){
                *rs_ptr = 255 ;
            }
            else{
                *rs_ptr = 0 ;
            }
        }
    }

    return ret_val;
}

bool ImageGrey::region_identify(unsigned int max_regions,
                                unsigned int &actual_regions,
                                ImageGrey    &result)
{
    bool         ret_val = true ;
    int          next_label = 1 ;
    equivalence  *equiv_latest=NULL;
    equivalence  *equiv_ptr=NULL;
    int          *ob_ptr ;
    int          *rs_ptr ;
    int          *ne_ptr = NULL ;
    unsigned int cnt_x;
    unsigned int cnt_y;
    int          cnt;
    int          cnt2;
    unsigned int region_no;
    unsigned int region_freq;
    int          neighbours[4] ;
    int          *frequency;
    bool         is_significant;

    ob_ptr = brightness ;
    rs_ptr = result.brightness ;

    if( max_regions <= 0            ||
        result.get_width() != width  ||
        result.get_height() != height ){
    /* Inputs are invalid */
        ret_val = false ;
    }
    else{
        /*
         * First sweep set true pixels to :
         *      next_label if surrounded by bgd or not determined pixels.
         *   or lowest label of a neighbouring pixel.
         * If there is more than one option add an equivalence to the
         * equivalence linked list pointed to by equiv_start.
         */

        for( ob_ptr = brightness, rs_ptr = result.brightness, cnt_x=0,
             cnt_y=0;
             (unsigned long)(ob_ptr - brightness) <width*height;
             cnt_x++, ob_ptr++, rs_ptr++){
            /* Loop round all pixels, this is coded for speed not beauty! */
                *rs_ptr = 0;
                if(cnt_x==width){
                /* Reached end of line */
                    cnt_x=0;
                    cnt_y++;
                }
                else if( cnt_y>0           &&
                         cnt_y<height-1    &&
                         cnt_x>0           &&
                         cnt_x<width-1     &&
                         *ob_ptr           ){
                /* This is a true, non-border pixel */
                    ne_ptr = rs_ptr - width - 1 ;
                    neighbours[0] = *(ne_ptr++);
                    neighbours[1] = *(ne_ptr++);
                    neighbours[2] = *ne_ptr++;
                    neighbours[3] = *(rs_ptr-1);
                    qsort(neighbours , 4, sizeof(int), int_compare);

                    if(neighbours[3] == 0){
                    /* No processed neighbours are 'object' */
                        *rs_ptr = next_label++;
                    }
                    else{
                    /* Pixel has neighbours that are 'object' */
                        for(cnt=3; cnt>=0 && ret_val; cnt--){
                            if(neighbours[cnt]){
                                *rs_ptr=neighbours[cnt];
                                if(neighbours[cnt] != neighbours[3]){
                                    ret_val = new_equivalence(neighbours[cnt],
                                                              neighbours[3],
                                                              &equiv_latest);
                                }
                            }
                        }
                    }
                }
         }
        /*
         * Second sweep sets all equivalences to the lowest. We start with the
         * highest equivalence.
         */
        frequency = new int[next_label];

        if(frequency==NULL) ret_val = false ;
        for(cnt=0; cnt<next_label && ret_val; cnt++) frequency[cnt]=0;

        for( ob_ptr = brightness, rs_ptr = result.brightness, cnt_x=0,
             cnt_y=0;
             (unsigned long)(ob_ptr - brightness) <width*height && ret_val;
             cnt_x++, ob_ptr++, rs_ptr++){
            /* Loop round all pixels, this is coded for speed not beauty! */
                if(cnt_x==width){
                /* Reached end of line */
                    cnt_x=0;
                    cnt_y++;
                }
                else if( cnt_y>0           &&
                         cnt_y<height-1    &&
                         cnt_x>0           &&
                         cnt_x<width-1     &&
                         *rs_ptr           ){
                /* This is a true, non-border pixel */
                    for(cnt=next_label-1; cnt>=0; cnt--){
                        equiv_ptr = equiv_latest ;
                        while(equiv_ptr != NULL){
                        /* Loop through equivalences */
                            if(*rs_ptr == cnt && equiv_ptr->eq2 == cnt){
                                *rs_ptr = equiv_ptr->eq1;
                            }
                            equiv_ptr = equiv_ptr->next ;
                        }
                    }
                    frequency[*rs_ptr]++;
                }
        }

    /*
     * At this point regions have a unique region numbers however they are not
     * in frequency order and there are gaps.
     */
        free_equivalence(&equiv_latest);

        actual_regions = 0 ;

        for(cnt=0; cnt<(int)max_regions; cnt++){
            region_freq=frequency[1];
            region_no = 1;
            for(cnt2=1; cnt2<next_label-1; cnt2++){
                if(frequency[cnt2] > (int)region_freq){
                    region_freq = frequency[cnt2];
                    region_no   = cnt2;
                }
            }
            if(region_freq){
                frequency[region_no] = 0;
                ret_val = new_equivalence(cnt+1,
                                          region_no,
                                          &equiv_latest);
                actual_regions++;
            }
        }

   /*
    * Third loop sets result to ordered value an removes smaller objects.
    */
        for( ob_ptr = brightness, rs_ptr = result.brightness, cnt_x=0,
             cnt_y=0;
             (unsigned long)(ob_ptr - brightness) <width*height && ret_val;
             cnt_x++, ob_ptr++, rs_ptr++){
            /* Loop round all pixels, this is coded for speed not beauty! */
                if(cnt_x==width){
                /* Reached end of line */
                    cnt_x=0;
                    cnt_y++;
                }
                else if( cnt_y>0           &&
                         cnt_y<height-1    &&
                         cnt_x>0           &&
                         cnt_x<width-1     &&
                         *rs_ptr           ){
                /* This is a true, non-border pixel */
                        equiv_ptr = equiv_latest ;
                        is_significant=false;
                        while(equiv_ptr != NULL && !is_significant){
                        /* Loop through equivalences */
                            if(*rs_ptr == equiv_ptr->eq2){
                                *rs_ptr = equiv_ptr->eq1;
                                is_significant=true;
                            }
                            equiv_ptr = equiv_ptr->next ;
                        }
                        if(!is_significant){
                        /* Remove small object */
                            *rs_ptr = 0;
                        }
                }
        }
    }

    return ret_val;

}

static int int_compare(const void *i1, const void *i2)
{
/*
 * Function to compare two integers for qsort().
 */

    int ret_val ;

    if( *((int *)i1) < *((int *)i2) ){
        ret_val = -1 ;
    }
    else{
        ret_val = 1 ;
    }
    return ret_val ;
}

bool ImageGrey::new_equivalence(int low, int high, equivalence **latest)
{
/*
 * This method adds a new equivalence to the linked list who's last element is
 * pointed to by *latest. It returns false if memory allocation is unsuccessful.
 */

    bool        ret_val = true;
    equivalence *new_eq=NULL;
    equivalence *eq_ptr = *latest ;
    bool        exists = false ;

    /* First search for an equivalence with the same upper number */

    while(eq_ptr != NULL && !exists){
        if(eq_ptr->eq2 == high) {exists = true;}
        else                    {eq_ptr = eq_ptr->next;}
    }

    if(!exists){
    /* Create an equivalence */
        new_eq = new equivalence ;
        if(new_eq == NULL) ret_val = false;
    }

    if(ret_val && exists && low<eq_ptr->eq1){
    /* Modify existing equivalence and create new one */
        ret_val = new_equivalence(low,eq_ptr->eq1,latest);
        eq_ptr->eq1 = low;
    }
    else if(ret_val && exists && low>eq_ptr->eq1){
        ret_val = new_equivalence(eq_ptr->eq1, low, latest);
    }

    if (ret_val && !exists){

        new_eq->eq1 = low;
        new_eq->eq2 = high;
        new_eq->next = *latest;
        *latest = new_eq ;
    }

    return ret_val;
}

void ImageGrey::free_equivalence(equivalence **eq)
{
/*
 * Function to free memory allocated to equivalence linked list.
 */

    equivalence *eq_ptr1 = *eq;
    equivalence *eq_ptr2;

    while(eq_ptr1 != NULL){
        eq_ptr2 = eq_ptr1->next ;
        delete eq_ptr1;
        eq_ptr1= eq_ptr2 ;
    }

    *eq = NULL;
}

void ImageGrey::normalise_intensity(float mean,
                                    float sd)
{
    unsigned int cnt ;
    unsigned int sz ;
    int          *ptr ;
    float        c_mean = 0 ;
    float        c_sd = 0 ;
    float        nval ;
    float        sf ;

    sz = width*height ;

    // First calculate actual mean / sd
    for(cnt=0, ptr=brightness ; cnt<sz ; cnt++, ptr++){
        c_mean += *ptr ;
        c_sd += *ptr * *ptr ;
    }
    c_mean /= sz ;
    c_sd /= sz ;
    c_sd -= c_mean*c_mean ;
    c_sd = sqrt(c_sd) ;

    sf = sd / c_sd ;

    // Normalise 
     for(cnt=0, ptr=brightness ; cnt<sz ; cnt++, ptr++){
        nval = *ptr - c_mean ;
        nval *= sf ;
        nval += mean ;
        *ptr = (int)nval ;
    }
}

void ImageGrey::normalise_intensity_range() 
{
    unsigned int cnt ;
    unsigned int sz ;
    int          *ptr ;
    int          min = INT_MAX ;
    int          max = -INT_MAX ;
    float        sf ;

    sz = width*height ;

    // Calculate min / max
    for(cnt=0, ptr=brightness ; cnt<sz ; cnt++, ptr++){
        if(*ptr < min) min = *ptr ;
        if(*ptr > max) max = *ptr ;
    }

    if(max > min){
        sf = 255.0 / (float)(max - min) ;
    }
    else{
        return ;
    }

    // Normalise
    for(cnt=0, ptr=brightness ; cnt<sz ; cnt++, ptr++){
        *ptr -= min ;
        *ptr *= sf ;
    }
}

bool ImageGrey:: median(unsigned int n, unsigned int m, 
                        ImageGrey &result)
{
    bool         ret_val = true ;
    unsigned int border_x;
    unsigned int border_y;
    unsigned int image_x ;
    unsigned int image_y ;
    unsigned int matrix_x ;
    unsigned int matrix_y ;
    int          *neighbourhood;
    unsigned int n_count;

    border_x = n / 2 ;
    border_y = m/ 2 ;

    /* Create an array for neighbourhood data for use by qsort() */
    neighbourhood = new int[n*m] ;

    for(image_y=0; image_y<height; image_y++){
        for(image_x=0; image_x<width; image_x++){
            if(image_y<border_y ||
               image_x<border_x ||
               image_y>=height-border_y ||
               image_x>=width-border_x){
            /* This is a border pixel, set it to black */
                result.brightness[image_y*width + image_x] = 0;
            }
            else{
                n_count = 0 ;
                for(matrix_x=0; matrix_x<n; matrix_x++){
                    for(matrix_y=0; matrix_y<m; matrix_y++){
                        neighbourhood[n_count] = 
                             brightness[(image_y+matrix_y-border_y)*width
                                         + image_x + matrix_x - border_x] ;
                        n_count++;
                    }
                }
                qsort(neighbourhood , n*m, sizeof(int), int_compare);
                result.brightness[image_y*width + image_x] =
                                          neighbourhood[n*m/2+1] ;
            }
        }
    }

    delete [] neighbourhood ;

    return ret_val ;
}

bool ImageGrey::remove_point_noise(ImageGrey &result)
{
    unsigned int image_x ;
    unsigned int image_y ;
    int          *ptr ;
    int          *rptr ;

    ptr = brightness ;
    rptr = result.brightness ;
    for(image_y=0; image_y<height; image_y++){
        for(image_x=0; image_x<width; image_x++, ptr++, rptr++){
            if(*ptr &&
               ((image_x==0 && image_y==0) || *(ptr-width-1)==0) &&
               (image_y==0 || *(ptr-width)==0) &&
               ((image_x==width-1 && image_y==0) || *(ptr-width+1)==0) &&
               (image_x==0 || *(ptr-1)==0) &&
               (image_x==width-1 || *(ptr+1)==0) &&
               ((image_x==0 && image_y==height-1) || *(ptr+width-1)==0) &&
               (image_y==height-1 || *(ptr+width)==0) &&
               ((image_x==width-1 && image_y==height-1) || *(ptr+width+1)==0)){
            // True, but no neighbours
                *rptr = 0 ;
            }
            else{
                *rptr = *ptr ;
            }
        }
    } 

    return true ;
}

---