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

Generated on 2 Sep 2015 for Hugintrunk by  doxygen 1.4.7