1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkRasterReprojectionFilter.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 #include "vtkRasterReprojectionFilter.h"
16
17 // VTK includes
18 #include "vtkCellData.h"
19 #include "vtkDataArray.h"
20 #include "vtkDataObject.h"
21 #include "vtkGDAL.h"
22 #include "vtkGDALRasterConverter.h"
23 #include "vtkGDALRasterReprojection.h"
24 #include "vtkImageData.h"
25 #include "vtkInformation.h"
26 #include "vtkInformationVector.h"
27 #include "vtkMath.h"
28 #include "vtkNew.h"
29 #include "vtkObjectFactory.h"
30 #include "vtkPointData.h"
31 #include "vtkStreamingDemandDrivenPipeline.h"
32 #include "vtkUniformGrid.h"
33
34 // GDAL includes
35 #include <gdal_priv.h>
36
37 // STL includes
38 #include <algorithm>
39
40 vtkStandardNewMacro(vtkRasterReprojectionFilter);
41
42 //----------------------------------------------------------------------------
43 class vtkRasterReprojectionFilter::vtkRasterReprojectionFilterInternal
44 {
45 public:
46 vtkRasterReprojectionFilterInternal();
47 ~vtkRasterReprojectionFilterInternal();
48
49 vtkGDALRasterConverter* GDALConverter;
50 vtkGDALRasterReprojection* GDALReprojection;
51
52 // Data saved during RequestInformation()
53 int InputImageExtent[6];
54 double OutputImageGeoTransform[6];
55 };
56
57 //----------------------------------------------------------------------------
58 vtkRasterReprojectionFilter::vtkRasterReprojectionFilterInternal::
vtkRasterReprojectionFilterInternal()59 vtkRasterReprojectionFilterInternal()
60 {
61 this->GDALConverter = vtkGDALRasterConverter::New();
62 this->GDALReprojection = vtkGDALRasterReprojection::New();
63 std::fill(this->InputImageExtent, this->InputImageExtent + 6, 0);
64 std::fill(
65 this->OutputImageGeoTransform, this->OutputImageGeoTransform + 6, 0.0);
66 }
67
68 //----------------------------------------------------------------------------
69 vtkRasterReprojectionFilter::vtkRasterReprojectionFilterInternal::
~vtkRasterReprojectionFilterInternal()70 ~vtkRasterReprojectionFilterInternal()
71 {
72 this->GDALConverter->Delete();
73 this->GDALReprojection->Delete();
74 }
75
76 //----------------------------------------------------------------------------
vtkRasterReprojectionFilter()77 vtkRasterReprojectionFilter::vtkRasterReprojectionFilter()
78 {
79 this->Internal = new vtkRasterReprojectionFilterInternal;
80 this->InputProjection = NULL;
81 this->OutputProjection = NULL;
82 this->OutputDimensions[0] = this->OutputDimensions[1] = 0;
83 this->NoDataValue = vtkMath::Nan();
84 this->MaxError = 0.0;
85 this->ResamplingAlgorithm = 0;
86
87 // Enable all the drivers.
88 GDALAllRegister();
89 }
90
91 //----------------------------------------------------------------------------
~vtkRasterReprojectionFilter()92 vtkRasterReprojectionFilter::~vtkRasterReprojectionFilter()
93 {
94 if (this->InputProjection)
95 {
96 delete[] this->InputProjection;
97 }
98 if (this->OutputProjection)
99 {
100 delete[] this->OutputProjection;
101 }
102 delete this->Internal;
103 }
104
105 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)106 void vtkRasterReprojectionFilter::PrintSelf(ostream& os, vtkIndent indent)
107 {
108 this->Superclass::PrintSelf(os, indent);
109
110 os << indent << "InputProjection: ";
111 if (this->InputProjection)
112 {
113 os << this->InputProjection;
114 }
115 else
116 {
117 os << "(not specified)";
118 }
119 os << "\n";
120
121 os << indent << "OutputProjection: ";
122 if (this->OutputProjection)
123 {
124 os << this->OutputProjection;
125 }
126 else
127 {
128 os << "(not specified)";
129 }
130 os << "\n";
131
132 os << indent << "OutputDimensions: " << OutputDimensions[0] << ", "
133 << OutputDimensions[1] << "\n"
134 << indent << "NoDataValue: " << this->NoDataValue << "\n"
135 << indent << "MaxError: " << this->MaxError << "\n"
136 << indent << "ResamplingAlgorithm: " << this->ResamplingAlgorithm << "\n"
137 << std::endl;
138 }
139
140 //-----------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)141 int vtkRasterReprojectionFilter::RequestData(
142 vtkInformation* vtkNotUsed(request),
143 vtkInformationVector** inputVector,
144 vtkInformationVector* outputVector)
145 {
146 // Get the input image data
147 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
148 if (!inInfo)
149 {
150 return 0;
151 }
152
153 vtkDataObject* inDataObject = inInfo->Get(vtkDataObject::DATA_OBJECT());
154 if (!inDataObject)
155 {
156 return 0;
157 }
158
159 vtkImageData* inImageData = vtkImageData::SafeDownCast(inDataObject);
160 if (!inImageData)
161 {
162 return 0;
163 }
164 // std::cout << "RequestData() has image to reproject!" << std::endl;
165
166 // Get the output image
167 vtkInformation* outInfo = outputVector->GetInformationObject(0);
168 if (!outInfo)
169 {
170 return 0;
171 }
172
173 // Convert input image to GDALDataset
174 GDALDataset* inputGDAL = this->Internal->GDALConverter->CreateGDALDataset(
175 inImageData, this->InputProjection);
176
177 if (this->Debug)
178 {
179 std::string tifFileName = "inputGDAL.tif";
180 this->Internal->GDALConverter->WriteTifFile(inputGDAL, tifFileName.c_str());
181 std::cout << "Wrote " << tifFileName << std::endl;
182
183 double minValue, maxValue;
184 this->Internal->GDALConverter->FindDataRange(
185 inputGDAL, 1, &minValue, &maxValue);
186 std::cout << "Min: " << minValue << " Max: " << maxValue << std::endl;
187 }
188
189 // Construct GDAL dataset for output image
190 vtkDataArray* array = inImageData->GetCellData()->GetScalars();
191 int vtkDataType = array->GetDataType();
192 int rasterCount = array->GetNumberOfComponents();
193 GDALDataset* outputGDAL =
194 this->Internal->GDALConverter->CreateGDALDataset(this->OutputDimensions[0],
195 this->OutputDimensions[1],
196 vtkDataType,
197 rasterCount);
198 this->Internal->GDALConverter->CopyBandInfo(inputGDAL, outputGDAL);
199 this->Internal->GDALConverter->SetGDALProjection(outputGDAL,
200 this->OutputProjection);
201 outputGDAL->SetGeoTransform(this->Internal->OutputImageGeoTransform);
202 this->Internal->GDALConverter->CopyNoDataValues(inputGDAL, outputGDAL);
203
204 // Apply the reprojection
205 this->Internal->GDALReprojection->SetMaxError(this->MaxError);
206 this->Internal->GDALReprojection->SetResamplingAlgorithm(
207 this->ResamplingAlgorithm);
208 this->Internal->GDALReprojection->Reproject(inputGDAL, outputGDAL);
209
210 if (this->Debug)
211 {
212 std::string tifFileName = "reprojectGDAL.tif";
213 this->Internal->GDALConverter->WriteTifFile(outputGDAL,
214 tifFileName.c_str());
215 std::cout << "Wrote " << tifFileName << std::endl;
216 double minValue, maxValue;
217 this->Internal->GDALConverter->FindDataRange(
218 outputGDAL, 1, &minValue, &maxValue);
219 std::cout << "Min: " << minValue << " Max: " << maxValue << std::endl;
220 }
221
222 // Done with input GDAL dataset
223 GDALClose(inputGDAL);
224
225 // Convert output dataset to vtkUniformGrid
226 vtkUniformGrid* reprojectedImage =
227 this->Internal->GDALConverter->CreateVTKUniformGrid(outputGDAL);
228
229 // Done with output GDAL dataset
230 GDALClose(outputGDAL);
231
232 // Update pipeline output instance
233 vtkUniformGrid* output = vtkUniformGrid::GetData(outInfo);
234 output->ShallowCopy(reprojectedImage);
235
236 reprojectedImage->Delete();
237 return VTK_OK;
238 }
239
240 //-----------------------------------------------------------------------------
RequestUpdateExtent(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * vtkNotUsed (outputVector))241 int vtkRasterReprojectionFilter::RequestUpdateExtent(
242 vtkInformation* vtkNotUsed(request),
243 vtkInformationVector** inputVector,
244 vtkInformationVector* vtkNotUsed(outputVector))
245 {
246 // Set input extent to values saved in last RequestInformation() call
247 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
248 inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),
249 this->Internal->InputImageExtent,
250 6);
251 return VTK_OK;
252 }
253
254 //-----------------------------------------------------------------------------
RequestInformation(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)255 int vtkRasterReprojectionFilter::RequestInformation(
256 vtkInformation* vtkNotUsed(request),
257 vtkInformationVector** inputVector,
258 vtkInformationVector* outputVector)
259 {
260 // Get the info objects
261 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
262 if (!inInfo->Has(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()) ||
263 !inInfo->Has(vtkDataObject::SPACING()) ||
264 !inInfo->Has(vtkDataObject::ORIGIN()))
265 {
266 vtkErrorMacro("Input information missing");
267 return VTK_ERROR;
268 }
269 int* inputDataExtent =
270 inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT());
271 std::copy(
272 inputDataExtent, inputDataExtent + 6, this->Internal->InputImageExtent);
273
274 double* inputOrigin = inInfo->Get(vtkDataObject::ORIGIN());
275 double* inputSpacing = inInfo->Get(vtkDataObject::SPACING());
276 // std::cout << "Whole extent: " << inputDataExtent[0]
277 // << ", " << inputDataExtent[1]
278 // << ", " << inputDataExtent[2]
279 // << ", " << inputDataExtent[3]
280 // << ", " << inputDataExtent[4]
281 // << ", " << inputDataExtent[5]
282 // << std::endl;
283 // std::cout << "Input spacing: " << inputSpacing[0]
284 // << ", " << inputSpacing[1]
285 // << ", " << inputSpacing[2] << std::endl;
286 // std::cout << "Input origin: " << inputOrigin[0]
287 // << ", " << inputOrigin[1]
288 // << ", " << inputOrigin[2] << std::endl;
289
290 // InputProjection can be overridden, so only get from pipeline if needed
291 if (!this->InputProjection)
292 {
293 if (!inInfo->Has(vtkGDAL::MAP_PROJECTION()))
294 {
295 vtkErrorMacro("No map-projection for input image");
296 return VTK_ERROR;
297 }
298 this->SetInputProjection(inInfo->Get(vtkGDAL::MAP_PROJECTION()));
299 }
300
301 vtkInformation* outInfo = outputVector->GetInformationObject(0);
302 if (!outInfo)
303 {
304 vtkErrorMacro("Invalid output information object");
305 return VTK_ERROR;
306 }
307
308 // Validate current settings
309 if (!this->OutputProjection)
310 {
311 vtkErrorMacro("No output projection specified");
312 return VTK_ERROR;
313 }
314
315 // Create GDALDataset to compute suggested output
316 int xDim = inputDataExtent[1] - inputDataExtent[0] + 1;
317 int yDim = inputDataExtent[3] - inputDataExtent[2] + 1;
318 GDALDataset* gdalDataset = this->Internal->GDALConverter->CreateGDALDataset(
319 xDim, yDim, VTK_UNSIGNED_CHAR, 1);
320 this->Internal->GDALConverter->SetGDALProjection(gdalDataset,
321 this->InputProjection);
322 this->Internal->GDALConverter->SetGDALGeoTransform(
323 gdalDataset, inputOrigin, inputSpacing);
324
325 int nPixels = 0;
326 int nLines = 0;
327 this->Internal->GDALReprojection->SuggestOutputDimensions(
328 gdalDataset,
329 this->OutputProjection,
330 this->Internal->OutputImageGeoTransform,
331 &nPixels,
332 &nLines);
333 GDALClose(gdalDataset);
334
335 if ((this->OutputDimensions[0] < 1) || (this->OutputDimensions[1] < 1))
336 {
337 this->OutputDimensions[0] = nPixels;
338 this->OutputDimensions[1] = nLines;
339 }
340
341 // Set output info
342 int outputDataExtent[6] = {};
343 outputDataExtent[1] = this->OutputDimensions[0] - 1;
344 outputDataExtent[3] = this->OutputDimensions[1] - 1;
345 outInfo->Set(
346 vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), outputDataExtent, 6);
347
348 double outputImageOrigin[3] = {};
349 outputImageOrigin[0] = this->Internal->OutputImageGeoTransform[0];
350 outputImageOrigin[1] = this->Internal->OutputImageGeoTransform[3];
351 outputImageOrigin[2] = 1.0;
352 outInfo->Set(vtkDataObject::SPACING(), outputImageOrigin, 3);
353
354 double outputImageSpacing[3] = {};
355 outputImageSpacing[0] = this->Internal->OutputImageGeoTransform[1];
356 outputImageSpacing[1] = -this->Internal->OutputImageGeoTransform[5];
357 outputImageSpacing[2] = 1.0;
358 outInfo->Set(vtkDataObject::ORIGIN(), outputImageSpacing, 3);
359
360 return VTK_OK;
361 }
362
363 //-----------------------------------------------------------------------------
FillInputPortInformation(int port,vtkInformation * info)364 int vtkRasterReprojectionFilter::FillInputPortInformation(int port,
365 vtkInformation* info)
366 {
367 this->Superclass::FillInputPortInformation(port, info);
368 if (port == 0)
369 {
370 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkImageData");
371 return VTK_OK;
372 }
373 else
374 {
375 vtkErrorMacro("Input port: " << port << " is not a valid port");
376 return VTK_ERROR;
377 }
378 return VTK_OK;
379 }
380
381 //-----------------------------------------------------------------------------
FillOutputPortInformation(int port,vtkInformation * info)382 int vtkRasterReprojectionFilter::FillOutputPortInformation(int port,
383 vtkInformation* info)
384 {
385 if (port == 0)
386 {
387 info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkUniformGrid");
388 return VTK_OK;
389 }
390 else
391 {
392 vtkErrorMacro("Output port: " << port << " is not a valid port");
393 return VTK_ERROR;
394 }
395 }
396