lsgrmSplitter.txx 7.58 KB
Newer Older
remicres's avatar
remicres committed
1
#include "lsgrmSplitter.h"
2 3 4

namespace lsgrm
{
remicres's avatar
remicres committed
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240

template <class TInputImage>
void SplitOTBImage(TInputImage * imagePtr, // input image
    unsigned int& tileWidth, // width of the tile
    unsigned int& tileHeight, // height of the tile
    const unsigned int margin, // stability margin
    const unsigned int niter,
    unsigned int &nbTilesX,
    unsigned int &nbTilesY,
    std::vector<ProcessingTile> &m_Tiles)
{

  // Image size
  auto imageWidth = imagePtr->GetLargestPossibleRegion().GetSize()[0];
  auto imageHeight = imagePtr->GetLargestPossibleRegion().GetSize()[1];

  // Find the largest tile with the lowest compactness
  long long unsigned int maximumNbOfPixels = tileWidth * tileHeight;
  unsigned int k = 1;
  unsigned int lowestCompactnessNCol = 1;
  unsigned int layoutNCol(1), layoutNRow(1);
  while(imageWidth*imageHeight / ((float) k) > maximumNbOfPixels)
    {
    // Get the multiples of k. For each one, compute the
    // layout which has the minimum tile compactness
    float lowestCompactness = itk::NumericTraits<float>::max();
    for (layoutNCol = 1; layoutNCol<=k; layoutNCol++)
      {
      if (k % layoutNCol == 0)      // Is it a multiple of k?
        {
        // Compute tile compactness
        layoutNRow = k / layoutNCol;
        tileWidth = imageWidth / layoutNCol;
        tileHeight = imageHeight / layoutNRow;
        float perimeter = tileWidth + tileHeight;
        float surface = tileWidth * tileHeight;
        float compactness = perimeter / surface;

        // Update minimum compactness
        if (lowestCompactness > compactness)
          {
          lowestCompactness = compactness;
          lowestCompactnessNCol = layoutNCol;
          }
        }
      } // for each multiple of k
    k++;
    } // while

  tileWidth = imageWidth / lowestCompactnessNCol;
  tileHeight = imageHeight / (k / lowestCompactnessNCol);
  std::cout << "New tiling scheme:"
      << "\n\t\size (px): " << tileWidth << " x " << tileHeight
      << "\n\t\layout : " << lowestCompactnessNCol << " x " << (k / lowestCompactnessNCol) << std::endl;

  //  if(tileWidth + 2*margin >= imageWidth || tileHeight + 2* margin >= imageHeight)
  //    {
  //    if(tileWidth + 2*margin >= imageWidth && tileHeight + 2* margin >= imageHeight)
  //      {
  //      std::cout << "Both dimensions of the tile is greater than the image, classic GRM should be used" << std::endl;
  //      exit(1);
  //      }
  //    long long maximumNumberOfPixels = (tileWidth + 2*margin)*(tileHeight + 2* margin);
  //    if(tileWidth + 2*margin >= imageWidth)
  //      {
  //      tileWidth = imageWidth;
  //      tileHeight = std::min( (float) imageHeight, (float) std::floor(maximumNumberOfPixels / tileWidth));
  //      }
  //    else if (tileHeight + 2* margin >= imageHeight)
  //      {
  //      tileHeight = imageHeight;
  //      tileWidth = std::min( (float) imageWidth,  (float) std::floor(maximumNumberOfPixels / tileHeight));
  //      }
  //    std::cout << "Changing tiles size to " << tileWidth << " x " << tileHeight << std::endl;
  //    }

  //  // Determine final size of the tiles according to the dimension of the image
  //  unsigned int k = std::ceil(imageWidth / tileWidth);
  //  std::cout << "k = " << k << std::endl;
  //  while(imageWidth % k)
  //    k++;
  //  tileWidth = imageWidth / k;
  //  assert(tileWidth > 0);
  //  k = std::ceil(imageHeight / tileHeight);
  //  std::cout << "k = " << k << std::endl;
  //  while(imageWidth % k)
  //    k++;
  //  tileHeight = imageHeight / k;
  //  assert(tileHeight > 0);
  //
  //  std::cout << "Final tiles size is " << tileWidth << " x " << tileHeight << std::endl;

  /* Determine the number of tiles according to the tile dimension given by the user. */
  nbTilesX = imageWidth / tileWidth + ( imageWidth % tileWidth > 0 ? 1 : 0);
  nbTilesY = imageHeight / tileHeight + ( imageHeight % tileHeight > 0 ? 1 : 0);

  /* Local variables for the next loop */
  unsigned int startX, startY; // Upper left coordinates of the tile.
  unsigned int sizeX, sizeY; // Size of the tiles.

  /* Loop over the tiles*/
  m_Tiles.assign(nbTilesX * nbTilesY, ProcessingTile());
  unsigned int i = 0;
  for(unsigned int row = 0; row < nbTilesY; ++row)
    {
    for(unsigned int col = 0; col < nbTilesX; ++col)
      {

      startX = col * tileWidth;
      startY = row * tileHeight;
      sizeX = std::min(tileWidth, static_cast<unsigned int>(imageWidth - startX));
      sizeY = std::min(tileHeight, static_cast<unsigned int>(imageHeight - startY));

      /* Margin at the top ? */
      if( row > 0 )
        {
        startY -= margin;
        sizeY += margin;
        }

      /* Margin at the right */
      if( col < nbTilesX - 1 )
        {
        sizeX += margin;
        }

      /* Margin at the bottom */
      if( row < nbTilesY - 1)
        {
        sizeY += margin;
        }

      /* Margin at the left */
      if( col > 0 )
        {
        startX -= margin;
        sizeX += margin;
        }

      /* Store the tile region */
      typename TInputImage::IndexType index;
      index[0] = startX;
      index[1] = startY;
      typename TInputImage::SizeType size;
      size[0] = sizeX;
      size[1] = sizeY;
      typename TInputImage::RegionType region(index, size);
      m_Tiles[i].region = region;

      /* Margin at the top ? */
      if( row > 0 )
        {
        m_Tiles[i].rows[0] = row * tileHeight;
        m_Tiles[i].tileNeighbors[0] = i - nbTilesX;
        m_Tiles[i].margin[0] = true;
        }
      else
        {
        m_Tiles[i].rows[0] = 0;
        m_Tiles[i].tileNeighbors[0] = -1;
        m_Tiles[i].margin[0] = false;
        }


      /* Margin at the right */
      if( col < nbTilesX - 1 )
        {
        m_Tiles[i].columns[1] = col * tileWidth + sizeX - 1; //sizeX
        m_Tiles[i].tileNeighbors[2] = i+1;
        m_Tiles[i].margin[1] = true;
        }
      else
        {
        m_Tiles[i].columns[1] = imageWidth - 1;
        m_Tiles[i].tileNeighbors[2] = -1;
        m_Tiles[i].margin[1] = false;
        }

      /* Margin at the bottom */
      if( row < nbTilesY - 1)
        {
        m_Tiles[i].rows[1] = row * tileHeight + sizeY - 1; // sizeY
        m_Tiles[i].tileNeighbors[4] = i + nbTilesX;
        m_Tiles[i].margin[2] = true;
        }
      else
        {
        m_Tiles[i].rows[1] = imageHeight - 1;
        m_Tiles[i].tileNeighbors[4] = -1;
        m_Tiles[i].margin[2] = false;
        }

      /* Margin at the left */
      if( col > 0 )
        {
        m_Tiles[i].columns[0] = col * tileWidth;
        m_Tiles[i].tileNeighbors[6] = i-1;
        m_Tiles[i].margin[3] = true;
        }
      else
        {
        m_Tiles[i].columns[0] = 0;
        m_Tiles[i].tileNeighbors[6] = -1;
        m_Tiles[i].margin[3] = false;
        }

      /* Is there a neighbor at the rop right */
      if(row > 0 && col < nbTilesX - 1)
        m_Tiles[i].tileNeighbors[1] = i - nbTilesX + 1;
      else
        m_Tiles[i].tileNeighbors[1] = -1;

      /* Is there a neighbor at the bottom right */
      if(col < nbTilesX - 1 && row < nbTilesY - 1)
        m_Tiles[i].tileNeighbors[3] = i + nbTilesX + 1;
      else
        m_Tiles[i].tileNeighbors[3] = -1;

      /* Is there a neighbor at the bottom left */
      if(row < nbTilesY - 1 && col > 0)
        m_Tiles[i].tileNeighbors[5] = i + nbTilesX - 1;
      else
        m_Tiles[i].tileNeighbors[5] = -1;

      /* Is there a neighbor at the top left */
      if(col > 0 && row > 0)
        m_Tiles[i].tileNeighbors[7] = i - nbTilesX - 1;
      else
        m_Tiles[i].tileNeighbors[7] = -1;

      i++;
      } // end for(unsigned int col = 0; col < nbTilesX; ++col)

    } // for(unsigned int row = 0; row < nbTilesY; ++row)

}
241
} // end of namespace lsgrm