evField.h

Go to the documentation of this file.

/*
 *
 * EVITA: Efficient Visualization of Terascale Datasets
 * Copyright (C) 2000-2016  Team Evita
 * 
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 * 
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 * 
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 675 Mass Ave, Cambridge,
 * MA 02139, USA.
 *
 */


#ifndef __evField_h
#define __evField_h


#include "evVector.h"
#include "evAddress.h"
#include "evFieldScaleInfo.h"


template <class T>
class evField : public evDataObject
{
public:
  evField(void)
    {
      this->Initialize();
    }
  virtual ~evField(void)
    { }

  virtual const evString GetClassName(void) const
  {
    return(evString("evField"));
  }
  
  bool IsAllocated(void) const;
  bool IsSameSize(const evField &field) const;

  int Copy(const evField &field);

  virtual int SetSize(const evAddress &size,
                      const bool allocate = true);
  virtual int Resize(const evAddress &size,
                     const bool linear = 0);

  void SetNumScales(const unsigned int num_scales);

  unsigned int GetNumDimensions(void) const;
  unsigned int GetSize(const unsigned int index) const;
  evAddress GetSize(void) const;
  unsigned int GetTotalSize(void) const;
  unsigned int GetNumScales(void) const;

  T GetField(const unsigned int index) const;
  T GetField(const evAddress &address) const;

  T &operator [] (const unsigned int index);
  const T &operator [] (const unsigned int index) const;
  T &operator [] (const evAddress &address);
  const T &operator [] (const evAddress &address) const;

  void PutField(const T &value);
  void PutField(const T &value, const unsigned int index);
  void PutField(const T &value, const evAddress &address);

  void ZeroField(void);

  int SubtractMean(T &mean);
  int AddMean(const T mean);

  double DistortionMSE(const evField &field) const;

  static unsigned int GetNumSubbands(const unsigned int num_scales,
                                     const unsigned int num_dimensions);
  static unsigned int NumSubbandsToNumScales(const unsigned int num_subbands,
                                             const unsigned int
                                             num_dimensions);
  static unsigned int CalcScaleFromSubband(const unsigned int subband,
                                           const unsigned int num_scales,
                                           const unsigned int num_dimensions);
  static int SubbandInfo(const unsigned int subband,
                         const unsigned int num_scales,
                         const evAddress &field_size,
                         const evAddress &field_origin,
                         evAddress &subband_size,
                         evAddress &subband_origin,
                         unsigned int &subband_scale);
  
  virtual int RasterScan(evField<T> &field,
                         const evAddress &field_origin) const;
  virtual int InverseRasterScan(evField<T> &field,
                                const evAddress &field_origin) const;
  virtual int InverseRasterScan(evField<T> &field,
                                const evAddress &field_origin,
                                const evAddress &full_size,
                                const unsigned int max_num_scales) const;

  virtual int GetScaleInfo(evFieldScaleInfo &scale_info,
                           const evAddress &field_origin) const;
  
  virtual int ExtractSubband(const unsigned int subband, evField<T> &field,
                             const evAddress &field_origin) const;
  int PreserveBaseband(const evAddress &field_origin);

  virtual void PrintSelf(const int tabbing) const;

protected:
  void Deallocate(void);
  virtual void Initialize(void);
  unsigned int NumDimensions;
  evAddress Size;
  std::vector<T> Field;

private:
  evField(const evField &field);

  int NullSubband(const evAddress &subband_size) const;

  unsigned int NumScales;
  unsigned int NumSubbands;
};


typedef evField<int> evMapField;
typedef evField<unsigned char> evFlagField;


template <class T>
inline void evField<T>::Deallocate(void)
{
  this->Field = std::vector<T>();
}


void Deallocate(void);
template <class T>
inline void evField<T>::Initialize(void)
{
  this->Deallocate();
  this->NumDimensions = 0;
  this->NumScales = 0;
  this->NumSubbands = 1;
}


template <class T>
inline bool evField<T>::IsAllocated(void) const
{
  return(!this->Field.empty());
}


template <class T>
inline bool evField<T>::IsSameSize(const evField &field) const
{
  return(this->Size == field.Size);
}


template <class T>
int evField<T>::Copy(const evField<T> &field)
{
  if (!this->IsAllocated())
    {
      if (this->SetSize(field.Size))
        {
          ErrorHandler.ErrorMessage("evField::Copy(): Error calling evField::SetSize()");
          return(1);
        }
    }
  else
    if (!this->IsSameSize(field))
      {
        ErrorHandler.ErrorMessage("evField::Copy(): Fields not the same size");
        return(1);
      }

  std::copy(field.Field.begin(),
            field.Field.end(),
            this->Field.begin());
  
  this->NumScales = field.NumScales;
  this->NumSubbands = field.NumSubbands;

  return(0);
}


template <class T>
int evField<T>::SetSize(const evAddress &size,
                        const bool allocate)
{
  this->Size = size;
  
  this->NumDimensions = this->Size.GetNumDimensions();

  unsigned int total_size = size.Product();

  T value = 0;

  this->Deallocate();

  if (allocate)
    this->Field.insert(this->Field.begin(), total_size, value);

  return(0);
}


template <class T>
int evField<T>::Resize(const evAddress &size,
                       const bool linear)
{
  if (this->Size == size)
    return(0);

  evField<T> old_field;
  old_field.Copy(*this);

  this->SetSize(size);

  if (linear)
    std::copy(old_field.Field.begin(),
              old_field.Field.begin() + 
              std::min(this->Field.size(), old_field.Field.size()),
              this->Field.begin());
  else
    {
      evAddress address(this->Size);
      address = 0;
      do
        {
          bool inside = true;
          for (unsigned int dimension = 0; dimension < address.GetNumDimensions();
               dimension++)
            if (address[dimension] >= old_field.Size[dimension])
              inside = false;
          if (inside)
            (*this)[address] = old_field[address];
        }
      while(!address.Increment(this->Size));
    }
  
  return(0);
}


template <class T>
inline void evField<T>::SetNumScales(const unsigned int num_scales)
{
  this->NumScales = num_scales;
  this->NumSubbands = this->GetNumSubbands(this->NumScales,
                                           this->NumDimensions);
}


template <class T>
inline unsigned int evField<T>::GetNumScales(void) const
{
  return(this->NumScales);
}


template <class T>
inline unsigned int evField<T>::GetNumDimensions(void) const
{
  return(this->NumDimensions);
}


template <class T>
inline unsigned int evField<T>::GetSize(const unsigned int index) const
{
  return(this->Size.Get(index));
}


template <class T>
inline evAddress evField<T>::GetSize(void) const
{
  return(this->Size);
}


template <class T>
inline unsigned int evField<T>::GetTotalSize(void) const
{
  return(this->Field.size());
}


template <class T>
inline T evField<T>::GetField(const unsigned int index) const
{
  return(this->Field[index]);
}


template <class T>
inline T evField<T>::GetField(const evAddress &address) const
{
  unsigned int dimension;
  unsigned int index;

  index = 0;
  for (dimension = this->NumDimensions - 1; dimension >= 1; dimension--)
    index = (index + address.Get(dimension)) * this->Size.Get(dimension - 1);
  index += address.Get(0);

  return(this->Field[index]);
}


template <class T>
inline T &evField<T>::operator [] (const unsigned int index)
{
  return(this->Field[index]);
}


template <class T>
inline const T &evField<T>::operator [] (const unsigned int index) const
{
  return(this->Field[index]);
}


template <class T>
inline T &evField<T>::operator [] (const evAddress &address)
{
  unsigned int dimension;
  unsigned int index;

  index = 0;
  for (dimension = this->NumDimensions - 1; dimension >= 1; dimension--)
    index = (index + address.Get(dimension)) * this->Size.Get(dimension - 1);
  index += address.Get(0);

  return(this->Field[index]);
}


template <class T>
inline const T &evField<T>::operator [] (const evAddress &address) const
{
  unsigned int dimension;
  unsigned int index;

  index = 0;
  for (dimension = this->NumDimensions - 1; dimension >= 1; dimension--)
    index = (index + address.Get(dimension)) * this->Size.Get(dimension - 1);
  index += address.Get(0);

  return(this->Field[index]);
}


template <class T>
inline void evField<T>::PutField(const T &value)
{
  std::fill(this->Field.begin(),
            this->Field.end(),
            value);
}


template <class T>
inline void evField<T>::PutField(const T &value, const unsigned int index)
{
  this->Field[index] = value;
}


template <class T>
inline void evField<T>::PutField(const T &value, const evAddress &address)
{
  unsigned int dimension;
  unsigned int index;

  for (index = 0, dimension = this->NumDimensions - 1;
       dimension >= 1;
       dimension--)
    index = (index + address.Get(dimension)) * this->Size.Get(dimension - 1);
  index += address.Get(0);

  this->Field[index] = value;
}


template <class T>
inline void evField<T>::ZeroField(void)
{
  if (!this->IsAllocated())
    return;

  std::fill(this->Field.begin(),
            this->Field.end(),
            0);
}


template <class T>
int evField<T>::SubtractMean(T &mean)
{
  int index;
  
  if (this->GetFieldMean(mean))
    {
      ErrorHandler.ErrorMessage("evField::SubtractMean(): Error calling evField::GetFieldMean()");
      return(1);
    }
  
  for (index = 0; index < this->GetTotalSize(); index++)
    this->Field[index] = this->Field[index] - mean;
  
  return(0);
}


template <class T>
int evField<T>::AddMean(const T mean)
{
  int index;
  
  if (!this->IsAllocated())
    {
      ErrorHandler.ErrorMessage("evField::AddMean(): Field not allocated");
      return(1);
    }
  
  for (index = 0; index < this->GetTotalSize(); index++)
    this->Field[index] = this->Field[index] + mean;
  
  return(0);
}


template <class T>
double evField<T>::DistortionMSE(const evField<T> &field) const
{
  int size;
  int index;
  T difference;
  double distortion = 0.0;
  
  if (!this->IsSameSize(field))
    {
      ErrorHandler.ErrorMessage("evField::DistortionMSE(): Data fields must be of same size and dimension");
      return(-1.0);
    }
  
  size = this->GetTotalSize();
  
  if (!size)
    return(0.0);
  
  for (index = 0; index < size; index++)
    {
      difference = this->Field[index] - field.Field[index];
      distortion += difference * difference;
    }
  
  distortion /= size;
  
  return(distortion);
}


template <class T>
unsigned int evField<T>::GetNumSubbands(const unsigned int num_scales,
                                        const unsigned int num_dimensions)
{
  int subband_constant;
  
  if (num_dimensions == 0)
    return(0);
  
  subband_constant = (1 << num_dimensions) - 1;
  
  return(num_scales * subband_constant + 1);
}


template <class T>
unsigned int evField<T>::NumSubbandsToNumScales(const unsigned int
                                                num_subbands,
                                                const unsigned int
                                                num_dimensions)
{
  int subband_constant;
  
  if (num_dimensions == 0)
    return(0);
  
  if (num_subbands == 0)
    return(0);
  
  subband_constant = (1 << num_dimensions) - 1;
  
  if ((num_subbands - 1) % subband_constant)
    return(0);
  
  return((num_subbands - 1) / subband_constant);
}


template <class T>
unsigned int evField<T>::CalcScaleFromSubband(const unsigned int subband,
                                              const unsigned int num_scales,
                                              const unsigned int
                                              num_dimensions)
{
  int subband_constant;
  
  if (num_dimensions == 0)
    return(0);
  
  subband_constant = (1 << num_dimensions) - 1;
  
  if (!num_scales)
    return(0);
  
  if (!subband)
    return(num_scales);
  
  return(num_scales - (subband - 1)/subband_constant);
}


template <class T>
int evField<T>::SubbandInfo(const unsigned int subband,
                            const unsigned int num_scales,
                            const evAddress &field_size,
                            const evAddress &field_origin,
                            evAddress &subband_size,
                            evAddress &subband_offset,
                            unsigned int &subband_scale)
{
  unsigned int dimension;
  unsigned int subband_constant;
  unsigned int subband_code;
  unsigned int num_dimensions = field_size.GetNumDimensions();
  evAddress baseband_size(num_dimensions);
  
  baseband_size = 0;
  subband_size = baseband_size;
  subband_offset = baseband_size;

  if (num_dimensions == 0)
    return(0);
  
  subband_constant = (1 << num_dimensions) - 1;
  
  subband_scale =
    evField<T>::CalcScaleFromSubband(subband,
                                     num_scales,
                                     num_dimensions);
  
  for (dimension = 0; dimension < num_dimensions; dimension++)
    baseband_size[dimension] =
      QccWAVWaveletDWTSubbandLength(field_size[dimension],
                                    subband_scale, 0,
                                    field_origin[dimension],
                                    0);
  
  if (!subband)
    {
      for (dimension = 0; dimension < num_dimensions; dimension++)
        {
          subband_offset[dimension] = 0;
          subband_size[dimension] = baseband_size[dimension];
        }
      return(0);
    }
  
  subband_code = ((subband - 1) % subband_constant) + 1;
  for (dimension = 0; dimension < num_dimensions; dimension++)
    {
      if (subband_code & 1)
        {
          subband_size[dimension] =
            QccWAVWaveletDWTSubbandLength(field_size[dimension],
                                          subband_scale, 1,
                                          field_origin[dimension],
                                          0);
          subband_offset[dimension] = baseband_size[dimension];
        }
      else
        {
          subband_size[dimension] =
            QccWAVWaveletDWTSubbandLength(field_size[dimension],
                                          subband_scale, 0,
                                          field_origin[dimension],
                                          0);
          subband_offset[dimension] = 0;
        }
      subband_code >>= 1;
    }
  
  return(0);
}


template <class T>
inline int evField<T>::NullSubband(const evAddress &subband_size) const
{
  unsigned int dimension;

  for (dimension = 0; dimension < subband_size.GetNumDimensions(); dimension++)
    if (subband_size.Get(dimension) <= 0)
      return(1);

  return(0);
}


template <class T>
int evField<T>::RasterScan(evField<T> &field,
                           const evAddress &field_origin) const
{
  evAddress subband_size(this->NumDimensions);
  evAddress subband_origin(this->NumDimensions);
  evAddress counter(this->NumDimensions);
  evAddress address(this->NumDimensions);
  unsigned int subband_scale;
  
  if (!this->IsAllocated())
    return(0);

  if (field.GetNumDimensions() != field_origin.GetNumDimensions())
    {
      ErrorHandler.ErrorMessage("evField::RasterScan(): field_origin does not mathc field");
      return(1);
    }

  if (field.SetSize(this->Size))
    {
      ErrorHandler.ErrorMessage("evField::RasterScan(): Error calling evField::SetSize()");
      return(1);
    }

  unsigned int index;
  unsigned int subband;
  for (index = 0, subband = 0;
       subband < this->NumSubbands; subband++)
    {
      if (this->SubbandInfo(subband,
                            this->NumScales,
                            this->Size,
                            field_origin,
                            subband_size,
                            subband_origin,
                            subband_scale))
        {
          ErrorHandler.ErrorMessage("evField::RasterScan(): Error calling evField::SubbandInfo()");
          return(1);
        }
      
      counter = 0;
      do
        {
          address = counter + subband_origin;
          
          if (!this->NullSubband(subband_size))
            field[index++] = (*this)[address];
        }
      while (!counter.Increment(subband_size));
    }
  
  return(0);
}


template <class T>
int evField<T>::InverseRasterScan(evField<T> &field,
                                  const evAddress &field_origin) const
{
  if (this->InverseRasterScan(field,
                              field_origin,
                              field.Size,
                              field.NumScales))
    {
      ErrorHandler.ErrorMessage("evField::InverseRasterScan(): Error calling evField::InverseRasterScan()");
      return(1);
    }

  return(0);
}


template <class T>
int evField<T>::InverseRasterScan(evField<T> &field,
                                  const evAddress &field_origin,
                                  const evAddress &full_size,
                                  const unsigned int max_num_scales) const
{
  unsigned int subband;
  unsigned int index;
  evAddress subband_size(field.NumDimensions);
  evAddress subband_origin(field.NumDimensions);
  evAddress counter(field.NumDimensions);
  evAddress address(field.NumDimensions);
  unsigned int subband_scale;
  
  if (!this->IsAllocated())
    return(0);

  if (field.GetNumDimensions() != field_origin.GetNumDimensions())
    {
      ErrorHandler.ErrorMessage("evField::InverseRasterScan(): field_origin does not match field");
      return(1);
    }
  
  for (index = 0, subband = 0; subband < field.NumSubbands; subband++)
    {
      if (this->SubbandInfo(subband,
                            max_num_scales,
                            full_size,
                            field_origin,
                            subband_size,
                            subband_origin,
                            subband_scale))
        {
          ErrorHandler.ErrorMessage("evField::InverseRasterScan(): Error calling evField::SubbandInfo()");
          return(1);
        }
      
      counter = 0;
      do
        {
          address = counter + subband_origin;
          
          if (!field.NullSubband(subband_size))
            field[address] = (*this)[index++];
        }
      while (!counter.Increment(subband_size));
    }
  
  return(0);
}


template <class T>
int evField<T>::GetScaleInfo(evFieldScaleInfo &scale_info,
                             const evAddress &field_origin) const
{
  evAddress subband_size(this->NumDimensions);
  evAddress subband_origin(this->NumDimensions);
  evAddress counter(this->NumDimensions);
  evAddress address(this->NumDimensions);
  unsigned int subband_scale;
  unsigned int subband;
  unsigned int start;
  unsigned int length;
  int scale = 0;
  int previous_scale = -1;

  for (subband = 0, start = 0, length = 0;
       subband < this->NumSubbands; subband++)
    {
      if (this->SubbandInfo(subband,
                            this->NumScales,
                            this->Size,
                            field_origin,
                            subband_size,
                            subband_origin,
                            subband_scale))
        {
          ErrorHandler.ErrorMessage("evField::RasterScan(): Error calling evField::SubbandInfo()");
          return(1);
        }

      scale = this->CalcScaleFromSubband(subband,
                                         this->NumScales,
                                         this->Size.GetNumDimensions());

      if (!subband)
        {
          length = subband_size.Product();
          scale_info.Append(scale + 1, start, length);
          start = length;
          length = 0;
          previous_scale = scale;
        }
      else
        if (scale != previous_scale)
          {
            scale_info.Append(previous_scale, start, length);
            
            start += length;
            length = subband_size.Product();
            previous_scale = scale;
          }
        else
          length += subband_size.Product();
    }
  
  scale_info.Append(scale, start, length);
  
  return(0);
}


template <class T>
int evField<T>::ExtractSubband(const unsigned int subband, evField<T> &field,
                               const evAddress &field_origin)
  const
{
  evAddress subband_size(this->NumDimensions);
  evAddress subband_origin(this->NumDimensions);
  evAddress counter(this->NumDimensions);
  evAddress address(this->NumDimensions);
  unsigned int subband_scale;

  if (!this->IsAllocated())
    return(0);

  if (this->GetNumDimensions() != field_origin.GetNumDimensions())
    {
      ErrorHandler.ErrorMessage("evField::ExtractSubband(): field_origin does not match field");
      return(1);
    }

  if (this->SubbandInfo(subband,
                        this->NumScales,
                        this->Size,
                        field_origin,
                        subband_size,
                        subband_origin,
                        subband_scale))
    {
      ErrorHandler.ErrorMessage("evField::ExtractSubband(): Error calling evField::SubbandInfo()");
      return(1);
    }
      
  if (this->NullSubband(subband_size))
    return(0);

  if (field.SetSize(subband_size))
    {
      ErrorHandler.ErrorMessage("evField::ExtractSubband(): Error calling evField::SetSize()");
      return(1);
    }

  unsigned int index = 0;
  counter = 0;
  do
    {
      address = counter + subband_origin;
      
      field[index++] = (*this)[address];
    }
  while (!counter.Increment(subband_size));
  
  field.NumScales = 0;
  field.NumSubbands = 1;

  return(0);
}


template <class T>
int evField<T>::PreserveBaseband(const evAddress &field_origin)
{
  unsigned int baseband_scale;

  if (!this->NumScales)
    return(0);

  if (this->GetNumDimensions() != field_origin.GetNumDimensions())
    {
      ErrorHandler.ErrorMessage("evField::PreserveBaseband(): field_origin does not mathc field");
      return(1);
    }

  evAddress baseband_size;
  evAddress baseband_origin;

  if (this->SubbandInfo(0,
                        this->NumScales,
                        this->Size,
                        field_origin,
                        baseband_size,
                        baseband_origin,
                        baseband_scale))
    {
      ErrorHandler.ErrorMessage("evField::PreserveBaseband(): Error calling evField::SubbandInfo()");
      return(1);
    }

  evField<T> baseband;
  if (this->ExtractSubband(0, baseband))
    {
      ErrorHandler.ErrorMessage("evField::PreserveBaseband(): Error calling evField::ExtractSubband()");
      return(1);
    }

  this->SetSize(baseband.Size);

  this->Copy(baseband);

  this->NumScales = 0;
  this->NumSubbands = 1;

  return(0);
}


template <class T>
void evField<T>::PrintSelf(const int tabbing) const
{
  std::ostream_iterator<T> out(std::cout, " ");

  if (tabbing != EVPRINT_NOFORMAT)
    {
      evPrintWithTab(tabbing + 2, "Size: ");
      this->Size.Print(EVPRINT_NOFORMAT);
    }

  if (tabbing != EVPRINT_NOFORMAT)
    evPrintWithTab(tabbing + 2, "Field: ");
  if (this->Field.empty())
    std::cout << "<NULL>" << std::endl;
  else
    std::copy(this->Field.begin(),
              this->Field.end(),
              out);

  if (tabbing != EVPRINT_NOFORMAT)
    std::cout << std::endl;
}


#endif