vtkGeoAlignedImageSource.cxx 13.9 KB
Newer Older
Ken Martin's avatar
Ken Martin committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
/*=============================================================================

  Program:   Visualization Toolkit
  Module:    vtkGeoAlignedImageSource.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=============================================================================*/
/*-------------------------------------------------------------------------
  Copyright 2008 Sandia Corporation.
  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
  the U.S. Government retains certain rights in this software.
-------------------------------------------------------------------------*/

#include "vtkGeoAlignedImageSource.h"

#include "vtkCommand.h"
24 25
#include "vtkGeoImageNode.h"
#include "vtkImageData.h"
Ken Martin's avatar
Ken Martin committed
26
#include "vtkImageShrink3D.h"
27
#include "vtkMultiBlockDataSet.h"
Ken Martin's avatar
Ken Martin committed
28
#include "vtkObjectFactory.h"
29 30
#include "vtkSmartPointer.h"
#include "vtkTexture.h"
Ken Martin's avatar
Ken Martin committed
31
#include "vtkTimerLog.h"
32
#include "vtkTransform.h"
Ken Martin's avatar
Ken Martin committed
33

34
#include <cassert>
Ken Martin's avatar
Ken Martin committed
35

36 37
vtkStandardNewMacro(vtkGeoAlignedImageSource);
vtkCxxSetObjectMacro(vtkGeoAlignedImageSource, Image, vtkImageData);
Ken Martin's avatar
Ken Martin committed
38 39 40 41 42 43 44

class vtkGeoAlignedImageSource::vtkProgressObserver : public vtkCommand
{
public:
  static vtkProgressObserver* New()
    { return new vtkProgressObserver(); }

45
  void Execute(vtkObject *, unsigned long eventId,
46
    void *callData) override
47
  {
Ken Martin's avatar
Ken Martin committed
48
    if (eventId == vtkCommand::ProgressEvent)
49
    {
Ken Martin's avatar
Ken Martin committed
50 51 52
      double progress = *reinterpret_cast<double*>(callData);
      progress = this->Offset + this->Scale * progress;
      if (this->Target)
53
      {
Ken Martin's avatar
Ken Martin committed
54 55 56
        this->Target->InvokeEvent(vtkCommand::ProgressEvent, &progress);
      }
    }
57
  }
Ken Martin's avatar
Ken Martin committed
58 59

  void SetTarget(vtkObject* t)
60
  {
Ken Martin's avatar
Ken Martin committed
61
    this->Target = t;
62
  }
Ken Martin's avatar
Ken Martin committed
63 64 65 66 67 68

  double Offset;
  double Scale;

private:
  vtkProgressObserver()
69
  {
70
    this->Target = nullptr;
Ken Martin's avatar
Ken Martin committed
71 72
    this->Offset = 0.0;
    this->Scale = 1.0;
73
  }
Ken Martin's avatar
Ken Martin committed
74 75 76 77
  vtkObject* Target;
};

//----------------------------------------------------------------------------
78
vtkGeoAlignedImageSource::vtkGeoAlignedImageSource()
Ken Martin's avatar
Ken Martin committed
79
{
Sankhesh Jhaveri's avatar
Sankhesh Jhaveri committed
80 81
  VTK_LEGACY_BODY(vtkGeoAlignedImageSource::vtkGeoAlignedImageSource,
                  "VTK 8.2");
82
  this->Image = nullptr;
83 84 85 86 87
  this->LevelImages = vtkMultiBlockDataSet::New();
  this->LatitudeRange[0] = -90;
  this->LatitudeRange[1] = 90;
  this->LongitudeRange[0] = -180;
  this->LongitudeRange[1] = 180;
Ken Martin's avatar
Ken Martin committed
88 89
  this->ProgressObserver = vtkProgressObserver::New();
  this->ProgressObserver->SetTarget(this);
90 91
  this->PowerOfTwoSize = true;
  this->Overlap = 0.0;
Ken Martin's avatar
Ken Martin committed
92 93 94
}

//-----------------------------------------------------------------------------
95 96
vtkGeoAlignedImageSource::~vtkGeoAlignedImageSource()
{
97
  this->SetImage(nullptr);
98 99
  this->LevelImages->Delete();

100
  this->ProgressObserver->SetTarget(nullptr);
Ken Martin's avatar
Ken Martin committed
101
  this->ProgressObserver->Delete();
102
  this->ProgressObserver = nullptr;
Ken Martin's avatar
Ken Martin committed
103 104 105 106 107
}

//-----------------------------------------------------------------------------
void vtkGeoAlignedImageSource::PrintSelf(ostream& os, vtkIndent indent)
{
108 109 110
  this->Superclass::PrintSelf(os, indent);
  os << indent << "Image: " << (this->Image ? "" : "(null)") << endl;
  if (this->Image)
111
  {
112
    this->Image->PrintSelf(os, indent.GetNextIndent());
113
  }
114 115
  os << indent << "LatitudeRange: " << this->LatitudeRange[0] << "," << this->LatitudeRange[1] << endl;
  os << indent << "LongitudeRange: " << this->LongitudeRange[0] << "," << this->LongitudeRange[1] << endl;
116 117
  os << indent << "PowerOfTwoSize: " << (this->PowerOfTwoSize ? "On" : "Off") << endl;
  os << indent << "Overlap: " << this->Overlap << endl;
Ken Martin's avatar
Ken Martin committed
118 119 120
}

//-----------------------------------------------------------------------------
121
bool vtkGeoAlignedImageSource::FetchRoot(vtkGeoTreeNode* r)
Ken Martin's avatar
Ken Martin committed
122
{
123
  vtkGeoImageNode* root = nullptr;
124
  if (!(root = vtkGeoImageNode::SafeDownCast(r)))
125
  {
126 127
    vtkErrorMacro(<< "Node must be an image node for this source.");
    return false;
128
  }
Ken Martin's avatar
Ken Martin committed
129
  int imageDims[3];
130
  this->Image->GetDimensions(imageDims);
Ken Martin's avatar
Ken Martin committed
131 132 133 134 135 136 137 138 139 140 141

  // I am ignoring the geometry of the image, and assuming the scalars
  // are cell data.  The normal shrink should not shift the image by half
  // a pixel.  I believe texture maps will preserve the image bounds.
  vtkSmartPointer<vtkImageShrink3D> shrink = vtkSmartPointer<vtkImageShrink3D>::New();
  shrink->SetShrinkFactors(2,2,1);
  shrink->AveragingOn();
  shrink->AddObserver(vtkCommand::ProgressEvent, this->ProgressObserver);

  // We count the number of times vtkImageShrink3D will be executed so that we
  // can report progress correctly.
142
  int numLevels = 0;
Ken Martin's avatar
Ken Martin committed
143
  while (imageDims[0] > 300 || imageDims[1] > 300)
144
  {
145
    imageDims[0] = static_cast<int>(floor(imageDims[0] /
Ken Martin's avatar
Ken Martin committed
146
        static_cast<double>(shrink->GetShrinkFactors()[0])));
147
    imageDims[1] = static_cast<int>(floor(imageDims[1] /
Ken Martin's avatar
Ken Martin committed
148
        static_cast<double>(shrink->GetShrinkFactors()[1])));
149
    numLevels++;
150
  }
151
  this->Image->GetDimensions(imageDims);
Ken Martin's avatar
Ken Martin committed
152 153 154 155 156 157 158 159 160

  // Nothing says that the images cannot overlap and be larger than
  // the terain pathces.  Nothing says that the images have to be
  // a same size for all nodes either.

  // The easiest this to do to get multiple resolutions is to reduce
  // the image size before traversing.  This way we can avoid issues
  // with the bottom up approach.  Specifically, we do not need
  // to combine tiles, or worry about seams from smoothing.
161

Ken Martin's avatar
Ken Martin committed
162 163
  // This is not the best termination condition, but it will do.
  // This should also work for images that do not cover the whole globe.
164 165 166 167 168 169 170
  vtkSmartPointer<vtkImageData> image = vtkSmartPointer<vtkImageData>::New();
  image->ShallowCopy(this->Image);
  vtkSmartPointer<vtkImageData> fullImage = vtkSmartPointer<vtkImageData>::New();
  fullImage->ShallowCopy(this->Image);
  vtkSmartPointer<vtkMultiBlockDataSet> tempBlocks = vtkSmartPointer<vtkMultiBlockDataSet>::New();
  tempBlocks->SetBlock(0, fullImage);
  for (unsigned int curIter=0; imageDims[0] > 300 || imageDims[1] > 300; ++curIter)
171
  {
172 173 174
    this->ProgressObserver->Offset = curIter * 1.0/numLevels;
    this->ProgressObserver->Scale = 1.0/numLevels;

Ken Martin's avatar
Ken Martin committed
175
    // Shrink image for the next level.
176
    shrink->SetInputData(image);
Ken Martin's avatar
Ken Martin committed
177 178
    shrink->Update();
    image->ShallowCopy(shrink->GetOutput());
179
    shrink->SetInputData(nullptr);
Ken Martin's avatar
Ken Martin committed
180
    image->GetDimensions(imageDims);
181 182 183 184 185 186 187

    // Store the image for the level.
    vtkSmartPointer<vtkImageData> block = vtkSmartPointer<vtkImageData>::New();
    block->ShallowCopy(shrink->GetOutput());
    block->SetOrigin(-180, -90, 0);
    block->SetSpacing(180, 90, 0);
    tempBlocks->SetBlock(curIter+1, block);
188
  }
189 190 191

  // Reverse the coarsened images so they are in order by level.
  for (unsigned int block = 0; block < tempBlocks->GetNumberOfBlocks(); ++block)
192
  {
193 194
    this->LevelImages->SetBlock(tempBlocks->GetNumberOfBlocks() - 1 - block,
      tempBlocks->GetBlock(block));
195
  }
196
  vtkSmartPointer<vtkTexture> texture = vtkSmartPointer<vtkTexture>::New();
197
  texture->SetInputData(this->LevelImages->GetBlock(0));
198 199 200 201 202 203 204 205 206
  vtkSmartPointer<vtkTransform> texTrans = vtkSmartPointer<vtkTransform>::New();
  // Start with (lat,lon)
  texTrans->PostMultiply();
  texTrans->RotateZ(90.0); // (-lon,lat)
  texTrans->Scale(-1.0, 1.0, 1.0); // (lon,lat)
  texTrans->Translate(-(-180.0), -(-90.0), 0.0); // to origin
  texTrans->Scale(1.0/360.0, 1.0/180.0, 1.0); // to [0,1]
  texture->SetTransform(texTrans);
  texture->InterpolateOn();
207 208
  texture->RepeatOff();
  texture->EdgeClampOn();
209 210 211 212 213 214

  root->SetLevel(-1);
  root->SetLatitudeRange(-270, 90);
  root->SetLongitudeRange(-180, 180);
  root->SetTexture(texture);
  return true;
Ken Martin's avatar
Ken Martin committed
215 216
}

217 218
//------------------------------------------------------------------------------
bool vtkGeoAlignedImageSource::FetchChild(vtkGeoTreeNode* p, int index, vtkGeoTreeNode* c)
Ken Martin's avatar
Ken Martin committed
219
{
220
  vtkGeoImageNode* parent = nullptr;
221
  if (!(parent = vtkGeoImageNode::SafeDownCast(p)))
222
  {
223 224
    vtkErrorMacro(<< "Node must be an image node for this source.");
    return false;
225
  }
226
  vtkGeoImageNode* child = nullptr;
227
  if (!(child = vtkGeoImageNode::SafeDownCast(c)))
228
  {
229 230
    vtkErrorMacro(<< "Node must be an image node for this source.");
    return false;
231
  }
232 233 234
  int level = parent->GetLevel() + 1;
  unsigned int blockLevel = level + 1;
  if (blockLevel >= this->LevelImages->GetNumberOfBlocks())
235
  {
236 237
    vtkDebugMacro(<< "Reached max number of blocks (" << this->LevelImages->GetNumberOfBlocks() << ")");
    return false;
238
  }
Ken Martin's avatar
Ken Martin committed
239

240 241 242 243 244 245 246 247 248 249
  double lonRange[2];
  double latRange[2];
  double center[2];
  parent->GetLongitudeRange(lonRange);
  parent->GetLatitudeRange(latRange);
  center[0] = (lonRange[1] + lonRange[0])/2.0;
  center[1] = (latRange[1] + latRange[0])/2.0;

  child->SetLevel(level);
  if (index / 2)
250
  {
251
    child->SetLatitudeRange(center[1], latRange[1]);
252
  }
253
  else
254
  {
255
    child->SetLatitudeRange(latRange[0], center[1]);
256
  }
257
  if (index % 2)
258
  {
259
    child->SetLongitudeRange(center[0], lonRange[1]);
260
  }
261
  else
262
  {
263
    child->SetLongitudeRange(lonRange[0], center[0]);
264
  }
265 266 267

  int id = 0;
  if (level == 0)
268
  {
269 270 271 272
    // Special case: in the first level, the western hemisphere has id 0, and
    // the eastern hemisphere has id 1. This is to be compatible with the old
    // tile database format.
    if (index == 2)
273
    {
274
      id = 0;
275
    }
276
    else if (index == 3)
277
    {
278
      id = 1;
279
    }
280
    else if (index == 0)
281
    {
282 283 284
      vtkSmartPointer<vtkImageData> dummyImageWest = vtkSmartPointer<vtkImageData>::New();
      dummyImageWest->SetOrigin(-180.0, -270.0, 0.0);
      dummyImageWest->SetSpacing(0.0, -90.0, 0.0);
285
      child->GetTexture()->SetInputData(dummyImageWest);
286 287 288 289
      child->SetLatitudeRange(-270, -90);
      child->SetLongitudeRange(-180, 0);
      child->SetId(2);
      return true;
290
    }
291
    else if (index == 1)
292
    {
293 294 295
      vtkSmartPointer<vtkImageData> dummyImageEast = vtkSmartPointer<vtkImageData>::New();
      dummyImageEast->SetOrigin(0.0, -270.0, 0.0);
      dummyImageEast->SetSpacing(180.0, -90.0, 0.0);
296
      child->GetTexture()->SetInputData(dummyImageEast);
297 298 299 300
      child->SetLatitudeRange(-270, -90);
      child->SetLongitudeRange(0, 180);
      child->SetId(3);
      return true;
Ken Martin's avatar
Ken Martin committed
301
    }
302
  }
303
  else
304
  {
305
    id = parent->GetId() | (index << (2*level - 1));
306
  }
307 308 309 310 311 312
  child->SetId(id);

  // Crop and save the image.
  // Overwrite an image if it already exists.
  this->CropImageForNode(child, vtkImageData::SafeDownCast(this->LevelImages->GetBlock(blockLevel)));
  return true;
Ken Martin's avatar
Ken Martin committed
313 314 315
}

//------------------------------------------------------------------------------
316
void vtkGeoAlignedImageSource::CropImageForNode(vtkGeoImageNode* node, vtkImageData* image)
Ken Martin's avatar
Ken Martin committed
317
{
318 319 320 321 322 323 324 325 326
  int ext[6];
  int wholeExt[6];

  // I am keeping this all external to the imageData object because
  // I consider pixels are cells not points.
  image->GetExtent(ext);
  image->GetExtent(wholeExt);
  double origin[2];
  double spacing[2];
Francois Bertel's avatar
Francois Bertel committed
327 328 329 330
  spacing[0] = (this->LongitudeRange[1]-this->LongitudeRange[0])/static_cast<double>(ext[1]-ext[0]+1);
  spacing[1] = (this->LatitudeRange[1]-this->LatitudeRange[0])/static_cast<double>(ext[3]-ext[2]+1);
  origin[0] = this->LongitudeRange[0] - ext[0]*spacing[0];
  origin[1] = this->LatitudeRange[0] - ext[2]*spacing[1];
331 332

  // Compute the minimum extent that covers the terrain patch.
333 334 335 336 337 338 339
  double overlapDist[2];
  overlapDist[0] = this->Overlap*(node->GetLongitudeRange()[1]-node->GetLongitudeRange()[0]);
  overlapDist[1] = this->Overlap*(node->GetLatitudeRange()[1]-node->GetLatitudeRange()[0]);
  ext[0] = static_cast<int>(floor((node->GetLongitudeRange()[0]-overlapDist[0]-origin[0])/spacing[0]));
  ext[1] = static_cast<int>(ceil((node->GetLongitudeRange()[1]+overlapDist[0]-origin[0])/spacing[0]));
  ext[2] = static_cast<int>(floor((node->GetLatitudeRange()[0]-overlapDist[1]-origin[1])/spacing[1]));
  ext[3] = static_cast<int>(ceil((node->GetLatitudeRange()[1]+overlapDist[1]-origin[1])/spacing[1]));
340 341

  int dims[2];
342
  if (this->PowerOfTwoSize)
343
  {
344 345 346 347
    dims[0] = this->PowerOfTwo(ext[1]-ext[0]+1);
    dims[1] = this->PowerOfTwo(ext[3]-ext[2]+1);
    ext[1] = ext[0] + dims[0] - 1;
    ext[3] = ext[2] + dims[1] - 1;
348
  }
349
  else
350
  {
351 352
    dims[0] = ext[1]-ext[0]+1;
    dims[1] = ext[3]-ext[2]+1;
353
  }
354

355
  if (ext[1] > wholeExt[1])
356
  {
357
    ext[1] = wholeExt[1];
358
  }
359
  if (ext[3] > wholeExt[3])
360
  {
361
    ext[3] = wholeExt[3];
362
  }
363 364 365
  ext[0] = ext[1] - dims[0] + 1;
  ext[2] = ext[3] - dims[1] + 1;
  if (ext[0] < wholeExt[0])
366
  {
367
    ext[0] = wholeExt[0];
368
  }
369
  if (ext[2] < wholeExt[2])
370
  {
371
    ext[2] = wholeExt[2];
372
  }
373 374 375

  vtkSmartPointer<vtkImageData> cropped = vtkSmartPointer<vtkImageData>::New();
  cropped->ShallowCopy(image);
376
  cropped->Crop(ext);
377 378 379 380

  // Now set the longitude and latitude range based on the actual image size.
  double lonRange[2];
  double latRange[2];
Francois Bertel's avatar
Francois Bertel committed
381 382 383 384
  lonRange[0] = origin[0] + ext[0]*spacing[0];
  lonRange[1] = origin[0] + (ext[1]+1)*spacing[0];
  latRange[0] = origin[1] + ext[2]*spacing[1];
  latRange[1] = origin[1] + (ext[3]+1)*spacing[1];
385 386 387 388 389 390 391 392 393 394 395 396 397 398
  cropped->SetOrigin(lonRange[0], latRange[0], 0);
  cropped->SetSpacing(lonRange[1], latRange[1], 0);
  //assert(lonRange[1] >= lonRange[0]);
  //assert(latRange[1] >= latRange[0]);

  vtkSmartPointer<vtkTexture> tex = vtkSmartPointer<vtkTexture>::New();
  vtkSmartPointer<vtkTransform> texTrans = vtkSmartPointer<vtkTransform>::New();
  // Start with (lat,lon)
  texTrans->PostMultiply();
  texTrans->RotateZ(90.0); // (-lon,lat)
  texTrans->Scale(-1.0, 1.0, 1.0); // (lon,lat)
  texTrans->Translate(-lonRange[0], -latRange[0], 0.0); // to origin
  texTrans->Scale(1.0/(lonRange[1] - lonRange[0]), 1.0/(latRange[1] - latRange[0]), 1.0); // to [0,1]
  tex->SetTransform(texTrans);
399
  tex->SetInputData(cropped);
400
  tex->InterpolateOn();
401 402
  tex->RepeatOff();
  tex->EdgeClampOn();
403 404 405 406 407 408 409 410 411 412 413 414

  node->SetTexture(tex);
}

//-----------------------------------------------------------------------------
int vtkGeoAlignedImageSource::PowerOfTwo(int val)
{
  // Pick the next highest power of two.
  int tmp;
  bool nextHigherFlag = false;
  tmp = 1;
  while (val)
415
  {
416
    if ((val & 1) && val > 1)
417
    {
418
      nextHigherFlag = true;
419
    }
420 421
    val = val >> 1;
    tmp = tmp << 1;
422
  }
Ken Martin's avatar
Ken Martin committed
423

424
  if ( ! nextHigherFlag)
425
  {
426
    tmp = tmp >> 1;
427
  }
428
  return tmp;
Ken Martin's avatar
Ken Martin committed
429
}