vil/file_formats/vil_j2k_image.cxx

#include <algorithm>
#include <cmath>
#include <limits>
#include <cstdlib>
#include "vil_j2k_image.h"
//:
// \file
// vil_j2k: Written by Rob Radtke (rob@) and Harry Voorhees (hlv@) of
// Stellar Science Ltd. Co. (stellarscience.com) for
// Air Force Research Laboratory, 2005.
// Write capability added by J. Mundy, April 2009
// Do not remove the following notice
// Modifications approved for public Release, distribution unlimited
// DISTAR Case 14074
//

#include <NCSFile.h>

#include <cassert>
#ifdef _MSC_VER
#  include "vcl_msvc_warnings.h"
#endif
#include "vil/vil_memory_chunk.h"
#include "vil/vil_image_view.h"
#include "vil/vil_load.h"
#include "vil/vil_open.h"
#include "vil/vil_new.h"
#include "NCSJPCVilIOStream.h"
#include <NCSTypes.h>

// Fix me - need to add UINT64 and INT64 - JLM
// Also note, float and double are defined but not handled by SDK 3.3
// they are reinterpreted as INT32 and INT64
NCSFileBandInfo
bandInfo(const vil_pixel_format & vilType)
{
  NCSFileBandInfo info;
  switch (vil_pixel_format_component_format(vilType))
  {
    case VIL_PIXEL_FORMAT_UINT_32: {
      info.nBits = sizeof(vxl_uint_32) * 8;
      info.bSigned = std::numeric_limits<vxl_uint_32>::is_signed;
      info.szDesc = 0;
      return info;
    }
    case VIL_PIXEL_FORMAT_INT_32: {
      info.nBits = sizeof(vxl_int_32) * 8;
      info.bSigned = std::numeric_limits<vxl_int_32>::is_signed;
      info.szDesc = 0;
      return info;
    }
    case VIL_PIXEL_FORMAT_UINT_16: {
      info.nBits = sizeof(vxl_uint_16) * 8;
      info.bSigned = std::numeric_limits<vxl_uint_16>::is_signed;
      info.szDesc = 0;
      return info;
    }
    case VIL_PIXEL_FORMAT_INT_16: {
      info.nBits = sizeof(vxl_int_16) * 8;
      info.bSigned = std::numeric_limits<vxl_int_16>::is_signed;
      info.szDesc = 0;
      return info;
    }
    case VIL_PIXEL_FORMAT_BYTE: {
      info.nBits = sizeof(vxl_byte) * 8;
      info.bSigned = std::numeric_limits<vxl_byte>::is_signed;
      info.szDesc = 0;
      return info;
    }
    case VIL_PIXEL_FORMAT_SBYTE: {
      info.nBits = sizeof(vxl_sbyte) * 8;
      info.bSigned = std::numeric_limits<vxl_sbyte>::is_signed;
      info.szDesc = 0;
      return info;
    }
    case VIL_PIXEL_FORMAT_FLOAT: {
      info.nBits = sizeof(float) * 8;
      info.bSigned = std::numeric_limits<float>::is_signed;
      info.szDesc = 0;
      return info;
    }
    case VIL_PIXEL_FORMAT_DOUBLE: {
      info.nBits = sizeof(double) * 8;
      info.bSigned = std::numeric_limits<double>::is_signed;
      info.szDesc = 0;
      return info;
    }
    case VIL_PIXEL_FORMAT_BOOL:
    case VIL_PIXEL_FORMAT_COMPLEX_FLOAT:
    case VIL_PIXEL_FORMAT_COMPLEX_DOUBLE:
    case VIL_PIXEL_FORMAT_UINT_64:
    case VIL_PIXEL_FORMAT_INT_64:
    case VIL_PIXEL_FORMAT_UNKNOWN:
    default: {
      assert(0);
      info.nBits = 0;
      info.bSigned = false;
      info.szDesc = 0;
      return info;
    }
  }
}

//--------------------------------------------------------------------------------
// class vil_j2k_file_format

static char const j2k_string[] = "j2k";

char const *
vil_j2k_file_format::tag() const
{
  return j2k_string;
}

vil_image_resource_sptr
vil_j2k_file_format::make_input_image(vil_stream * vs)
{
  vil_j2k_image * im = new vil_j2k_image(vs);
  if (!im->is_valid())
  {
    delete im;
    im = 0;
  }
  return im;
}


vil_image_resource_sptr
vil_j2k_file_format::make_output_image(vil_stream * vs,
                                       unsigned ni,
                                       unsigned nj,
                                       unsigned nplanes,
                                       enum vil_pixel_format format)
{
  vil_j2k_image * j2k_img = new vil_j2k_image(vs, ni, nj, nplanes, format, compression_ratio_);
  if (j2k_img->is_valid())
    return static_cast<vil_image_resource *>(j2k_img);
  return 0;
}

NCSEcwCellType
convertType(const vil_pixel_format & vilType)
{
  switch (vil_pixel_format_component_format(vilType))
  {
    case VIL_PIXEL_FORMAT_UINT_64:
      return NCSCT_UINT64;
    case VIL_PIXEL_FORMAT_INT_64:
      return NCSCT_INT64;
    case VIL_PIXEL_FORMAT_UINT_32:
      return NCSCT_UINT32;
    case VIL_PIXEL_FORMAT_INT_32:
      return NCSCT_INT32;
    case VIL_PIXEL_FORMAT_UINT_16:
      return NCSCT_UINT16;
    case VIL_PIXEL_FORMAT_INT_16:
      return NCSCT_INT16;
    case VIL_PIXEL_FORMAT_BYTE:
      return NCSCT_UINT8;
    case VIL_PIXEL_FORMAT_SBYTE:
      return NCSCT_INT8;
    case VIL_PIXEL_FORMAT_FLOAT:
      return NCSCT_IEEE4;
    case VIL_PIXEL_FORMAT_DOUBLE:
      return NCSCT_IEEE8;
    case VIL_PIXEL_FORMAT_BOOL:
    case VIL_PIXEL_FORMAT_COMPLEX_FLOAT:
    case VIL_PIXEL_FORMAT_COMPLEX_DOUBLE:
    case VIL_PIXEL_FORMAT_UNKNOWN:
    default:
      assert(0);
      return NCSCT_UINT8;
  }
}

// Note the J2K SDK defines IEEE4 and IEEE8 but they are not
// handled as float or double, they are reinterpreted as INT32 and INT64
vil_pixel_format
convertType(const NCSEcwCellType & ecwType)
{
  switch (ecwType)
  {
    case NCSCT_UINT64:
      return VIL_PIXEL_FORMAT_UINT_64;
    case NCSCT_INT64:
      return VIL_PIXEL_FORMAT_INT_64;
    case NCSCT_UINT32:
      return VIL_PIXEL_FORMAT_UINT_32;
    case NCSCT_INT32:
      return VIL_PIXEL_FORMAT_INT_32;
    case NCSCT_UINT16:
      return VIL_PIXEL_FORMAT_UINT_16;
    case NCSCT_INT16:
      return VIL_PIXEL_FORMAT_INT_16;
    case NCSCT_UINT8:
      return VIL_PIXEL_FORMAT_BYTE;
    case NCSCT_INT8:
      return VIL_PIXEL_FORMAT_SBYTE;
    case NCSCT_IEEE4:
      return VIL_PIXEL_FORMAT_FLOAT;
    case NCSCT_IEEE8:
      return VIL_PIXEL_FORMAT_DOUBLE;
    default:
      assert(0);
      return VIL_PIXEL_FORMAT_UNKNOWN;
  }
}

////////////////////////////////////////////////////
//                vil_j2k_image
////////////////////////////////////////////////////

vil_j2k_image::vil_j2k_image(const std::string & fileOrUrl)
  : vil_image_resource()
  , mFileResource(new CNCSFile())
  , mStr(0)
  , mMaxLocalDimension(5000)
  , // default value
  mMaxRemoteDimension(640)
  , // default value
  mRemoteFile(false)
  , mFinfo(0)
  , mBandinfo(0)
  , line_index_(0)
{
  if (mFileResource->Open((char *)fileOrUrl.c_str(), false, false) != NCS_SUCCESS)
  {
    mFileResource = 0;
    return;
  }
  if (fileOrUrl.substr(0, 7) == "ecwp://" || fileOrUrl.substr(0, 7) == "ECWP://")
  {
    mRemoteFile = true;
  }
}

vil_j2k_image::vil_j2k_image(vil_stream * is)
  : vil_image_resource()
  , mFileResource(new CNCSFile())
  , mStr(new CNCSJPCVilIOStream())
  , mMaxLocalDimension(5000)
  , // default value
  mMaxRemoteDimension(640)
  , // default value
  mRemoteFile(false)
  , mFinfo(0)
  , mBandinfo(0)
  , line_index_(0)
{
  mStr->Open(is);

  if ((static_cast<CNCSJP2FileView *>(mFileResource))->Open(mStr) != NCS_SUCCESS)
  {
    mFileResource = 0;
    return;
  }
}

vil_j2k_image::vil_j2k_image(vil_stream * vs,
                             unsigned ni,
                             unsigned nj,
                             unsigned nplanes,
                             enum vil_pixel_format format,
                             unsigned compression_ratio)
  : vil_image_resource()
  , mFileResource(new CNCSFile())
  , mStr(new CNCSJPCVilIOStream())
  , mMaxLocalDimension(5000)
  , // default value
  mMaxRemoteDimension(640)
  , // default value
  mRemoteFile(false)
  , mFinfo(new NCSFileViewFileInfoEx())
  , line_index_(0)
{
  mBandinfo = new NCSFileBandInfo[nplanes];
  CNCSError Error;
  for (unsigned i = 0; i < nplanes; ++i)
    mBandinfo[i] = bandInfo(format);

  // String names for each band, should specialize according to color space
  mBandinfo[0].szDesc = "grey";
  if (nplanes == 3)
  {
    mBandinfo[0].szDesc = "Red";
    mBandinfo[1].szDesc = "Green";
    mBandinfo[2].szDesc = "Blue";
  }
  NCSEcwCellType t = convertType(format);
  NCSFileViewFileInfoEx finfo = *mFinfo;
  finfo.pBands = mBandinfo;
  finfo.nSizeX = ni;
  finfo.nSizeY = nj;
  finfo.nBands = nplanes;
  finfo.eCellType = t;
  finfo.nCompressionRate = compression_ratio;
  finfo.eCellSizeUnits = ECW_CELL_UNITS_METERS;
  finfo.fCellIncrementX = 1.0;
  finfo.fCellIncrementY = 1.0;
  finfo.fOriginX = 0.0;
  finfo.fOriginY = 0.0;
  finfo.szDatum = "RAW";
  finfo.szProjection = "RAW";
  finfo.fCWRotationDegrees = 0.0;
  if (nplanes == 1)
    finfo.eColorSpace = NCSCS_GREYSCALE;
  else if (nplanes == 3)
    finfo.eColorSpace = NCSCS_sRGB;
  else
  {
    delete mFileResource;
    mFileResource = 0;
  }
  Error = mFileResource->SetFileInfo(finfo);
  if (Error != NCS_SUCCESS)
  {
    if (mFileResource)
      delete mFileResource;
    mFileResource = 0;
  }
  Error = mStr->Open(vs, true);
  if (Error != NCS_SUCCESS)
  {
    if (mFileResource)
      delete mFileResource;
    mFileResource = 0;
  }
  CNCSJP2FileView * fview = static_cast<CNCSJP2FileView *>(mFileResource);
  Error = fview->Open(static_cast<CNCSJPCIOStream *>(mStr));
  if (Error != NCS_SUCCESS)
  {
    if (mFileResource)
      delete mFileResource;
    mFileResource = 0;
    return;
  }
}

vil_j2k_image::~vil_j2k_image()
{
  if (mFileResource)
  {
    mFileResource->Close(true);
  }
  if (mStr)
    delete mStr;
  if (mBandinfo)
    delete mBandinfo;
  if (mFinfo)
    delete mFinfo;
}

unsigned
vil_j2k_image::nplanes() const
{
  assert(mFileResource);
  return mFileResource->GetFileInfo()->nBands;
}

unsigned
vil_j2k_image::ni() const
{
  assert(mFileResource);
  return mFileResource->GetFileInfo()->nSizeX;
}

unsigned
vil_j2k_image::nj() const
{
  assert(mFileResource);
  return mFileResource->GetFileInfo()->nSizeY;
}

enum vil_pixel_format
vil_j2k_image::pixel_format() const
{
  assert(mFileResource);
  return convertType(mFileResource->GetFileInfo()->eCellType);
}

char const *
vil_j2k_image::file_format() const
{
  return "j2k";
}

vil_image_view_base_sptr
vil_j2k_image::get_copy_view_decimated(unsigned sample0,
                                       unsigned num_samples,
                                       unsigned line0,
                                       unsigned numLines,
                                       double i_factor,
                                       double j_factor) const
{
  return get_copy_view_decimated_by_size(sample0,
                                         num_samples,
                                         line0,
                                         numLines,
                                         (unsigned int)(((double)num_samples) / i_factor),
                                         (unsigned int)(((double)numLines) / j_factor));
}

vil_image_view_base_sptr
vil_j2k_image::get_copy_view_decimated_by_size(unsigned sample0,
                                               unsigned num_samples,
                                               unsigned line0,
                                               unsigned numLines,
                                               unsigned int output_width,
                                               unsigned int output_height) const
{
  if (!(mFileResource) || !((sample0 + num_samples - 1) < ni() && (line0 + numLines - 1) < nj()))
  {
    return 0;
  }

  // we want all bands mapped in the same order as they come in the input file
  // eg. bandMap = {0,1,2,3...nBands}
  INT32 nBands = nplanes();
  INT32 * bandMap = (INT32 *)std::malloc(sizeof(UINT32) * nBands);
  for (int i = 0; i < nBands; i++)
  {
    bandMap[i] = i;
  }

  // this guards us from returning an image that is too big for the computer's memory
  //(or would take too long to download in the remote case).
  // We don't want infinite hangs or application crashes.
  unsigned int maxDim = mRemoteFile ? mMaxRemoteDimension : mMaxLocalDimension;
  if (output_width > maxDim || output_height > maxDim)
  {
    unsigned int biggestDim = (std::max)(output_width, output_height);
    double zoomFactor = ((double)maxDim) / ((double)biggestDim);
    output_width = (unsigned int)(((double)output_width) * zoomFactor);
    output_height = (unsigned int)(((double)output_height) * zoomFactor);
  }

  // set the view to be that specified by the function's input parameters
  // note that we don't want ECW to do any scaling for us.  That's why
  // the box created by (sample0,line0) and (sample0+num_samples-1,line0+numLines-1) is made to be exactly
  // num_samplesXnumLines.

  NCSError setViewError = mFileResource->SetView(nBands,
                                                 bandMap,
                                                 output_width,
                                                 output_height,
                                                 (INT32)sample0,
                                                 (INT32)line0,
                                                 (INT32)(sample0 + num_samples - 1),
                                                 (INT32)(line0 + numLines - 1));
  if (setViewError != NCS_SUCCESS)
  {
    free(bandMap);
    return 0;
  }

  // number of samples times the bytes per sample in each band
  double bitsPerSample = mFileResource->GetFileInfo()->pBands[0].nBits;
  unsigned int bytesPerSample = (unsigned int)std::ceil(bitsPerSample / 8.0);
  unsigned int singleBandLineSizeBytes = output_width * bytesPerSample;
  unsigned int allBandLineSizeBytes = singleBandLineSizeBytes * nBands;
  unsigned int dataPtrSizeBytes = allBandLineSizeBytes * output_height;
  // void* data_ptr = std::malloc( dataPtrSizeBytes );
  vil_memory_chunk_sptr data_ptr =
    new vil_memory_chunk(dataPtrSizeBytes, convertType(mFileResource->GetFileInfo()->eCellType));
  void ** linePtrPtr =
    (void **)std::malloc(nBands * sizeof(int * /*all pointers have same size, so eg char* would work too*/));
  // now read all the lines that we want
  for (unsigned int currLine = 0; currLine < output_height; currLine++)
  {
    for (int currBand = 0; currBand < nBands; currBand++)
    {
      linePtrPtr[currBand] =
        (void *)(((char *)data_ptr->data()) + currLine * allBandLineSizeBytes + currBand * singleBandLineSizeBytes);
    }
    NCSEcwReadStatus readStatus =
      mFileResource->ReadLineBIL(mFileResource->GetFileInfo()->eCellType, nBands, linePtrPtr, 0);
    if (readStatus != NCSECW_READ_OK)
    {
      free(bandMap);
      free(linePtrPtr);
      return 0;
    }
  }

  // free our temp resources
  free(bandMap);
  free(linePtrPtr);

  vil_image_view_base_sptr view = 0;
  // now create our vil_image_view
  // note that float and double are defaulted since the J2K SDK doesn't
  // implement these types properly
  switch (vil_pixel_format_component_format(data_ptr->pixel_format()))
  {
#define macro(F, T)                                                                                                    \
  case F:                                                                                                              \
    view = new vil_image_view<T>(data_ptr,                                                                             \
                                 reinterpret_cast<T *>(data_ptr->data()),                                              \
                                 output_width,                                                                         \
                                 output_height,                                                                        \
                                 nBands,                                                                               \
                                 1,                                                                                    \
                                 output_width * nBands,                                                                \
                                 output_width);                                                                        \
    break
    macro(VIL_PIXEL_FORMAT_BYTE, vxl_byte);
    macro(VIL_PIXEL_FORMAT_SBYTE, vxl_sbyte);
    macro(VIL_PIXEL_FORMAT_UINT_64, vxl_uint_64);
    macro(VIL_PIXEL_FORMAT_INT_64, vxl_int_64);
    macro(VIL_PIXEL_FORMAT_UINT_32, vxl_uint_32);
    macro(VIL_PIXEL_FORMAT_INT_32, vxl_int_32);
    macro(VIL_PIXEL_FORMAT_UINT_16, vxl_uint_16);
    macro(VIL_PIXEL_FORMAT_INT_16, vxl_int_16);
    macro(VIL_PIXEL_FORMAT_BOOL, bool);
#undef macro
    default:
      std::cerr << "Pixel format not supported by ERMapper SDK\n";
      assert(0);
      break;
  }

  return view;
}

vil_image_view_base_sptr
vil_j2k_image::get_copy_view(unsigned sample0, unsigned num_samples, unsigned line0, unsigned numLines) const
{
  return get_copy_view_decimated(sample0, num_samples, line0, numLines, 1.0, 1.0);
}

void
vil_j2k_image::unsetMaxImageDimension(bool remote)
{
#undef max
  setMaxImageDimension(std::numeric_limits<unsigned int>::max(), remote);
}

void
vil_j2k_image::setMaxImageDimension(unsigned int widthOrHeight, bool remote)
{
  if (remote)
  {
    mMaxRemoteDimension = widthOrHeight;
  }
  else
  {
    mMaxLocalDimension = widthOrHeight;
  }
}

vil_image_view_base_sptr
vil_j2k_image::s_decode_jpeg_2000(vil_stream * vs,
                                  unsigned i0,
                                  unsigned ni,
                                  unsigned j0,
                                  unsigned nj,
                                  double i_factor,
                                  double j_factor)
{
  vil_j2k_image * j2k_image = new vil_j2k_image(vs);
  // remove limit by default, since vil is not typically used remotely
  // but more commonly with large image files - JLM Jan 07, 2012
  j2k_image->unsetMaxImageDimension();
  vil_image_view_base_sptr view = j2k_image->get_copy_view_decimated(i0, ni, j0, nj, i_factor, j_factor);
  delete j2k_image;
  return view;
}

vil_image_view_base_sptr
vil_j2k_image::s_decode_jpeg_2000_by_size(vil_stream * vs,
                                          unsigned i0,
                                          unsigned ni,
                                          unsigned j0,
                                          unsigned nj,
                                          unsigned int output_width,
                                          unsigned int output_height)
{
  vil_j2k_image * j2k_image = new vil_j2k_image(vs);
  vil_image_view_base_sptr view =
    j2k_image->get_copy_view_decimated_by_size(i0, ni, j0, nj, output_width, output_height);
  delete j2k_image;
  return view;
}

template <class T>
static bool
write_line_BIL(vil_memory_chunk_sptr & chunk,
               unsigned ni,
               unsigned nplanes,
               unsigned istep,
               unsigned planestep,
               unsigned bytes_per_pixel,
               CNCSFile * f_resource,
               NCSEcwCellType t)
{
  T * cdata = reinterpret_cast<T *>(chunk->data());
  T ** line_ptr = new T *[nplanes];
  for (unsigned p = 0; p < nplanes; ++p)
    line_ptr[p] = new T[ni * bytes_per_pixel];

  for (unsigned p = 0; p < nplanes; ++p)
  {
    T * wline = line_ptr[p];
    for (unsigned i = 0; i < ni; ++i)
    {
      *(wline + i) = *(cdata + i * istep + p * planestep);
    }
  }
  bool good = true;
  void ** outbuf = reinterpret_cast<void **>(line_ptr);
  CNCSError writeError = f_resource->WriteLineBIL(t, nplanes, outbuf);
  if (writeError != NCS_SUCCESS)
    good = false;

  for (unsigned p = 0; p < nplanes; ++p)
    delete[] line_ptr[p];
  delete[] line_ptr;
  return good;
}


//: JPEG2K compress by inserting an image row (line) at a time
//  When the full image has been inserted, call put_line with
// image_row == nj(). This call causes the resource to be closed
// and is no longer valid. The lines must be inserted in strict row order.
bool
vil_j2k_image::put_line(const vil_image_view_base & im)
{
  if (!mFileResource)
    return false;
  vil_pixel_format format = this->pixel_format();
  unsigned ni = this->ni(), nj = this->nj(), nplanes = this->nplanes();
  if (line_index_ >= nj)
  {
    mFileResource->Close(true);
    if (mFileResource)
      delete mFileResource;
    mFileResource = 0;
    return true;
  }
  unsigned bytes_per_pixel = 0;
  NCSEcwCellType t = convertType(format);
  vil_memory_chunk_sptr chunk;
  // now write out the image line
  // note that float and double are defaulted since the J2K SDK doesn't
  // implement these types properly
  switch (vil_pixel_format_component_format(format))
  {
#define macro(F, T)                                                                                                    \
  case F: {                                                                                                            \
    bytes_per_pixel = sizeof(T);                                                                                       \
    const vil_image_view<T> & view = static_cast<const vil_image_view<T> &>(im);                                       \
    chunk = view.memory_chunk();                                                                                       \
    if (!write_line_BIL<T>(chunk, ni, nplanes, view.istep(), view.planestep(), bytes_per_pixel, mFileResource, t))     \
      return false;                                                                                                    \
  }                                                                                                                    \
  break
    macro(VIL_PIXEL_FORMAT_BYTE, vxl_byte);
    macro(VIL_PIXEL_FORMAT_SBYTE, vxl_sbyte);
    macro(VIL_PIXEL_FORMAT_UINT_32, vxl_uint_32);
    macro(VIL_PIXEL_FORMAT_INT_32, vxl_int_32);
    macro(VIL_PIXEL_FORMAT_UINT_16, vxl_uint_16);
    macro(VIL_PIXEL_FORMAT_INT_16, vxl_int_16);
    macro(VIL_PIXEL_FORMAT_BOOL, bool);
#undef macro
    default:
      std::cerr << "Pixel format not supported by ERMapper SDK\n";
      assert(0);
      break;
  }
  ++line_index_;
  return true;
}

//: JPEG2K compress the view and of the full image and insert in resource
//  The file is closed after putting the view into the resource
//  and becomes an invalid resource.
bool
vil_j2k_image::put_view(const vil_image_view_base & im)
{
  if (!this->view_fits(im, 0, 0))
    return false;
  if (!mFileResource)
    return false;
  unsigned ni = im.ni(), nj = im.nj();
  vil_image_resource_sptr mem_res = vil_new_image_resource_of_view(im);
  vil_image_view_base_sptr view;
  for (unsigned j = 0; j < nj; ++j)
  {
    view = mem_res->get_copy_view(0, ni, j, 1);
    if (!this->put_line(*view))
      return false;
  }
  return true;
}

//: Check that a view will fit into the data at the given offset.
bool
vil_j2k_image::view_fits(const vil_image_view_base & im, unsigned i0, unsigned j0)
{
  unsigned ni_view = im.ni(), nj_view = im.nj(), nplanes_view = im.nplanes();
  return i0 + 1 < ni_view && j0 + 1 < nj_view && ni_view <= this->ni() && nj_view <= this->nj() &&
         nplanes_view <= this->nplanes();
}

//:
//  Encode an entire image by loading the input resource from stream
//  and compressing the input line by line by extracting an image view
//  of a block of lines at a time, thus works for arbitrarily large images.
//  The num_lines_block parameter is the number of image rows in the
//  block which is read into memory from the resource
bool
vil_j2k_image::s_encode_jpeg2000(vil_stream * vs,
                                 const char * out_filename,
                                 unsigned compression_ratio,
                                 unsigned num_lines_block,
                                 bool verbose)
{
  vil_image_resource_sptr in_res = vil_load_image_resource_raw(vs);
  if (!in_res)
    return false;
  unsigned ni = in_res->ni(), nj = in_res->nj(), nplanes = in_res->nplanes();
  vil_pixel_format format = in_res->pixel_format();
  vil_stream * os = vil_open(out_filename, "w");
  if (!vs)
    return false;
  vil_j2k_file_format fmt;
  fmt.set_compression_ratio(compression_ratio);
  vil_image_resource_sptr res = fmt.make_output_image(os, ni, nj, nplanes, format);
  if (!res)
    return false;
  vil_j2k_image * j2k_img = reinterpret_cast<vil_j2k_image *>(res.ptr());

  // number of full blocks in image height
  unsigned n_blocks = nj / num_lines_block;
  unsigned jb = 0;
  for (unsigned b = 0; b < n_blocks; b++, jb += num_lines_block)
  {
    // read a block from the file: width = ni, height = num_lines_block
    vil_image_view_base_sptr block_view = in_res->get_view(0, ni, jb, num_lines_block);
    if (!block_view)
      return false;

    // wrap the view in a memory resident resource
    vil_image_resource_sptr block_res = vil_new_image_resource_of_view(*block_view);

    // compress the block, line by line
    for (unsigned j = 0; j < num_lines_block; ++j)
    {
      vil_image_view_base_sptr line_view = block_res->get_copy_view(0, ni, j, 1);
      if (!j2k_img->put_line(*line_view))
        return false;
      if (verbose)
        if (j % 100 == 0) // output a dot every 100 lines
          std::cout << '.';
    }
  }
  // output the remaining lines left over after loading block-sized chunks
  unsigned remaining_lines = nj - jb;
  if (remaining_lines)
  {
    vil_image_view_base_sptr residual_view = in_res->get_view(0, ni, jb, remaining_lines);
    vil_image_resource_sptr residual_res = vil_new_image_resource_of_view(*residual_view);
    vil_image_view_base_sptr view;
    for (unsigned j = 0; j < remaining_lines; ++j)
    {
      view = residual_res->get_copy_view(0, ni, j, 1);
      if (!j2k_img->put_line(*view))
        return false;
      if (verbose)
        if (j % 100 == 0) // output a dot every 100 lines
          std::cout << '.';
    }
  }
  if (verbose)
    std::cout << '\n';
  return true;
}

bool
vil_j2k_image::s_encode_jpeg2000(const char * in_filename,
                                 const char * out_filename,
                                 unsigned compression_ratio,
                                 unsigned num_lines_block,
                                 bool verbose)
{
  vil_stream * vs = vil_open(in_filename);
  vs->ref();
  bool success = vil_j2k_image::s_encode_jpeg2000(vs, out_filename, compression_ratio, num_lines_block, verbose);
  vs->unref();
  return success;
}