Adaptation of Selective Gaussian Blur Noise Reduction (SGBNR) that works with PPM images

selective_gaussian_blur.c

/** CENG466 HW1 - 2006/2007 - 1st Semester
*
* - Team -
* Cagatay Calli
* Felix Wolfsteller
*
* - Summary -
* Implemented an original adaptation of Selective Gaussian Blur (SGBNR)
*(given name: Varying Multi-pass Isotropic Selective Gaussian Blur, VMISGB), 
* using a Gaussian mask having double type values. 
* This Gaussian mask is generated by mask_gaussian() function automatically,
* according to given sigma and mask size.
* Image is filtered with a median filter before SGBNR.
*
**/

#include <string>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <vector>
#include <math.h>

#define IND(i,j,w)  ((i)*(w)+(j))

#define GAUSS_SCALER 1
#define GAUSS_SIGMA 30

using namespace std;

typedef struct pixel {
  unsigned short r,g,b;
   
  void operator=(struct pixel p){
     r = p.r;
     g = p.g;
     b = p.b;
  }
} PIXEL;

class Image {
public:
          short type;
          short num_bytes;
          int width,height;
          int maxval;
          PIXEL *data;
          
          
          // constructors, destructor and overloaded operators
          
          Image() {
              num_bytes = 2;
              data = NULL;        
          }

	  
	  Image(const Image &copy){
	      type = copy.type;
	      num_bytes = copy.num_bytes;
	      width = copy.width;
	      height = copy.height;
	      maxval = copy.maxval;
	      
	      data = (PIXEL *) malloc(sizeof(PIXEL)*width*height);

	      for(int i=0;i<width*height;i++){
	        data[i] = copy.data[i];
	      }
	  }
	  
	  Image& operator=(const Image &b){
         type = b.type;
         num_bytes = b.num_bytes;
         width = b.width;
         height = b.height;
         maxval = b.maxval;
         
         data = (PIXEL *) malloc(sizeof(PIXEL)*width*height);
         
  	      for(int i=0;i<width*height;i++){
	        data[i] = b.data[i];
	      }
	      
	      return *this;
                
      }
	  
      ~Image(){
         if(data != NULL) free(data);         
      }

      // end of constructors,destructor and overloaded operators

	  void emptyCopy(const Image &copy){
	     type = copy.type;
	     num_bytes = copy.num_bytes;
	     width = copy.width;
	     height = copy.height;
	     maxval = copy.maxval;

	     data = (PIXEL *) malloc(sizeof(PIXEL)*width*height);
	  }

	 void prepareOPImage(const Image &copy){
	     type = copy.type;
	     num_bytes = copy.num_bytes;
	     width = copy.width*2;
	     height = copy.height;
	     maxval = copy.maxval;

	     data = (PIXEL *) malloc(sizeof(PIXEL)*2*width*height);
             for(int i = 0; i<width*height/2; i++){
                  data[(width/2+ i/(width/2)*width )+i%(width/2)] = copy.data[i];
             }
	  }
	  
          
          void readppm(char * filename){
                   int i,j;
                   ifstream infile;
                   
                   infile.open(filename,ios::in | ios::binary);
          
                   if(infile.is_open()){
                            infile.ignore(); // read "P" character       
                            infile >> type; // read second character as file type
                            infile.ignore();
       
                            while(infile.peek() =='#') infile.ignore(500,'\n');
                            infile >> width;
       
                            while(infile.peek() =='#') infile.ignore(500,'\n');
                            infile >> height;
       
                            while(infile.peek() =='#') infile.ignore(500,'\n');
                            infile >> maxval;
                            infile.ignore();
       
                            if(maxval < 256){
                               num_bytes = 1;          
                            }
       
       
                            // create an array dynamically
                            data = (PIXEL *) malloc(sizeof(PIXEL)*width*height);
                            // array created
       
                            switch(type){
                                 case 6:
                                    if(num_bytes == 1){
                                      char rgb[3];
                     
                                      for(i=0;i<height;i++){
                                         for(j=0;j<width;j++){
                                              infile.read(rgb,3);
                             
                                              data[IND(i,j,width)].r = ((unsigned char) rgb[0]); // for red
                                              data[IND(i,j,width)].g = ((unsigned char) rgb[1]); // for green
                                              data[IND(i,j,width)].b = ((unsigned char) rgb[2]); // for blue
                                         }                      
                                     }
                            
                                    }
                            
                                    else if(num_bytes == 2){
                                       char rgb[6];
                                       unsigned short temp = 0;
                     
                                       for(i=0;i<height;i++){
                                           for(j=0;j<width;j++){
                                               infile.read(rgb,6);
                             
                                               temp = ((unsigned char) rgb[0]);
                                               temp = temp << 8;
                                               temp = temp | ((unsigned char) rgb[1]);
                                              data[IND(i,j,width)].r = temp; // for red
                                              
                                              temp = ((unsigned char) rgb[2]);
                                               temp = temp << 8;
                                               temp = temp | ((unsigned char) rgb[3]);
                                              data[IND(i,j,width)].g = temp; // for green
                                              
                                              temp = ((unsigned char) rgb[4]);
                                               temp = temp << 8;
                                               temp = temp | ((unsigned char) rgb[5]);
                                              data[IND(i,j,width)].b = temp; // for blue
                                           }                      
                                       }
                            
                                   }
               
                                   break;
                                   
                                default:
                                   break;
                                }
                   }
       
                   // image read
                   infile.close();          
          }
          
          void writeppm(char * filename){
                  int i,j;
                  ofstream outfile;
                  outfile.open (filename, ios::out | ios::trunc | ios::binary);
    

       
                  // start output
                  outfile << 'P' << type << ' ';
                  outfile << width << ' ' << height << ' ';
                  outfile << maxval << ' ';
       
       
                  for(i=0;i<height;i++){
                     for(j=0;j<width;j++){
                           // r1,g1,b1 are first(MSB) byte of 2-byte long values
                           PIXEL curpixel = data[IND(i,j,width)];
                           char r1,g1,b1;
                           
                           if(num_bytes == 2){
                               r1 = curpixel.r >> 8;
                               g1 = curpixel.g >> 8;
                               b1 = curpixel.b >> 8;
                           }
                           
                           if(num_bytes == 2) outfile.put(r1);
                           outfile.put(curpixel.r);
                           if(num_bytes == 2) outfile.put(g1);
                           outfile.put(curpixel.g); 
                           if(num_bytes == 2) outfile.put(b1);
                           outfile.put(curpixel.b);  
                     }                      
                  }

                  outfile.close();      
          }
          
};

/**
* sets r,g and b value of a pixel, specified through x and y in an image img
* to a value. The resulting color will appear grey.
**/
void setGrayValue(Image &img, int x, int y,int grayvalue){
	img.data[IND(x,y,img.width)].r = grayvalue;
	img.data[IND(x,y,img.width)].g = grayvalue;
	img.data[IND(x,y,img.width)].b = grayvalue;
}

/**
* median() finds the median of given integer vector.
*
**/
int median(const vector<int> &vec){
    bool even_sized = ( ((vec.size()+1)/2) == (vec.size()/2) );
    int middle = vec.size()/2;
    int result;
    vector<int> copy(vec);
    sort(copy.begin(), copy.end());
    
    if(even_sized){
       // for even sized vector:: median = middle + (middle-1) / 2
       result = (copy[middle-1] + copy[middle]) / 2;
    }
    else {
       result = copy[middle];
    }
    
    return result;
}


void median_filter(Image &input,Image &output,const vector<int> &med_mask){
     int mask_size = (int) sqrt((double)med_mask.size());
     int dist = (mask_size-1) / 2;
     int addnum;
     int cur_pix;
     int z;
     //output.emptyCopy(input);
     vector<int> vec;
     
     // for all pixels
     for(int i=0;i<input.height;i++){
        for(int j=0;j<input.width;j++){
            vec.clear();
                
            // build the vector to be sorted    
            for(int u=(-dist);u<=dist;u++){
               for(int t=(-dist);t<=dist;t++){
			      if( ((i+u) >=0) && ((i+u) < input.height) && ((j+t) >=0) && ((j+t) < input.width)){
                      cur_pix = input.data[IND(i+u,j+t,input.width)].r;
                      addnum = med_mask[IND(dist+u,dist+t,mask_size)];
                      for(z=0;z<addnum;z++){
                         vec.push_back(cur_pix);
                      }
			      }

			   }
		    }
		    
		    int med = median(vec);
		    setGrayValue(output, i,j,med);
        }
     }
     
}

/** mask_gaussian()
* Creates a Gaussian mask as big as the given vector size, with a given sigma.
**/
void mask_gaussian(vector<double> &matrix,double sigma){
    int mask_dim = (int) sqrt((double)matrix.size());
	int dist = (mask_dim-1)/2;
	int rsquare;
	double gauss_val;

	// 2*(sigma^2) constants
	double sigma2sqr = sigma*sigma*2;

	// gaussian multipliers
	double mult_gauss = 1 / (sigma2sqr * M_PI);
	
	// exponentials
	double exp_gauss;
	
	// calculate only for 1 quadrant and write all quadrants
	for(int i=(-dist);i<=0;i++){
          for(int j=(-dist);j<=0;j++){
              rsquare = (i*i+j*j); 
              exp_gauss = -rsquare / sigma2sqr;
              
              gauss_val = GAUSS_SCALER * (mult_gauss * exp(exp_gauss));
              //cout << gauss_val << endl;
              
              matrix[IND((i+dist),(j+dist),mask_dim)] = 
              matrix[IND((-i+dist),(j+dist),mask_dim)] =
              matrix[IND((i+dist),(-j+dist),mask_dim)] =
              matrix[IND((-i+dist),(-j+dist),mask_dim)] = gauss_val;
              
         }
    }
	
} 


/** smooth()
* Takes an a (double) mask and a pixel difference threshold (maxdelta) as input.
**/
void sel_gauss(Image &input, Image &output, vector<double> mask){
     double pixelcalc;
     double divider;
     int cur_pix,tmp_pix;
     //int mn,dev;
      double mult;
      int diff;
     
     // mask dimension (NxN) = (maskdim x maskdim)
     int maskdim = (int) sqrt((double)mask.size()); 

     int dist = (maskdim-1) / 2; // distance from center pixel to edge of mask
     int maskiter;
     
     //output.emptyCopy(input);

     // for all pixels
     for(int i=0;i<input.height;i++){
        for(int j=0;j<input.width;j++){
                 maskiter = 0;
                 divider = 0;
                 pixelcalc = 0;
                 
                 //mn = get_mean(input,i,j,dist);                 
                 cur_pix = input.data[IND(i,j,input.width)].r;

                 for(int u=(-dist);u<=dist;u++){
	                     for(int t=(-dist);t<=dist;t++){
                            // leave a one pixel boarder, due to missing information
		                    if( ((i+u) >=0) && ((i+u) < input.height) && ((j+t) >=0) && ((j+t) < input.width)){
           
                               tmp_pix = input.data[IND(i+u,j+t,input.width)].r;
                               
                               diff = abs(tmp_pix-cur_pix);
                               mult = ((double) diff / input.maxval);
                               mult = 1 - mult;
                               
                               if( 0.9 < mult ) mult = 1;
                               else if( (0.80 < mult) && (mult <= 0.9) ) mult = mult * 0.256;
                               else if( (0.70 < mult) && (mult <= 0.8) ) mult = mult * 0.016;
                               else if( (0.60 < mult) && (mult <= 0.7) ) mult = mult * 0.0020;
                               else if( (0.50 < mult) && (mult <= 0.6) ) mult = mult * 0.0005;
                               else if( (0.30 < mult) && (mult <= 0.50) ) mult = mult * 0.00025;
                               else if( mult <= 0.30) mult = 0;
                               
                               if(mult > 1) mult = 1;
                               
                               
                               if(mult > 0) { 
		                          pixelcalc += tmp_pix * mask[maskiter] * mult; 
		                          divider += mask[maskiter] * mult;  
                               }
                               maskiter++; 
                            }
                        }
                 }
                 
                 if(divider) pixelcalc = pixelcalc / divider;
                 else pixelcalc = cur_pix;

		    setGrayValue(output,i,j,(int)pixelcalc);
        }        
     } 
     
}

/**
* printUsage()
* prints help for the user, how to use this piece of art.
**/
void printUsage(){
     cout << "usage: hw1 input.ppm [output.ppm]\n"
          << "\tif output.ppm is dropped, output will be written in "
	  << "a file named after the pattern resultinput.ppm\n";
}


int main(int argc, char *argv[])
{
    string outfilename;
    int blur_radius = 1;
    int median_radius = 1;

    // basic user error handling

    switch(argc){
      case 2:
           outfilename = "result";
     	   outfilename.append(argv[1]);
           break;
      case 3:
           outfilename = argv[2];
           break;
      case 4:
           outfilename = argv[2];
           median_radius = atoi(argv[3]);
           break;
      case 5:
           outfilename = argv[2];
           median_radius = atoi(argv[3]);
           blur_radius = atoi(argv[4]);
           break;      
      default:
           printUsage();
	       exit(0);                 
	       break;
    }

    /* read PPM image file */
    Image my_image,temp_image,median_image,gauss_image;
    my_image.readppm(argv[1]);

    /* create a 3x3 mask vector
    *  for MEDIAN FILTER */
    // int med_dim = 2*median_radius+1;
    int med_a[] = {0,1,0,
                   1,1,1,
                   0,1,0 };
    vector<int> med(med_a,med_a+9);
    /* This code piece is for generating square-shaped median mask. 
    
      for(int w=0;w<med_dim*med_dim;w++){
      med[w] = 1;
    } */
    
    
    /* create an NxN Gaussian mask vector
    * for GAUSSIAN BLUR using sigma=30 */
    
    int mask_dim = 2*blur_radius + 1;
    
    vector<double> gauss(mask_dim * mask_dim);
    mask_gaussian(gauss,GAUSS_SIGMA);
    /* This code piece is for generating averaging blur.
    
       for(int w=0;w<blur_radius*blur_radius;w++){
       gauss[w] = 1;
    } */

     /* MEDIAN FILTER */
    /* anti salt n pepper */
    median_image.emptyCopy(my_image);
    median_filter(my_image,median_image,med);
    /** GAUSS and... **/
    gauss_image.emptyCopy(my_image);
    
    // do the main job.. selective gaussian blur.
    temp_image = median_image;
    
    for(int w=0;w<4;w++){
      sel_gauss(temp_image,gauss_image,gauss);
      temp_image = gauss_image;
      
      blur_radius++;
      mask_dim = 2*blur_radius + 1;
    
      gauss.resize(mask_dim * mask_dim);
      mask_gaussian(gauss,GAUSS_SIGMA);
    }
    
    /* 
    cout << "median_radius: " << median_radius << endl;
    cout << "blur_radius: " << blur_radius << endl;
    */
    
    gauss_image.writeppm((char*)outfilename.data());
           
    //system("PAUSE");
    return EXIT_SUCCESS;
}