<commit_before>// __BEGIN_LICENSE__
// 
// Copyright (C) 2008 United States Government as represented by the
// Administrator of the National Aeronautics and Space Administration
// (NASA).  All Rights Reserved.
// 
// Copyright 2008 <NAME>. All rights reserved.
// 
// This software is distributed under the NASA Open Source Agreement
// (NOSA), version 1.3.  The NOSA has been approved by the Open Source
// Initiative.  See the file COPYING at the top of the distribution
// directory tree for the complete NOSA document.
// 
// THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY
// KIND, EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT
// LIMITED TO, ANY WARRANTY THAT THE SUBJECT SOFTWARE WILL CONFORM TO
// SPECIFICATIONS, ANY IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR
// A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY WARRANTY THAT
// THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT
// DOCUMENTATION, IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE.
//
// __END_LICENSE__

/// \file StereoSessionIsis.cc
///

#include <boost/shared_ptr.hpp>

#include <vw/Image.h>
#include <vw/Math.h>
#include <vw/FileIO.h>
#include <vw/InterestPoint.h>
#include <vw/Stereo.h>
#include <vw/Cartography.h>
#include "Isis/StereoSessionIsis.h"
#include "Isis/IsisCameraModel.h"
#include "StereoSettings.h"

// Boost
#include "boost/filesystem.hpp"   
using namespace boost::filesystem; 

using namespace vw;
using namespace vw::camera;
using namespace vw::cartography;
using namespace vw::stereo;
using namespace vw::ip;


// Allows FileIO to correctly read/write these pixel types
namespace vw {
  template<> struct PixelFormatID<Vector3>   { static const PixelFormatEnum value = VW_PIXEL_GENERIC_3_CHANNEL; };
  template<> struct PixelFormatID<PixelDisparity<float> >   { static const PixelFormatEnum value = VW_PIXEL_GENERIC_3_CHANNEL; };
}


static std::string prefix_from_filename(std::string const& filename) {
  std::string result = filename;
  int index = result.rfind(".");
  if (index != -1) 
    result.erase(index, result.size());
  return result;
}

// Duplicate matches for any given interest point probably indicate a
// poor match, so we cull those out here.
static void remove_duplicates(std::vector<Vector3> &ip1, std::vector<Vector3> &ip2) {
  std::vector<Vector3> new_ip1, new_ip2;
  
  for (unsigned i = 0; i < ip1.size(); ++i) {
    bool bad_entry = false;
    for (unsigned j = 0; j < ip1.size(); ++j) {
      if (i != j && 
          (ip1[i] == ip1[j] || ip2[i] == ip2[j])) {
        bad_entry = true;
      }
    }
    if (!bad_entry) {
      new_ip1.push_back(ip1[i]);
      new_ip2.push_back(ip2[i]);
    }
  }
  
  ip1 = new_ip1;
  ip2 = new_ip2;
}

vw::math::Matrix<double> StereoSessionIsis::determine_image_alignment(std::string const& input_file1, std::string const& input_file2, float lo, float hi) {
  std::string left_align_image_file(input_file1), right_align_image_file(input_file2);

  // Load the two images
  DiskImageView<PixelGray<float> > left_disk_image(input_file1);
  DiskImageView<PixelGray<float> > right_disk_image(input_file2);

  // Image Alignment
  //
  // Images are aligned by computing interest points, matching
  // them using a standard 2-Norm nearest-neighor metric, and then
  // rejecting outliers by fitting a similarity between the
  // putative matches using RANSAC.  

  // Interest points are matched in image chunk of <= 2048x2048
  // pixels to conserve memory.
  vw_out(InfoMessage) << "\nInterest Point Detection:\n";
  InterestPointList ip1, ip2;
  
  // If the interest points are cached in a file, read that file in here.
  if ( boost::filesystem::exists( prefix_from_filename(input_file1) + ".vwip" ) && 
       boost::filesystem::exists( prefix_from_filename(input_file2) + ".vwip" )) {
    std::cout << "Found cached interest point files: \n";
    std::cout << "\t" << (prefix_from_filename(input_file1) + ".vwip") << "\n";
    std::cout << "\t" << (prefix_from_filename(input_file2) + ".vwip") << "\n";
    std::cout << "Skipping interest point detection step.\n";
    
  } else {

    ImageViewRef<PixelGray<float> > left_image = normalize(remove_isis_special_pixels(left_disk_image, lo), lo, hi, 0, 1.0);
    ImageViewRef<PixelGray<float> > right_image = normalize(remove_isis_special_pixels(right_disk_image, lo), lo, hi, 0, 1.0);
    
    // Interest Point module detector code.
    LogInterestOperator log_detector;
    ScaledInterestPointDetector<LogInterestOperator> detector(log_detector, 500);
    std::cout << "Processing " << input_file1 << "\n";
    ip1 = detect_interest_points(left_image, detector);
    std::cout << "\nProcessing " << input_file2 << "\n";
    ip2 = detect_interest_points(right_image, detector);
    vw_out(InfoMessage) << "\nLeft image: " << ip1.size() << " points.  Right image: " << ip2.size() << "\n"; 

    // Generate descriptors for interest points.
    // TODO: Switch to a more sophisticated descriptor
    vw_out(InfoMessage) << "\tGenerating descriptors... " << std::flush;
    PatchDescriptorGenerator descriptor;
    descriptor(left_image, ip1);
    descriptor(right_image, ip2);
    vw_out(InfoMessage) << "done.\n";
  
    // Write out the results
    vw_out(InfoMessage) << "\tWriting interest points to disk... " << std::flush;
    write_binary_ip_file(prefix_from_filename(input_file1)+".vwip", ip1);
    write_binary_ip_file(prefix_from_filename(input_file2)+".vwip", ip2);
    vw_out(InfoMessage) << "done.\n";
  }

  // The basic interest point matcher does not impose any
  // constraints on the matched interest points.
  vw_out(InfoMessage) << "\nInterest Point Matching:\n";
  double matcher_threshold = 0.8;

  // RANSAC needs the matches as a vector, and so does the matcher.
  // this is messy, but for now we simply make a copy.
  std::vector<InterestPoint> ip1_copy, ip2_copy;
  ip1_copy = read_binary_ip_file(prefix_from_filename(input_file1)+".vwip");
  ip2_copy = read_binary_ip_file(prefix_from_filename(input_file2)+".vwip");

  InterestPointMatcher<L2NormMetric,NullConstraint> matcher(matcher_threshold);
  std::vector<InterestPoint> matched_ip1, matched_ip2;
  matcher(ip1_copy, ip2_copy, matched_ip1, matched_ip2, false, TerminalProgressCallback());
  vw_out(InfoMessage) << "Found " << matched_ip1.size() << " putative matches.\n";
       
  std::vector<Vector3> ransac_ip1(matched_ip1.size());
  std::vector<Vector3> ransac_ip2(matched_ip2.size());
  for (unsigned i = 0; i < matched_ip1.size();++i ) {
    ransac_ip1[i] = Vector3(matched_ip1[i].x, matched_ip1[i].y,1);
    ransac_ip2[i] = Vector3(matched_ip2[i].x, matched_ip2[i].y,1);
  }
  
  remove_duplicates(ransac_ip1, ransac_ip2);

  // RANSAC is used to fit a similarity transform between the
  // matched sets of points  
  RandomSampleConsensus<math::HomographyFittingFunctor, InterestPointErrorMetric> ransac( vw::math::HomographyFittingFunctor(),
                                                                                          InterestPointErrorMetric(), 
                                                                                          10 ); // inlier_threshold

  std::vector<Vector3> result_ip1;
  std::vector<Vector3> result_ip2;
  vw_out(InfoMessage) << "\nRunning RANSAC:\n";
  Matrix<double> T;
  try {
    T = ransac(ransac_ip2,ransac_ip1);
  } catch (vw::ip::RANSACErr &e) {
    std::cout << "\n*************************************************************\n";
    std::cout << "WARNING: Automatic Alignment Failed!  Proceed with caution...\n";
    std::cout << "*************************************************************\n\n";
        T.set_size(3,3);
    T.set_identity();
  }
  return T;
}

// Pre-align the ISIS images.  If the ISIS images are map projected,
// we can perform pre-alignment by transforming them both into a
// common map projection.  Otherwise, we resort to feature-based image
// matching techniques to align the right image to the left image.
void StereoSessionIsis::pre_preprocessing_hook(std::string const& input_file1, std::string const& input_file2,
                                               std::string & output_file1, std::string & output_file2) {

  DiskImageView<PixelGray<float> > left_disk_image(input_file1);
  DiskImageView<PixelGray<float> > right_disk_image(input_file2);
  output_file1 = m_out_prefix + "-L.tif";
  output_file2 = m_out_prefix + "-R.tif";

  GeoReference input_georef1, input_georef2;
  try {
    // Read georeferencing information (if it exists...)
    DiskImageResourceGDAL file_resource1( input_file1 );
    DiskImageResourceGDAL file_resource2( input_file2 );
    read_georeference( input_georef1, file_resource1 );
    read_georeference( input_georef2, file_resource2 );
  } catch (ArgumentErr &e) {
    std::cout << "Warning: Couldn't read georeference data from input images using the GDAL driver.\n";
  }

  // Make sure the images are normalized
  vw_out(InfoMessage) << "Computing min/max values for normalization.  " << std::flush;
  float left_lo, left_hi, right_lo, right_hi;
  isis_min_max_channel_values(left_disk_image, left_lo, left_hi);
  vw_out(InfoMessage) << "Left: [" << left_lo << " " << left_hi << "]    " << std::flush;
  isis_min_max_channel_values(right_disk_image, right_lo, right_hi);
  vw_out(InfoMessage) << "Right: [" << right_lo << " " << right_hi << "]\n\n";
  float lo = std::min (left_lo, right_lo);
  float hi = std::min (left_hi, right_hi);

  // If this is a map projected cube, we skip the step of aligning the
  // images, because the map projected images are probable very nearly
  // aligned already.  For unprojected cubes, we align in the "usual"
  // way using interest points.
  if (input_georef1.transform() != math::identity_matrix<3>() && 
      input_georef2.transform() != math::identity_matrix<3>() ) {
    std::cout << "\nMap projected ISIS cubes detected.  Placing both images into the same map projection.\n\n";
    
    // If we are using map-projected cubes, we need to put them into a
    // common projection.  We adopt the projection of the first image.
    GeoReference common_georef = input_georef1;
    BBox2i common_bbox(0,0,left_disk_image.cols(), left_disk_image.rows());

    // Create the geotransform objects and determine the common size
    // of the output images and apply it to the right image.
    GeoTransform trans2(input_georef2, common_georef);
    ImageViewRef<PixelGray<float> > Limg = normalize(remove_isis_special_pixels(left_disk_image, lo), lo, hi, 0, 1.0);
    ImageViewRef<PixelGray<float> > Rimg = crop(transform(normalize(remove_isis_special_pixels(right_disk_image, lo),lo,hi,0.0,1.0),trans2),common_bbox);

    // Write the results to disk.
    write_image(output_file1, channel_cast_rescale<uint8>(Limg), TerminalProgressCallback());
    write_image(output_file2, channel_cast_rescale<uint8>(Rimg), TerminalProgressCallback()); 
  
  } else {
    // For unprojected ISIS images, we resort to the "old style" of
    // image alignment: determine the alignment matrix using keypoint
    // matching techniques.
    std::cout << "\nUnprojected ISIS cubes detected.  Aligning images using feature-based matching techniques.\n\n";

    Matrix<double> align_matrix(3,3);
    align_matrix.set_identity();
    if (stereo_settings().keypoint_alignment)
      align_matrix = determine_image_alignment(input_file1, input_file2, lo, hi);
    ::write_matrix(m_out_prefix + "-align.exr", align_matrix);

    // Apply the alignment transformation to the right image.
    ImageViewRef<PixelGray<float> > Limg = normalize(remove_isis_special_pixels(left_disk_image, lo),lo,hi,0.0,1.0);
    ImageViewRef<PixelGray<float> > Rimg = transform(normalize(remove_isis_special_pixels(right_disk_image,lo),lo,hi,0.0,1.0), 
                                                     HomographyTransform(align_matrix),
                                                     left_disk_image.cols(), left_disk_image.rows());

    // Write the results to disk.
    write_image(output_file1, channel_cast_rescale<uint8>(Limg), TerminalProgressCallback());
    write_image(output_file2, channel_cast_rescale<uint8>(Rimg), TerminalProgressCallback()); 
  }
}

// Stage 2: Correlation
//
// Pre file is a pair of grayscale images.  ( ImageView<PixelGray<float> > )
// Post file is a disparity map.            ( ImageView<PixelDisparity<float> > )
void StereoSessionIsis::pre_filtering_hook(std::string const& input_file, std::string & output_file) {
  std::cout << "******** Call to post ISIS3 correlation hook: " << input_file << "! *************\n";
  output_file = m_out_prefix + "-D-masked.exr";
  
  ImageView<uint8> Lmask, Rmask;
  DiskImageView<PixelGray<float> > left_disk_image(m_left_image_file);
  DiskImageView<PixelGray<float> > right_disk_image(m_right_image_file);
  
  read_image(Lmask,m_out_prefix + "-lMask.tif");
  read_image(Rmask,m_out_prefix + "-rMask.tif");
  disparity::mask_black_pixels(clamp(left_disk_image,0,1e6), Lmask);
  disparity::mask_black_pixels(clamp(right_disk_image,0,1e6), Rmask);
  write_image(m_out_prefix + "-lMaskDebug.tif", Lmask);
  write_image(m_out_prefix + "-rMaskDebug.tif", Rmask);
  
  DiskImageView<PixelDisparity<float> > disparity_disk_image(input_file);
  ImageViewRef<PixelDisparity<float> > disparity_map = disparity::mask(disparity_disk_image, Lmask, Rmask);

  DiskImageResourceOpenEXR disparity_map_rsrc(output_file, disparity_map.format() );
  disparity_map_rsrc.set_tiled_write(std::min(2048,disparity_map.cols()),std::min(2048, disparity_map.rows()));
  block_write_image( disparity_map_rsrc, disparity_map, TerminalProgressCallback() );
}

// Reverse any pre-alignment that was done to the images.
void StereoSessionIsis::pre_pointcloud_hook(std::string const& input_file, std::string & output_file) {

  DiskImageView<PixelDisparity<float> > disparity_map(input_file);
  output_file = m_out_prefix + "-F-corrected.exr";<commit_msg>Removed (by default, at least) the code that detects dark portions of the ISIS camera frames, as this causes trouble when processing other isis images.  It's a power user feature for now.<commit_after>// __BEGIN_LICENSE__
// 
// Copyright (C) 2008 United States Government as represented by the
// Administrator of the National Aeronautics and Space Administration
// (NASA).  All Rights Reserved.
// 
// Copyright 2008 <NAME>. All rights reserved.
// 
// This software is distributed under the NASA Open Source Agreement
// (NOSA), version 1.3.  The NOSA has been approved by the Open Source
// Initiative.  See the file COPYING at the top of the distribution
// directory tree for the complete NOSA document.
// 
// THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY
// KIND, EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT
// LIMITED TO, ANY WARRANTY THAT THE SUBJECT SOFTWARE WILL CONFORM TO
// SPECIFICATIONS, ANY IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR
// A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY WARRANTY THAT
// THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT
// DOCUMENTATION, IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE.
//
// __END_LICENSE__

/// \file StereoSessionIsis.cc
///

#include <boost/shared_ptr.hpp>

#include <vw/Image.h>
#include <vw/Math.h>
#include <vw/FileIO.h>
#include <vw/InterestPoint.h>
#include <vw/Stereo.h>
#include <vw/Cartography.h>
#include "Isis/StereoSessionIsis.h"
#include "Isis/IsisCameraModel.h"
#include "StereoSettings.h"

// Boost
#include "boost/filesystem.hpp"   
using namespace boost::filesystem; 

using namespace vw;
using namespace vw::camera;
using namespace vw::cartography;
using namespace vw::stereo;
using namespace vw::ip;


// Allows FileIO to correctly read/write these pixel types
namespace vw {
  template<> struct PixelFormatID<Vector3>   { static const PixelFormatEnum value = VW_PIXEL_GENERIC_3_CHANNEL; };
  template<> struct PixelFormatID<PixelDisparity<float> >   { static const PixelFormatEnum value = VW_PIXEL_GENERIC_3_CHANNEL; };
}


static std::string prefix_from_filename(std::string const& filename) {
  std::string result = filename;
  int index = result.rfind(".");
  if (index != -1) 
    result.erase(index, result.size());
  return result;
}

// Duplicate matches for any given interest point probably indicate a
// poor match, so we cull those out here.
static void remove_duplicates(std::vector<Vector3> &ip1, std::vector<Vector3> &ip2) {
  std::vector<Vector3> new_ip1, new_ip2;
  
  for (unsigned i = 0; i < ip1.size(); ++i) {
    bool bad_entry = false;
    for (unsigned j = 0; j < ip1.size(); ++j) {
      if (i != j && 
          (ip1[i] == ip1[j] || ip2[i] == ip2[j])) {
        bad_entry = true;
      }
    }
    if (!bad_entry) {
      new_ip1.push_back(ip1[i]);
      new_ip2.push_back(ip2[i]);
    }
  }
  
  ip1 = new_ip1;
  ip2 = new_ip2;
}

vw::math::Matrix<double> StereoSessionIsis::determine_image_alignment(std::string const& input_file1, std::string const& input_file2, float lo, float hi) {
  std::string left_align_image_file(input_file1), right_align_image_file(input_file2);

  // Load the two images
  DiskImageView<PixelGray<float> > left_disk_image(input_file1);
  DiskImageView<PixelGray<float> > right_disk_image(input_file2);

  // Image Alignment
  //
  // Images are aligned by computing interest points, matching
  // them using a standard 2-Norm nearest-neighor metric, and then
  // rejecting outliers by fitting a similarity between the
  // putative matches using RANSAC.  

  // Interest points are matched in image chunk of <= 2048x2048
  // pixels to conserve memory.
  vw_out(InfoMessage) << "\nInterest Point Detection:\n";
  InterestPointList ip1, ip2;
  
  // If the interest points are cached in a file, read that file in here.
  if ( boost::filesystem::exists( prefix_from_filename(input_file1) + ".vwip" ) && 
       boost::filesystem::exists( prefix_from_filename(input_file2) + ".vwip" )) {
    std::cout << "Found cached interest point files: \n";
    std::cout << "\t" << (prefix_from_filename(input_file1) + ".vwip") << "\n";
    std::cout << "\t" << (prefix_from_filename(input_file2) + ".vwip") << "\n";
    std::cout << "Skipping interest point detection step.\n";
    
  } else {

    ImageViewRef<PixelGray<float> > left_image = normalize(remove_isis_special_pixels(left_disk_image, lo), lo, hi, 0, 1.0);
    ImageViewRef<PixelGray<float> > right_image = normalize(remove_isis_special_pixels(right_disk_image, lo), lo, hi, 0, 1.0);
    
    // Interest Point module detector code.
    LogInterestOperator log_detector;
    ScaledInterestPointDetector<LogInterestOperator> detector(log_detector, 500);
    std::cout << "Processing " << input_file1 << "\n";
    ip1 = detect_interest_points(left_image, detector);
    std::cout << "\nProcessing " << input_file2 << "\n";
    ip2 = detect_interest_points(right_image, detector);
    vw_out(InfoMessage) << "\nLeft image: " << ip1.size() << " points.  Right image: " << ip2.size() << "\n"; 

    // Generate descriptors for interest points.
    // TODO: Switch to a more sophisticated descriptor
    vw_out(InfoMessage) << "\tGenerating descriptors... " << std::flush;
    PatchDescriptorGenerator descriptor;
    descriptor(left_image, ip1);
    descriptor(right_image, ip2);
    vw_out(InfoMessage) << "done.\n";
  
    // Write out the results
    vw_out(InfoMessage) << "\tWriting interest points to disk... " << std::flush;
    write_binary_ip_file(prefix_from_filename(input_file1)+".vwip", ip1);
    write_binary_ip_file(prefix_from_filename(input_file2)+".vwip", ip2);
    vw_out(InfoMessage) << "done.\n";
  }

  // The basic interest point matcher does not impose any
  // constraints on the matched interest points.
  vw_out(InfoMessage) << "\nInterest Point Matching:\n";
  double matcher_threshold = 0.8;

  // RANSAC needs the matches as a vector, and so does the matcher.
  // this is messy, but for now we simply make a copy.
  std::vector<InterestPoint> ip1_copy, ip2_copy;
  ip1_copy = read_binary_ip_file(prefix_from_filename(input_file1)+".vwip");
  ip2_copy = read_binary_ip_file(prefix_from_filename(input_file2)+".vwip");

  InterestPointMatcher<L2NormMetric,NullConstraint> matcher(matcher_threshold);
  std::vector<InterestPoint> matched_ip1, matched_ip2;
  matcher(ip1_copy, ip2_copy, matched_ip1, matched_ip2, false, TerminalProgressCallback());
  vw_out(InfoMessage) << "Found " << matched_ip1.size() << " putative matches.\n";
       
  std::vector<Vector3> ransac_ip1(matched_ip1.size());
  std::vector<Vector3> ransac_ip2(matched_ip2.size());
  for (unsigned i = 0; i < matched_ip1.size();++i ) {
    ransac_ip1[i] = Vector3(matched_ip1[i].x, matched_ip1[i].y,1);
    ransac_ip2[i] = Vector3(matched_ip2[i].x, matched_ip2[i].y,1);
  }
  
  remove_duplicates(ransac_ip1, ransac_ip2);

  // RANSAC is used to fit a similarity transform between the
  // matched sets of points  
  RandomSampleConsensus<math::HomographyFittingFunctor, InterestPointErrorMetric> ransac( vw::math::HomographyFittingFunctor(),
                                                                                          InterestPointErrorMetric(), 
                                                                                          10 ); // inlier_threshold

  std::vector<Vector3> result_ip1;
  std::vector<Vector3> result_ip2;
  vw_out(InfoMessage) << "\nRunning RANSAC:\n";
  Matrix<double> T;
  try {
    T = ransac(ransac_ip2,ransac_ip1);
  } catch (vw::ip::RANSACErr &e) {
    std::cout << "\n*************************************************************\n";
    std::cout << "WARNING: Automatic Alignment Failed!  Proceed with caution...\n";
    std::cout << "*************************************************************\n\n";
        T.set_size(3,3);
    T.set_identity();
  }
  return T;
}

// Pre-align the ISIS images.  If the ISIS images are map projected,
// we can perform pre-alignment by transforming them both into a
// common map projection.  Otherwise, we resort to feature-based image
// matching techniques to align the right image to the left image.
void StereoSessionIsis::pre_preprocessing_hook(std::string const& input_file1, std::string const& input_file2,
                                               std::string & output_file1, std::string & output_file2) {

  DiskImageView<PixelGray<float> > left_disk_image(input_file1);
  DiskImageView<PixelGray<float> > right_disk_image(input_file2);
  output_file1 = m_out_prefix + "-L.tif";
  output_file2 = m_out_prefix + "-R.tif";

  GeoReference input_georef1, input_georef2;
  try {
    // Read georeferencing information (if it exists...)
    DiskImageResourceGDAL file_resource1( input_file1 );
    DiskImageResourceGDAL file_resource2( input_file2 );
    read_georeference( input_georef1, file_resource1 );
    read_georeference( input_georef2, file_resource2 );
  } catch (ArgumentErr &e) {
    std::cout << "Warning: Couldn't read georeference data from input images using the GDAL driver.\n";
  }

  // Make sure the images are normalized
  vw_out(InfoMessage) << "Computing min/max values for normalization.  " << std::flush;
  float left_lo, left_hi, right_lo, right_hi;
  isis_min_max_channel_values(left_disk_image, left_lo, left_hi);
  vw_out(InfoMessage) << "Left: [" << left_lo << " " << left_hi << "]    " << std::flush;
  isis_min_max_channel_values(right_disk_image, right_lo, right_hi);
  vw_out(InfoMessage) << "Right: [" << right_lo << " " << right_hi << "]\n\n";
  float lo = std::min (left_lo, right_lo);
  float hi = std::min (left_hi, right_hi);

  // If this is a map projected cube, we skip the step of aligning the
  // images, because the map projected images are probable very nearly
  // aligned already.  For unprojected cubes, we align in the "usual"
  // way using interest points.
  if (input_georef1.transform() != math::identity_matrix<3>() && 
      input_georef2.transform() != math::identity_matrix<3>() ) {
    std::cout << "\nMap projected ISIS cubes detected.  Placing both images into the same map projection.\n\n";
    
    // If we are using map-projected cubes, we need to put them into a
    // common projection.  We adopt the projection of the first image.
    GeoReference common_georef = input_georef1;
    BBox2i common_bbox(0,0,left_disk_image.cols(), left_disk_image.rows());

    // Create the geotransform objects and determine the common size
    // of the output images and apply it to the right image.
    GeoTransform trans2(input_georef2, common_georef);
    ImageViewRef<PixelGray<float> > Limg = normalize(remove_isis_special_pixels(left_disk_image, lo), lo, hi, 0, 1.0);
    ImageViewRef<PixelGray<float> > Rimg = crop(transform(normalize(remove_isis_special_pixels(right_disk_image, lo),lo,hi,0.0,1.0),trans2),common_bbox);

    // Write the results to disk.
    write_image(output_file1, channel_cast_rescale<uint8>(Limg), TerminalProgressCallback());
    write_image(output_file2, channel_cast_rescale<uint8>(Rimg), TerminalProgressCallback()); 
  
  } else {
    // For unprojected ISIS images, we resort to the "old style" of
    // image alignment: determine the alignment matrix using keypoint
    // matching techniques.
    std::cout << "\nUnprojected ISIS cubes detected.  Aligning images using feature-based matching techniques.\n\n";

    Matrix<double> align_matrix(3,3);
    align_matrix.set_identity();
    if (stereo_settings().keypoint_alignment)
      align_matrix = determine_image_alignment(input_file1, input_file2, lo, hi);
    ::write_matrix(m_out_prefix + "-align.exr", align_matrix);

    // Apply the alignment transformation to the right image.
    ImageViewRef<PixelGray<float> > Limg = normalize(remove_isis_special_pixels(left_disk_image, lo),lo,hi,0.0,1.0);
    ImageViewRef<PixelGray<float> > Rimg = transform(normalize(remove_isis_special_pixels(right_disk_image,lo),lo,hi,0.0,1.0), 
                                                     HomographyTransform(align_matrix),
                                                     left_disk_image.cols(), left_disk_image.rows());

    // Write the results to disk.
    write_image(output_file1, channel_cast_rescale<uint8>(Limg), TerminalProgressCallback());
    write_image(output_file2, channel_cast_rescale<uint8>(Rimg), TerminalProgressCallback()); 
  }
}

// Stage 2: Correlation
//
// Pre file is a pair of grayscale images.  ( ImageView<PixelGray<float> > )
// Post file is a disparity map.            ( ImageView<PixelDisparity<float> > )
void StereoSessionIsis::pre_filtering_hook(std::string const& input_file, std::string & output_file) {
  output_file = input_file;
  // ****************************************************
  // The following code is for Apollo Metric Camera ONLY!
  // (uncomment at your own risk)

//   output_file = m_out_prefix + "-D-masked.exr";
  
//   ImageView<uint8> Lmask, Rmask;
//   DiskImageView<PixelGray<float> > left_disk_image(m_left_image_file);
//   DiskImageView<PixelGray<float> > right_disk_image(m_right_image_file);
  
//   read_image(Lmask,m_out_prefix + "-lMask.tif");
//   read_image(Rmask,m_out_prefix + "-rMask.tif");
//   disparity::mask_black_pixels(clamp(left_disk_image,0,1e6), Lmask);
//   disparity::mask_black_pixels(clamp(right_disk_image,0,1e6), Rmask);
//   write_image(m_out_prefix + "-lMaskDebug.tif", Lmask);
//   write_image(m_out_prefix + "-rMaskDebug.tif", Rmask);
  
//   DiskImageView<PixelDisparity<float> > disparity_disk_image(input_file);
//   ImageViewRef<PixelDisparity<float> > disparity_map = disparity::mask(disparity_disk_image, Lmask, Rmask);

//   DiskImageResourceOpenEXR disparity_map_rsrc(output_file, disparity_map.format() );
//   disparity_map_rsrc.set_tiled_write(std::min(2048,disparity_map.cols()),std::min(2048, disparity_map.rows()));
//   block_write_image( disparity_map_rsrc, disparity_map, TerminalProgressCallback() );
}

// Reverse any pre-alignment that was done to the images.
void StereoSessionIsis::pre_pointcloud_hook(std::string const& input_file, std::string & output_file) {

  DiskImageView<PixelDisparity<float> > disparity_map(input_file);
  output_file = m_out_prefix + "-F-corrected.exr";
{'fraction_non_alphanumeric': 0.1307210477614065, 'fraction_numerical': 0.014057940066908676, 'mean_word_length': 4.033858831959871, 'pattern_counts': {'":': 0, '<': 355, '<?xml version=': 0, '>': 141, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '98221', 'n_tokens_mistral': 9380, 'n_tokens_neox': 8927, 'n_words': 2197}