CelesteTrain.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2008 by Tim Nugent                                      *
 *   timnugent@gmail.com                                                   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>. *
 ***************************************************************************/

#include <iostream>
#include "vigra/stdimage.hxx"
#include "vigra/resizeimage.hxx"
#include "vigra/impex.hxx"
#include "vigra/colorconversions.hxx"
#include <sys/types.h>
#include <sys/stat.h>
#include <stdlib.h>
#include "Gabor.h"
#include "Utilities.h"
#include "CelesteGlobals.h"

using namespace vigra; 
using namespace std; 

void usage (char* filename){

        cout << "Usage: " << filename << " <image file> <mask file>" << endl;
        cout << "or for a positive example: " << filename << " -p <image file>" << endl;
        cout << "or for a negative example: " << filename << " -n <image file>" << endl;
        cout << "(supported formats: " << vigra::impexListFormats() << ")" << endl;
}

int main(int argc, char ** argv){

        unsigned int a = 1;
        unsigned int pos = 0;
        unsigned int neg = 0;
        string image_file,mask_file;
        vigra::FRGBImage mask_in;
        vigra::FRGBImage::iterator mask_iter;

        // Deal with arguments
        if(argc < 3 || argc >= 4){
        
                usage(argv[0]); exit(1);
        }

        while( a < argc ){

                if( argv[a][0] == '-'){
    
                        a++;
                        switch(argv[a-1][1]){

                                // Put some more argument options in here later
                                case 'h' : {usage(argv[0]); exit(1);}
                                case 'p' : {pos = 1;
                                                cerr << "Ignoring mask. Image is a positive example." << endl;
                                                image_file = argv[argc-1];
                                                break;
                                            }
                                case 'n' : {neg = 1;
                                                cerr << "Ignoring mask. Image is a negative example." << endl;
                                                image_file = argv[argc-1];
                                                break;
                                                }
                                default  : {usage(argv[0]); exit(1);}
                        }

                }

                a++;
        }

        if (pos == 0 && neg == 0){

                image_file = argv[argc-2];
                mask_file = argv[argc-1];
                
                //cout << "Image is " << argv[argc-2] << endl;
                //cout << "Mask is " << argv[argc-1] << endl;
                                
        }
    
        try{

                cerr << "Opening image file: " << image_file << endl;

                // Read image given as first argument
                // File type is determined automatically
                vigra::ImageImportInfo info(image_file.c_str());

                // Create RGB images of appropriate size        
                vigra::FRGBImage in(info.width(), info.height());

                // Import the image
                importImage(info, destImage(in));

                // Deal with mask unless -p || -n flags are on
                if (pos == 0 && neg == 0){
                
                        cerr << "Opening mask file:  " << mask_file << endl;

                        // Read mask given as second argument
                        vigra::ImageImportInfo mask_info(mask_file.c_str());
                
                        // Check dimensions are the same
                        if (info.width() != mask_info.width() || info.height() != mask_info.height()){

                                cerr << "Error - image and mask files are different dimension:" << endl;
                                cerr << "Image:\t" << info.width() << " x " << info.height() << endl;
                                cerr << "Mask:\t" << mask_info.width() << " x " << mask_info.height() << endl << endl;
                                return 1;

                        }
                
                        // Create RGB images of appropriate size        
                        vigra::FRGBImage tmp_mask_in(info.width(), info.height());
                        // Import the mask
                        importImage(mask_info, destImage(tmp_mask_in));
                        mask_in = tmp_mask_in;
                
                }
                

                // Max dimension
                const int resize_dimension = 800;
                
                // Re-size to 800 max dimension
                if (info.width() >= resize_dimension || info.height() >= resize_dimension){

                        // Method (0 - pixel repetition, 1 - linear, 2 - spline
                        // Pixel repetition seems to be the fastest
                        int method = 0;
                        double sizefactor;
        
                        if (info.width() >= info.height()){
                                sizefactor = (double)resize_dimension/info.width();
                        }else{
                                sizefactor = (double)resize_dimension/info.height();
                        }
        
                        // calculate new image size
                        int nw = (int)(sizefactor*info.width());
                        int nh = (int)(sizefactor*info.height());
 
                        cerr << "Image dimensions:\t" << info.width() << " x " << info.height() << endl;
                        cerr << "Scaling by:\t\t" << sizefactor << endl;
                        cerr << "New dimensions:\t\t" << nw << " x " << nh << endl;
        
                        // create a RGB image of appropriate size               
                        vigra::FRGBImage out(nw, nh);
                                    
                        switch(method){
              
                                case 0:
                                // resize the image, using a bi-cubic spline algorithms
                                resizeImageNoInterpolation(srcImageRange(in),destImageRange(out));
                                break;
                
                                case 1:
                                // resize the image, using a bi-cubic spline algorithms
                                resizeImageLinearInterpolation(srcImageRange(in),destImageRange(out));
                                break;
                
                                default:
                                // resize the image, using a bi-cubic spline algorithms
                                resizeImageSplineInterpolation(srcImageRange(in),destImageRange(out));
                        }

                        in = out;

                        // Deal with mask unless -p || -n flags are on
                        if (pos == 0 && neg == 0){

                                vigra::FRGBImage mask_out(nw, nh);

                                switch(method){
              
                                        case 0:
                                        // resize the image, using a bi-cubic spline algorithms
                                        resizeImageNoInterpolation(srcImageRange(mask_in),destImageRange(mask_out));
                                        break;
                
                                        case 1:
                                        // resize the image, using a bi-cubic spline algorithms
                                        resizeImageLinearInterpolation(srcImageRange(mask_in),destImageRange(mask_out));
                                        break;
                
                                        default:
                                        // resize the image, using a bi-cubic spline algorithms
                                        resizeImageSplineInterpolation(srcImageRange(mask_in),destImageRange(mask_out));
                                        
                                }
                                
                                mask_in = mask_out;
                        }

                }else{
                        cerr << "Image dimensions:\t" << info.width() << " x " << info.height() << endl;
                }

                cerr << "Total pixels:\t\t" << in.width()*in.height() << endl;
        
                // Convert to LUV colour space
                cerr << "Converting to LUV colourspace..." << endl;             
                
                FRGBImage luv(in.width(),in.height());
                transformImage(srcImageRange(in), destImage(luv), RGBPrime2LuvFunctor<double>() );
                
                // Create grid of fiducial points
                cerr << "Generating fiducial points..." << endl;

                for (int i = gRadius; i < in.height() - gRadius; i += spacing ){
                        for (int j = gRadius; j < in.width() - gRadius; j += spacing ){
                                gNumLocs++;
                        }
                }

                // Create the storage matrix
                gLocations = CreateMatrix( (int)0, gNumLocs, 2);
                gNumLocs = 0;
                for (int i = gRadius; i < in.height() - gRadius; i += spacing ){
                        for (int j = gRadius; j < in.width() - gRadius; j += spacing ){
                        
                                gLocations[gNumLocs][0] = j;
                                gLocations[gNumLocs][1] = i;
                                
                                //cerr << "fPoint " << gNumLocs << ":\t" << i << " " << j << endl;
                                
                                gNumLocs++;
                        }
                }
                cerr << "Total fiducial points: " << gNumLocs << endl;

                // Prepare Gabor API array
                float *frameBuf = new float[in.width()*in.height()];
                float *u_values = new float[in.width()*in.height()];
                float *v_values = new float[in.width()*in.height()];
                
                
                float** pixels = CreateMatrix( (float)0, in.height(), in.width() );

                // Do something with each pixel...
                int counter = 0;
                vigra::FRGBImage::iterator img_iter(luv.begin()),end(luv.end());
                for(; img_iter != end; ++img_iter) {

                        // [0] is L, [1] is U, [2] is V
                        // We only want L for Gabor filtering
                        frameBuf[counter] = (*img_iter)[0];
                        
                        u_values[counter] = (*img_iter)[1];
                        v_values[counter] = (*img_iter)[2];
                
                        //cout << "Pixel " << counter << " - L: " << (*img_iter)[0] << endl;
                        //cout << "Pixel " << counter << " - U: " << (*img_iter)[1] << endl;
                        //cout << "Pixel " << counter << " - V: " << (*img_iter)[2] << endl;
                        counter++;
                }


                int k = 0;
                for (int i = 0; i < in.height(); i++ ){
                        for (int j = 0; j < in.width(); j++ ){
                                pixels[i][j] = frameBuf[k];
                                //cout << i << " " << j << " = " << k << " - " << frameBuf[k] << endl;
                                k++;
                        }
                }

                char basename[] = "gabor_filters/celeste"; 

                float *response = NULL;
                int len = 0;

                response = ProcessChannel( pixels, in.height(), in.width(), response, &len, basename );

                // Deal with mask unless -p || -n flags are on
                if (pos == 0 && neg == 0){
                
                        // Iterator for mask
                        vigra::FRGBImage::iterator tmp_mask_iter(mask_in.begin()),mask_end(mask_in.end());
                        mask_iter = tmp_mask_iter;
                        mask_end = NULL;
                }
                
                cerr << "Generating feature vector by Gabor filtering..." << endl;
                
                // Do something nice with response
                int vector_length = (int)len/gNumLocs;
                for ( int j = 0; j < gNumLocs; j++ ){
                        
                        //cerr << gLocations[j][0] << "," << gLocations[j][1] << endl;
                        
                        int pixel_number = gLocations[j][0] + (in.width() * (gLocations[j][1] - 1)) - 1;

                        // Apply label                  
                        if (pos){
                                cout << "+1 ";
                        }else if (neg){
                                cout << "-1 ";
                        }else{
                                // In the mask, clouds are black, non-clouds are white
                                if (mask_iter[pixel_number][0] != 0 || mask_iter[pixel_number][1] != 0 || mask_iter[pixel_number][2] != 0){
                                        // So this must be non-cloud = -ve example
                                        cout << "-1 ";
                                }else{
                                        // So this must be cloud = +ve example
                                        cout << "+1 ";
                                }                       
                        }

                        int feature = 1;
                        for ( int v = (j * vector_length); v < ((j + 1) * vector_length); v++){
                                cout << feature << ":" << response[v] << " ";
                                feature++;
                        }

                        // Work out average colour - U + V channels
                        
                        float u_sum = 0, v_sum = 0;
                        int pixels_in_region = (gRadius * 2)*(gRadius * 2);
                        
                        // height
                        //cout << endl << "Pixel " << pixel_number << " x:" << gLocations[j][0] << " y:" <<
                        //gLocations[j][1] << " width: " <<  in.width() << endl;
                        
                        for (int t = 1 - gRadius; t <= gRadius; t++){                   
                        
                                int this_y_pixel = pixel_number + (t * in.width());
                                //cout << "t = " << t << " pixel= " << pixel_number << " + " << (t * in.width()) << endl;
                                // width
                                for (int r = 1 - gRadius; r <= gRadius; r++){                   
                        
                                        int this_x_pixel = this_y_pixel + r;
                                        
                                        //if (this_x_pixel == pixel_number){
                                        //      cout << "*" << this_x_pixel << "* ";
                                        //}else{
                                        //      cout << this_x_pixel << " ";
                                        //}
                                        
                                        u_sum += u_values[this_x_pixel];
                                        v_sum += v_values[this_x_pixel];
                                        
                                        //cout << "u = " << u_values[this_x_pixel] << "(" << u_sum << ")" << endl;
                                        //cout << "v = " << v_values[this_x_pixel] << "(" << v_sum << ")" << endl;
                        
                                }
                                //cout << endl;
                        }
                        //cout << endl;

                        float u_ave = (float)u_sum/pixels_in_region;
                        float v_ave = (float)v_sum/pixels_in_region;

                        // Now work out standard deviation for U and V channels
                        u_sum = 0, v_sum = 0;

                        for (int t = 1 - gRadius; t <= gRadius; t++){                   
                                int this_y_pixel = pixel_number + (t * in.width());
                                
                                for (int r = 1 - gRadius; r <= gRadius; r++){                   
                        
                                        int this_x_pixel = this_y_pixel + r;
                                        
                                        u_sum +=  pow(u_values[this_x_pixel]-u_ave,2);
                                        v_sum +=  pow(v_values[this_x_pixel]-v_ave,2);

                                }
                        }
                        
                        
                        float std_u = sqrt(u_sum/(pixels_in_region-1));
                        float std_v = sqrt(v_sum/(pixels_in_region-1)); 
                        
                        //cerr << "U sum = " << u_sum << " U ave = " << u_ave << " U std = " << std_u << endl;
                        //cerr << "V sum = " << v_sum << " V ave = " << v_ave << " V std = " << std_v << endl << endl;
                        
                        cout << feature << ":" << u_ave << " ";
                        feature++;
                        cout << feature << ":" << std_u << " ";
                        feature++;
                        cout << feature << ":" << v_ave << " ";
                        feature++;
                        cout << feature << ":" << std_v << " ";
                        feature++;                              

                        // Add the U/V values for this actual pixel

                        cout << feature << ":" << u_values[pixel_number] << " ";
                        feature++;
                        cout << feature << ":" << v_values[pixel_number] << " ";
                        feature++;                                      
                                
                                
                        //cout << "# FP_" << j << "_Pixel_" << pixel_number << "_Pos_" << gLocations[j][0] << "-" << gLocations[j][1] <<
                        //"_Mask_RGB_" << mask_iter[pixel_number][0] << "," << mask_iter[pixel_number][1] << "," <<
                        //mask_iter[pixel_number][2] << endl;
                        cout << endl;
                }

                delete[] response;
                delete[] frameBuf;
        delete[] u_values;
        delete[] v_values;

                // Export images
                //exportImage(srcImageRange(in), ImageExportInfo("output_images/image_reduced.jpg"));
                //exportImage(srcImageRange(luv), ImageExportInfo("output_images/luv_reduced.jpg"));
                //
                //if (pos == 0 && neg == 0){
                //      exportImage(srcImageRange(mask_in), ImageExportInfo("output_images/mask_reduced.jpg"));
                //}
                        
                //cout << "Done." << endl;
        
        }catch (vigra::StdException & e){
        
                // catch any errors that might have occurred and print their reason
                std::cout << e.what() << std::endl;
                return 1;
                
        }
    
        return 0;
        
}