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

Generated on 27 Sep 2016 for Hugintrunk by  doxygen 1.4.7