landsat.cpp

/* vim: set expandtab shiftwidth=2 softtabstop=2 tw=70: */

#include <Rcpp.h>
using namespace Rcpp;

// Cross-reference work:
// 1. update ../src/registerDynamicSymbol.c with an item for this
// 2. main code should use the autogenerated wrapper in ../R/RcppExports.R
//
// NB. breaking work up into two functions may slow things down but
// it makes it simpler to check.

// #define DEBUG // uncomment this for debugging info


// [[Rcpp::export]]
RawMatrix do_landsat_transpose_flip(RawMatrix m)
{
#ifdef DEBUG
  Rprintf("do_landsat_transpose_flip() start\n");
#endif
  int nrow = m.nrow();
  int ncol = m.ncol();
#ifdef DEBUG
  Rprintf("nrow=%d ncol=%d\n", nrow, ncol);
#endif
  RawMatrix res(ncol, nrow); // this is the *transpose* dimension, so ncol then nrow
  int nrow_res = ncol;
  int ncol_res = nrow;
  // Transpose
#ifdef DEBUG
  Rprintf("start transpose ...\n");
#endif
  for (int i = 0; i < nrow; i++) {
    for (int j = 0; j < ncol; j++) {
      res(j, i) = m(i, j);
    }
  }
#ifdef DEBUG
  Rprintf("... transpose done\n");
#endif
  // Flip in second dimension
  // ncol_res_half is to avoid undoing the flip by repeating it.
#ifdef DEBUG
  Rprintf("start flip ...\n");
#endif
  int ncol_res_half = (int)floor(ncol_res / 2.0);
  for (int i = 0; i < nrow_res; i++) {
    for (int j = 0; j < ncol_res_half; j++) {
      unsigned char tmp = res(i, j);
      res(i, j) = res(i, ncol_res-j-1);
      res(i, ncol_res-j-1) = tmp;
    }
  }
#ifdef DEBUG
  Rprintf("... end flip\n");
#endif
  return(res);
}

// [[Rcpp::export]]
List do_landsat_numeric_to_bytes(NumericMatrix m, IntegerVector bits)
{
#ifdef DEBUG
  Rprintf("do_landsat_numeric_to_bytes() start\n");
#endif

  int nrow = m.nrow();
  int ncol = m.ncol();
  int two_byte = bits[0] > 8;

  // The most-significant matrix is just 1x1 for 8-bit data
  RawMatrix lsb(nrow, ncol);
  RawMatrix msb(two_byte?nrow:1, two_byte?ncol:1);

  // Check endianness
  unsigned int x = 1;
  char *c = (char*) &x;
  int little_endian = (int)*c;
#ifdef DEBUG
  Rprintf("little_endian: %d\n", little_endian);
#endif
  // fill up arrays
  // No need to index by i and j here; this will speed up
  int n = nrow * ncol;
  if (two_byte) {
    if (little_endian) {
      unsigned int mij_int;
      unsigned char ms, ls;
#ifdef DEBUG
      Rprintf("little-endian two-byte nrow=%d ncol=%d\n", nrow, ncol);
#endif
      for (int i = 0; i < n; i++) {
        mij_int = (unsigned int)(65535*m[i]);
        ms = (mij_int & 0xFF00) >> 8;
        ls = mij_int & 0x00FF;
        //#ifdef DEBUG
        //        Rprintf("i %d, m: %f -> %d -> msb 0x%02x lsb 0x%02x (little endian two-byte)\n", i, m[i], mij_int, ms, ls);
        //#endif
        lsb[i] = ls;
        msb[i] = ms;
      }
    } else {
#ifdef DEBUG
      Rprintf("big-endian two-byte nrow=%d ncol=%d\n", nrow, ncol);
#endif

      // big endian below
      for (int i = 0; i < n; i++) {
        double mij = m[i];
        unsigned int mij_int = (unsigned int)(65535*mij);
        unsigned char ls = (mij_int & 0xFF00) >> 8;
        unsigned char ms = mij_int & 0x00FF;
#ifdef DEBUG
        Rprintf("i %d, m: %f -> %d -> msb 0x%02x lsb 0x%02x (big endian two-byte)\n", i, mij, mij_int, ms, ls);
#endif
        lsb[i] = ls;
        msb[i] = ms;
      }
    }
  } else {
    // !two_byte: msb is scalar 0; just fill up lsb
    if (little_endian) {
#ifdef DEBUG
      Rprintf("little-endian one-byte nrow=%d ncol=%d\n", nrow, ncol);
#endif
      for (int i = 0; i < n; i++) {
        double mij = m[i];
        unsigned int mij_int = (unsigned int)(255*mij);
        unsigned char ls = mij_int; // & 0xFF00;
        //#ifdef DEBUG
        //        Rprintf("i=%d   %f -> %d -> lsb 0x%02x (little endian one-byte)\n", i, mij, mij_int, ls);
        //#endif
        lsb[i] = ls;
      }
    } else {
      // big endian
#ifdef DEBUG
      Rprintf("big-endian one-byte nrow=%d ncol=%d\n", nrow, ncol);
#endif
      unsigned int mij_int;
      unsigned char ls;
      for (int i = 0; i < n; i++) {
        mij_int = (unsigned int)(255*m[i]);
        ls = mij_int; // & 0x00FF;
        //#ifdef DEBUG
        //        Rprintf("i=%d   %f -> %d -> lsb 0x%02x (big endian one-byte)\n", i, m[i], mij_int, ls);
        //#endif
        lsb[i] = ls;
      }
    }
  }
#ifdef DEBUG
  Rprintf("do_landsat_numeric_to_bytes() done\n");
#endif
  return(List::create(Named("lsb")=lsb, Named("msb")=msb));
}