Celeste.cpp

Go to the documentation of this file.
00001 /***************************************************************************
00002  *   Copyright (C) 2008 by Tim Nugent                                      *
00003  *   timnugent@gmail.com                                                   *
00004  *                                                                         *
00005  *   This program is free software; you can redistribute it and/or modify  *
00006  *   it under the terms of the GNU General Public License as published by  *
00007  *   the Free Software Foundation; either version 2 of the License, or     *
00008  *   (at your option) any later version.                                   *
00009  *                                                                         *
00010  *   This program is distributed in the hope that it will be useful,       *
00011  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00012  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00013  *   GNU General Public License for more details.                          *
00014  *                                                                         *
00015  *   You should have received a copy of the GNU General Public License     *
00016  *   along with this program; if not, write to the                         *
00017  *   Free Software Foundation, Inc.,                                       *
00018  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
00019  ***************************************************************************/
00020 
00021 #include <iostream>
00022 #include <vigra/stdimage.hxx>
00023 #include <vigra/resizeimage.hxx>
00024 #include <vigra/inspectimage.hxx>
00025 #include <vigra/copyimage.hxx>
00026 #include <vigra/transformimage.hxx>
00027 #include <vigra/initimage.hxx>
00028 #include "vigra/colorconversions.hxx"
00029 #include <sys/types.h>
00030 #include <sys/stat.h>
00031 #include <stdlib.h>
00032 #include <string>
00033 #include <vector>
00034 #include "Gabor.h"
00035 #include "Utilities.h"
00036 #include "CelesteGlobals.h"
00037 #include "Celeste.h"
00038 #include "svm.h"
00039 #include <stdio.h>
00040 
00041 using namespace vigra; 
00042 using namespace std; 
00043 
00044 namespace celeste
00045 {
00046 
00047 typedef vigra::BRGBImage::PixelType RGB;
00048 
00049 // load SVM model
00050 bool loadSVMmodel(struct svm_model*& model, string& model_file)
00051 {
00052     if((model = svm_load_model(model_file.c_str())) == 0)
00053     {
00054         cout << "Couldn't load model file '" << model_file << "'" << endl << endl;
00055         return false;
00056     }
00057     else
00058     {
00059         cout << "Loaded model file:\t" << model_file << endl;
00060         return true;
00061     }
00062 };
00063 
00064 // destroy SVM model
00065 void destroySVMmodel(struct svm_model*& model)
00066 {
00067     svm_destroy_model(model);
00068 };
00069 
00070 // prepare image for use with celeste (downscale, converting in Luv)
00071 void prepareCelesteImage(vigra::UInt16RGBImage& rgb,vigra::UInt16RGBImage& luv, int& resize_dimension, double& sizefactor,bool verbose)
00072 {
00073     // Max dimension
00074     sizefactor = 1;
00075     int nw=rgb.width();
00076     int nh=rgb.height();
00077 
00078     if(nw >= nh)
00079     {
00080         if (resize_dimension >= nw )
00081         {
00082             resize_dimension = nw;
00083         }
00084     }
00085     else
00086     {
00087         if (resize_dimension >= nh)
00088         {
00089             resize_dimension = nh;
00090         }
00091     }
00092     //cout << "Re-size dimenstion:\t" << resize_dimension << endl;
00093     if(verbose)
00094         cout << "Image dimensions:\t" << rgb.width() << " x " << rgb.height() << endl;
00095 
00096     // Re-size to max dimension
00097     vigra::UInt16RGBImage temp;
00098 
00099     if (rgb.width() > resize_dimension || rgb.height() > resize_dimension)
00100     {
00101         if (rgb.width() >= rgb.height())
00102         {
00103             sizefactor = (double)resize_dimension/rgb.width();
00104             // calculate new image size
00105             nw = resize_dimension;
00106             nh = static_cast<int>(0.5 + (sizefactor*rgb.height()));
00107         }
00108         else
00109         {
00110             sizefactor = (double)resize_dimension/rgb.height();
00111             // calculate new image size
00112             nw = static_cast<int>(0.5 + (sizefactor*rgb.width()));
00113             nh = resize_dimension;
00114         }
00115 
00116         if(verbose)
00117         {
00118             cout << "Scaling by:\t\t" << sizefactor << endl;
00119             cout << "New dimensions:\t\t" << nw << " x " << nh << endl;
00120         };
00121 
00122         // create a RGB image of appropriate size
00123         temp.resize(nw,nh);
00124 
00125         // resize the image, using a bi-cubic spline algorithm
00126         vigra::resizeImageNoInterpolation(srcImageRange(rgb),destImageRange(temp));
00127         //exportImage(srcImageRange(out), ImageExportInfo("test.tif").setPixelType("UINT16"));
00128     }
00129     else
00130     {
00131         temp.resize(nw,nh);
00132         vigra::copyImage(srcImageRange(rgb),destImage(temp));
00133     };
00134 
00135     // Convert to Luv colour space
00136     luv.resize(temp.width(),temp.height());
00137     vigra::transformImage(srcImageRange(temp), destImage(luv), RGB2LuvFunctor<double>() );
00138     //exportImage(srcImageRange(luv), ImageExportInfo("test_luv.tif").setPixelType("UINT16"));
00139     temp.resize(0,0);
00140 };
00141 
00142 // converts the given Luv image into arrays for Gabors filtering
00143 void prepareGaborImage(vigra::UInt16RGBImage& luv, float**& pixels)
00144 {
00145     pixels = CreateMatrix( (float)0, luv.height(), luv.width() );
00146     // Prepare framebuf for Gabor API, we need only L channel
00147     for (int i = 0; i < luv.height(); i++ )
00148     {
00149         for (int j = 0; j < luv.width(); j++ )
00150         {
00151             pixels[i][j] = luv(j,i)[0];
00152             //cout << i << " " << j << " = " << k << " - " << frameBuf[k] << endl;
00153         }
00154     }
00155 };
00156 
00157 //classify the points with SVM
00158 vector<double> classifySVM(struct svm_model* model, int gNumLocs,int**& gLocations,int width,int height,int vector_length, float*& response,int gRadius,vigra::UInt16RGBImage& luv,bool needsMoreIndex=false)
00159 {
00160     std::vector<double> svm_response;
00161     // Integers and containers for libsvm
00162     int max_nr_attr = 56;
00163     int nr_class=svm_get_nr_class(model);
00164     struct svm_node *gabor_responses = (struct svm_node *) malloc(max_nr_attr*sizeof(struct svm_node));
00165     double *prob_estimates = (double *) malloc(nr_class*sizeof(double));
00166 
00167     for (int j = 0; j < gNumLocs; j++)
00168     {
00169         unsigned int feature = 1;
00170 
00171         if(needsMoreIndex)
00172         {
00173             if(j >= max_nr_attr - 1)
00174             {
00175                 max_nr_attr *= 2;
00176                 gabor_responses = (struct svm_node *) realloc(gabor_responses,max_nr_attr*sizeof(struct svm_node));
00177             }
00178         };
00179 
00180         for ( int v = (j * vector_length); v < ((j + 1) * vector_length); v++)
00181         {
00182             gabor_responses[feature-1].index = feature;
00183             gabor_responses[feature-1].value = response[v];                                     
00184             //cout << feature << ":" << response[v] << " ";
00185             feature++;
00186         }
00187 
00188         // Work out average colour and variance
00189         vigra::FindAverageAndVariance<vigra::UInt16RGBImage::PixelType> average;
00190         vigra::inspectImage(srcIterRange(
00191             luv.upperLeft()+vigra::Diff2D(gLocations[j][0]-gRadius,gLocations[j][1]-gRadius),
00192             luv.upperLeft()+vigra::Diff2D(gLocations[j][0]+gRadius,gLocations[j][1]+gRadius)
00193             ),average);
00194         // Add these colour features to feature vector                                                  
00195 
00196         gabor_responses[feature-1].index = feature;
00197         gabor_responses[feature-1].value = average.average()[1];
00198         //cout << feature << ":" << u_ave << " ";
00199         feature++;
00200         gabor_responses[feature-1].index = feature;
00201         gabor_responses[feature-1].value = sqrt(average.variance()[1]);
00202         //cout << feature << ":" << std_u << " ";
00203         feature++;
00204         gabor_responses[feature-1].index = feature;
00205         gabor_responses[feature-1].value = average.average()[2];
00206         //cout << feature << ":" << v_ave << " ";
00207         feature++;
00208         gabor_responses[feature-1].index = feature;
00209         gabor_responses[feature-1].value = sqrt(average.variance()[2]);
00210         //cout << feature << ":" << std_v << " ";
00211         feature++;
00212         gabor_responses[feature-1].index = feature;
00213         gabor_responses[feature-1].value = luv(gLocations[j][0],gLocations[j][1])[1];
00214         //cout << feature << ":" << u_values[pixel_number] << " ";
00215         feature++;
00216         gabor_responses[feature-1].index = feature;
00217         gabor_responses[feature-1].value = luv(gLocations[j][0],gLocations[j][1])[2];
00218         //cout << feature << ":" << v_values[pixel_number] << " " << endl;
00219         gabor_responses[feature].index = -1;
00220 
00221         svm_predict_probability(model,gabor_responses,prob_estimates);  
00222         svm_response.push_back(prob_estimates[0]);
00223     }
00224     // Free up libsvm stuff
00225     free(gabor_responses);
00226     free(prob_estimates);
00227     return svm_response;
00228 };
00229 
00230 // create a grid of control points for creating masks
00231 void createGrid(int& gNumLocs,int**& gLocations,int gRadius,int width, int height)
00232 {
00233     int spacing=(gRadius*2)+1;
00234     for (int i = gRadius; i < height - gRadius; i += spacing )
00235     {
00236         for (int j = gRadius; j < width - gRadius; j += spacing )
00237         {
00238             gNumLocs++;
00239         }
00240         // Add extra FP at the end of each row in case width % gRadius
00241         gNumLocs++;
00242     }
00243 
00244     // Add extra FP at the end of each row in case nh % gRadius 
00245     for (int j = gRadius; j < width - gRadius; j += spacing )
00246     {
00247         gNumLocs++;
00248     }
00249 
00250     // Create the storage matrix
00251     gLocations = CreateMatrix( (int)0, gNumLocs, 2);
00252     gNumLocs = 0;
00253     for (int i = gRadius; i < height - gRadius; i += spacing )
00254     {
00255         for (int j = gRadius; j < width - gRadius; j += spacing )
00256         {
00257             gLocations[gNumLocs][0] = j;
00258             gLocations[gNumLocs][1] = i;
00259             //cout << "fPoint " << gNumLocs << ":\t" << i << " " << j << endl;
00260             gNumLocs++;
00261         }
00262 
00263         // Add extra FP at the end of each row in case width % spacing
00264         if (width % spacing)
00265         {
00266             gLocations[gNumLocs][0] = width - gRadius - 1;
00267             gLocations[gNumLocs][1] = i;
00268             //cout << "efPoint " << gNumLocs << ":\t" << i << " " << nw - gRadius - 1 << endl;
00269             gNumLocs++;
00270         }
00271     }
00272 
00273     // Add extra FP at the end of each row in case height % spacing
00274     if (height % spacing)
00275     {
00276         for (int j = gRadius; j < width - gRadius; j += spacing )
00277         {
00278             gLocations[gNumLocs][0] = j;
00279             gLocations[gNumLocs][1] = height - gRadius - 1;
00280             //cout << "efPoint " << gNumLocs << ":\t" << nh - gRadius - 1 << " " << j << endl;
00281             gNumLocs++;
00282         }
00283     }
00284 };
00285 
00286 //generates the celeste mask on base of given responses and locations
00287 void generateMask(vigra::BImage& mask,int& gNumLocs, int**& gLocations,std::vector<double> svm_responses,int gRadius, double threshold)
00288 {
00289     for ( int j = 0; j < gNumLocs; j++ )
00290     {
00291         if (svm_responses[j] >= threshold)
00292         {
00293             unsigned int sub_x0 = gLocations[j][0] - gRadius;
00294             unsigned int sub_y0 = gLocations[j][1] - gRadius;
00295             unsigned int sub_x1 = gLocations[j][0] + gRadius + 1;
00296             unsigned int sub_y1 = gLocations[j][1] + gRadius + 1;
00297             //cout << sub_x0 << ","<< sub_y0 << " - " << sub_x1 << "," << sub_y1 << endl;
00298 
00299             // Set region to black
00300             vigra::initImage(srcIterRange(mask.upperLeft() + vigra::Diff2D(sub_x0, sub_y0),
00301                         mask.upperLeft() + vigra::Diff2D(sub_x1, sub_y1)), 0);                          
00302         }
00303         else
00304         {
00305             //cout << "Non-cloud\t(score " << prob_estimates[0] << " <= " << threshold << ")" << endl;  
00306         }
00307     }
00308 };
00309 
00310 vigra::BImage getCelesteMask(struct svm_model* model, vigra::UInt16RGBImage& input, int radius, float threshold, int resize_dimension,bool adaptThreshold,bool verbose)
00311 {
00312     vigra::UInt16RGBImage luv;
00313     double sizefactor=1.0;
00314     prepareCelesteImage(input,luv,resize_dimension,sizefactor,verbose);
00315 
00316     // Prepare Gabor API array
00317     float** pixels=NULL;
00318     prepareGaborImage(luv,pixels);
00319 
00320     int** gLocations = NULL;
00321     int gNumLocs = 0;
00322 
00323     // Create grid of fiducial points
00324     createGrid(gNumLocs,gLocations,radius,luv.width(),luv.height());
00325 
00326     int len = 0;
00327     float* mask_response=NULL;
00328     mask_response = ProcessChannel(pixels,luv.width(),luv.height(),gNumLocs,gLocations,radius,mask_response,&len);
00329     // Turn the response into SVM vector, and add colour features
00330     vector<double> svm_responses=classifySVM(model,gNumLocs,gLocations,luv.width(),luv.height(),(int)len/gNumLocs,mask_response,radius,luv);
00331     delete [] mask_response;
00332 
00333     if(adaptThreshold)
00334     {
00335         double minVal=1;
00336         for(unsigned int i=0;i<svm_responses.size();i++)
00337         {
00338             if(svm_responses[i]<minVal)
00339             {
00340                 minVal=svm_responses[i];
00341             };
00342         };
00343         if(threshold<minVal)
00344         {
00345             threshold=min(minVal+0.1,1.0);
00346         };
00347     };
00348     // Create mask of same dimensions
00349     vigra::BImage mask_out(luv.width(), luv.height(),255);
00350     generateMask(mask_out,gNumLocs,gLocations,svm_responses,radius,threshold);
00351     // Re-size mask to match original image
00352     vigra::BImage mask_resize(input.width(),input.height());                       
00353     resizeImageNoInterpolation(srcImageRange(mask_out),destImageRange(mask_resize));            
00354     DisposeMatrix(pixels,luv.height());
00355     DisposeMatrix(gLocations,gNumLocs);
00356     mask_out.resize(0,0);
00357     return mask_resize;
00358 };
00359 
00360 HuginBase::UIntSet getCelesteControlPoints(struct svm_model* model, vigra::UInt16RGBImage& input, HuginBase::CPointVector cps, int radius, float threshold, int resize_dimension,bool verbose)
00361 {
00362     HuginBase::UIntSet cloudCP;
00363     vigra::UInt16RGBImage luv;
00364     double sizefactor=1.0;
00365     prepareCelesteImage(input,luv,resize_dimension,sizefactor,verbose);
00366 
00367     // Prepare Gabor API array
00368     float** pixels=NULL;
00369     prepareGaborImage(luv,pixels);
00370 
00371     int gNumLocs = cps.size();
00372     int** gLocations = CreateMatrix( (int)0, gNumLocs, 2);
00373     for(unsigned int j=0;j<cps.size();j++)
00374     {
00375         HuginBase::ControlPoint cp=cps[j].second;
00376         gLocations[j][0] = int(cp.x1 * sizefactor);
00377         gLocations[j][1] = int(cp.y1 * sizefactor);
00378         // Move CPs to border if the filter radius is out of bounds
00379         if (gLocations[j][0] <= radius)
00380         {
00381             gLocations[j][0] = radius + 1;
00382         }
00383         if (gLocations[j][1] <= radius)
00384         {
00385             gLocations[j][1] = radius + 1;
00386         }
00387         if (gLocations[j][0] >= luv.width() - radius)
00388         {
00389             gLocations[j][0] = luv.width() - radius - 1;
00390         }
00391         if (gLocations[j][1] >= luv.height() - radius)
00392         {
00393             gLocations[j][1] = luv.height() - radius - 1;
00394         }
00395     };
00396 
00397     int len = 0;
00398     float* response=NULL;
00399     response = ProcessChannel(pixels,luv.width(),luv.height(),gNumLocs,gLocations,radius,response,&len);
00400     // Turn the response into SVM vector, and add colour features
00401     vector<double> svm_responses=classifySVM(model,gNumLocs,gLocations,luv.width(),luv.height(),(int)len/gNumLocs,response,radius,luv);
00402     delete [] response;
00403 
00404     for(unsigned int i=0;i<svm_responses.size();i++)
00405     {
00406         if(svm_responses[i]>=threshold)
00407         {
00408             cloudCP.insert(cps[i].first);
00409         };
00410     };
00411     DisposeMatrix(pixels,luv.height());
00412     DisposeMatrix(gLocations,gNumLocs);
00413 
00414     return cloudCP;
00415 };
00416 
00417 } // end of namespace

Generated on 23 Oct 2014 for Hugintrunk by  doxygen 1.4.7