geo/tcChannelFile.cpp

   1 /***************************************************************************
   2  *   This file is part of Tecorrec.                                        *
   3  *   Copyright 2008 James Hogan <james@albanarts.com>                      *
   4  *                                                                         *
   5  *   Tecorrec is free software: you can redistribute it and/or modify      *
   6  *   it under the terms of the GNU General Public License as published by  *
   7  *   the Free Software Foundation, either version 2 of the License, or     *
   8  *   (at your option) any later version.                                   *
   9  *                                                                         *
  10  *   Tecorrec is distributed in the hope that it will be useful,           *
  11  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  12  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
  13  *   GNU General Public License for more details.                          *
  14  *                                                                         *
  15  *   You should have received a copy of the GNU General Public License     *
  16  *   along with Tecorrec.  If not, write to the Free Software Foundation,  *
  17  *   Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.   *
  18  ***************************************************************************/
  19
  20 /**
  21  * @file tcChannelFile.cpp
  22  * @brief An image channel from an image file.
  23  */
  24
  25 #include "tcChannelFile.h"
  26 #include "tcPixelData.h"
  27
  28 #include <gdal_priv.h>
  29 #include <ogr_spatialref.h>
  30
  31 #include <QtDebug>
  32
  33 #include <GL/glu.h>
  34
  35 /*
  36  * Constructors + destructor
  37  */
  38
  39 /// Primary constructor.
  40 tcChannelFile::tcChannelFile(const QString& filename, const QString& name, const QString& description)
  41 : tcChannel(name, description)
  42 , m_width(0)
  43 , m_height(0)
  44 , m_dataset(0)
  45 , m_thumbnail(0)
  46 , m_projToGeo(0)
  47 , m_geoToProj(0)
  48 , m_pixelsToProj()
  49 , m_projToPixels()
  50 , m_texToGeo()
  51 , m_geoToTex()
  52 {
  53   // Open the file
  54   /// @todo this probably isn't the best place to set it up
  55   GDALAllRegister();
  56
  57   m_dataset = (GDALDataset*)GDALOpen((const char*)filename.toAscii(), GA_ReadOnly);
  58   if (0 != m_dataset)
  59   {
  60     GDALRasterBand* band = m_dataset->GetRasterBand(1);
  61     Q_ASSERT(0 != band);
  62     m_width = band->GetXSize();
  63     m_height = band->GetYSize();
  64
  65     OGRCoordinateTransformation* result = 0;
  66     if (0 != m_dataset->GetProjectionRef())
  67     {
  68       OGRSpatialReference srs(m_dataset->GetProjectionRef());
  69       OGRSpatialReference* geog = srs.CloneGeogCS();
  70       m_projToGeo = OGRCreateCoordinateTransformation(&srs, geog);
  71       m_geoToProj = OGRCreateCoordinateTransformation(geog, &srs);
  72       delete geog;
  73       m_dataset->GetGeoTransform(m_pixelsToProj);
  74       GDALInvGeoTransform(m_pixelsToProj, m_projToPixels);
  75
  76       // Transform each corner
  77       maths::Vector<2,double> corners[4] = {
  78         maths::Vector<2,double>( 0.0,       0.0        ),
  79         maths::Vector<2,double>( 0.0,       m_height-1 ),
  80         maths::Vector<2,double>( m_width-1, 0.0        ),
  81         maths::Vector<2,double>( m_width-1, m_height-1 )
  82       };
  83       maths::Vector<2,double> min, max;
  84       for (int i = 0; i < 4; ++i)
  85       {
  86         maths::Vector<2,double> lonlat = m_pixelsToProj * corners[i];
  87         m_projToGeo->Transform(1, &lonlat[0], &lonlat[1]);
  88         // Keep track of extents in geographical space
  89         if (i > 0)
  90         {
  91           for (int j = 0; j < 2; ++j)
  92           {
  93             if (lonlat[j] < min[j])
  94             {
  95               min[i] = lonlat[j];
  96             }
  97             else if (lonlat[j] > max[j])
  98             {
  99               max[j] = lonlat[j];
 100             }
 101           }
 102         }
 103         else
 104         {
 105           for (int j = 0; j < 2; ++j)
 106           {
 107             min[j] = max[j] = lonlat[j];
 108           }
 109         }
 110       }
 111       // Convert extents into radians
 112       min *= M_PI/180;
 113       max *= M_PI/180;
 114       // Texture is mirrored in y
 115       qSwap(min[1], max[1]);
 116       // Get affine transformation matrices
 117       m_texToGeo = tcAffineTransform2<double>(min, max-min);
 118       m_geoToTex = m_texToGeo.inverse();
 119     }
 120   }
 121 }
 122
 123 /// Destructor.
 124 tcChannelFile::~tcChannelFile()
 125 {
 126   delete m_dataset;
 127   if (0 != m_thumbnail)
 128   {
 129     glDeleteTextures(1, &m_thumbnail);
 130   }
 131
 132   delete m_projToGeo;
 133   delete m_geoToProj;
 134 }
 135
 136 /*
 137  * Accessors
 138  */
 139
 140 /// Load the transformation from projection to texture.
 141 const tcAffineTransform2<double> tcChannelFile::geoToTex() const
 142 {
 143   return m_geoToTex;
 144 }
 145
 146 /// Load the transformation from projection to texture.
 147 const tcAffineTransform2<double> tcChannelFile::texToGeo() const
 148 {
 149   return m_texToGeo;
 150 }
 151
 152 /*
 153  * Main image interface
 154  */
 155
 156 // Reimplemented
 157 GLuint tcChannelFile::thumbnailTexture()
 158 {
 159   if (0 == m_thumbnail && 0 != m_dataset)
 160   {
 161     GDALRasterBand* band = m_dataset->GetRasterBand(1);
 162     Q_ASSERT(0 != band);
 163
 164     // Get an overview
 165     GDALRasterBand* thumbnailBand = 0;
 166     int numOverviews = band->GetOverviewCount();
 167     if (0 == numOverviews)
 168     {
 169
 170       // we want to reduce to around 1000000 pixels
 171       /// @todo Do we need to be careful of overflow here?
 172       int pixels = m_width*m_height;
 173       int overview = (int)sqrt(pixels >> 22);
 174       if (overview > 1)
 175       {
 176         if (CE_None != m_dataset->BuildOverviews("NEAREST", 1, &overview, 0, 0,
 177               GDALDummyProgress, 0))
 178         {
 179           thumbnailBand = band->GetOverview(0);
 180         }
 181       }
 182     }
 183     else
 184     {
 185       /// @todo Pick an appropriate overview if more than one
 186       thumbnailBand = band->GetOverview(0);
 187     }
 188
 189     if (thumbnailBand)
 190     {
 191       // Turn into a GL texture
 192       int width  = thumbnailBand->GetXSize();
 193       int height = thumbnailBand->GetYSize();
 194       int tnWidth = 2048;
 195       int tnHeight = 2048;
 196       tcPixelData<GLubyte> data(tnWidth, tnHeight);
 197       tcPixelData<GLubyte> geoData(tnWidth, tnHeight);
 198       if (CE_None == thumbnailBand->RasterIO(GF_Read, 0, 0, width, height,
 199                                              data.buffer(), tnWidth, tnHeight, GDT_Byte, 0, 0))
 200       {
 201         int index = 0;
 202         for (int j = 0; j < tnHeight; ++j)
 203         {
 204           for (int i = 0; i < tnWidth; ++i)
 205           {
 206             // Geo coordinates
 207             maths::Vector<2,double> tex((double)i/(tnWidth-1), (double)j/(tnHeight-1));
 208             maths::Vector<2,double> geo = (m_texToGeo * tex) * 180.0f/M_PI;
 209             // Local coordinates
 210             m_geoToProj->Transform(1, &geo[0], &geo[1]);
 211             // To pixel coordinates
 212             maths::Vector<2,double> pix = m_projToPixels*(maths::Vector<2,double>)geo;
 213             geoData.buffer()[index] = (GLubyte)data.tcPixelData<GLubyte>::sampleFloat(pix[0]/(m_width-1), pix[1]/(m_height-1));
 214             ++index;
 215           }
 216         }
 217         glGenTextures(1, &m_thumbnail);
 218         glBindTexture(GL_TEXTURE_2D, m_thumbnail);
 219         if (0 == gluBuild2DMipmaps(GL_TEXTURE_2D, 1, tnWidth, tnHeight, GL_LUMINANCE, GL_UNSIGNED_BYTE, geoData.buffer()))
 220         {
 221           glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_NEAREST);
 222           glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
 223           glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
 224           glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
 225           glBindTexture(GL_TEXTURE_2D, 0);
 226         }
 227         else
 228         {
 229           glBindTexture(GL_TEXTURE_2D, 0);
 230           glDeleteTextures(1, &m_thumbnail);
 231           m_thumbnail = 0;
 232         }
 233       }
 234     }
 235   }
 236   return m_thumbnail;
 237 }
 238
 239 /*
 240  * Interface for derived class to implement
 241  */
 242
 243 void tcChannelFile::roundPortion(double* x1, double* y1, double* x2, double* y2)
 244 {
 245   // We're resampling the texture anyway so no need to adjust range
 246 }
 247
 248 tcAbstractPixelData* tcChannelFile::loadPortion(double x1, double y1, double x2, double y2)
 249 {
 250   tcPixelData<GLubyte>* buffer = 0;
 251   if (0 != m_dataset)
 252   {
 253     GDALRasterBand* band = m_dataset->GetRasterBand(1);
 254     Q_ASSERT(0 != band);
 255
 256     // Transform each corner
 257     tcAffineTransform2<double> portionToTex(maths::Vector<2,double>(x1,y1), maths::Vector<2,double>(x2-x1, y2-y1));
 258     tcAffineTransform2<double> portionToGeo = m_texToGeo * portionToTex;
 259     maths::Vector<2,double> corners[4] = {
 260       maths::Vector<2,double>( x1, y1 ),
 261       maths::Vector<2,double>( x1, y2 ),
 262       maths::Vector<2,double>( x2, y1 ),
 263       maths::Vector<2,double>( x2, y2 )
 264     };
 265     maths::Vector<2,double> min, max;
 266     for (int i = 0; i < 4; ++i)
 267     {
 268       maths::Vector<2,double> lonlat = m_texToGeo * corners[i];
 269       lonlat *= 180.0f/M_PI;
 270       m_geoToProj->Transform(1, &lonlat[0], &lonlat[1]);
 271       maths::Vector<2,double> pix = m_projToPixels * lonlat;
 272       // Keep track of extents in pixel space
 273       if (i > 0)
 274       {
 275         for (int j = 0; j < 2; ++j)
 276         {
 277           if (pix[j] < min[j])
 278           {
 279             min[j] = pix[j];
 280           }
 281           else if (pix[j] > max[j])
 282           {
 283             max[j] = pix[j];
 284           }
 285         }
 286       }
 287       else
 288       {
 289         min = max = pix;
 290       }
 291     }
 292     // Round min down and max up
 293     min[0] = floor(min[0]);
 294     min[1] = floor(min[1]);
 295     max[0] = ceil(max[0]);
 296     max[1] = ceil(max[1]);
 297     // Keep within limits of original texture
 298     if (min[0] < 0.0)        { min[0] = 0.0;        }
 299     if (min[1] < 0.0)        { min[1] = 0.0;        }
 300     if (max[0] > m_width-1)  { max[0] = m_width-1;  }
 301     if (max[1] > m_height-1) { max[1] = m_height-1; }
 302     // Get affine transformation matrices
 303     tcAffineTransform2<double> outOfWindow = tcAffineTransform2<double>(min, max-min);
 304     tcAffineTransform2<double> inToWindow = outOfWindow.inverse();
 305
 306     int s1 = (int)min[0];
 307     int t1 = (int)min[1];
 308     int partWidth  = (int)max[0] - s1;
 309     int partHeight = (int)max[1] - t1;
 310
 311     int resw = partWidth;
 312     int resh = partHeight;
 313     // ensure resw and resh are powers of 2
 314     for (int i = 1; i < 20; ++i)
 315     {
 316       if (0 == (resw & ~((1<<i)-1)))
 317       {
 318         resw = 1<<i;
 319         break;
 320       }
 321     }
 322     for (int i = 1; i < 20; ++i)
 323     {
 324       if (0 == (resh & ~((1<<i)-1)))
 325       {
 326         resh = 1<<i;
 327         break;
 328       }
 329     }
 330     qDebug() << "Extract resolution from" << partWidth << "x" << partHeight << "to" << resw << "x" << resh;
 331     tcPixelData<GLubyte> data(resw, resh);
 332     if (CE_None == band->RasterIO(GF_Read, s1, t1, partWidth, partHeight,
 333           data.buffer(), resw, resh, GDT_Byte, 0, 0))
 334     {
 335
 336       buffer = new tcPixelData<GLubyte>(resw, resh);
 337       int index = 0;
 338       for (int j = 0; j < resh; ++j)
 339       {
 340         for (int i = 0; i < resw; ++i)
 341         {
 342           // Geo coordinates
 343           maths::Vector<2,double> portion((double)i/(resw-1), (double)j/(resh-1));
 344           maths::Vector<2,double> tex = portionToTex * portion;
 345           maths::Vector<2,double> geo = (m_texToGeo * tex) * 180.0f/M_PI;
 346           // Local coordinates
 347           m_geoToProj->Transform(1, &geo[0], &geo[1]);
 348           // To pixel coordinates
 349           maths::Vector<2,double> pix = m_projToPixels * geo;
 350           maths::Vector<2,double> window = inToWindow * pix;
 351           buffer->buffer()[index] = (GLubyte)data.tcPixelData<GLubyte>::sampleFloat(window[0], window[1]);
 352           //buffer->buffer()[index] = 255.0f*portion[0];//data.buffer()[index];
 353           ++index;
 354         }
 355       }
 356     }
 357   }
 358   return buffer;
 359 }