1 /***************************************************************************
2 * This file is part of Tecorrec. *
3 * Copyright 2008 James Hogan <james@albanarts.com> *
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. *
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. *
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 ***************************************************************************/
21 * @file tcChannelFile.cpp
22 * @brief An image channel from an image file.
25 #include "tcChannelFile.h"
26 #include "tcPixelData.h"
28 #include <gdal_priv.h>
29 #include <ogr_spatialref.h>
36 * Constructors + destructor
39 /// Primary constructor.
40 tcChannelFile::tcChannelFile(const QString
& filename
, const QString
& name
, const QString
& description
, float radianceOffset
, float radianceBias
)
41 : tcChannel(name
, description
)
53 setRadianceTransform(radianceOffset
, radianceBias
);
56 /// @todo this probably isn't the best place to set it up
59 m_dataset
= (GDALDataset
*)GDALOpen((const char*)filename
.toAscii(), GA_ReadOnly
);
62 GDALRasterBand
* band
= m_dataset
->GetRasterBand(1);
64 m_width
= band
->GetXSize();
65 m_height
= band
->GetYSize();
67 OGRCoordinateTransformation
* result
= 0;
68 if (0 != m_dataset
->GetProjectionRef())
70 OGRSpatialReference
srs(m_dataset
->GetProjectionRef());
71 OGRSpatialReference
* geog
= srs
.CloneGeogCS();
72 m_projToGeo
= OGRCreateCoordinateTransformation(&srs
, geog
);
73 m_geoToProj
= OGRCreateCoordinateTransformation(geog
, &srs
);
75 m_dataset
->GetGeoTransform(m_pixelsToProj
);
76 GDALInvGeoTransform(m_pixelsToProj
, m_projToPixels
);
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 )
85 maths::Vector
<2,double> min
, max
;
86 for (int i
= 0; i
< 4; ++i
)
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
93 for (int j
= 0; j
< 2; ++j
)
95 if (lonlat
[j
] < min
[j
])
99 else if (lonlat
[j
] > max
[j
])
107 for (int j
= 0; j
< 2; ++j
)
109 min
[j
] = max
[j
] = lonlat
[j
];
113 // Convert extents into radians
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();
126 tcChannelFile::~tcChannelFile()
129 if (0 != m_thumbnail
)
131 glDeleteTextures(1, &m_thumbnail
);
142 /// Load the transformation from projection to texture.
143 const tcAffineTransform2
<double> tcChannelFile::geoToTex() const
148 /// Load the transformation from projection to texture.
149 const tcAffineTransform2
<double> tcChannelFile::texToGeo() const
155 * Main image interface
159 GLuint
tcChannelFile::thumbnailTexture()
161 if (0 == m_thumbnail
&& 0 != m_dataset
)
163 GDALRasterBand
* band
= m_dataset
->GetRasterBand(1);
167 GDALRasterBand
* thumbnailBand
= 0;
168 int numOverviews
= band
->GetOverviewCount();
169 if (0 == numOverviews
)
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);
178 if (CE_None
!= m_dataset
->BuildOverviews("NEAREST", 1, &overview
, 0, 0,
179 GDALDummyProgress
, 0))
181 thumbnailBand
= band
->GetOverview(0);
187 /// @todo Pick an appropriate overview if more than one
188 thumbnailBand
= band
->GetOverview(0);
193 // Turn into a GL texture
194 int width
= thumbnailBand
->GetXSize();
195 int height
= thumbnailBand
->GetYSize();
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))
204 for (int j
= 0; j
< tnHeight
; ++j
)
206 for (int i
= 0; i
< tnWidth
; ++i
)
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
;
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));
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()))
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);
231 glBindTexture(GL_TEXTURE_2D
, 0);
232 glDeleteTextures(1, &m_thumbnail
);
242 * Interface for derived class to implement
245 void tcChannelFile::roundPortion(double* x1
, double* y1
, double* x2
, double* y2
)
247 // We're resampling the texture anyway so no need to adjust range
250 tcAbstractPixelData
* tcChannelFile::loadPortion(double x1
, double y1
, double x2
, double y2
, bool changed
)
252 tcPixelData
<GLubyte
>* buffer
= 0;
255 GDALRasterBand
* band
= m_dataset
->GetRasterBand(1);
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
)
267 maths::Vector
<2,double> min
, max
;
268 for (int i
= 0; i
< 4; ++i
)
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
277 for (int j
= 0; j
< 2; ++j
)
283 else if (pix
[j
] > max
[j
])
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();
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
;
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
)
318 if (0 == (resw
& ~((1<<i
)-1)))
324 for (int i
= 1; i
< 20; ++i
)
326 if (0 == (resh
& ~((1<<i
)-1)))
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))
338 buffer
= new tcPixelData
<GLubyte
>(resw
, resh
);
340 for (int j
= 0; j
< resh
; ++j
)
342 for (int i
= 0; i
< resw
; ++i
)
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
;
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];