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, float radianceOffset, float radianceBias)
  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   setRadianceTransform(radianceOffset, radianceBias);
  54
  55   // Open the file
  56   /// @todo this probably isn't the best place to set it up
  57   GDALAllRegister();
  58
  59   m_dataset = (GDALDataset*)GDALOpen((const char*)filename.toAscii(), GA_ReadOnly);
  60   if (0 != m_dataset)
  61   {
  62     GDALRasterBand* band = m_dataset->GetRasterBand(1);
  63     Q_ASSERT(0 != band);
  64     m_width = band->GetXSize();
  65     m_height = band->GetYSize();
  66
  67     OGRCoordinateTransformation* result = 0;
  68     if (0 != m_dataset->GetProjectionRef())
  69     {
  70       OGRSpatialReference srs(m_dataset->GetProjectionRef());
  71       OGRSpatialReference* geog = srs.CloneGeogCS();
  72       m_projToGeo = OGRCreateCoordinateTransformation(&srs, geog);
  73       m_geoToProj = OGRCreateCoordinateTransformation(geog, &srs);
  74       delete geog;
  75       m_dataset->GetGeoTransform(m_pixelsToProj);
  76       GDALInvGeoTransform(m_pixelsToProj, m_projToPixels);
  77
  78       // Transform each corner
  79       maths::Vector<2,double> corners[4] = {
  80         maths::Vector<2,double>( 0.0,       0.0        ),
  81         maths::Vector<2,double>( 0.0,       m_height-1 ),
  82         maths::Vector<2,double>( m_width-1, 0.0        ),
  83         maths::Vector<2,double>( m_width-1, m_height-1 )
  84       };
  85       maths::Vector<2,double> min, max;
  86       for (int i = 0; i < 4; ++i)
  87       {
  88         maths::Vector<2,double> lonlat = m_pixelsToProj * corners[i];
  89         m_projToGeo->Transform(1, &lonlat[0], &lonlat[1]);
  90         // Keep track of extents in geographical space
  91         if (i > 0)
  92         {
  93           for (int j = 0; j < 2; ++j)
  94           {
  95             if (lonlat[j] < min[j])
  96             {
  97               min[i] = lonlat[j];
  98             }
  99             else if (lonlat[j] > max[j])
 100             {
 101               max[j] = lonlat[j];
 102             }
 103           }
 104         }
 105         else
 106         {
 107           for (int j = 0; j < 2; ++j)
 108           {
 109             min[j] = max[j] = lonlat[j];
 110           }
 111         }
 112       }
 113       // Convert extents into radians
 114       min *= M_PI/180;
 115       max *= M_PI/180;
 116       // Texture is mirrored in y
 117       qSwap(min[1], max[1]);
 118       // Get affine transformation matrices
 119       m_texToGeo = tcAffineTransform2<double>(min, max-min);
 120       m_geoToTex = m_texToGeo.inverse();
 121     }
 122   }
 123 }
 124
 125 /// Destructor.
 126 tcChannelFile::~tcChannelFile()
 127 {
 128   delete m_dataset;
 129   if (0 != m_thumbnail)
 130   {
 131     glDeleteTextures(1, &m_thumbnail);
 132   }
 133
 134   delete m_projToGeo;
 135   delete m_geoToProj;
 136 }
 137
 138 /*
 139  * Accessors
 140  */
 141
 142 /// Load the transformation from projection to texture.
 143 const tcAffineTransform2<double> tcChannelFile::geoToTex() const
 144 {
 145   return m_geoToTex;
 146 }
 147
 148 /// Load the transformation from projection to texture.
 149 const tcAffineTransform2<double> tcChannelFile::texToGeo() const
 150 {
 151   return m_texToGeo;
 152 }
 153
 154 /*
 155  * Main image interface
 156  */
 157
 158 // Reimplemented
 159 GLuint tcChannelFile::thumbnailTexture()
 160 {
 161   if (0 == m_thumbnail && 0 != m_dataset)
 162   {
 163     GDALRasterBand* band = m_dataset->GetRasterBand(1);
 164     Q_ASSERT(0 != band);
 165
 166     // Get an overview
 167     GDALRasterBand* thumbnailBand = 0;
 168     int numOverviews = band->GetOverviewCount();
 169     if (0 == numOverviews)
 170     {
 171
 172       // we want to reduce to around 1000000 pixels
 173       /// @todo Do we need to be careful of overflow here?
 174       int pixels = m_width*m_height;
 175       int overview = (int)sqrt(pixels >> 22);
 176       if (overview > 1)
 177       {
 178         if (CE_None != m_dataset->BuildOverviews("NEAREST", 1, &overview, 0, 0,
 179               GDALDummyProgress, 0))
 180         {
 181           thumbnailBand = band->GetOverview(0);
 182         }
 183       }
 184     }
 185     else
 186     {
 187       /// @todo Pick an appropriate overview if more than one
 188       thumbnailBand = band->GetOverview(0);
 189     }
 190
 191     if (thumbnailBand)
 192     {
 193       // Turn into a GL texture
 194       int width  = thumbnailBand->GetXSize();
 195       int height = thumbnailBand->GetYSize();
 196       int tnWidth = 2048;
 197       int tnHeight = 2048;
 198       tcPixelData<GLubyte> data(tnWidth, tnHeight);
 199       tcPixelData<GLubyte> geoData(tnWidth, tnHeight);
 200       if (CE_None == thumbnailBand->RasterIO(GF_Read, 0, 0, width, height,
 201                                              data.buffer(), tnWidth, tnHeight, GDT_Byte, 0, 0))
 202       {
 203         int index = 0;
 204         for (int j = 0; j < tnHeight; ++j)
 205         {
 206           for (int i = 0; i < tnWidth; ++i)
 207           {
 208             // Geo coordinates
 209             maths::Vector<2,double> tex((double)i/(tnWidth-1), (double)j/(tnHeight-1));
 210             maths::Vector<2,double> geo = (m_texToGeo * tex) * 180.0f/M_PI;
 211             // Local coordinates
 212             m_geoToProj->Transform(1, &geo[0], &geo[1]);
 213             // To pixel coordinates
 214             maths::Vector<2,double> pix = m_projToPixels*(maths::Vector<2,double>)geo;
 215             geoData.buffer()[index] = (GLubyte)data.tcPixelData<GLubyte>::sampleFloat(pix[0]/(m_width-1), pix[1]/(m_height-1));
 216             ++index;
 217           }
 218         }
 219         glGenTextures(1, &m_thumbnail);
 220         glBindTexture(GL_TEXTURE_2D, m_thumbnail);
 221         if (0 == gluBuild2DMipmaps(GL_TEXTURE_2D, 1, tnWidth, tnHeight, GL_LUMINANCE, GL_UNSIGNED_BYTE, geoData.buffer()))
 222         {
 223           glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_NEAREST);
 224           glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
 225           glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
 226           glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
 227           glBindTexture(GL_TEXTURE_2D, 0);
 228         }
 229         else
 230         {
 231           glBindTexture(GL_TEXTURE_2D, 0);
 232           glDeleteTextures(1, &m_thumbnail);
 233           m_thumbnail = 0;
 234         }
 235       }
 236     }
 237   }
 238   return m_thumbnail;
 239 }
 240
 241 /*
 242  * Interface for derived class to implement
 243  */
 244
 245 void tcChannelFile::roundPortion(double* x1, double* y1, double* x2, double* y2)
 246 {
 247   // We're resampling the texture anyway so no need to adjust range
 248 }
 249
 250 tcAbstractPixelData* tcChannelFile::loadPortion(double x1, double y1, double x2, double y2, bool changed)
 251 {
 252   tcPixelData<GLubyte>* buffer = 0;
 253   if (0 != m_dataset)
 254   {
 255     GDALRasterBand* band = m_dataset->GetRasterBand(1);
 256     Q_ASSERT(0 != band);
 257
 258     // Transform each corner
 259     tcAffineTransform2<double> portionToTex(maths::Vector<2,double>(x1,y1), maths::Vector<2,double>(x2-x1, y2-y1));
 260     tcAffineTransform2<double> portionToGeo = m_texToGeo * portionToTex;
 261     maths::Vector<2,double> corners[4] = {
 262       maths::Vector<2,double>( x1, y1 ),
 263       maths::Vector<2,double>( x1, y2 ),
 264       maths::Vector<2,double>( x2, y1 ),
 265       maths::Vector<2,double>( x2, y2 )
 266     };
 267     maths::Vector<2,double> min, max;
 268     for (int i = 0; i < 4; ++i)
 269     {
 270       maths::Vector<2,double> lonlat = m_texToGeo * corners[i];
 271       lonlat *= 180.0f/M_PI;
 272       m_geoToProj->Transform(1, &lonlat[0], &lonlat[1]);
 273       maths::Vector<2,double> pix = m_projToPixels * lonlat;
 274       // Keep track of extents in pixel space
 275       if (i > 0)
 276       {
 277         for (int j = 0; j < 2; ++j)
 278         {
 279           if (pix[j] < min[j])
 280           {
 281             min[j] = pix[j];
 282           }
 283           else if (pix[j] > max[j])
 284           {
 285             max[j] = pix[j];
 286           }
 287         }
 288       }
 289       else
 290       {
 291         min = max = pix;
 292       }
 293     }
 294     // Round min down and max up
 295     min[0] = floor(min[0]);
 296     min[1] = floor(min[1]);
 297     max[0] = ceil(max[0]);
 298     max[1] = ceil(max[1]);
 299     // Keep within limits of original texture
 300     if (min[0] < 0.0)        { min[0] = 0.0;        }
 301     if (min[1] < 0.0)        { min[1] = 0.0;        }
 302     if (max[0] > m_width-1)  { max[0] = m_width-1;  }
 303     if (max[1] > m_height-1) { max[1] = m_height-1; }
 304     // Get affine transformation matrices
 305     tcAffineTransform2<double> outOfWindow = tcAffineTransform2<double>(min, max-min);
 306     tcAffineTransform2<double> inToWindow = outOfWindow.inverse();
 307
 308     int s1 = (int)min[0];
 309     int t1 = (int)min[1];
 310     int partWidth  = (int)max[0] - s1;
 311     int partHeight = (int)max[1] - t1;
 312
 313     int resw = partWidth;
 314     int resh = partHeight;
 315     // ensure resw and resh are powers of 2
 316     for (int i = 1; i < 20; ++i)
 317     {
 318       if (0 == (resw & ~((1<<i)-1)))
 319       {
 320         resw = 1<<i;
 321         break;
 322       }
 323     }
 324     for (int i = 1; i < 20; ++i)
 325     {
 326       if (0 == (resh & ~((1<<i)-1)))
 327       {
 328         resh = 1<<i;
 329         break;
 330       }
 331     }
 332     qDebug() << "Extract resolution from" << partWidth << "x" << partHeight << "to" << resw << "x" << resh;
 333     tcPixelData<GLubyte> data(resw, resh);
 334     if (CE_None == band->RasterIO(GF_Read, s1, t1, partWidth, partHeight,
 335           data.buffer(), resw, resh, GDT_Byte, 0, 0))
 336     {
 337
 338       buffer = new tcPixelData<GLubyte>(resw, resh);
 339       int index = 0;
 340       for (int j = 0; j < resh; ++j)
 341       {
 342         for (int i = 0; i < resw; ++i)
 343         {
 344           // Geo coordinates
 345           maths::Vector<2,double> portion((double)i/(resw-1), (double)j/(resh-1));
 346           maths::Vector<2,double> tex = portionToTex * portion;
 347           maths::Vector<2,double> geo = (m_texToGeo * tex) * 180.0f/M_PI;
 348           // Local coordinates
 349           m_geoToProj->Transform(1, &geo[0], &geo[1]);
 350           // To pixel coordinates
 351           maths::Vector<2,double> pix = m_projToPixels * geo;
 352           maths::Vector<2,double> window = inToWindow * pix;
 353           buffer->buffer()[index] = (GLubyte)data.tcPixelData<GLubyte>::sampleFloat(window[0], window[1]);
 354           //buffer->buffer()[index] = 255.0f*portion[0];//data.buffer()[index];
 355           ++index;
 356         }
 357       }
 358     }
 359   }
 360   return buffer;
 361 }