<|start_filename|>sys.mk<|end_filename|>
HOST_SYS =
HOST_CPU =

ifdef COMSPEC
  HOST_SYS = windows
  HOST_CPU = $(shell uname -m)
else
  UNAME = $(shell uname -s)
ifeq ($(UNAME), Linux)
  HOST_SYS = linux
  HOST_CPU = $(shell uname -m)
else
ifeq ($(UNAME), Darwin)
  HOST_SYS = osx
  HOST_CPU = universal
endif
endif
endif

ifndef HOST_SYS
  HOST_SYS = generic
  $(warning WARNING: Could not determine host system, assuming $(HOST_SYS))
endif

ifndef HOST_CPU
  HOST_CPU = generic
  $(warning WARNING: Could not determine host CPU, assuming $(HOST_CPU))
endif

HOST = $(HOST_SYS)-$(HOST_CPU)

ifndef TARGET
  TARGET = $(HOST)
endif
TARGET_SYS = $(shell echo $(TARGET) | sed "s/-.*$$//")
TARGET_CPU = $(shell echo $(TARGET) | sed "s/^$(TARGET_SYS)-//")

CC =
CXX =
CMAKE_TOOLCHAIN =
LIPO =
SDKFLAGS =
STRIP = strip
INSTFLAGS =

ifeq ($(HOST), $(TARGET))
  CC = gcc
  CXX = g++
ifeq ($(HOST_SYS), osx)
  LIPO = lipo
endif
else				# Cross compilation
ifeq ($(HOST_SYS), linux)
ifeq ($(TARGET), windows-x86_64)
  CC = x86_64-w64-mingw32-gcc
  CXX = x86_64-w64-mingw32-g++
  CMAKE_TOOLCHAIN = x86_64-w64-mingw32.cmake
else
ifeq ($(TARGET), windows-i686)
  CC = i686-w64-mingw32-gcc
  CXX = i686-w64-mingw32-g++
  CMAKE_TOOLCHAIN = i686-w64-mingw32.cmake
else
ifeq ($(TARGET), osx-x86_64)
  CC = x86_64-apple-darwin11-gcc
  CXX = x86_64-apple-darwin11-g++
  CMAKE_TOOLCHAIN = x86_64-apple-darwin11.cmake
  SDKFLAGS = -mmacosx-version-min=10.7 -arch x86_64
  STRIP = x86_64-apple-darwin11-strip
else
ifeq ($(TARGET), osx-i386)
  CC = i386-apple-darwin11-gcc
  CXX = i386-apple-darwin11-g++
  CMAKE_TOOLCHAIN = i386-apple-darwin11.cmake
  SDKFLAGS = -mmacosx-version-min=10.7 -arch i386
  STRIP = i386-apple-darwin11-strip
else
ifeq ($(TARGET), osx-universal)
  LIPO = x86_64-apple-darwin11-lipo
  STRIP = x86_64-apple-darwin11-strip
endif
endif
endif
endif
endif
else
ifeq ($(HOST_SYS), osx)
ifeq ($(TARGET_SYS), osx)
  CC = gcc
  CXX = g++
  SDKFLAGS = -isysroot /Developer/SDKs/MacOSX10.7.sdk -mmacosx-version-min=10.7 -arch $(TARGET_CPU)
endif
endif
endif
endif

ifneq ($(STRIP), strip)
  INSTFLAGS = --strip-program=$(STRIP)
endif

ifeq ($(TARGET), osx-universal)

ifndef LIPO
  $(error Building universal binaries for target $(TARGET) on host $(HOST) is unsupported)
endif

else

ifndef CC
  $(error C compilation for target $(TARGET) on host $(HOST) is unsupported)
endif

ifndef CXX
  $(error C++ compilation for target $(TARGET) on host $(HOST) is unsupported)
endif

endif


<|start_filename|>src/x3f_denoise_aniso.cpp<|end_filename|>
/* X3F_DENOISE_ANISO.CPP
 *
 * Library for anisotropic denoising of X3F image data.
 *
 * Copyright 2015 - <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#include <iostream>
#include "x3f_denoise_aniso.h"
#include "x3f_printf.h"
#include <algorithm> //for std::sort for the median filter

static inline float determine_pixel_difference(const float& v1_1, const float& v1_2, const float& v1_3,
                                               const float& v2_1, const float& v2_2, const float& v2_3)
{ //assumes 3 channels for now, may be a bad assumption
  //also doing the exp version rather than the other technique
  //since the diff is symmetric, can precalc these and cut down on processing time by a factor of 2
  const float c1 = (v1_1 - v2_1);
  const float c2 = (v1_2 - v2_2);
  const float c3 = (v1_3 - v2_3);
  const float K = 0.0000050f;
  const float l2norm = (c1 * c1 + c2 * c2 + c3 * c3); //should be sqrt'd to be mathy, but it's immediately squared
  return -1.0f * exp(-1.0f * l2norm/K); //not the lack of squaring of l2
  //return -1.0f * exp(-1.0f * l2norm*l2norm/K); //expensive but better
  //return -1.0f/(1.0f - l2norm*l2norm/K);
}

//need a median filter to remove the pepper noise
//there are faster ways to do this, but this variant is immediately implementable in gl
//probably a variant in opencv for this
void median_filter(x3f_area16_t *image)
{
  uint32_t x, y, i, j;
  int dx, dy; //needs to be signed, derp
  const uint32_t rows = image->rows;
  const uint32_t cols = image->columns;
  const uint32_t edgeless_cols = cols - 1;
  const uint32_t edgeless_rows = rows - 1;
  const uint32_t jump = image->channels;
  const uint32_t stride = image->row_stride;

  uint16_t* current_set = new uint16_t[9];
  memset(current_set, 0, sizeof(uint16_t)*9);
  uint16_t* current_ptr;
  float median, current_val; //so we can have negatives
  uint16_t thresh = 1000;
  uint16_t background_thresh = 100;  //entirely arbitrarily chosen
  //gonna do median filtering by channel for now, to start with.
  uint16_t* scratch = new uint16_t[rows * cols * jump];
  float avg;

  memcpy(scratch, image->data, rows*cols*jump*sizeof(uint16_t));

  for (y = 1; y < edgeless_rows; y++){
    for (x = 1; x < edgeless_cols; x++) {
      for (i = 0; i < jump; i++){
        current_ptr = current_set;
        current_val = image->data[y*stride + x*jump + i];
        for (dy = -1; dy <= 1; dy++){
          for (dx = -1; dx <= 1; dx++, current_ptr++){
            *current_ptr = image->data[(y + dy)*stride + (x + dx)*jump + i];
          }
        }
	std::sort(current_set, current_set + 9);
        //if the current pixel is very very different from the median
        //then replace.  'very very different' means 'more than 50% of the dynamic
        //range of the image away' or 'zero'
        median = current_set[4];
        if (fabs(current_val - median) > thresh || current_val < background_thresh)
	  {
	    //replace with the local average in this channel
	    //scratch[y*stride + x*jump + i] = median;
	    for (j = 0; j < jump; j++){
	      avg = 0;
	      for (dy = -1; dy <= 1; dy++){
		for (dx = -1; dx <= 1; dx++, current_ptr++){
		  if (dx == 0 && dy == 0) continue;
		  avg += image->data[(y + dy)*stride + (x + dx)*jump + j];
		}
	      }
	      scratch[y*stride + x*jump + j] = avg/8.0f; //should be int cast when put back in
	    }
	  }
      }
    }
  }
  //copy from scratch back into the image
  memcpy(image->data, scratch, rows*cols*jump*sizeof(uint16_t));


  delete [] scratch;
  delete [] current_set;
}

//crude denoising.  Ignores boundary conditions by being super lazy.
void denoise_aniso_float(const uint32_t& in_rows, const uint32_t& in_columns, float* image)
{
  uint32_t x, y, i;
  const uint32_t jump = 3;
  const uint32_t num_rows = in_rows;
  const uint32_t row_stride = in_columns * jump;
  const uint32_t edgeless_rows = num_rows - 1;
  const uint32_t edgeless_cols = in_columns - 1;  //used for boundaries
  float* in_ptr;
  float* out_ptr;
  float* out_image = new float[in_rows * in_columns * jump];
  float coeff_fwd, coeff_bkwd, coeff_up, coeff_down, lambda1, lambda2;
  //do the memcpy to keep the boundaries as they were before processing
  //not ideal, because there will be noise there, but it should avoid the green boundary problem.
  //there are a number of ways to solve boundary problems, but none of them are efficient.
  memcpy(out_image, image, in_rows * in_columns * jump * sizeof(float));

  //diffusion for the first channel is different, to prevent color bleed
  float* h = new float[jump];
  h[0] = 5.0f;
  h[1] = 1.0f;
  h[2] = 1.0f;

  for (y = 1; y < edgeless_rows; y++){
    for (x = 1,
	   out_ptr = (out_image + y * row_stride + jump),
	   in_ptr = (image + y * row_stride + jump);
         x < edgeless_cols;
         x++, out_ptr += jump, in_ptr += jump){

      coeff_fwd = determine_pixel_difference(*in_ptr, *(in_ptr+1), *(in_ptr+2),
                                             *(in_ptr + jump), *(in_ptr + jump + 1), *(in_ptr + jump + 2));
      coeff_bkwd = determine_pixel_difference(*in_ptr, *(in_ptr+1), *(in_ptr+2),
					      *(in_ptr - jump), *(in_ptr - jump + 1), *(in_ptr - jump + 2));
      coeff_down = determine_pixel_difference(*in_ptr, *(in_ptr+1), *(in_ptr+2),
					      *(in_ptr + row_stride), *(in_ptr + row_stride + 1), *(in_ptr + row_stride + 2));
      coeff_up = determine_pixel_difference(*in_ptr, *(in_ptr+1), *(in_ptr+2),
                                            *(in_ptr - row_stride), *(in_ptr - row_stride + 1), *(in_ptr - row_stride + 2));

      //uncomment out this line to do isotropic noise reduction
      //coeff_fwd = coeff_up = coeff_down = coeff_bkwd = -1.0f;
      lambda1 = coeff_fwd + coeff_bkwd + coeff_down + coeff_up;
      lambda2 = fabs(lambda1);

      for (i = 0; i < jump; i++)
	{
	  *(out_ptr + i) = *(in_ptr + i) - (lambda2 * *(in_ptr + i) + coeff_fwd * *(in_ptr + jump + i)
					    + coeff_bkwd * *(in_ptr - jump + i)
					    + coeff_down * *(in_ptr + row_stride + i)
					    + coeff_up * *(in_ptr - row_stride + i))/(h[i] *lambda2);

	}
    }
  }

  memcpy(image, out_image, in_rows * in_columns * jump * sizeof(float));
  delete [] out_image;
  delete [] h;
}

//puts the aniso denoising into float format, does the denoising, then
//converts back to shorts.  Along the way, ditches the luminosity channel in the back conversion.
void denoise_aniso(x3f_area16_t *image, const int& in_iterations)
{
  median_filter(image);
  float* float_image = convert_to_float_image(image);
  for (int i = 0; i < in_iterations; i++){
    denoise_aniso_float(image->rows, image->columns, float_image);
    x3f_printf(DEBUG, "iteration: %d\n", i);
  }
  convert_from_float_image(image, float_image);
  delete [] float_image;
}


//even cruder denoising.  Ignores boundary conditions by being super lazy.
//also, is basically a gaussian blur.
void denoise_iso_float(const uint32_t& in_rows, const uint32_t& in_columns, float* image)
{
  uint32_t x, y, i;
  const uint32_t jump = 3;
  const uint32_t num_rows = in_rows;
  const uint32_t row_stride = in_columns * jump;
  const uint32_t edgeless_rows = num_rows - 1;
  const uint32_t edgeless_cols = in_columns - 1;  //used for boundaries
  float* in_ptr;
  float* out_ptr;
  float* out_image = new float[in_rows * in_columns * jump];
  float coeff_fwd, coeff_bkwd, coeff_up, coeff_down, lambda1, lambda2;
  //do the memcpy to keep the boundaries as they were before processing
  //not ideal, because there will be noise there, but it should avoid the green boundary problem.
  //there are a number of ways to solve boundary problems, but none of them are efficient.
  memcpy(out_image, image, in_rows * in_columns * jump * sizeof(float));

  //diffusion for the first channel is different, to prevent color bleed
  float* h = new float[jump];
  h[0] = 5.0f;
  h[1] = 1.0f;
  h[2] = 1.0f;

  coeff_fwd = coeff_up = coeff_down = coeff_bkwd = -1.0f;
  lambda1 = coeff_fwd + coeff_bkwd + coeff_down + coeff_up;
  lambda2 = fabs(lambda1);

  for (y = 1; y < edgeless_rows; y++){
    for (x = 1,
         out_ptr = (out_image + y * row_stride + jump),
         in_ptr = (image + y * row_stride + jump);
         x < edgeless_cols;
         x++, out_ptr += jump, in_ptr += jump){


      for (i = 0; i < jump; i++)
      {
        *(out_ptr + i) = *(in_ptr + i) - (lambda2 * *(in_ptr + i) + coeff_fwd * *(in_ptr + jump + i)
                                          + coeff_bkwd * *(in_ptr - jump + i)
                                          + coeff_down * *(in_ptr + row_stride + i)
                                          + coeff_up * *(in_ptr - row_stride + i))/(h[i] *lambda2);

      }
    }
  }

  memcpy(image, out_image, in_rows * in_columns * jump * sizeof(float));
  delete [] out_image;
  delete [] h;
}

//puts the iso denoising into float format, does the denoising, then
//converts back to shorts.  Along the way, ditches the luminosity channel in the back conversion.
void denoise_iso(x3f_area16_t *image, const int& in_iterations)
{
  median_filter(image);
  float* float_image = convert_to_float_image(image);
  for (int i = 0; i < in_iterations; i++){
    denoise_iso_float(image->rows, image->columns, float_image);
    x3f_printf(DEBUG, "iteration: %d\n", i);
  }
  convert_from_float_image(image, float_image);
  delete [] float_image;
}

//use these two functions as function pointers for the dilate/erosion operations
inline uint16_t give_max(const uint16_t& left, const uint16_t& right)
{
  return left > right ? left : right;
}

inline uint16_t give_min(const uint16_t& left, const uint16_t& right)
{
  return left < right ? left : right;
}

//denoising using the morphological operations that compose the 'splotchify' algorithm
//ideally directly implementable in opencv, come to think.
void morphological_op(x3f_area16_t *image, const int& in_radius, const bool& erode)
{
  const int ysize = image->rows;  //forcing a cast to signed for boundary conditions
  const int xsize = image->columns;

  uint16_t (*determinant)(const uint16_t&, const uint16_t&);
  determinant = &give_max;
  uint16_t seed = 0;
  if (erode){
    determinant = &give_min;
    seed = 65535;
  }

  const int the_edge = in_radius*2 + 1;
  bool** theMask = new bool*[the_edge];
  const int sqrRad = in_radius*in_radius;

  int x, y;
  for (y = -in_radius; y <= in_radius; y++){
    theMask[y + in_radius] = new bool[the_edge];
    for (x = -in_radius; x <= in_radius; x++){
      if (x*x + y*y <= sqrRad){
        theMask[y+in_radius][x+in_radius] = true;
      } else {
        theMask[y+in_radius][x+in_radius] = false;
      }
    }
  }

  int dx, dy, i;
  const int jump = image->channels;
  uint16_t* outPixelData = new uint16_t[ysize*xsize*jump];
  memcpy(outPixelData, image->data, ysize*xsize*jump*sizeof(uint16_t));
  uint16_t* in_pxl_ptr = image->data;
  uint16_t* out_pxl_ptr = outPixelData;
  uint16_t* maxes = new uint16_t[jump];
  for (y = 0; y < ysize; y++){
    for (x = 0; x < xsize; x++, in_pxl_ptr+=jump, out_pxl_ptr+=jump){
      for (i = 0; i < jump; i++){
        *(maxes+i) = seed;
      }
      for (dy = -in_radius; dy <= in_radius; dy++){
        if (y + dy < 0 || y + dy >= ysize) continue;
        for (dx = -in_radius; dx <= in_radius; dx++){
          if (x + dx < 0 || x + dx >= xsize) continue;
          if (!theMask[dy+in_radius][dx+in_radius]) continue;
          for (i = 0; i < jump; i++){
            *(maxes + i) = determinant(*(outPixelData + i), *(maxes + i));
          }
        }
      }
      for (i = 0; i < jump; i++){
        *(in_pxl_ptr + i) = *(maxes + i);
      }
    }
  }

  //memcpy(image->data, outPixelData, ysize*xsize*jump*sizeof(uint16_t));
  delete [] outPixelData;
  delete [] maxes;

  for (y = -in_radius; y <= in_radius; y++){
    delete [] theMask[y + in_radius];
  }
  delete [] theMask;
}

//uses morphological operations for every aggressive noise reduction for high ISO images.
//the effect is somewhat artistic (and probably pretty useless) at low ISOs.
void denoise_splotchify(x3f_area16_t *image, const int& in_radius)
{
  size_t i;
  size_t x, y;
  const size_t xsize = image->columns;
  const size_t ysize = image->rows;
  const size_t channel_size = xsize*ysize;
  uint16_t* channel = new uint16_t[channel_size];
  uint16_t* chan_ptr = channel;
  uint16_t* image_ptr;

  int morph_elem = 2; //ellipse
  int morph_size = in_radius;

  cv::Mat element = cv::getStructuringElement( morph_elem, cv::Size( 2*morph_size + 1, 2*morph_size+1 ), cv::Point( morph_size, morph_size ) );
  for (i = 0; i < image->channels; i++){
    //first, get the data, one channel at a time.
    for (y = 0; y < ysize; y++){
      for (x = 0, chan_ptr = &(channel[y*xsize]), image_ptr = &(image->data[y*image->row_stride + i]);
           x < xsize; x++, chan_ptr++, image_ptr+=3){
        *chan_ptr = *image_ptr;
      }
    }

    if (i>0){//skipping the first channel
      cv::Mat img(image->rows, image->columns, CV_16U, //apparently doesn't work on color images
              channel, sizeof(uint16_t)*xsize);

      cv::UMat out;
      cv::morphologyEx( img, out, 3, element );
      cv::morphologyEx( out, img, 2, element );


      for (y = 0; y < ysize; y++){
        for (x = 0, chan_ptr = &(channel[y*xsize]), image_ptr = &(image->data[y*image->row_stride + i]);
             x < xsize; x++, chan_ptr++, image_ptr+=3){
          *image_ptr = *chan_ptr;
        }
      }
    }
  }
  delete [] channel;

}


<|start_filename|>src/x3f_denoise_utils.h<|end_filename|>
/* X3F_DENOISE_UTILS.H
 *
 * Library for support functions for noise reduction algorithms,
 * including color conversions and conversions to/from float type images
 *
 * Copyright 2015 - <NAME> and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#ifndef DENOISE_UTILS
#define DENOISE_UTILS

#ifndef __cplusplus
#error This file can only be included from C++
#endif

#include "x3f_io.h"

#include <opencv2/core.hpp>
#include <opencv2/core/ocl.hpp>
#include <opencv2/photo.hpp>
#include <opencv2/imgproc.hpp>

typedef void (*conv_t)(x3f_area16_t *image);


typedef struct {
  float h;
  conv_t BMT_to_YUV;
  conv_t YUV_to_BMT;
} denoise_desc_t;



void BMT_to_YUV_YisT(x3f_area16_t *image);
void YUV_to_BMT_YisT(x3f_area16_t *image);
void BMT_to_YUV_STD(x3f_area16_t *image);
void YUV_to_BMT_STD(x3f_area16_t *image);
void BMT_to_YUV_Yis4T(x3f_area16_t *image);
void YUV_to_BMT_Yis4T(x3f_area16_t *image);


float* convert_to_float_image(x3f_area16_t *image);
void convert_from_float_image(x3f_area16_t *image, float* in_float_image);


const denoise_desc_t denoise_types[] = {
  {100.0, BMT_to_YUV_STD, YUV_to_BMT_STD},
  {70.0, BMT_to_YUV_YisT, YUV_to_BMT_YisT},
  {300.0, BMT_to_YUV_Yis4T, YUV_to_BMT_Yis4T},
};

#endif //DENOISE_UTILS


<|start_filename|>src/mingw_dowildcard.c<|end_filename|>
int _dowildcard = 1;


<|start_filename|>src/x3f_denoise.cpp<|end_filename|>
/* X3F_DENOISE.CPP
 *
 * Library for denoising of X3F image data.
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#include <iostream>
#include <inttypes.h>

#include <opencv2/core.hpp>
#include <opencv2/core/ocl.hpp>
#include <opencv2/photo.hpp>
#include <opencv2/imgproc.hpp>

#include "x3f_denoise_utils.h"
#include "x3f_denoise.h"
#include "x3f_io.h"
#include "x3f_printf.h"

using namespace cv;

static void denoise_nlm(Mat& img, float h)
{
  UMat out, sub, sub_dn, sub_res, res;
  float h1[3] = {0.0, h, h}, h2[3] = {0.0, h/8, h/4};

  x3f_printf(DEBUG, "BEGIN denoising\n");
  fastNlMeansDenoising(img, out, std::vector<float>(h1, h1+3),
		       3, 11, NORM_L1);
  x3f_printf(DEBUG, "END denoising\n");

  x3f_printf(DEBUG, "BEGIN V median filtering\n");
  UMat V(out.size(), CV_16U);
  int get_V[2] = { 2,0 }, set_V[2] = { 0,2 };
  mixChannels(std::vector<UMat>(1, out), std::vector<UMat>(2, V), get_V, 1);
  medianBlur(V, V, 3);
  mixChannels(std::vector<UMat>(1, V), std::vector<UMat>(2, out), set_V, 1);
  x3f_printf(DEBUG, "END V median filtering\n");

  x3f_printf(DEBUG, "BEGIN low-frequency denoising\n");
  resize(out, sub, Size(), 1.0/4, 1.0/4, INTER_AREA);
  fastNlMeansDenoising(sub, sub_dn, std::vector<float>(h2, h2+3),
		       3, 21, NORM_L1);
  subtract(sub, sub_dn, sub_res, noArray(), CV_16S);
  resize(sub_res, res, out.size(), 0.0, 0.0, INTER_CUBIC);
  subtract(out, res, out, noArray(), CV_16U);
  x3f_printf(DEBUG, "END low-frequency denoising\n");

  out.copyTo(img);
}

void x3f_denoise(x3f_area16_t *image, x3f_denoise_type_t type)
{
  assert(image->channels == 3);
  assert(type < sizeof(denoise_types)/sizeof(denoise_desc_t));
  const denoise_desc_t *d = &denoise_types[type];

  d->BMT_to_YUV(image);

  Mat img(image->rows, image->columns, CV_16UC3,
	 image->data, sizeof(uint16_t)*image->row_stride);
  denoise_nlm(img, d->h);

  d->YUV_to_BMT(image);
}

// NOTE: active has to be a subaera of image, i.e. they have to share
//       the same data area.
// NOTE: image, active and qtop will be destructively modified in place.
void x3f_expand_quattro(x3f_area16_t *image, x3f_area16_t *active,
			x3f_area16_t *qtop,
			x3f_area16_t *expanded, x3f_area16_t *active_exp)
{
  assert(image->channels == 3);
  assert(qtop->channels == 1);
  assert(X3F_DENOISE_F23 < sizeof(denoise_types)/sizeof(denoise_desc_t));
  const denoise_desc_t *d = &denoise_types[X3F_DENOISE_F23];

  d->BMT_to_YUV(image);

  Mat img(image->rows, image->columns, CV_16UC3,
	  image->data, sizeof(uint16_t)*image->row_stride);
  Mat qt(qtop->rows, qtop->columns, CV_16U,
	 qtop->data, sizeof(uint16_t)*qtop->row_stride);
  Mat exp(expanded->rows, expanded->columns, CV_16UC3,
	  expanded->data, sizeof(uint16_t)*expanded->row_stride);

  assert(qt.size() == exp.size());

  if (active) {
    assert(active->channels == 3);
    Mat act(active->rows, active->columns, CV_16UC3,
	    active->data, sizeof(uint16_t)*active->row_stride);
    denoise_nlm(act, d->h);
  }

  resize(img, exp, exp.size(), 0.0, 0.0, INTER_CUBIC);
  qt *= 4;
  int from_to[] = { 0,0 };
  mixChannels(&qt, 1, &exp, 1, from_to, 1);

  if (active_exp) {
    assert(active_exp->channels == 3);
    Mat act_exp(active_exp->rows, active_exp->columns, CV_16UC3,
		active_exp->data, sizeof(uint16_t)*active_exp->row_stride);
    UMat out;
    float h[3] = {0.0, d->h, d->h*2};

    x3f_printf(DEBUG, "BEGIN Quattro full-resolution denoising\n");
    fastNlMeansDenoising(act_exp, out, std::vector<float>(h, h+3),
			 3, 11, NORM_L1);
    x3f_printf(DEBUG, "END Quattro full-resolution denoising\n");

    out.copyTo(act_exp);
  }

  d->YUV_to_BMT(expanded);
}

void x3f_set_use_opencl(int flag)
{
  ocl::setUseOpenCL(flag);

  if (flag) {
    if (ocl::useOpenCL()) {
      ocl::Device dev = ocl::Device::getDefault();
      x3f_printf(INFO, "OpenCL device name: %s\n", dev.name().c_str());
      x3f_printf(INFO, "OpenCL device version: %s\n", dev.version().c_str());
    }
    else x3f_printf(WARN, "OpenCL is not available\n");
  }
  else x3f_printf(DEBUG, "OpenCL is disabled\n");
}


<|start_filename|>i386-apple-darwin11.cmake<|end_filename|>
SET(CMAKE_SYSTEM_NAME Darwin-GNU)
SET(CMAKE_C_COMPILER i386-apple-darwin11-gcc)
SET(CMAKE_CXX_COMPILER i386-apple-darwin11-g++)
SET(_CMAKE_TOOLCHAIN_PREFIX i386-apple-darwin11-)
EXECUTE_PROCESS(COMMAND i386-apple-darwin11-gcc -print-sysroot OUTPUT_VARIABLE CMAKE_FIND_ROOT_PATH OUTPUT_STRIP_TRAILING_WHITESPACE)
SET(CMAKE_FIND_ROOT_PATH_MODE_PROGRAM NEVER)
SET(CMAKE_FIND_ROOT_PATH_MODE_LIBRARY ONLY)
SET(CMAKE_FIND_ROOT_PATH_MODE_INCLUDE ONLY)


<|start_filename|>src/x3f_matrix_test.c<|end_filename|>
#include <stdio.h>
#include "x3f_matrix.h"

double a[9] = {11.0, 12.0, 13.0,
	       21.0, 22.0, 23.0,
	       31.0, 32.0, 33.0};
double b[9] = {0.0, 0.0, 1.0,
	       0.0, 1.0, 0.0,
	       1.0, 0.0, 0.0};
double x[3] = {1.0,
	       2.0,
	       3.0};


int main(int argc, char *argv[])
{
  double c[9];
  double y[3];
  double diag[9];
  double argb[9], argb_inv[9], identity[9];
  double srgb[9];
  double pprgb[9];

  x3f_3x3_3x3_mul(a, b, c);
  x3f_3x3_3x1_mul(c, x, y);
  x3f_3x3_diag(x, diag);
  x3f_XYZ_to_AdobeRGB(argb);
  x3f_3x3_inverse(argb, argb_inv);
  x3f_3x3_3x3_mul(argb_inv, argb, identity);
  x3f_XYZ_to_sRGB(srgb);
  x3f_XYZ_to_ProPhotoRGB(pprgb);

  printf("a\n");
  x3f_3x3_print(INFO, a);
  printf("b\n");
  x3f_3x3_print(INFO, b);
  printf("c\n");
  x3f_3x3_print(INFO, c);

  printf("x\n");
  x3f_3x1_print(INFO, x);
  printf("y\n");
  x3f_3x1_print(INFO, y);

  printf("diag\n");
  x3f_3x3_print(INFO, diag);

  printf("Adobe RGB\n");
  x3f_3x3_print(INFO, argb);
  printf("Adobe RGB inverse\n");
  x3f_3x3_print(INFO, argb_inv);
  printf("Identity\n");
  x3f_3x3_print(INFO, identity);
  printf("sRGB\n");
  x3f_3x3_print(INFO, srgb);
  printf("Pro Photo RGB\n");
  x3f_3x3_print(INFO, pprgb);

  return 0;
}


<|start_filename|>src/x3f_io.c<|end_filename|>
/* X3F_IO.C
 *
 * Library for accessing X3F Files.
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#include "x3f_io.h"
#include "x3f_printf.h"

#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>

#if defined(_WIN32) || defined (_WIN64)
#include <windows.h>
#else
#include <iconv.h>
#endif

/* --------------------------------------------------------------------- */
/* Hacky external flags                                                 */
/* --------------------------------------------------------------------- */

/* extern */ int legacy_offset = 0;
/* extern */ bool_t auto_legacy_offset = 1;

/* --------------------------------------------------------------------- */
/* Huffman Decode Macros                                                 */
/* --------------------------------------------------------------------- */

#define HUF_TREE_MAX_LENGTH 27
#define HUF_TREE_MAX_NODES(_leaves) ((HUF_TREE_MAX_LENGTH+1)*(_leaves))
#define HUF_TREE_GET_LENGTH(_v) (((_v)>>27)&0x1f)
#define HUF_TREE_GET_CODE(_v) ((_v)&0x07ffffff)

/* --------------------------------------------------------------------- */
/* Reading and writing - assuming little endian in the file              */
/* --------------------------------------------------------------------- */

static int x3f_get1(FILE *f)
{
  /* Little endian file */
  return getc(f);
}

static int x3f_get2(FILE *f)
{
  /* Little endian file */
  return (getc(f)<<0) + (getc(f)<<8);
}

static int x3f_get4(FILE *f)
{
  /* Little endian file */
  return (getc(f)<<0) + (getc(f)<<8) + (getc(f)<<16) + (getc(f)<<24);
}

#define FREE(P) do { free(P); (P) = NULL; } while (0)

#define PUT_GET_N(_buffer,_size,_file,_func)			\
  do								\
    {								\
      int _left = _size;					\
      while (_left != 0) {					\
	int _cur = _func(_buffer,1,_left,_file);		\
	if (_cur == 0) {					\
	  x3f_printf(ERR, "Failure to access file\n");		\
	  exit(1);						\
	}							\
	_left -= _cur;						\
      }								\
    } while(0)

#define GET1(_v) do {(_v) = x3f_get1(I->input.file);} while (0)
#define GET2(_v) do {(_v) = x3f_get2(I->input.file);} while (0)
#define GET4(_v) do {(_v) = x3f_get4(I->input.file);} while (0)
#define GET4F(_v)				\
  do {						\
    union {int32_t i; float f;} _tmp;		\
    _tmp.i = x3f_get4(I->input.file);		\
    (_v) = _tmp.f;				\
  } while (0)
#define GETN(_v,_s) PUT_GET_N(_v,_s,I->input.file,fread)

#define GET_TABLE(_T, _GETX, _NUM)					\
  do {									\
    int _i;								\
    (_T).size = (_NUM);							\
    (_T).element = (void *)realloc((_T).element,			\
				   (_NUM)*sizeof((_T).element[0]));	\
    for (_i = 0; _i < (_T).size; _i++)					\
      _GETX((_T).element[_i]);						\
  } while (0)

#define GET_PROPERTY_TABLE(_T, _NUM)					\
  do {									\
    int _i;								\
    (_T).size = (_NUM);							\
    (_T).element = (void *)realloc((_T).element,			\
				   (_NUM)*sizeof((_T).element[0]));	\
    for (_i = 0; _i < (_T).size; _i++) {				\
      GET4((_T).element[_i].name_offset);				\
      GET4((_T).element[_i].value_offset);				\
    }									\
  } while (0)

#define GET_TRUE_HUFF_TABLE(_T)						\
  do {									\
    int _i;								\
    (_T).element = NULL;						\
    for (_i = 0; ; _i++) {						\
      (_T).size = _i + 1;						\
      (_T).element = (void *)realloc((_T).element,			\
				     (_i + 1)*sizeof((_T).element[0]));	\
      GET1((_T).element[_i].code_size);					\
      GET1((_T).element[_i].code);					\
      if ((_T).element[_i].code_size == 0) break;			\
    }									\
  } while (0)

/* --------------------------------------------------------------------- */
/* Allocating Huffman tree help data                                   */
/* --------------------------------------------------------------------- */

static void cleanup_huffman_tree(x3f_hufftree_t *HTP)
{
  free(HTP->nodes);
}

static void new_huffman_tree(x3f_hufftree_t *HTP, int bits)
{
  int leaves = 1<<bits;

  HTP->free_node_index = 0;
  HTP->nodes = (x3f_huffnode_t *)
    calloc(1, HUF_TREE_MAX_NODES(leaves)*sizeof(x3f_huffnode_t));
}

/* --------------------------------------------------------------------- */
/* Allocating TRUE engine RAW help data                                  */
/* --------------------------------------------------------------------- */

static void cleanup_true(x3f_true_t **TRUP)
{
  x3f_true_t *TRU = *TRUP;

  if (TRU == NULL) return;

  x3f_printf(DEBUG, "Cleanup TRUE data\n");

  FREE(TRU->table.element);
  FREE(TRU->plane_size.element);
  cleanup_huffman_tree(&TRU->tree);
  FREE(TRU->x3rgb16.buf);

  FREE(TRU);

  *TRUP = NULL;
}

static x3f_true_t *new_true(x3f_true_t **TRUP)
{
  x3f_true_t *TRU = (x3f_true_t *)calloc(1, sizeof(x3f_true_t));

  cleanup_true(TRUP);

  TRU->table.size = 0;
  TRU->table.element = NULL;
  TRU->plane_size.size = 0;
  TRU->plane_size.element = NULL;
  TRU->tree.nodes = NULL;
  TRU->x3rgb16.data = NULL;
  TRU->x3rgb16.buf = NULL;

  *TRUP = TRU;

  return TRU;
}

static void cleanup_quattro(x3f_quattro_t **QP)
{
  x3f_quattro_t *Q = *QP;

  if (Q == NULL) return;

  x3f_printf(DEBUG, "Cleanup Quattro\n");

  FREE(Q->top16.buf);
  FREE(Q);

  *QP = NULL;
}

static x3f_quattro_t *new_quattro(x3f_quattro_t **QP)
{
  x3f_quattro_t *Q = (x3f_quattro_t *)calloc(1, sizeof(x3f_quattro_t));
  int i;

  cleanup_quattro(QP);

  for (i=0; i<TRUE_PLANES; i++) {
    Q->plane[i].columns = 0;
    Q->plane[i].rows = 0;
  }

  Q->unknown = 0;

  Q->top16.data = NULL;
  Q->top16.buf = NULL;

  *QP = Q;

  return Q;
}

/* --------------------------------------------------------------------- */
/* Allocating Huffman engine help data                                   */
/* --------------------------------------------------------------------- */

static void cleanup_huffman(x3f_huffman_t **HUFP)
{
  x3f_huffman_t *HUF = *HUFP;

  if (HUF == NULL) return;

  x3f_printf(DEBUG, "Cleanup Huffman\n");

  FREE(HUF->mapping.element);
  FREE(HUF->table.element);
  cleanup_huffman_tree(&HUF->tree);
  FREE(HUF->row_offsets.element);
  FREE(HUF->rgb8.buf);
  FREE(HUF->x3rgb16.buf);
  FREE(HUF);

  *HUFP = NULL;
}

static x3f_huffman_t *new_huffman(x3f_huffman_t **HUFP)
{
  x3f_huffman_t *HUF = (x3f_huffman_t *)calloc(1, sizeof(x3f_huffman_t));

  cleanup_huffman(HUFP);

  /* Set all not read data block pointers to NULL */
  HUF->mapping.size = 0;
  HUF->mapping.element = NULL;
  HUF->table.size = 0;
  HUF->table.element = NULL;
  HUF->tree.nodes = NULL;
  HUF->row_offsets.size = 0;
  HUF->row_offsets.element = NULL;
  HUF->rgb8.data = NULL;
  HUF->rgb8.buf = NULL;
  HUF->x3rgb16.data = NULL;
  HUF->x3rgb16.buf = NULL;

  *HUFP = HUF;

  return HUF;
}

/* --------------------------------------------------------------------- */
/* Creating a new x3f structure from file                                */
/* --------------------------------------------------------------------- */

/* extern */ x3f_t *x3f_new_from_file(FILE *infile)
{
  x3f_t *x3f = (x3f_t *)calloc(1, sizeof(x3f_t));
  x3f_info_t *I = NULL;
  x3f_header_t *H = NULL;
  x3f_directory_section_t *DS = NULL;
  int i, d;

  I = &x3f->info;
  I->error = NULL;
  I->input.file = infile;
  I->output.file = NULL;

  if (infile == NULL) {
    I->error = "No infile";
    return x3f;
  }

  /* Read file header */
  H = &x3f->header;
  fseek(infile, 0, SEEK_SET);
  GET4(H->identifier);

  if (H->identifier != X3F_FOVb) {
    x3f_printf(ERR, "Faulty file type\n");
    x3f_delete(x3f);
    return NULL;
  }

  GET4(H->version);
  GETN(H->unique_identifier, SIZE_UNIQUE_IDENTIFIER);
  /* TODO: the meaning of the rest of the header for version >= 4.0
           (Quattro) is unknown */
  if (H->version < X3F_VERSION_4_0) {
    GET4(H->mark_bits);
    GET4(H->columns);
    GET4(H->rows);
    GET4(H->rotation);
    if (H->version >= X3F_VERSION_2_1) {
      int num_ext_data =
	H->version >= X3F_VERSION_3_0 ? NUM_EXT_DATA_3_0 : NUM_EXT_DATA_2_1;

      GETN(H->white_balance, SIZE_WHITE_BALANCE);
      if (H->version >= X3F_VERSION_2_3)
	GETN(H->color_mode, SIZE_COLOR_MODE);
      GETN(H->extended_types, num_ext_data);
      for (i=0; i<num_ext_data; i++)
	GET4F(H->extended_data[i]);
    }
  }

  /* Go to the beginning of the directory */
  fseek(infile, -4, SEEK_END);
  fseek(infile, x3f_get4(infile), SEEK_SET);

  /* Read the directory header */
  DS = &x3f->directory_section;
  GET4(DS->identifier);
  GET4(DS->version);
  GET4(DS->num_directory_entries);

  if (DS->num_directory_entries > 0) {
    size_t size = DS->num_directory_entries * sizeof(x3f_directory_entry_t);
    DS->directory_entry = (x3f_directory_entry_t *)calloc(1, size);
  }

  /* Traverse the directory */
  for (d=0; d<DS->num_directory_entries; d++) {
    x3f_directory_entry_t *DE = &DS->directory_entry[d];
    x3f_directory_entry_header_t *DEH = &DE->header;
    uint32_t save_dir_pos;

    /* Read the directory entry info */
    GET4(DE->input.offset);
    GET4(DE->input.size);

    DE->output.offset = 0;
    DE->output.size = 0;

    GET4(DE->type);

    /* Save current pos and go to the entry */
    save_dir_pos = ftell(infile);
    fseek(infile, DE->input.offset, SEEK_SET);

    /* Read the type independent part of the entry header */
    DEH = &DE->header;
    GET4(DEH->identifier);
    GET4(DEH->version);

    /* NOTE - the tests below could be made on DE->type instead */

    if (DEH->identifier == X3F_SECp) {
      x3f_property_list_t *PL = &DEH->data_subsection.property_list;

      /* Read the property part of the header */
      GET4(PL->num_properties);
      GET4(PL->character_format);
      GET4(PL->reserved);
      GET4(PL->total_length);

      /* Set all not read data block pointers to NULL */
      PL->data = NULL;
      PL->data_size = 0;
    }

    if (DEH->identifier == X3F_SECi) {
      x3f_image_data_t *ID = &DEH->data_subsection.image_data;

      /* Read the image part of the header */
      GET4(ID->type);
      GET4(ID->format);
      ID->type_format = (ID->type << 16) + (ID->format);
      GET4(ID->columns);
      GET4(ID->rows);
      GET4(ID->row_stride);

      /* Set all not read data block pointers to NULL */
      ID->huffman = NULL;

      ID->data = NULL;
      ID->data_size = 0;
    }

    if (DEH->identifier == X3F_SECc) {
      x3f_camf_t *CAMF = &DEH->data_subsection.camf;

      /* Read the CAMF part of the header */
      GET4(CAMF->type);
      GET4(CAMF->tN.val0);
      GET4(CAMF->tN.val1);
      GET4(CAMF->tN.val2);
      GET4(CAMF->tN.val3);

      /* Set all not read data block pointers to NULL */
      CAMF->data = NULL;
      CAMF->data_size = 0;

      /* Set all not allocated help pointers to NULL */
      CAMF->table.element = NULL;
      CAMF->table.size = 0;
      CAMF->tree.nodes = NULL;
      CAMF->decoded_data = NULL;
      CAMF->decoded_data_size = 0;
      CAMF->entry_table.element = NULL;
      CAMF->entry_table.size = 0;
    }

    /* Reset the file pointer back to the directory */
    fseek(infile, save_dir_pos, SEEK_SET);
  }

  return x3f;
}

/* --------------------------------------------------------------------- */
/* Clean up an x3f structure                                             */
/* --------------------------------------------------------------------- */

static void free_camf_entry(camf_entry_t *entry)
{
  FREE(entry->property_name);
  FREE(entry->property_value);
  FREE(entry->matrix_decoded);
  FREE(entry->matrix_dim_entry);
}

/* extern */ x3f_return_t x3f_delete(x3f_t *x3f)
{
  x3f_directory_section_t *DS;
  int d;

  if (x3f == NULL)
    return X3F_ARGUMENT_ERROR;

  x3f_printf(DEBUG, "X3F Delete\n");

  DS = &x3f->directory_section;

  for (d=0; d<DS->num_directory_entries; d++) {
    x3f_directory_entry_t *DE = &DS->directory_entry[d];
    x3f_directory_entry_header_t *DEH = &DE->header;

    if (DEH->identifier == X3F_SECp) {
      x3f_property_list_t *PL = &DEH->data_subsection.property_list;
      int i;

      for (i=0; i<PL->property_table.size; i++) {
	FREE(PL->property_table.element[i].name_utf8);
	FREE(PL->property_table.element[i].value_utf8);
      }

      FREE(PL->property_table.element);
      FREE(PL->data);
    }

    if (DEH->identifier == X3F_SECi) {
      x3f_image_data_t *ID = &DEH->data_subsection.image_data;

      cleanup_huffman(&ID->huffman);

      cleanup_true(&ID->tru);

      cleanup_quattro(&ID->quattro);

      FREE(ID->data);
    }

    if (DEH->identifier == X3F_SECc) {
      x3f_camf_t *CAMF = &DEH->data_subsection.camf;
      int i;

      FREE(CAMF->data);
      FREE(CAMF->table.element);
      cleanup_huffman_tree(&CAMF->tree);
      FREE(CAMF->decoded_data);
      for (i=0; i < CAMF->entry_table.size; i++) {
	free_camf_entry(&CAMF->entry_table.element[i]);
      }
      FREE(CAMF->entry_table.element);
    }
  }

  FREE(DS->directory_entry);
  FREE(x3f);

  return X3F_OK;
}

/* --------------------------------------------------------------------- */
/* Getting a reference to a directory entry                              */
/* --------------------------------------------------------------------- */

/* TODO: all those only get the first instance */

static x3f_directory_entry_t *x3f_get(x3f_t *x3f,
                                      uint32_t type,
                                      uint32_t image_type)
{
  x3f_directory_section_t *DS;
  int d;

  if (x3f == NULL) return NULL;

  DS = &x3f->directory_section;

  for (d=0; d<DS->num_directory_entries; d++) {
    x3f_directory_entry_t *DE = &DS->directory_entry[d];
    x3f_directory_entry_header_t *DEH = &DE->header;

    if (DEH->identifier == type) {
      switch (DEH->identifier) {
      case X3F_SECi:
        {
          x3f_image_data_t *ID = &DEH->data_subsection.image_data;

          if (ID->type_format == image_type)
            return DE;
        }
        break;
      default:
        return DE;
      }
    }
  }

  return NULL;
}

/* extern */ x3f_directory_entry_t *x3f_get_raw(x3f_t *x3f)
{
  x3f_directory_entry_t *DE;

  if ((DE = x3f_get(x3f, X3F_SECi, X3F_IMAGE_RAW_HUFFMAN_X530)) != NULL)
    return DE;

  if ((DE = x3f_get(x3f, X3F_SECi, X3F_IMAGE_RAW_HUFFMAN_10BIT)) != NULL)
    return DE;

  if ((DE = x3f_get(x3f, X3F_SECi, X3F_IMAGE_RAW_TRUE)) != NULL)
    return DE;

  if ((DE = x3f_get(x3f, X3F_SECi, X3F_IMAGE_RAW_MERRILL)) != NULL)
    return DE;

  if ((DE = x3f_get(x3f, X3F_SECi, X3F_IMAGE_RAW_QUATTRO)) != NULL)
    return DE;

  if ((DE = x3f_get(x3f, X3F_SECi, X3F_IMAGE_RAW_SDQ)) != NULL)
    return DE;

  if ((DE = x3f_get(x3f, X3F_SECi, X3F_IMAGE_RAW_SDQH)) != NULL)
    return DE;

  return NULL;
}

/* extern */ x3f_directory_entry_t *x3f_get_thumb_plain(x3f_t *x3f)
{
  return x3f_get(x3f, X3F_SECi, X3F_IMAGE_THUMB_PLAIN);
}

/* extern */ x3f_directory_entry_t *x3f_get_thumb_huffman(x3f_t *x3f)
{
  return x3f_get(x3f, X3F_SECi, X3F_IMAGE_THUMB_HUFFMAN);
}

/* extern */ x3f_directory_entry_t *x3f_get_thumb_jpeg(x3f_t *x3f)
{
  return x3f_get(x3f, X3F_SECi, X3F_IMAGE_THUMB_JPEG);
}

/* extern */ x3f_directory_entry_t *x3f_get_camf(x3f_t *x3f)
{
  return x3f_get(x3f, X3F_SECc, 0);
}

/* extern */ x3f_directory_entry_t *x3f_get_prop(x3f_t *x3f)
{
  return x3f_get(x3f, X3F_SECp, 0);
}

/* For some obscure reason, the bit numbering is weird. It is
   generally some kind of "big endian" style - e.g. the bit 7 is the
   first in a byte and bit 31 first in a 4 byte int. For patterns in
   the huffman pattern table, bit 27 is the first bit and bit 26 the
   next one. */

#define PATTERN_BIT_POS(_len, _bit) ((_len) - (_bit) - 1)
#define MEMORY_BIT_POS(_bit) PATTERN_BIT_POS(8, _bit)

/* --------------------------------------------------------------------- */
/* Huffman Decode                                                        */
/* --------------------------------------------------------------------- */

/* Make the huffman tree */

#ifdef DBG_PRNT
static char *display_code(int length, uint32_t code, char *buffer)
{
  int i;

  for (i=0; i<length; i++) {
    int pos = PATTERN_BIT_POS(length, i);
    buffer[i] = ((code>>pos)&1) == 0 ? '0' : '1';
  }

  buffer[i] = 0;

  return buffer;
}
#endif

static x3f_huffnode_t *new_node(x3f_hufftree_t *tree)
{
  x3f_huffnode_t *t = &tree->nodes[tree->free_node_index];

  t->branch[0] = NULL;
  t->branch[1] = NULL;
  t->leaf = UNDEFINED_LEAF;

  tree->free_node_index++;

  return t;
}

static void add_code_to_tree(x3f_hufftree_t *tree,
                             int length, uint32_t code, uint32_t value)
{
  int i;

  x3f_huffnode_t *t = tree->nodes;

  for (i=0; i<length; i++) {
    int pos = PATTERN_BIT_POS(length, i);
    int bit = (code>>pos)&1;
    x3f_huffnode_t *t_next = t->branch[bit];

    if (t_next == NULL)
      t_next = t->branch[bit] = new_node(tree);

    t = t_next;
  }

  t->leaf = value;
}

static void populate_true_huffman_tree(x3f_hufftree_t *tree,
				       x3f_true_huffman_t *table)
{
  int i;

  new_node(tree);

  for (i=0; i<table->size; i++) {
    x3f_true_huffman_element_t *element = &table->element[i];
    uint32_t length = element->code_size;

    if (length != 0) {
      /* add_code_to_tree wants the code right adjusted */
      uint32_t code = ((element->code) >> (8 - length)) & 0xff;
      uint32_t value = i;

      add_code_to_tree(tree, length, code, value);

#ifdef DBG_PRNT
      {
        char buffer[100];

        x3f_printf(DEBUG, "H %5d : %5x : %5d : %02x %08x (%08x) (%s)\n",
		   i, i, value, length, code, value,
		   display_code(length, code, buffer));
      }
#endif
    }
  }
}

static void populate_huffman_tree(x3f_hufftree_t *tree,
				  x3f_table32_t *table,
				  x3f_table16_t *mapping)
{
  int i;

  new_node(tree);

  for (i=0; i<table->size; i++) {
    uint32_t element = table->element[i];

    if (element != 0) {
      uint32_t length = HUF_TREE_GET_LENGTH(element);
      uint32_t code = HUF_TREE_GET_CODE(element);
      uint32_t value;

      /* If we have a valid mapping table - then the value from the
         mapping table shall be used. Otherwise we use the current
         index in the table as value. */
      if (table->size == mapping->size)
        value = mapping->element[i];
      else
        value = i;

      add_code_to_tree(tree, length, code, value);

#ifdef DBG_PRNT
      {
        char buffer[100];

        x3f_printf(DEBUG, "H %5d : %5x : %5d : %02x %08x (%08x) (%s)\n",
		   i, i, value, length, code, element,
		   display_code(length, code, buffer));
      }
#endif
    }
  }
}

#ifdef DBG_PRNT
static void print_huffman_tree(x3f_huffnode_t *t, int length, uint32_t code)
{
  char buf1[100];
  char buf2[100];

  x3f_printf(DEBUG, "%*s (%s,%s) %s (%s)\n",
	     length, length < 1 ? "-" : (code&1) ? "1" : "0",
	     t->branch[0]==NULL ? "-" : "0",
	     t->branch[1]==NULL ? "-" : "1",
	     t->leaf==UNDEFINED_LEAF ? "-" : (sprintf(buf1, "%x", t->leaf),buf1),
	     display_code(length, code, buf2));

  code = code << 1;
  if (t->branch[0]) print_huffman_tree(t->branch[0], length+1, code+0);
  if (t->branch[1]) print_huffman_tree(t->branch[1], length+1, code+1);
}
#endif

/* Help machinery for reading bits in a memory */

typedef struct bit_state_s {
  uint8_t *next_address;
  uint8_t bit_offset;
  uint8_t bits[8];
} bit_state_t;

static void set_bit_state(bit_state_t *BS, uint8_t *address)
{
  BS->next_address = address;
  BS->bit_offset = 8;
}

static uint8_t get_bit(bit_state_t *BS)
{
  if (BS->bit_offset == 8) {
    uint8_t byte = *BS->next_address;
    int i;

    for (i=7; i>= 0; i--) {
      BS->bits[i] = byte&1;
      byte = byte >> 1;
    }
    BS->next_address++;
    BS->bit_offset = 0;
  }

  return BS->bits[BS->bit_offset++];
}

/* Decode use the TRUE algorithm */

static int32_t get_true_diff(bit_state_t *BS, x3f_hufftree_t *HTP)
{
  int32_t diff;
  x3f_huffnode_t *node = &HTP->nodes[0];
  uint8_t bits;

  while (node->branch[0] != NULL || node->branch[1] != NULL) {
    uint8_t bit = get_bit(BS);
    x3f_huffnode_t *new_node = node->branch[bit];

    node = new_node;
    if (node == NULL) {
      /* TODO: Shouldn't this be treated as a fatal error? */
      x3f_printf(ERR, "Huffman coding got unexpected bit\n");
      return 0;
    }
  }

  bits = node->leaf;

  if (bits == 0)
    diff = 0;
  else {
    uint8_t first_bit = get_bit(BS);
    int i;

    diff = first_bit;

    for (i=1; i<bits; i++)
      diff = (diff << 1) + get_bit(BS);

    if (first_bit == 0)
      diff -= (1<<bits) - 1;
  }

  return diff;
}

/* This code (that decodes one of the X3F color planes, really is a
   decoding of a compression algorithm suited for Bayer CFA data. In
   Bayer CFA the data is divided into 2x2 squares that represents
   (R,G1,G2,B) data. Those four positions are (in this compression)
   treated as one data stream each, where you store the differences to
   previous data in the stream. The reason for this is, of course,
   that the date is more often than not near to the next data in a
   stream that represents the same color. */

/* TODO: write more about the compression */

static void true_decode_one_color(x3f_image_data_t *ID, int color)
{
  x3f_true_t *TRU = ID->tru;
  x3f_quattro_t *Q = ID->quattro;
  uint32_t seed = TRU->seed[color]; /* TODO : Is this correct ? */
  int row;

  x3f_hufftree_t *tree = &TRU->tree;
  bit_state_t BS;

  int32_t row_start_acc[2][2];
  uint32_t rows = ID->rows;
  uint32_t cols = ID->columns;
  x3f_area16_t *area = &TRU->x3rgb16;
  uint16_t *dst = area->data + color;

  set_bit_state(&BS, TRU->plane_address[color]);

  row_start_acc[0][0] = seed;
  row_start_acc[0][1] = seed;
  row_start_acc[1][0] = seed;
  row_start_acc[1][1] = seed;

  if (ID->type_format == X3F_IMAGE_RAW_QUATTRO ||
      ID->type_format == X3F_IMAGE_RAW_SDQ ||
      ID->type_format == X3F_IMAGE_RAW_SDQH) {
    rows = Q->plane[color].rows;
    cols = Q->plane[color].columns;

    if (Q->quattro_layout && color == 2) {
      area = &Q->top16;
      dst = area->data;
    }
    x3f_printf(DEBUG, "Quattro decode one color (%d) rows=%d cols=%d\n",
	       color, rows, cols);
  } else {
    x3f_printf(DEBUG, "TRUE decode one color (%d) rows=%d cols=%d\n",
	       color, rows, cols);
  }

  assert(rows == area->rows && cols >= area->columns);

  for (row = 0; row < rows; row++) {
    int col;
    bool_t odd_row = row&1;
    int32_t acc[2];

    for (col = 0; col < cols; col++) {
      bool_t odd_col = col&1;
      int32_t diff = get_true_diff(&BS, tree);
      int32_t prev = col < 2 ?
	row_start_acc[odd_row][odd_col] :
	acc[odd_col];
      int32_t value = prev + diff;

      acc[odd_col] = value;
      if (col < 2)
	row_start_acc[odd_row][odd_col] = value;

      /* Discard additional data at the right for binned Quattro plane 2 */
      if (col >= area->columns) continue;

      *dst = value;
      dst += area->channels;
    }
  }
}

static void true_decode(x3f_info_t *I,
			x3f_directory_entry_t *DE)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_image_data_t *ID = &DEH->data_subsection.image_data;
  int color;

  for (color = 0; color < 3; color++) {
    true_decode_one_color(ID, color);
  }
}

/* Decode use the huffman tree */

static int32_t get_huffman_diff(bit_state_t *BS, x3f_hufftree_t *HTP)
{
  int32_t diff;
  x3f_huffnode_t *node = &HTP->nodes[0];

  while (node->branch[0] != NULL || node->branch[1] != NULL) {
    uint8_t bit = get_bit(BS);
    x3f_huffnode_t *new_node = node->branch[bit];

    node = new_node;
    if (node == NULL) {
      /* TODO: Shouldn't this be treated as a fatal error? */
      x3f_printf(ERR, "Huffman coding got unexpected bit\n");
      return 0;
    }
  }

  diff = node->leaf;

  return diff;
}

static void huffman_decode_row(x3f_info_t *I,
                               x3f_directory_entry_t *DE,
                               int bits,
                               int row,
                               int offset,
                               int *minimum)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_image_data_t *ID = &DEH->data_subsection.image_data;
  x3f_huffman_t *HUF = ID->huffman;

  int16_t c[3] = {offset,offset,offset};
  int col;
  bit_state_t BS;

  set_bit_state(&BS, ID->data + HUF->row_offsets.element[row]);

  for (col = 0; col < ID->columns; col++) {
    int color;

    for (color = 0; color < 3; color++) {
      uint16_t c_fix;

      c[color] += get_huffman_diff(&BS, &HUF->tree);
      if (c[color] < 0) {
        c_fix = 0;
        if (c[color] < *minimum)
          *minimum = c[color];
      } else {
        c_fix = c[color];
      }

      switch (ID->type_format) {
      case X3F_IMAGE_RAW_HUFFMAN_X530:
      case X3F_IMAGE_RAW_HUFFMAN_10BIT:
        HUF->x3rgb16.data[3*(row*ID->columns + col) + color] = (uint16_t)c_fix;
        break;
      case X3F_IMAGE_THUMB_HUFFMAN:
        HUF->rgb8.data[3*(row*ID->columns + col) + color] = (uint8_t)c_fix;
        break;
      default:
	/* TODO: Shouldn't this be treated as a fatal error? */
        x3f_printf(ERR, "Unknown huffman image type\n");
      }
    }
  }
}

static void huffman_decode(x3f_info_t *I,
                           x3f_directory_entry_t *DE,
                           int bits)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_image_data_t *ID = &DEH->data_subsection.image_data;

  int row;
  int minimum = 0;
  int offset = legacy_offset;

  x3f_printf(DEBUG, "Huffman decode with offset: %d\n", offset);
  for (row = 0; row < ID->rows; row++)
    huffman_decode_row(I, DE, bits, row, offset, &minimum);

  if (auto_legacy_offset && minimum < 0) {
    offset = -minimum;
    x3f_printf(DEBUG, "Redo with offset: %d\n", offset);
    for (row = 0; row < ID->rows; row++)
      huffman_decode_row(I, DE, bits, row, offset, &minimum);
  }
}

static int32_t get_simple_diff(x3f_huffman_t *HUF, uint16_t index)
{
  if (HUF->mapping.size == 0)
    return index;
  else
    return HUF->mapping.element[index];
}

static void simple_decode_row(x3f_info_t *I,
                              x3f_directory_entry_t *DE,
                              int bits,
                              int row,
                              int row_stride)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_image_data_t *ID = &DEH->data_subsection.image_data;
  x3f_huffman_t *HUF = ID->huffman;

  uint32_t *data = (uint32_t *)(ID->data + row*row_stride);

  uint16_t c[3] = {0,0,0};
  int col;

  uint32_t mask = 0;

  switch (bits) {
  case 8:
    mask = 0x0ff;
    break;
  case 9:
    mask = 0x1ff;
    break;
  case 10:
    mask = 0x3ff;
    break;
  case 11:
    mask = 0x7ff;
    break;
  case 12:
    mask = 0xfff;
    break;
  default:
    /* TODO: Shouldn't this be treated as a fatal error? */
    x3f_printf(ERR, "Unknown number of bits: %d\n", bits);
    mask = 0;
    break;
  }

  for (col = 0; col < ID->columns; col++) {
    int color;
    uint32_t val = data[col];

    for (color = 0; color < 3; color++) {
      uint16_t c_fix;
      c[color] += get_simple_diff(HUF, (val>>(color*bits))&mask);

      switch (ID->type_format) {
      case X3F_IMAGE_RAW_HUFFMAN_X530:
      case X3F_IMAGE_RAW_HUFFMAN_10BIT:
        c_fix = (int16_t)c[color] > 0 ? c[color] : 0;

        HUF->x3rgb16.data[3*(row*ID->columns + col) + color] = c_fix;
        break;
      case X3F_IMAGE_THUMB_HUFFMAN:
        c_fix = (int8_t)c[color] > 0 ? c[color] : 0;

        HUF->rgb8.data[3*(row*ID->columns + col) + color] = c_fix;
        break;
      default:
	/* TODO: Shouldn't this be treated as a fatal error? */
        x3f_printf(ERR, "Unknown huffman image type\n");
      }
    }
  }
}

static void simple_decode(x3f_info_t *I,
                          x3f_directory_entry_t *DE,
                          int bits,
                          int row_stride)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_image_data_t *ID = &DEH->data_subsection.image_data;

  int row;

  for (row = 0; row < ID->rows; row++)
    simple_decode_row(I, DE, bits, row, row_stride);
}

/* --------------------------------------------------------------------- */
/* Loading the data in a directory entry                                 */
/* --------------------------------------------------------------------- */

/* First you set the offset to where to start reading the data ... */

static void read_data_set_offset(x3f_info_t *I,
                                 x3f_directory_entry_t *DE,
                                 uint32_t header_size)
{
  uint32_t i_off = DE->input.offset + header_size;

  fseek(I->input.file, i_off, SEEK_SET);
}

/* ... then you read the data, block for block */

static uint32_t read_data_block(void **data,
                                x3f_info_t *I,
                                x3f_directory_entry_t *DE,
                                uint32_t footer)
{
  uint32_t size =
    DE->input.size + DE->input.offset - ftell(I->input.file) - footer;

  *data = (void *)malloc(size);

  GETN(*data, size);

  return size;
}

static void x3f_load_image_verbatim(x3f_info_t *I, x3f_directory_entry_t *DE)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_image_data_t *ID = &DEH->data_subsection.image_data;

  x3f_printf(DEBUG, "Load image verbatim\n");

  ID->data_size = read_data_block(&ID->data, I, DE, 0);
}

#if defined(_WIN32) || defined (_WIN64)
static char *utf16le_to_utf8(utf16_t *str)
{
  size_t osize = WideCharToMultiByte(CP_UTF8, 0, str, -1, NULL, 0, NULL, NULL);
  char *buf = malloc(osize);

  WideCharToMultiByte(CP_UTF8, 0, str, -1, buf, osize, NULL, NULL);

  return buf;
}
#else
static char *utf16le_to_utf8(utf16_t *str)
{
  iconv_t ic = iconv_open("UTF-8", "UTF-16LE");
  size_t isize, osize;
  char *buf, *ibuf, *obuf;

  assert(ic != (iconv_t)-1);

  for (isize=0; str[isize]; isize++);
  isize *= 2;			/* Size in bytes */
  osize = 2*isize;		/* Worst case scenario */

  buf = malloc(osize+1);
  ibuf = (char *)str;
  obuf = buf;

  assert(iconv(ic, &ibuf, &isize, &obuf, &osize) != -1);
  *obuf = 0;

  iconv_close(ic);

  return realloc(buf, obuf-buf+1);
}
#endif

static void x3f_load_property_list(x3f_info_t *I, x3f_directory_entry_t *DE)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_property_list_t *PL = &DEH->data_subsection.property_list;
  int i;

  read_data_set_offset(I, DE, X3F_PROPERTY_LIST_HEADER_SIZE);

  GET_PROPERTY_TABLE(PL->property_table, PL->num_properties);

  PL->data_size = read_data_block(&PL->data, I, DE, 0);

  for (i=0; i<PL->num_properties; i++) {
    x3f_property_t *P = &PL->property_table.element[i];

    P->name = ((utf16_t *)PL->data + P->name_offset);
    P->value = ((utf16_t *)PL->data + P->value_offset);
    P->name_utf8 = utf16le_to_utf8(P->name);
    P->value_utf8 = utf16le_to_utf8(P->value);
  }
}

static void x3f_load_true(x3f_info_t *I,
			  x3f_directory_entry_t *DE)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_image_data_t *ID = &DEH->data_subsection.image_data;
  x3f_true_t *TRU = new_true(&ID->tru);
  x3f_quattro_t *Q = NULL;
  int i;

  if (ID->type_format == X3F_IMAGE_RAW_QUATTRO ||
      ID->type_format == X3F_IMAGE_RAW_SDQ ||
      ID->type_format == X3F_IMAGE_RAW_SDQH) {
    x3f_printf(DEBUG, "Load Quattro extra info\n");

    Q = new_quattro(&ID->quattro);

    for (i=0; i<TRUE_PLANES; i++) {
      GET2(Q->plane[i].columns);
      GET2(Q->plane[i].rows);
    }

    if (Q->plane[0].rows == ID->rows/2) {
      x3f_printf(DEBUG, "Quattro layout\n");
      Q->quattro_layout = 1;
    } else if (Q->plane[0].rows == ID->rows) {
      x3f_printf(DEBUG, "Binned Quattro\n");
      Q->quattro_layout = 0;
    } else {
      x3f_printf(ERR, "Quattro file with unknown layer size\n");
      assert(0);
    }
  }

  x3f_printf(DEBUG, "Load TRUE\n");

  /* Read TRUE header data */
  GET2(TRU->seed[0]);
  GET2(TRU->seed[1]);
  GET2(TRU->seed[2]);
  GET2(TRU->unknown);
  GET_TRUE_HUFF_TABLE(TRU->table);

  if (ID->type_format == X3F_IMAGE_RAW_QUATTRO ||
      ID->type_format == X3F_IMAGE_RAW_SDQ ||
      ID->type_format == X3F_IMAGE_RAW_SDQH) {

    x3f_printf(DEBUG, "Load Quattro extra info 2\n");

    GET4(Q->unknown);
  }

  GET_TABLE(TRU->plane_size, GET4, TRUE_PLANES);

  /* Read image data */
  ID->data_size = read_data_block(&ID->data, I, DE, 0);

  /* TODO: can it be fewer than 8 bits? Maybe taken from TRU->table? */
  new_huffman_tree(&TRU->tree, 8);

  populate_true_huffman_tree(&TRU->tree, &TRU->table);

#ifdef DBG_PRNT
  print_huffman_tree(TRU->tree.nodes, 0, 0);
#endif

  TRU->plane_address[0] = ID->data;
  for (i=1; i<TRUE_PLANES; i++)
    TRU->plane_address[i] =
      TRU->plane_address[i-1] +
      (((TRU->plane_size.element[i-1] + 15) / 16) * 16);

  if ( (ID->type_format == X3F_IMAGE_RAW_QUATTRO ||
	ID->type_format == X3F_IMAGE_RAW_SDQ ||
	ID->type_format == X3F_IMAGE_RAW_SDQH ) &&
       Q->quattro_layout) {
    uint32_t columns = Q->plane[0].columns;
    uint32_t rows = Q->plane[0].rows;
    uint32_t channels = 3;
    uint32_t size = columns * rows * channels;

    TRU->x3rgb16.columns = columns;
    TRU->x3rgb16.rows = rows;
    TRU->x3rgb16.channels = channels;
    TRU->x3rgb16.row_stride = columns * channels;
    TRU->x3rgb16.data = TRU->x3rgb16.buf =
      (uint16_t *)malloc(sizeof(uint16_t)*size);

    columns = Q->plane[2].columns;
    rows = Q->plane[2].rows;
    channels = 1;
    size = columns * rows * channels;

    Q->top16.columns = columns;
    Q->top16.rows = rows;
    Q->top16.channels = channels;
    Q->top16.row_stride = columns * channels;
    Q->top16.data = Q->top16.buf =
      (uint16_t *)malloc(sizeof(uint16_t)*size);
  } else {
    uint32_t size = ID->columns * ID->rows * 3;

    TRU->x3rgb16.columns = ID->columns;
    TRU->x3rgb16.rows = ID->rows;
    TRU->x3rgb16.channels = 3;
    TRU->x3rgb16.row_stride = ID->columns * 3;
    TRU->x3rgb16.data = TRU->x3rgb16.buf =
      (uint16_t *)malloc(sizeof(uint16_t)*size);
  }

  true_decode(I, DE);
}

static void x3f_load_huffman_compressed(x3f_info_t *I,
                                        x3f_directory_entry_t *DE,
                                        int bits,
                                        int use_map_table)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_image_data_t *ID = &DEH->data_subsection.image_data;
  x3f_huffman_t *HUF = ID->huffman;
  int table_size = 1<<bits;
  int row_offsets_size = ID->rows * sizeof(HUF->row_offsets.element[0]);

  x3f_printf(DEBUG, "Load huffman compressed\n");

  GET_TABLE(HUF->table, GET4, table_size);

  ID->data_size = read_data_block(&ID->data, I, DE, row_offsets_size);

  GET_TABLE(HUF->row_offsets, GET4, ID->rows);

  x3f_printf(DEBUG, "Make huffman tree ...\n");
  new_huffman_tree(&HUF->tree, bits);
  populate_huffman_tree(&HUF->tree, &HUF->table, &HUF->mapping);
  x3f_printf(DEBUG, "... DONE\n");

#ifdef DBG_PRNT
  print_huffman_tree(HUF->tree.nodes, 0, 0);
#endif

  huffman_decode(I, DE, bits);
}

static void x3f_load_huffman_not_compressed(x3f_info_t *I,
                                            x3f_directory_entry_t *DE,
                                            int bits,
                                            int use_map_table,
                                            int row_stride)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_image_data_t *ID = &DEH->data_subsection.image_data;

  x3f_printf(DEBUG, "Load huffman not compressed\n");

  ID->data_size = read_data_block(&ID->data, I, DE, 0);

  simple_decode(I, DE, bits, row_stride);
}

static void x3f_load_huffman(x3f_info_t *I,
                             x3f_directory_entry_t *DE,
                             int bits,
                             int use_map_table,
                             int row_stride)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_image_data_t *ID = &DEH->data_subsection.image_data;
  x3f_huffman_t *HUF = new_huffman(&ID->huffman);
  uint32_t size;

  if (use_map_table) {
    int table_size = 1<<bits;

    GET_TABLE(HUF->mapping, GET2, table_size);
  }

  switch (ID->type_format) {
  case X3F_IMAGE_RAW_HUFFMAN_X530:
  case X3F_IMAGE_RAW_HUFFMAN_10BIT:
    size = ID->columns * ID->rows * 3;
    HUF->x3rgb16.columns = ID->columns;
    HUF->x3rgb16.rows = ID->rows;
    HUF->x3rgb16.channels = 3;
    HUF->x3rgb16.row_stride = ID->columns * 3;
    HUF->x3rgb16.data = HUF->x3rgb16.buf =
      (uint16_t *)malloc(sizeof(uint16_t)*size);
    break;
  case X3F_IMAGE_THUMB_HUFFMAN:
    size = ID->columns * ID->rows * 3;
    HUF->rgb8.columns = ID->columns;
    HUF->rgb8.columns = ID->rows;
    HUF->rgb8.channels = 3;
    HUF->rgb8.row_stride = ID->columns * 3;
    HUF->rgb8.data = HUF->rgb8.buf =
      (uint8_t *)malloc(sizeof(uint8_t)*size);
    break;
  default:
    /* TODO: Shouldn't this be treated as a fatal error? */
    x3f_printf(ERR, "Unknown huffman image type\n");
  }

  if (row_stride == 0)
    return x3f_load_huffman_compressed(I, DE, bits, use_map_table);
  else
    return x3f_load_huffman_not_compressed(I, DE, bits, use_map_table, row_stride);
}

static void x3f_load_pixmap(x3f_info_t *I, x3f_directory_entry_t *DE)
{
  x3f_printf(DEBUG, "Load pixmap\n");
  x3f_load_image_verbatim(I, DE);
}

static void x3f_load_jpeg(x3f_info_t *I, x3f_directory_entry_t *DE)
{
  x3f_printf(DEBUG, "Load JPEG\n");
  x3f_load_image_verbatim(I, DE);
}

static void x3f_load_image(x3f_info_t *I, x3f_directory_entry_t *DE)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_image_data_t *ID = &DEH->data_subsection.image_data;

  read_data_set_offset(I, DE, X3F_IMAGE_HEADER_SIZE);

  switch (ID->type_format) {
  case X3F_IMAGE_RAW_TRUE:
  case X3F_IMAGE_RAW_MERRILL:
  case X3F_IMAGE_RAW_QUATTRO:
  case X3F_IMAGE_RAW_SDQ:
  case X3F_IMAGE_RAW_SDQH:
    x3f_load_true(I, DE);
    break;
  case X3F_IMAGE_RAW_HUFFMAN_X530:
  case X3F_IMAGE_RAW_HUFFMAN_10BIT:
    x3f_load_huffman(I, DE, 10, 1, ID->row_stride);
    break;
  case X3F_IMAGE_THUMB_PLAIN:
    x3f_load_pixmap(I, DE);
    break;
  case X3F_IMAGE_THUMB_HUFFMAN:
    x3f_load_huffman(I, DE, 8, 0, ID->row_stride);
    break;
  case X3F_IMAGE_THUMB_JPEG:
    x3f_load_jpeg(I, DE);
    break;
  default:
    /* TODO: Shouldn't this be treated as a fatal error? */
    x3f_printf(ERR, "Unknown image type\n");
  }
}

static void x3f_load_camf_decode_type2(x3f_camf_t *CAMF)
{
  uint32_t key = CAMF->t2.crypt_key;
  int i;

  CAMF->decoded_data_size = CAMF->data_size;
  CAMF->decoded_data = malloc(CAMF->decoded_data_size);

  for (i=0; i<CAMF->data_size; i++) {
    uint8_t old, new;
    uint32_t tmp;

    old = ((uint8_t *)CAMF->data)[i];
    key = (key * 1597 + 51749) % 244944;
    tmp = (uint32_t)(key * ((int64_t)301593171) >> 24);
    new = (uint8_t)(old ^ (uint8_t)(((((key << 8) - tmp) >> 1) + tmp) >> 17));
    ((uint8_t *)CAMF->decoded_data)[i] = new;
  }
}


/* NOTE: the unpacking in this code is in big respects identical to
   true_decode_one_color(). The difference is in the output you
   build. It might be possible to make some parts shared. NOTE ALSO:
   This means that the meta data is obfuscated using an image
   compression algorithm. */

static void camf_decode_type4(x3f_camf_t *CAMF)
{
  uint32_t seed = CAMF->t4.decode_bias;
  int row;

  uint8_t *dst;
  uint32_t dst_size = CAMF->t4.decoded_data_size;
  uint8_t *dst_end;

  bool_t odd_dst = 0;

  x3f_hufftree_t *tree = &CAMF->tree;
  bit_state_t BS;

  int32_t row_start_acc[2][2];
  uint32_t rows = CAMF->t4.block_count;
  uint32_t cols = CAMF->t4.block_size;

  CAMF->decoded_data_size = dst_size;

  CAMF->decoded_data = malloc(CAMF->decoded_data_size);
  memset(CAMF->decoded_data, 0, CAMF->decoded_data_size);

  dst = (uint8_t *)CAMF->decoded_data;
  dst_end = dst + dst_size;

  set_bit_state(&BS, CAMF->decoding_start);

  row_start_acc[0][0] = seed;
  row_start_acc[0][1] = seed;
  row_start_acc[1][0] = seed;
  row_start_acc[1][1] = seed;

  for (row = 0; row < rows; row++) {
    int col;
    bool_t odd_row = row&1;
    int32_t acc[2];

    /* We loop through all the columns and the rows. But the actual
       data is smaller than that, so we break the loop when reaching
       the end. */
    for (col = 0; col < cols; col++) {
      bool_t odd_col = col&1;
      int32_t diff = get_true_diff(&BS, tree);
      int32_t prev = col < 2 ?
	row_start_acc[odd_row][odd_col] :
	acc[odd_col];
      int32_t value = prev + diff;

      acc[odd_col] = value;
      if (col < 2)
	row_start_acc[odd_row][odd_col] = value;

      switch(odd_dst) {
      case 0:
	*dst++  = (uint8_t)((value>>4)&0xff);

	if (dst >= dst_end) {
	  goto ready;
	}

	*dst    = (uint8_t)((value<<4)&0xf0);
	break;
      case 1:
	*dst++ |= (uint8_t)((value>>8)&0x0f);

	if (dst >= dst_end) {
	  goto ready;
	}

	*dst++  = (uint8_t)((value<<0)&0xff);

	if (dst >= dst_end) {
	  goto ready;
	}

	break;
      }

      odd_dst = !odd_dst;
    } /* end col */
  } /* end row */

 ready:;
}

static void x3f_load_camf_decode_type4(x3f_camf_t *CAMF)
{
  int i;
  uint8_t *p;
  x3f_true_huffman_element_t *element = NULL;

  for (i=0, p = CAMF->data; *p != 0; i++) {
    /* TODO: Is this too expensive ??*/
    element =
      (x3f_true_huffman_element_t *)realloc(element, (i+1)*sizeof(*element));

    element[i].code_size = *p++;
    element[i].code = *p++;
  }

  CAMF->table.size = i;
  CAMF->table.element = element;

  /* TODO: where does the values 28 and 32 come from? */
#define CAMF_T4_DATA_SIZE_OFFSET 28
#define CAMF_T4_DATA_OFFSET 32
  CAMF->decoding_size = *(uint32_t *)(CAMF->data + CAMF_T4_DATA_SIZE_OFFSET);
  CAMF->decoding_start = (uint8_t *)CAMF->data + CAMF_T4_DATA_OFFSET;

  /* TODO: can it be fewer than 8 bits? Maybe taken from TRU->table? */
  new_huffman_tree(&CAMF->tree, 8);

  populate_true_huffman_tree(&CAMF->tree, &CAMF->table);

#ifdef DBG_PRNT
  print_huffman_tree(CAMF->tree.nodes, 0, 0);
#endif

  camf_decode_type4(CAMF);
}

static void camf_decode_type5(x3f_camf_t *CAMF)
{
  int32_t acc = CAMF->t5.decode_bias;

  uint8_t *dst;

  x3f_hufftree_t *tree = &CAMF->tree;
  bit_state_t BS;

  int32_t i;

  CAMF->decoded_data_size = CAMF->t5.decoded_data_size;
  CAMF->decoded_data = malloc(CAMF->decoded_data_size);

  dst = (uint8_t *)CAMF->decoded_data;

  set_bit_state(&BS, CAMF->decoding_start);

  for (i = 0; i < CAMF->decoded_data_size; i++) {
    int32_t diff = get_true_diff(&BS, tree);

    acc = acc + diff;
    *dst++ = (uint8_t)(acc & 0xff);
  }
}

static void x3f_load_camf_decode_type5(x3f_camf_t *CAMF)
{
  int i;
  uint8_t *p;
  x3f_true_huffman_element_t *element = NULL;

  for (i=0, p = CAMF->data; *p != 0; i++) {
    /* TODO: Is this too expensive ??*/
    element =
      (x3f_true_huffman_element_t *)realloc(element, (i+1)*sizeof(*element));

    element[i].code_size = *p++;
    element[i].code = *p++;
  }

  CAMF->table.size = i;
  CAMF->table.element = element;

  /* TODO: where does the values 28 and 32 come from? */
#define CAMF_T5_DATA_SIZE_OFFSET 28
#define CAMF_T5_DATA_OFFSET 32
  CAMF->decoding_size = *(uint32_t *)(CAMF->data + CAMF_T5_DATA_SIZE_OFFSET);
  CAMF->decoding_start = (uint8_t *)CAMF->data + CAMF_T5_DATA_OFFSET;

  /* TODO: can it be fewer than 8 bits? Maybe taken from TRU->table? */
  new_huffman_tree(&CAMF->tree, 8);

  populate_true_huffman_tree(&CAMF->tree, &CAMF->table);

#ifdef DBG_PRNT
  print_huffman_tree(CAMF->tree.nodes, 0, 0);
#endif

  camf_decode_type5(CAMF);
}

static void x3f_setup_camf_text_entry(camf_entry_t *entry)
{
  entry->text_size = *(uint32_t *)entry->value_address;
  entry->text = entry->value_address + 4;
}

static void x3f_setup_camf_property_entry(camf_entry_t *entry)
{
  int i;
  uint8_t *e =
    entry->entry;
  uint8_t *v =
    entry->value_address;
  uint32_t num =
    entry->property_num = *(uint32_t *)v;
  uint32_t off = *(uint32_t *)(v + 4);

  entry->property_name = (char **)malloc(num*sizeof(uint8_t*));
  entry->property_value = (uint8_t **)malloc(num*sizeof(uint8_t*));

  for (i=0; i<num; i++) {
    uint32_t name_off = off + *(uint32_t *)(v + 8 + 8*i);
    uint32_t value_off = off + *(uint32_t *)(v + 8 + 8*i + 4);

    entry->property_name[i] = (char *)(e + name_off);
    entry->property_value[i] = e + value_off;
  }
}

static void set_matrix_element_info(uint32_t type,
				    uint32_t *size,
				    matrix_type_t *decoded_type)
{
  switch (type) {
  case 0:
    *size = 2;
    *decoded_type = M_INT; /* known to be true */
    break;
  case 1:
    *size = 4;
    *decoded_type = M_UINT; /* TODO: unknown ???? */
    break;
  case 2:
    *size = 4;
    *decoded_type = M_UINT; /* TODO: unknown ???? */
    break;
  case 3:
    *size = 4;
    *decoded_type = M_FLOAT; /* known to be true */
    break;
  case 5:
    *size = 1;
    *decoded_type = M_UINT; /* TODO: unknown ???? */
    break;
  case 6:
    *size = 2;
    *decoded_type = M_UINT; /* TODO: unknown ???? */
    break;
  default:
    x3f_printf(ERR, "Unknown matrix type (%ud)\n", type);
    assert(0);
  }
}

static void get_matrix_copy(camf_entry_t *entry)
{
  uint32_t element_size = entry->matrix_element_size;
  uint32_t elements = entry->matrix_elements;
  int i, size = (entry->matrix_decoded_type==M_FLOAT ?
		 sizeof(double) :
		 sizeof(uint32_t)) * elements;

  entry->matrix_decoded = malloc(size);

  switch (element_size) {
  case 4:
    switch (entry->matrix_decoded_type) {
    case M_INT:
    case M_UINT:
      memcpy(entry->matrix_decoded, entry->matrix_data, size);
      break;
    case M_FLOAT:
      for (i=0; i<elements; i++)
	((double *)entry->matrix_decoded)[i] =
	  (double)((float *)entry->matrix_data)[i];
      break;
    default:
      x3f_printf(ERR, "Invalid matrix element type of size 4\n");
      assert(0);
    }
    break;
  case 2:
    switch (entry->matrix_decoded_type) {
    case M_INT:
      for (i=0; i<elements; i++)
	((int32_t *)entry->matrix_decoded)[i] =
	  (int32_t)((int16_t *)entry->matrix_data)[i];
      break;
    case M_UINT:
      for (i=0; i<elements; i++)
	((uint32_t *)entry->matrix_decoded)[i] =
	  (uint32_t)((uint16_t *)entry->matrix_data)[i];
      break;
    default:
      x3f_printf(ERR, "Invalid matrix element type of size 2\n");
      assert(0);
    }
    break;
  case 1:
    switch (entry->matrix_decoded_type) {
    case M_INT:
      for (i=0; i<elements; i++)
	((int32_t *)entry->matrix_decoded)[i] =
	  (int32_t)((int8_t *)entry->matrix_data)[i];
      break;
    case M_UINT:
      for (i=0; i<elements; i++)
	((uint32_t *)entry->matrix_decoded)[i] =
	  (uint32_t)((uint8_t *)entry->matrix_data)[i];
      break;
    default:
      x3f_printf(ERR, "Invalid matrix element type of size 1\n");
      assert(0);
    }
    break;
  default:
    x3f_printf(ERR, "Unknown size %d\n", element_size);
    assert(0);
  }
}

static void x3f_setup_camf_matrix_entry(camf_entry_t *entry)
{
  int i;
  int totalsize = 1;

  uint8_t *e =
    entry->entry;
  uint8_t *v =
    entry->value_address;
  uint32_t type =
    entry->matrix_type = *(uint32_t *)(v + 0);
  uint32_t dim =
    entry->matrix_dim = *(uint32_t *)(v + 4);
  uint32_t off =
    entry->matrix_data_off = *(uint32_t *)(v + 8);
  camf_dim_entry_t *dentry =
    entry->matrix_dim_entry =
    (camf_dim_entry_t*)malloc(dim*sizeof(camf_dim_entry_t));

  for (i=0; i<dim; i++) {
    uint32_t size =
      dentry[i].size = *(uint32_t *)(v + 12 + 12*i + 0);
    dentry[i].name_offset = *(uint32_t *)(v + 12 + 12*i + 4);
    dentry[i].n = *(uint32_t *)(v + 12 + 12*i + 8);
    dentry[i].name = (char *)(e + dentry[i].name_offset);

    if (dentry[i].n != i) {
      /* TODO: is something needed to be made in this case */
      x3f_printf(DEBUG,
		 "Matrix entry for %s/%s is out of order "
		 "(index/%d != order/%d)\n",
		 entry->name_address, dentry[i].name, dentry[i].n, i);
    }

    totalsize *= size;
  }

  set_matrix_element_info(type,
			  &entry->matrix_element_size,
			  &entry->matrix_decoded_type);
  entry->matrix_data = (void *)(e + off);

  entry->matrix_elements = totalsize;
  entry->matrix_used_space = entry->entry_size - off;

  /* This estimate only works for matrices above a certain size */
  entry->matrix_estimated_element_size = entry->matrix_used_space / totalsize;

  get_matrix_copy(entry);
}

static void x3f_setup_camf_entries(x3f_camf_t *CAMF)
{
  uint8_t *p = (uint8_t *)CAMF->decoded_data;
  uint8_t *end = p + CAMF->decoded_data_size;
  camf_entry_t *entry = NULL;
  int i;

  x3f_printf(DEBUG, "SETUP CAMF ENTRIES\n");

  for (i=0; p < end; i++) {
    uint32_t *p4 = (uint32_t *)p;

    switch (*p4) {
    case X3F_CMbP:
    case X3F_CMbT:
    case X3F_CMbM:
      break;
    default:
      x3f_printf(ERR, "Unknown CAMF entry %x @ %p\n", *p4, p4);
      x3f_printf(ERR, "  start = %p end = %p\n", CAMF->decoded_data, end);
      x3f_printf(ERR, "  left = %ld\n", (long)(end - p));
      x3f_printf(ERR, "Stop parsing CAMF\n");
      /* TODO: Shouldn't this be treated as a fatal error? */
      goto stop;
    }

    /* TODO: lots of realloc - may be inefficient */
    entry = (camf_entry_t *)realloc(entry, (i+1)*sizeof(camf_entry_t));

    /* Pointer */
    entry[i].entry = p;

    /* Header */
    entry[i].id = *p4++;
    entry[i].version = *p4++;
    entry[i].entry_size = *p4++;
    entry[i].name_offset = *p4++;
    entry[i].value_offset = *p4++;

    /* Compute adresses and sizes */
    entry[i].name_address = (char *)(p + entry[i].name_offset);
    entry[i].value_address = p + entry[i].value_offset;
    entry[i].name_size = entry[i].value_offset - entry[i].name_offset;
    entry[i].value_size = entry[i].entry_size - entry[i].value_offset;

    entry[i].text_size = 0;
    entry[i].text = NULL;
    entry[i].property_num = 0;
    entry[i].property_name = NULL;
    entry[i].property_value = NULL;
    entry[i].matrix_type = 0;
    entry[i].matrix_dim = 0;
    entry[i].matrix_data_off = 0;
    entry[i].matrix_data = NULL;
    entry[i].matrix_dim_entry = NULL;

    entry[i].matrix_decoded = NULL;

    switch (entry[i].id) {
    case X3F_CMbP:
      x3f_setup_camf_property_entry(&entry[i]);
      break;
    case X3F_CMbT:
      x3f_setup_camf_text_entry(&entry[i]);
      break;
    case X3F_CMbM:
      x3f_setup_camf_matrix_entry(&entry[i]);
      break;
    }

    p += entry[i].entry_size;
  }

 stop:

  CAMF->entry_table.size = i;
  CAMF->entry_table.element = entry;

  x3f_printf(DEBUG, "SETUP CAMF ENTRIES (READY) Found %d entries\n", i);
}

static void x3f_load_camf(x3f_info_t *I, x3f_directory_entry_t *DE)
{
  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_camf_t *CAMF = &DEH->data_subsection.camf;

  x3f_printf(DEBUG, "Loading CAMF of type %d\n", CAMF->type);

  read_data_set_offset(I, DE, X3F_CAMF_HEADER_SIZE);

  CAMF->data_size = read_data_block(&CAMF->data, I, DE, 0);

  switch (CAMF->type) {
  case 2:			/* Older SD9-SD14 */
    x3f_load_camf_decode_type2(CAMF);
    break;
  case 4:			/* TRUE ... Merrill */
    x3f_load_camf_decode_type4(CAMF);
    break;
  case 5:			/* Quattro ... */
    x3f_load_camf_decode_type5(CAMF);
    break;
  default:
    /* TODO: Shouldn't this be treated as a fatal error? */
    x3f_printf(ERR, "Unknown CAMF type\n");
  }

  if (CAMF->decoded_data != NULL)
    x3f_setup_camf_entries(CAMF);
  else
    /* TODO: Shouldn't this be treated as a fatal error? */
    x3f_printf(ERR, "No decoded CAMF data\n");
}

/* extern */ x3f_return_t x3f_load_data(x3f_t *x3f, x3f_directory_entry_t *DE)
{
  x3f_info_t *I = &x3f->info;

  if (DE == NULL)
    return X3F_ARGUMENT_ERROR;

  switch (DE->header.identifier) {
  case X3F_SECp:
    x3f_load_property_list(I, DE);
    break;
  case X3F_SECi:
    x3f_load_image(I, DE);
    break;
  case X3F_SECc:
    x3f_load_camf(I, DE);
    break;
  default:
    x3f_printf(ERR, "Unknown directory entry type\n");
    return X3F_INTERNAL_ERROR;
  }

  return X3F_OK;
}

/* extern */ x3f_return_t x3f_load_image_block(x3f_t *x3f, x3f_directory_entry_t *DE)
{
  x3f_info_t *I = &x3f->info;

  if (DE == NULL)
    return X3F_ARGUMENT_ERROR;

  x3f_printf(DEBUG, "Load image block\n");

  switch (DE->header.identifier) {
  case X3F_SECi:
    read_data_set_offset(I, DE, X3F_IMAGE_HEADER_SIZE);
    x3f_load_image_verbatim(I, DE);
    break;
  default:
    x3f_printf(ERR, "Unknown image directory entry type\n");
    return X3F_INTERNAL_ERROR;
  }

  return X3F_OK;
}

/* extern */ char *x3f_err(x3f_return_t err)
{
  switch (err) {
  case X3F_OK:             return "ok";
  case X3F_ARGUMENT_ERROR: return "argument error";
  case X3F_INFILE_ERROR:   return "infile error";
  case X3F_OUTFILE_ERROR:  return "outfile error";
  case X3F_INTERNAL_ERROR: return "internal error";
  default:                 return "undefined error";
  }
}

/* --------------------------------------------------------------------- */
/* The End                                                               */
/* --------------------------------------------------------------------- */


<|start_filename|>src/x3f_denoise.h<|end_filename|>
/* X3F_DENOISE.H
 *
 * Library for denoising of X3F image data.
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#ifndef X3F_DENOISE_H
#define X3F_DENOISE_H

#include "x3f_io.h"

#ifdef __cplusplus
extern "C" {
#endif

typedef enum {
  X3F_DENOISE_STD=0,
  X3F_DENOISE_F20=1,
  X3F_DENOISE_F23=2,
} x3f_denoise_type_t;

extern void x3f_denoise(x3f_area16_t *image, x3f_denoise_type_t type);
extern void x3f_expand_quattro(x3f_area16_t *image,
			       x3f_area16_t *active,
			       x3f_area16_t *qtop,
			       x3f_area16_t *expanded,
			       x3f_area16_t *active_exp);

extern void x3f_set_use_opencl(int flag);

#ifdef __cplusplus
}
#endif

#endif


<|start_filename|>src/x3f_print_meta.c<|end_filename|>
/* X3F_PRINT_META.C
 *
 * Library for printing meta data found in X3F files.
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#include "x3f_print_meta.h"
#include "x3f_io.h"

#include <stdio.h>

/* --------------------------------------------------------------------- */
/* Hacky external flags                                                 */
/* --------------------------------------------------------------------- */

/* extern */ uint32_t max_printed_matrix_elements = 100;

/* --------------------------------------------------------------------- */
/* Pretty print the x3f structure                                        */
/* --------------------------------------------------------------------- */

static char x3f_id_buf[5] = {0,0,0,0,0};

static char *x3f_id(uint32_t id)
{
  x3f_id_buf[0] = (id>>0) & 0xff;
  x3f_id_buf[1] = (id>>8) & 0xff;
  x3f_id_buf[2] = (id>>16) & 0xff;
  x3f_id_buf[3] = (id>>24) & 0xff;

  return x3f_id_buf;
}

static char *id_to_str(uint32_t id)
{
  char *idp = (char *)&id;
  static char buf[5];

  buf[0] = *(idp + 0);
  buf[1] = *(idp + 1);
  buf[2] = *(idp + 2);
  buf[3] = *(idp + 3);
  buf[4] = 0;

  return buf;
}

static void print_matrix_element(FILE *f_out, camf_entry_t *entry, uint32_t i)
{
  switch (entry->matrix_decoded_type) {
  case M_FLOAT:
    fprintf(f_out, "%12g ", ((double *)(entry->matrix_decoded))[i]);
    break;
  case M_INT:
    fprintf(f_out, "%12d ", ((int32_t *)(entry->matrix_decoded))[i]);
    break;
  case M_UINT:
    fprintf(f_out, "%12d ", ((uint32_t *)(entry->matrix_decoded))[i]);
    break;
  }
}

static void print_matrix(FILE *f_out, camf_entry_t *entry)
{
  uint32_t dim = entry->matrix_dim;
  uint32_t linesize = entry->matrix_dim_entry[dim-1].size;
  uint32_t blocksize = (uint32_t)(-1);
  uint32_t totalsize = entry->matrix_elements;
  int i;

  switch (entry->matrix_decoded_type) {
  case M_FLOAT:
    fprintf(f_out, "float ");
    break;
  case M_INT:
    fprintf(f_out, "integer ");
    break;
  case M_UINT:
    fprintf(f_out, "unsigned integer ");
    break;
  }

  switch (dim) {
  case 1:
    fprintf(f_out, "[%d]\n", entry->matrix_dim_entry[0].size);
    fprintf(f_out, "x: %s\n", entry->matrix_dim_entry[0].name);
    break;
  case 2:
    fprintf(f_out, "[%d][%d]\n",
	    entry->matrix_dim_entry[0].size,
	    entry->matrix_dim_entry[1].size);
    fprintf(f_out, "x: %s\n", entry->matrix_dim_entry[1].name);
    fprintf(f_out, "y: %s\n", entry->matrix_dim_entry[0].name);
    break;
  case 3:
    fprintf(f_out, "[%d][%d][%d]\n",
	    entry->matrix_dim_entry[0].size,
	    entry->matrix_dim_entry[1].size,
	    entry->matrix_dim_entry[2].size);
    fprintf(f_out, "x: %s\n", entry->matrix_dim_entry[2].name);
    fprintf(f_out, "y: %s\n", entry->matrix_dim_entry[1].name);
    fprintf(f_out, "z: %s (i.e. group)\n", entry->matrix_dim_entry[0].name);
    blocksize = linesize * entry->matrix_dim_entry[dim-2].size;
    break;
  default:
    fprintf(f_out, "\nNot support for higher than 3D in printout\n");
    fprintf(stderr, "Not support for higher than 3D in printout\n");
  }

  for (i=0; i<totalsize; i++) {
    print_matrix_element(f_out, entry, i);
    if ((i+1)%linesize == 0) fprintf(f_out, "\n");
    if ((i+1)%blocksize == 0) fprintf(f_out, "\n");
    if (i >= (max_printed_matrix_elements-1)) {
      fprintf(f_out, "\n... (%d skipped) ...\n", totalsize-i-1);
      break;
    }
  }
}

static void print_file_header_meta_data(FILE *f_out, x3f_t *x3f)
{
  x3f_header_t *H = NULL;

  fprintf(f_out, "BEGIN: file header meta data\n\n");

  H = &x3f->header;
  fprintf(f_out, "header.\n");
  fprintf(f_out, "  identifier        = %08x (%s)\n", H->identifier, x3f_id(H->identifier));
  fprintf(f_out, "  version           = %08x\n", H->version);
  /* TODO: the meaning of the rest of the header for version >= 4.0
           (Quattro) is unknown */
  if (H->version < X3F_VERSION_4_0) {
    fprintf(f_out, "  unique_identifier = %02x...\n", H->unique_identifier[0]);
    fprintf(f_out, "  mark_bits         = %08x\n", H->mark_bits);
    fprintf(f_out, "  columns           = %08x (%d)\n", H->columns, H->columns);
    fprintf(f_out, "  rows              = %08x (%d)\n", H->rows, H->rows);
    fprintf(f_out, "  rotation          = %08x (%d)\n", H->rotation, H->rotation);
    if (x3f->header.version >= X3F_VERSION_2_1) {
      int num_ext_data =
	H->version >= X3F_VERSION_3_0 ? NUM_EXT_DATA_3_0 : NUM_EXT_DATA_2_1;
      int i;

      fprintf(f_out, "  white_balance     = %s\n", H->white_balance);
      if (x3f->header.version >= X3F_VERSION_2_3)
	fprintf(f_out, "  color_mode        = %s\n", H->color_mode);

      fprintf(f_out, "  extended_types\n");
      for (i=0; i<num_ext_data; i++) {
	uint8_t type = H->extended_types[i];
	float data = H->extended_data[i];

	fprintf(f_out, "    %2d: %3d = %9f\n", i, type, data);
      }
    }
  }

  fprintf(f_out, "END: file header meta data\n\n");
}

static void print_camf_meta_data2(FILE *f_out, x3f_camf_t *CAMF)
{
  fprintf(f_out, "BEGIN: CAMF meta data\n\n");

  if (CAMF->entry_table.size != 0) {
    camf_entry_t *entry = CAMF->entry_table.element;
    int i;

    for (i=0; i<CAMF->entry_table.size; i++) {
      if (f_out == stdout) {
	fprintf(f_out, "          element[%d].name = \"%s\"\n",
		i, entry[i].name_address);
	fprintf(f_out, "            id = %x (%s)\n",
		entry[i].id, id_to_str(entry[i].id));
	fprintf(f_out, "            entry_size = %d\n",
		entry[i].entry_size);
	fprintf(f_out, "            name_size = %d\n",
		entry[i].name_size);
	fprintf(f_out, "            value_size = %d\n",
		entry[i].value_size);
      }

      /* Text CAMF */
      if (entry[i].text_size != 0) {
	fprintf(f_out, "BEGIN: CAMF text meta data (%s)\n",
		entry[i].name_address);
	fprintf(f_out, "\"%s\"\n", entry[i].text);
	fprintf(f_out, "END: CAMF text meta data\n\n");
      }

      /* Property CAMF */
      if (entry[i].property_num != 0) {
	int j;
	fprintf(f_out, "BEGIN: CAMF property meta data (%s)\n",
		entry[i].name_address);

	for (j=0; j<entry[i].property_num; j++) {
	  fprintf(f_out, "              \"%s\" = \"%s\"\n",
		  entry[i].property_name[j],
		  entry[i].property_value[j]);
	}

	fprintf(f_out, "END: CAMF property meta data\n\n");
      }

      /* Matrix CAMF */
      if (entry[i].matrix_dim != 0) {
	int j;
	camf_dim_entry_t *dentry = entry[i].matrix_dim_entry;

	fprintf(f_out, "BEGIN: CAMF matrix meta data (%s)\n",
		entry[i].name_address);

	if (f_out == stdout) {
	  fprintf(f_out, "            matrix_type = %d\n", entry[i].matrix_type);
	  fprintf(f_out, "            matrix_dim = %d\n", entry[i].matrix_dim);
	  fprintf(f_out, "            matrix_data_off = %d\n", entry[i].matrix_data_off);

	  for (j=0; j<entry[i].matrix_dim; j++) {
	    fprintf(f_out, "            %d\n", j);
	    fprintf(f_out, "              size = %d\n", dentry[j].size);
	    fprintf(f_out, "              name_offset = %d\n", dentry[j].name_offset);
	    fprintf(f_out, "              n = %d%s\n", dentry[j].n, j==dentry[j].n ? "" : " (out of order)");
	    fprintf(f_out, "              name = \"%s\"\n", dentry[j].name);
	  }

	  fprintf(f_out, "            matrix_element_size = %d\n", entry[i].matrix_element_size);
	  fprintf(f_out, "            matrix_elements = %d\n", entry[i].matrix_elements);
	  fprintf(f_out, "            matrix_estimated_element_size = %g\n", entry[i].matrix_estimated_element_size);
	}

	print_matrix(f_out, &entry[i]);

	fprintf(f_out, "END: CAMF matrix meta data\n\n");
      }
    }
  }

  fprintf(f_out, "END: CAMF meta data\n\n");
}

static void print_camf_meta_data(FILE *f_out, x3f_t *x3f)
{
  x3f_directory_entry_t *DE = x3f_get_camf(x3f);

  if (DE == NULL) {
    fprintf(f_out, "INFO: No CAMF meta data found\n\n");
    return;
  }

  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_camf_t *CAMF = &DEH->data_subsection.camf;

  print_camf_meta_data2(f_out, CAMF);
}

static void print_prop_meta_data2(FILE *f_out, x3f_property_list_t *PL)
{
  fprintf(f_out, "BEGIN: PROP meta data\n\n");

  if (PL->property_table.size != 0) {
    int i;
    x3f_property_t *P = PL->property_table.element;

    for (i=0; i<PL->num_properties; i++)
      fprintf(f_out, "          [%d] \"%s\" = \"%s\"\n",
	      i, P[i].name_utf8, P[i].value_utf8);
  }

  fprintf(f_out, "END: PROP meta data\n\n");
}

static void print_prop_meta_data(FILE *f_out, x3f_t *x3f)
{
  x3f_directory_entry_t *DE = x3f_get_prop(x3f);

  if (DE == NULL) {
    fprintf(f_out, "INFO: No PROP meta data found\n\n");
    return;
  }

  x3f_directory_entry_header_t *DEH = &DE->header;
  x3f_property_list_t *PL = &DEH->data_subsection.property_list;

  print_prop_meta_data2(f_out, PL);
}

/* extern */ void x3f_print_meta(x3f_t *x3f)
{
  int d;
  x3f_directory_section_t *DS = NULL;
  x3f_info_t *I = NULL;

  if (x3f == NULL) {
    printf("Null x3f\n");
    return;
  }

  I = &x3f->info;
  printf("info.\n");
  printf("  error = %s\n", I->error);
  printf("  input.\n");
  printf("    file = %p\n", I->input.file);
  printf("  output.\n");
  printf("    file = %p\n", I->output.file);

  print_file_header_meta_data(stdout, x3f);

  DS = &x3f->directory_section;
  printf("directory_section.\n");
  printf("  identifier            = %08x (%s)\n", DS->identifier, x3f_id(DS->identifier));
  printf("  version               = %08x\n", DS->version);
  printf("  num_directory_entries = %08x\n", DS->num_directory_entries);
  printf("  directory_entry       = %p\n", DS->directory_entry);

  for (d=0; d<DS->num_directory_entries; d++) {
    x3f_directory_entry_t *DE = &DS->directory_entry[d];
    x3f_directory_entry_header_t *DEH = &DE->header;

    printf("  directory_entry.\n");
    printf("    input.\n");
    printf("      offset = %08x\n", DE->input.offset);
    printf("      size   = %08x\n", DE->input.size);
    printf("    output.\n");
    printf("      offset = %08x\n", DE->output.offset);
    printf("      size   = %08x\n", DE->output.size);
    printf("    type     = %08x (%s)\n", DE->type, x3f_id(DE->type));
    DEH = &DE->header;
    printf("    header.\n");
    printf("      identifier = %08x (%s)\n", DEH->identifier, x3f_id(DEH->identifier));
    printf("      version    = %08x\n", DEH->version);

    /* NOTE - the tests below could be made on DE->type instead */

    if (DEH->identifier == X3F_SECp) {
      x3f_property_list_t *PL = &DEH->data_subsection.property_list;
      printf("      data_subsection.property_list.\n");
      printf("        num_properties   = %08x\n", PL->num_properties);
      printf("        character_format = %08x\n", PL->character_format);
      printf("        reserved         = %08x\n", PL->reserved);
      printf("        total_length     = %08x\n", PL->total_length);
      printf("    property_table       = %x %p\n",
	      PL->property_table.size, PL->property_table.element);
      printf("        data             = %p\n", PL->data);
      printf("        data_size        = %x\n", PL->data_size);

      print_prop_meta_data2(stdout, PL);
    }

    if (DEH->identifier == X3F_SECi) {
      x3f_image_data_t *ID = &DEH->data_subsection.image_data;
      x3f_huffman_t *HUF = ID->huffman;
      x3f_true_t *TRU = ID->tru;
      x3f_quattro_t *Q = ID->quattro;

      printf("      data_subsection.image_data.\n");
      printf("        type        = %08x\n", ID->type);
      printf("        format      = %08x\n", ID->format);
      printf("        type_format = %08x\n", ID->type_format);
      printf("        columns     = %08x (%d)\n",
             ID->columns, ID->columns);
      printf("        rows        = %08x (%d)\n",
             ID->rows, ID->rows);
      printf("        row_stride  = %08x (%d)\n",
             ID->row_stride, ID->row_stride);

      if (HUF == NULL) {
        printf("        huffman     = %p\n", HUF);
      } else {
        printf("        huffman->\n");
        printf("          mapping     = %x %p\n",
               HUF->mapping.size, HUF->mapping.element);
        printf("          table       = %x %p\n",
               HUF->table.size, HUF->table.element);
        printf("          tree        = %d %p\n",
	       HUF->tree.free_node_index, HUF->tree.nodes);
        printf("          row_offsets = %x %p\n",
               HUF->row_offsets.size, HUF->row_offsets.element);
        printf("          rgb8        = %x %x %p\n",
               HUF->rgb8.columns, HUF->rgb8.rows, HUF->rgb8.data);
        printf("          x3rgb16     = %x %x %p\n",
               HUF->x3rgb16.columns, HUF->x3rgb16.rows, HUF->x3rgb16.data);
      }

      if (TRU == NULL) {
        printf("        tru         = %p\n", TRU);
      } else {
	int i;

        printf("        tru->\n");
        printf("          seed[0]     = %x\n", TRU->seed[0]);
        printf("          seed[1]     = %x\n", TRU->seed[1]);
        printf("          seed[2]     = %x\n", TRU->seed[2]);
        printf("          unknown     = %x\n", TRU->unknown);
        printf("          table       = %x %p\n",
               TRU->table.size, TRU->table.element);
        printf("          plane_size  = %x %p (",
               TRU->plane_size.size, TRU->plane_size.element);
	for (i=0; i<TRU->plane_size.size; i++)
	  printf(" %d", TRU->plane_size.element[i]);
	printf(" )\n");
	printf("          plane_address (");
	for (i=0; i<TRUE_PLANES; i++)
	  printf(" %p", TRU->plane_address[i]);
	printf(" )\n");
	printf("          tree        = %d %p\n",
	       TRU->tree.free_node_index, TRU->tree.nodes);
        printf("          x3rgb16     = %x %x %p\n",
               TRU->x3rgb16.columns, TRU->x3rgb16.rows, TRU->x3rgb16.data);
      }

      if (Q == NULL) {
	printf("        quattro     = %p\n", Q);
      } else {
	int i;

	printf("        quattro->\n");
	printf("          planes (");
	for (i=0; i<TRUE_PLANES; i++) {
	  printf(" %d", Q->plane[i].columns);
	  printf("x%d", Q->plane[i].rows);
	}
	printf(" )\n");
        printf("          unknown     = %x %d\n",
               Q->unknown, Q->unknown);
      }

      printf("        data        = %p\n", ID->data);
    }

    if (DEH->identifier == X3F_SECc) {
      x3f_camf_t *CAMF = &DEH->data_subsection.camf;
      printf("      data_subsection.camf.\n");
      printf("        type             = %x\n", CAMF->type);

      switch (CAMF->type) {
      case 2:
        printf("        type2\n");
        printf("          reserved         = %x\n", CAMF->t2.reserved);
        printf("          infotype         = %08x (%s)\n", CAMF->t2.infotype, x3f_id(CAMF->t2.infotype));
        printf("          infotype_version = %08x\n", CAMF->t2.infotype_version);
        printf("          crypt_key        = %08x\n", CAMF->t2.crypt_key);
        break;
      case 4:
        printf("        type4\n");
        printf("          decoded_data_size= %x\n", CAMF->t4.decoded_data_size);
        printf("          decode_bias      = %x\n", CAMF->t4.decode_bias);
        printf("          block_size       = %x\n", CAMF->t4.block_size);
        printf("          block_count      = %x\n", CAMF->t4.block_count);
        break;
      case 5:
        printf("        type5\n");
        printf("          decoded_data_size= %x\n", CAMF->t5.decoded_data_size);
        printf("          decode_bias      = %x\n", CAMF->t5.decode_bias);
        printf("          unknown2         = %x\n", CAMF->t5.unknown2);
        printf("          unknown3         = %x\n", CAMF->t5.unknown3);
        break;
      default:
        printf("       (Unknown CAMF type)\n");
        printf("          val0             = %x\n", CAMF->tN.val0);
        printf("          val1             = %x\n", CAMF->tN.val1);
        printf("          val2             = %x\n", CAMF->tN.val2);
        printf("          val3             = %x\n", CAMF->tN.val3);
      }

      printf("        data             = %p\n", CAMF->data);
      printf("        data_size        = %x\n", CAMF->data_size);
      printf("        decoding_start   = %p\n", CAMF->decoding_start);
      printf("        decoding_size    = %x\n", CAMF->decoding_size);
      printf("        table            = %x %p\n",
	     CAMF->table.size, CAMF->table.element);
      printf("          tree           = %d %p\n",
	     CAMF->tree.free_node_index, CAMF->tree.nodes);
      printf("        decoded_data     = %p\n", CAMF->decoded_data);
      printf("        decoded_data_size= %x\n", CAMF->decoded_data_size);
      printf("        entry_table      = %x %p\n",
	     CAMF->entry_table.size, CAMF->entry_table.element);

      print_camf_meta_data2(stdout, CAMF);
    }
  }
}

/* extern */ x3f_return_t x3f_dump_meta_data(x3f_t *x3f, char *outfilename)
{
  FILE *f_out = fopen(outfilename, "wb");

  if (f_out == NULL) {
    return X3F_OUTFILE_ERROR;
  }

  print_file_header_meta_data(f_out, x3f);

  print_camf_meta_data(f_out, x3f);

  print_prop_meta_data(f_out, x3f);

  /* We assume that the JPEG meta data is not needed. Therefore JPEG
     is not loaded and no call to any EXIF extraction tool either
     called. */

  fclose(f_out);

  return X3F_OK;
}


<|start_filename|>src/x3f_output_ppm.h<|end_filename|>
/* X3F_OUTPUT_PPM.H
 *
 * Library for writing the image as PPM.
 *
 * Copyright 2015 - <NAME> <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#ifndef X3F_OUTPUT_PPM_H
#define X3F_OUTPUT_PPM_H

#include "x3f_io.h"
#include "x3f_process.h"

extern x3f_return_t x3f_dump_raw_data_as_ppm(x3f_t *x3f, char *outfilename,
                                             x3f_color_encoding_t encoding,
					     int crop,
					     int fix_bad,
					     int denoise,
					     int apply_sgain,
					     char *wb,
                                             int binary);

#endif


<|start_filename|>x86_64-w64-mingw32.cmake<|end_filename|>
SET(CMAKE_SYSTEM_NAME Windows-GNU)
SET(CMAKE_C_COMPILER x86_64-w64-mingw32-gcc)
SET(CMAKE_CXX_COMPILER x86_64-w64-mingw32-g++)
SET(CMAKE_RC_COMPILER x86_64-w64-mingw32-windres)
SET(CMAKE_FIND_ROOT_PATH /usr/x86_64-w64-mingw32)
SET(CMAKE_FIND_ROOT_PATH_MODE_PROGRAM NEVER)
SET(CMAKE_FIND_ROOT_PATH_MODE_LIBRARY ONLY)
SET(CMAKE_FIND_ROOT_PATH_MODE_INCLUDE ONLY)


<|start_filename|>src/x3f_version.c<|end_filename|>
/* X3F_VERSION.C
 *
 * Contains the version of X3F Tools
 *
 * Copyright 2016 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

char *version = VERSION;


<|start_filename|>src/x3f_matrix.h<|end_filename|>
/* X3F_MATRIX.H
 *
 * Library for manipulating matrices.
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#ifndef X3F_MATRIX_H
#define X3F_MATRIX_H

#include "x3f_printf.h"

#include <inttypes.h>

extern void x3f_3x1_invert(double *a, double *ainv);
extern void x3f_3x1_comp_mul(double *a, double *b, double *c);

extern void x3f_scalar_3x1_mul(double a, double *b, double *c);
extern void x3f_scalar_3x3_mul(double a, double *b, double *c);

extern void x3f_3x3_3x1_mul(double *a, double *b, double *c);
extern void x3f_3x3_3x3_mul(double *a, double *b, double *c);

extern void x3f_3x3_inverse(double *a, double *ainv);

extern void x3f_3x3_identity(double *a);
extern void x3f_3x3_diag(double *a, double *b);
extern void x3f_3x3_ones(double *a);

extern void x3f_3x1_print(x3f_verbosity_t level, double *a);
extern void x3f_3x3_print(x3f_verbosity_t level, double *a);

extern void x3f_XYZ_to_ProPhotoRGB(double *a);
extern void x3f_ProPhotoRGB_to_XYZ(double *a);

extern void x3f_XYZ_to_AdobeRGB(double *a);
extern void x3f_AdobeRGB_to_XYZ(double *a);

extern void x3f_XYZ_to_sRGB(double *a);
extern void x3f_sRGB_to_XYZ(double *a);

extern void x3f_CIERGB_to_XYZ(double *a);

extern void x3f_Bradford_D50_to_D65(double *a);
extern void x3f_Bradford_D65_to_D50(double *a);

extern void x3f_sRGB_LUT(double *lut, int size, uint16_t max);
extern void x3f_gamma_LUT(double *lut, int size, uint16_t max, double gamma);
extern uint16_t x3f_LUT_lookup(double *lut, int size, double val);

#endif	/* X3F_MATRIX_H */


<|start_filename|>src/x3f_dump.c<|end_filename|>
/* X3F_DUMP.C
 *
 * Library for writing unprocessed RAW and JPEG data.
 *
 * Copyright 2015 - <NAME> <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#include "x3f_dump.h"
#include "x3f_io.h"

#include <stdio.h>

/* extern */ x3f_return_t x3f_dump_raw_data(x3f_t *x3f,
                                            char *outfilename)
{
  x3f_directory_entry_t *DE = x3f_get_raw(x3f);

  if (DE == NULL) {
    return X3F_ARGUMENT_ERROR;
  } else {
    x3f_directory_entry_header_t *DEH = &DE->header;
    x3f_image_data_t *ID = &DEH->data_subsection.image_data;
    void *data = ID->data;

    if (data == NULL) {
      return X3F_INTERNAL_ERROR;
    } else {
      FILE *f_out = fopen(outfilename, "wb");

      if (f_out == NULL) {
        return X3F_OUTFILE_ERROR;
      } else {
        fwrite(data, 1, DE->input.size, f_out);
        fclose(f_out);
      }
    }
  }

  return X3F_OK;
}

/* extern */ x3f_return_t x3f_dump_jpeg(x3f_t *x3f, char *outfilename)
{
  x3f_directory_entry_t *DE = x3f_get_thumb_jpeg(x3f);

  if (DE == NULL) {
    return X3F_ARGUMENT_ERROR;
  } else {
    x3f_directory_entry_header_t *DEH = &DE->header;
    x3f_image_data_t *ID = &DEH->data_subsection.image_data;
    void *data = ID->data;

    if (data == NULL) {
      return X3F_INTERNAL_ERROR;
    } else {
      FILE *f_out = fopen(outfilename, "wb");

      if (f_out == NULL) {
        return X3F_OUTFILE_ERROR;
      } else {
        fwrite(data, 1, DE->input.size, f_out);
        fclose(f_out);
      }
    }
  }

  return X3F_OK;
}


<|start_filename|>x86_64-apple-darwin11.cmake<|end_filename|>
SET(CMAKE_SYSTEM_NAME Darwin-GNU)
SET(CMAKE_C_COMPILER x86_64-apple-darwin11-gcc)
SET(CMAKE_CXX_COMPILER x86_64-apple-darwin11-g++)
SET(_CMAKE_TOOLCHAIN_PREFIX x86_64-apple-darwin11-)
EXECUTE_PROCESS(COMMAND x86_64-apple-darwin11-gcc -print-sysroot OUTPUT_VARIABLE CMAKE_FIND_ROOT_PATH OUTPUT_STRIP_TRAILING_WHITESPACE)
SET(CMAKE_FIND_ROOT_PATH_MODE_PROGRAM NEVER)
SET(CMAKE_FIND_ROOT_PATH_MODE_LIBRARY ONLY)
SET(CMAKE_FIND_ROOT_PATH_MODE_INCLUDE ONLY)


<|start_filename|>src/x3f_denoise_utils.cpp<|end_filename|>
/* X3F_DENOISE_UTILS.CPP
 *
 * Library for support functions for noise reduction algorithms,
 * including color conversions and conversions to/from float type images
 *
 * Copyright 2015 - <NAME> and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#include <iostream>
#include <inttypes.h>

#include "x3f_denoise_utils.h"

using namespace cv;

static const int32_t O_UV = 32768; // To avoid clipping negative values in U,V

/* This file contains support functions for matrix transformations,
   etc, used in denoising images.  It is not the denoising functions themselves.
*/

// Matrix used to convert BMT to YUV:
//  0    0    1
//  2    0   -2
//  1   -2    1
void BMT_to_YUV_YisT(x3f_area16_t *image)
{
  for (uint32_t row=0; row < image->rows; row++)
    for (uint32_t col=0; col < image->columns; col++) {
      uint16_t *p = &image->data[row*image->row_stride + col*image->channels];

      int32_t B = (int32_t)p[0];
      int32_t M = (int32_t)p[1];
      int32_t T = (int32_t)p[2];

      int32_t Y =               +T;
      int32_t U =   +2*B      -2*T;
      int32_t V =     +B -2*M   +T;

      p[0] = saturate_cast<uint16_t>(Y);
      p[1] = saturate_cast<uint16_t>(U + O_UV);
      p[2] = saturate_cast<uint16_t>(V + O_UV);
    }
}

// Matrix used to convert YUV to BMT:
//  1    1/2  0
//  1    1/4 -1/2
//  1    0    0
void YUV_to_BMT_YisT(x3f_area16_t *image)
{
  for (uint32_t row=0; row < image->rows; row++)
    for (uint32_t col=0; col < image->columns; col++) {
      uint16_t *p = &image->data[row*image->row_stride + col*image->channels];

      int32_t Y = (int32_t)p[0];
      int32_t U = (int32_t)p[1] - O_UV;
      int32_t V = (int32_t)p[2] - O_UV;

      int32_t B = ( +2*Y   +U      ) / 2;
      int32_t M = ( +4*Y   +U -2*V ) / 4;
      int32_t T = (   +Y           )    ;

      p[0] = saturate_cast<uint16_t>(B);
      p[1] = saturate_cast<uint16_t>(M);
      p[2] = saturate_cast<uint16_t>(T);
    }
}

// Matrix used to convert BMT to YUV:
// 0 0 4
// 2 0 -2
// 1 -2 1
void BMT_to_YUV_Yis4T(x3f_area16_t *image)
{
  for (uint32_t row=0; row < image->rows; row++)
    for (uint32_t col=0; col < image->columns; col++) {
      uint16_t *p = &image->data[row*image->row_stride + col*image->channels];
      int32_t B = (int32_t)p[0];
      int32_t M = (int32_t)p[1];
      int32_t T = (int32_t)p[2];
      
      int32_t Y = +4*T;
      int32_t U = +2*B -2*T;
      int32_t V = +B -2*M +T;
      
      p[0] = saturate_cast<uint16_t>(Y);
      p[1] = saturate_cast<uint16_t>(U + O_UV);
      p[2] = saturate_cast<uint16_t>(V + O_UV);
    }
}

// Matrix used to convert YUV to BMT:
// 1/4 1/2 0
// 1/4 1/4 -1/2
// 1/4 0 0
void YUV_to_BMT_Yis4T(x3f_area16_t *image)
{
  for (uint32_t row=0; row < image->rows; row++)
    for (uint32_t col=0; col < image->columns; col++) {
      uint16_t *p = &image->data[row*image->row_stride + col*image->channels];
      int32_t Y = (int32_t)p[0];
      int32_t U = (int32_t)p[1] - O_UV;
      int32_t V = (int32_t)p[2] - O_UV;
      
      int32_t B = ( +Y +2*U + 2 ) / 4;
      int32_t M = ( +Y +U -2*V + 2 ) / 4;
      int32_t T = ( +Y + 2 ) / 4;
      
      p[0] = saturate_cast<uint16_t>(B);
      p[1] = saturate_cast<uint16_t>(M);
      p[2] = saturate_cast<uint16_t>(T);
    }
}

// Matrix used to convert BMT to YUV:
//  1/3  1/3  1/3
//  2    0   -2
//  1   -2    1
void BMT_to_YUV_STD(x3f_area16_t *image)
{
  for (uint32_t row=0; row < image->rows; row++)
    for (uint32_t col=0; col < image->columns; col++) {
      uint16_t *p = &image->data[row*image->row_stride + col*image->channels];

      int32_t B = (int32_t)p[0];
      int32_t M = (int32_t)p[1];
      int32_t T = (int32_t)p[2];

      int32_t Y = (   +B   +M   +T ) / 3;
      int32_t U =   +2*B      -2*T;
      int32_t V =     +B -2*M   +T;

      p[0] = saturate_cast<uint16_t>(Y);
      p[1] = saturate_cast<uint16_t>(U + O_UV);
      p[2] = saturate_cast<uint16_t>(V + O_UV);
    }
}

// Matrix used to convert YUV to BMT:
//  1    1/4  1/6
//  1    0   -1/3
//  1   -1/4  1/6
void YUV_to_BMT_STD(x3f_area16_t *image)
{
  for (uint32_t row=0; row < image->rows; row++)
    for (uint32_t col=0; col < image->columns; col++) {
      uint16_t *p = &image->data[row*image->row_stride + col*image->channels];

      int32_t Y = (int32_t)p[0];
      int32_t U = (int32_t)p[1] - O_UV;
      int32_t V = (int32_t)p[2] - O_UV;

      int32_t B = ( +12*Y +3*U +2*V ) / 12;
      int32_t M = (  +3*Y        -V ) /  3;
      int32_t T = ( +12*Y -3*U +2*V ) / 12;

      p[0] = saturate_cast<uint16_t>(B);
      p[1] = saturate_cast<uint16_t>(M);
      p[2] = saturate_cast<uint16_t>(T);
    }
}

float* convert_to_float_image(x3f_area16_t *image)
{
  uint32_t x, y, i;
  const uint32_t num_rows = image->rows;
  const uint32_t num_cols = image->columns;
  const uint32_t jump = image->channels;
  float* outImage = new float[num_cols * num_rows * jump];
  float* out_ptr = outImage;
  uint16_t* in_ptr;
  
  for (y = 0; y < num_rows; y++)
    {
      for (x = 0, in_ptr = (image->data + y * image->row_stride); x < num_cols; x++)
	{
	  for (i = 0; i < jump; i++, out_ptr++, in_ptr++)
	    {
	      *out_ptr = (float(*in_ptr)/65535.0f);
	    }
	}
    }
  return outImage;
}

void convert_from_float_image(x3f_area16_t *image, float* in_float_image)
{
  uint32_t x, y, i;
  const uint32_t num_rows = image->rows;
  const uint32_t num_cols = image->columns;
  const uint32_t jump = image->channels;
  float* in_float_ptr = in_float_image;
  uint16_t* in_ptr;
  float new_val;
  for (y = 0; y < num_rows; y++)
    {
      for (x = 0, in_ptr = (image->data + y * image->row_stride); x < num_cols; x++)
	{
	  //skipping the first channel, leaving that behind
	  in_ptr++;
	  in_float_ptr++;
	  for (i = 1; i < jump; i++, in_ptr++, in_float_ptr++)
	    {
	      new_val = *in_float_ptr;
	      if (new_val < 0) new_val = 0;
	      if (new_val > 1.0f) new_val = 1.0f; //clamping
	      *in_ptr = new_val * 65535.0f;
	    }
	}
    }
}


<|start_filename|>src/x3f_denoise_aniso.h<|end_filename|>
/* X3F_DENOISE_ANISO.H
 *
 * Library for anisotropic denoising of X3F image data.
 *
 * Copyright 2015 - <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#ifndef DENOISE_ANISO
#define DENOISE_ANISO

#include "x3f_denoise_utils.h"  //has the C++ error define in it

//but adding it, just in case it's not entirely clear
#ifndef __cplusplus
#error This file can only be included from C++
#endif

void median_filter(x3f_area16_t *image);
void denoise_aniso(x3f_area16_t *image, const int& in_iterations);
void denoise_iso(x3f_area16_t *image, const int& in_iterations);
void denoise_splotchify(x3f_area16_t *image, const int& in_radius);

#endif //DENOISE_ANISO


<|start_filename|>src/x3f_spatial_gain.h<|end_filename|>
/* X3F_SPATIAL_GAIN.H
 *
 * Library for adjusting for spatial gain in X3F images.
 *
 * Copyright 2015 - <NAME> <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#ifndef X3F_SPATIAL_GAIN_H
#define X3F_SPATIAL_GAIN_H

#include "x3f_io.h"

typedef struct {
  double weight;

  uint32_t *gain;
  double mingain, delta;
} x3f_spatial_gain_corr_merrill_t;

typedef struct {
  double *gain;		    /* final gain (interpolated if necessary) */
  int malloc;		    /* 1 if gain is allocated with malloc() */

  int rows, cols;
  int rowoff, coloff, rowpitch, colpitch;
  int chan, channels;

  x3f_spatial_gain_corr_merrill_t mgain[4]; /* raw Merrill-type gains */
  int mgain_num;
} x3f_spatial_gain_corr_t;

#define MAXCORR 6 /* Quattro HP: R, G, B0, B1, B2, B3 */

extern int x3f_get_merrill_type_spatial_gain(x3f_t *x3f, int hp_flag,
					     x3f_spatial_gain_corr_t *corr);
extern int x3f_get_interp_merrill_type_spatial_gain(x3f_t *x3f, int hp_flag,
						    x3f_spatial_gain_corr_t *corr);
extern int x3f_get_classic_spatial_gain(x3f_t *x3f, char *wb,
					x3f_spatial_gain_corr_t *corr);
extern int x3f_get_spatial_gain(x3f_t *x3f, char *wb,
				x3f_spatial_gain_corr_t *corr);
extern void x3f_cleanup_spatial_gain(x3f_spatial_gain_corr_t *corr,
				     int corr_num);
extern double x3f_calc_spatial_gain(x3f_spatial_gain_corr_t *corr,
				    int corr_num,
				    int row, int col, int chan,
				    int rows, int cols);

#endif


<|start_filename|>src/x3f_process.c<|end_filename|>
/* X3F_PROCESS.C
 *
 * Library for processing X3F data.
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#include "x3f_io.h"
#include "x3f_process.h"
#include "x3f_meta.h"
#include "x3f_image.h"
#include "x3f_matrix.h"
#include "x3f_denoise.h"
#include "x3f_spatial_gain.h"
#include "x3f_printf.h"

#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <assert.h>

static int sum_area(x3f_area16_t area, int colors,
		    uint64_t *sum)
{
  int row, col, color;

  for (color=0; color<colors; color++) sum[color] = 0;

  for (row = 0; row < area.rows; row++)
    for (col = 0; col < area.columns; col++)
      for (color = 0; color < colors; color++)
	sum[color] += (uint64_t)area.data[area.row_stride*row +
					  area.channels*col + color];

  return area.columns*area.rows;
}

static int sum_area_sqdev(x3f_area16_t area, int colors, double *mean,
			  double *sum)
{
  int row, col, color;

  for (color=0; color<colors; color++) sum[color] = 0.0;

  for (row = 0; row < area.rows; row++)
    for (col = 0; col < area.columns; col++)
      for (color = 0; color < colors; color++) {
	double dev = area.data[area.row_stride*row +
			       area.channels*col + color] - mean[color];
	sum[color] += dev*dev;
      }

  return area.columns*area.rows;
}

static int get_black_level(x3f_t *x3f,
			   x3f_area16_t *image, int rescale, int colors,
			   double *black_level, double *black_dev)
{
  uint64_t *black, *black_sum;
  double *black_sqdev, *black_sqdev_sum;
  int pixels_sum, i;

  col_side_t side[4] = {COL_SIDE_WRONG,
			COL_SIDE_WRONG,
			COL_SIDE_LEFT,
			COL_SIDE_RIGHT};
  char *name[4] = {"DarkShieldTop",
		   "DarkShieldBottom",
		   "Left", /* Only used in printout */
		   "Right" /* Only used in printout */
  };
  int use[4] = {1, 1, 1, 1};
  x3f_area16_t area[4];

  if (image->channels < colors) return 0;

#define BOTTOM 1

  /* Workaround for bug in DP2 firmware. DarkShieldBottom is specified
     incorrectly and thus ignored. */
  {
    char *cammodel;

    if (x3f_get_prop_entry(x3f, "CAMMODEL", &cammodel))
      if (!strcmp(cammodel, "SIGMA DP2"))
	use[BOTTOM] = 0;
  }

  /* Workaround for bug in sd Quattro H firmaware. DarkShieldBottom is
     specified incorrectly and thus ignored. */
  {
    uint32_t cameraid;

    if (x3f_get_camf_unsigned(x3f, "CAMERAID", &cameraid))
      if (cameraid == X3F_CAMERAID_SDQH)
	use[BOTTOM] = 0;
  }

  /* Real CAMF rects */
  for (i=0; i<2; i++)
    if (use[i])
      use[i] = x3f_crop_area_camf(x3f, name[i], image, rescale, &area[i]);

  /* Column based rects */
  for (i=2; i<4; i++)
    if (use[i])
      use[i] = x3f_crop_area_column(x3f, side[i], image, rescale, &area[i]);

  for (i=0; i<4; i++)
    if (use[i])
      x3f_printf(DEBUG, "Calculate black level for %s\n", name[i]);
    else
      x3f_printf(DEBUG, "Do not calculate black level for %s\n", name[i]);

  pixels_sum = 0;
  black = alloca(colors*sizeof(uint64_t));
  black_sum = alloca(colors*sizeof(uint64_t));
  for (i=0; i<colors; i++) black_sum[i] = 0;

  x3f_printf(DEBUG, "Dark level\n");

  for (i=0; i<4; i++)
    if (use[i]) {
      int color;
      int pixels = sum_area(area[i], colors, black);

      pixels_sum += pixels;

      x3f_printf(DEBUG, "  %s (%d)\n", name[i], pixels);

      for (color = 0; color < colors; color++) {
	x3f_printf(DEBUG, "    mean[%d] = %f\n",
		   color,
		   (double)black[color]/pixels);
	black_sum[color] += black[color];
      }
    }

  if (pixels_sum == 0) return 0;

  for (i=0; i<colors; i++)
    black_level[i] = (double)black_sum[i]/pixels_sum;

  pixels_sum = 0;
  black_sqdev = alloca(colors*sizeof(double));
  black_sqdev_sum = alloca(colors*sizeof(double));
  for (i=0; i<colors; i++) black_sqdev_sum[i] = 0.0;

  for (i=0; i<4; i++)
    if (use[i]) {
      int color;
      int pixels = sum_area_sqdev(area[i], colors, black_level, black_sqdev);

      pixels_sum += pixels;

      for (color = 0; color < colors; color++) {
	black_sqdev_sum[color] += black_sqdev[color]; 
      }
    }

  if (pixels_sum == 0) return 0;

  x3f_printf(DEBUG, "  SUM\n");

  for (i=0; i<colors; i++) {
    black_dev[i] = sqrt(black_sqdev_sum[i]/pixels_sum);
    x3f_printf(DEBUG, "    level[%d] = %f\n",
	       i,
	       black_level[i]);
    x3f_printf(DEBUG, "    dev[%d] = %g\n",
	       i,
	       black_dev[i]);
  }

  return 1;
}

static void get_raw_neutral(double *raw_to_xyz, double *raw_neutral)
{
  double d65_xyz[3] = {0.95047, 1.00000, 1.08883};
  double xyz_to_raw[9];

  x3f_3x3_inverse(raw_to_xyz, xyz_to_raw);
  x3f_3x3_3x1_mul(xyz_to_raw, d65_xyz, raw_neutral);
}

/* extern */ int x3f_get_gain(x3f_t *x3f, char *wb, double *gain)
{
  double cam_to_xyz[9], wb_correction[9], gain_fact[3];

  if (x3f_get_camf_matrix_for_wb(x3f, "WhiteBalanceGains", wb, 3, 0, gain) ||
      x3f_get_camf_matrix_for_wb(x3f, "DP1_WhiteBalanceGains", wb, 3, 0, gain));
  else if (x3f_get_camf_matrix_for_wb(x3f, "WhiteBalanceIlluminants", wb,
				      3, 3, cam_to_xyz) &&
	   x3f_get_camf_matrix_for_wb(x3f, "WhiteBalanceCorrections", wb,
				      3, 3, wb_correction)) {
    double raw_to_xyz[9], raw_neutral[3];

    x3f_3x3_3x3_mul(wb_correction, cam_to_xyz, raw_to_xyz);
    get_raw_neutral(raw_to_xyz, raw_neutral);
    x3f_3x1_invert(raw_neutral, gain);
  }
  else
    return 0;

  if (x3f_get_camf_float_vector(x3f, "SensorAdjustmentGainFact", gain_fact))
    x3f_3x1_comp_mul(gain_fact, gain, gain);

  if (x3f_get_camf_float_vector(x3f, "TempGainFact", gain_fact))
    x3f_3x1_comp_mul(gain_fact, gain, gain);

  if (x3f_get_camf_float_vector(x3f, "FNumberGainFact", gain_fact))
    x3f_3x1_comp_mul(gain_fact, gain, gain);

  x3f_printf(DEBUG, "gain\n");
  x3f_3x1_print(DEBUG, gain);

  return 1;
}

/* extern */ int x3f_get_bmt_to_xyz(x3f_t *x3f, char *wb, double *bmt_to_xyz)
{
  double cc_matrix[9], cam_to_xyz[9], wb_correction[9];

  if (x3f_get_camf_matrix_for_wb(x3f, "WhiteBalanceColorCorrections", wb,
				 3, 3, cc_matrix) ||
      x3f_get_camf_matrix_for_wb(x3f, "DP1_WhiteBalanceColorCorrections", wb,
				 3, 3, cc_matrix)) {
    double srgb_to_xyz[9];

    x3f_sRGB_to_XYZ(srgb_to_xyz);
    x3f_3x3_3x3_mul(srgb_to_xyz, cc_matrix, bmt_to_xyz);
  }
  else if (x3f_get_camf_matrix_for_wb(x3f, "WhiteBalanceIlluminants", wb,
				      3, 3, cam_to_xyz) &&
	   x3f_get_camf_matrix_for_wb(x3f, "WhiteBalanceCorrections", wb,
				      3, 3, wb_correction)) {
    double raw_to_xyz[9], raw_neutral[3], raw_neutral_mat[9];

    x3f_3x3_3x3_mul(wb_correction, cam_to_xyz, raw_to_xyz);
    get_raw_neutral(raw_to_xyz, raw_neutral);
    x3f_3x3_diag(raw_neutral, raw_neutral_mat);
    x3f_3x3_3x3_mul(raw_to_xyz, raw_neutral_mat, bmt_to_xyz);
  }
  else
    return 0;

  x3f_printf(DEBUG, "bmt_to_xyz\n");
  x3f_3x3_print(DEBUG, bmt_to_xyz);

  return 1;
}

/* extern */ int x3f_get_raw_to_xyz(x3f_t *x3f, char *wb, double *raw_to_xyz)
{
  double bmt_to_xyz[9], gain[9], gain_mat[9];

  if (!x3f_get_gain(x3f, wb, gain)) return 0;
  if (!x3f_get_bmt_to_xyz(x3f, wb, bmt_to_xyz)) return 0;

  x3f_3x3_diag(gain, gain_mat);
  x3f_3x3_3x3_mul(bmt_to_xyz, gain_mat, raw_to_xyz);

  x3f_printf(DEBUG, "raw_to_xyz\n");
  x3f_3x3_print(DEBUG, raw_to_xyz);

  return 1;
}

/* x3f_denoise expects a 14-bit image since rescaling by a factor of 4
   takes place internally. */
#define INTERMEDIATE_DEPTH 14
#define INTERMEDIATE_UNIT ((1<<INTERMEDIATE_DEPTH) - 1)
#define INTERMEDIATE_BIAS_FACTOR 4.0

static int get_max_intermediate(x3f_t *x3f, char *wb,
				double intermediate_bias,
				uint32_t *max_intermediate)
{
  double gain[3], maxgain = 0.0;
  int i;

  if (!x3f_get_gain(x3f, wb, gain)) return 0;

  /* Cap the gains to 1.0 to avoid clipping */
  for (i=0; i<3; i++)
    if (gain[i] > maxgain) maxgain = gain[i];
  for (i=0; i<3; i++)
    max_intermediate[i] =
      (int32_t)round(gain[i]*(INTERMEDIATE_UNIT - intermediate_bias)/maxgain +
		     intermediate_bias);

  return 1;
}

static int get_intermediate_bias(x3f_t *x3f, char *wb,
				 double *black_level, double *black_dev,
				 double *intermediate_bias)
{
  uint32_t max_raw[3], max_intermediate[3];
  int i;

  if (!x3f_get_max_raw(x3f, max_raw)) return 0;
  if (!get_max_intermediate(x3f, wb, 0, max_intermediate)) return 0;

  *intermediate_bias = 0.0;
  for (i=0; i<3; i++) {
    double bias = INTERMEDIATE_BIAS_FACTOR * black_dev[i] *
      max_intermediate[i] / (max_raw[i] - black_level[i]);
    if (bias > *intermediate_bias) *intermediate_bias = bias;
  }

  return 1;
}

typedef struct bad_pixel_s {
  int c, r;
  struct bad_pixel_s *prev, *next;
} bad_pixel_t;

typedef struct {
  /* c = column, r = row; i = intial, f = final, p = pitch, s = size */
  int ci, cf, cp, cs, ri, rf, rp, rs;
} grid_t;

/* Address pixel at column _c and row _r */
#define _PN(_c, _r, _cs) ((_r)*(_cs) + (_c))

/* Test if a pixel (_c,_r) is within a rectancle */
#define _INB(_c, _r, _cs, _rs)					\
  ((_c) >= 0 && (_c) < (_cs) && (_r) >= 0 && (_r) < (_rs))

/* Test if a pixel has been marked in the bad pixel vector */
#define TEST_PIX(_vec, _c, _r, _cs, _rs)				\
  (_INB((_c), (_r), (_cs), (_rs)) ?					\
   (_vec)[_PN((_c), (_r), (_cs)) >> 5] &				\
   1 << (_PN((_c), (_r), (_cs)) & 0x1f) : 1)

/* Mark the pixel, in the bad pixel vector and the bad pixel list */
#define MARK_PIX(_list, _vec, _c, _r, _cs, _rs)				\
  do {									\
    if (!TEST_PIX((_vec), (_c), (_r), (_cs), (_rs))) {			\
      bad_pixel_t *_p = malloc(sizeof(bad_pixel_t));			\
      _p->c = (_c);							\
      _p->r = (_r);							\
      _p->prev = NULL;							\
      _p->next = (_list);						\
      if (_list) (_list)->prev = (_p);					\
      (_list) = _p;							\
      (_vec)[_PN((_c), (_r), (_cs)) >> 5] |=				\
	1 << (_PN((_c), (_r), (_cs)) & 0x1f);				\
    }									\
    else if (!_INB((_c), (_r), (_cs), (_rs)))				\
      x3f_printf(WARN, "Bad pixel (%u,%u) out of bounds : (%u,%u)\n",   \
		 (_c), (_r), (_cs), (_rs));				\
  } while (0)

/* Clear the mark in the bad pixel vector */
#define CLEAR_PIX(_vec, _c, _r, _cs, _rs)				\
  do {									\
    assert(_INB((_c), (_r), (_cs), (_rs)));				\
    _vec[_PN((_c), (_r), (_cs)) >> 5] &=				\
      ~(1 << (_PN((_c), (_r), (_cs)) & 0x1f));				\
  } while (0)

static void interpolate_bad_pixels(x3f_t *x3f, x3f_area16_t *image, int colors)
{
  bad_pixel_t *bad_pixel_list = NULL;
  uint32_t *bad_pixel_vec = calloc((image->rows*image->columns + 31)/32,
				   sizeof(uint32_t));
  int row, col, color, i;
  uint32_t *bpf23, cameraid;
  int bpf23_len;
  int stat_pass = 0;		/* Statistics */
  int fix_corner = 0;		/* By default, do not accept corners */

  /* BEGIN - collecting bad pixels. This part reads meta data and
     collects all bad pixels both in the list 'bad_pixel_list' and the
     vector 'bad_pixel_vec' */

  if (colors == 3) {
    uint32_t keep[4], hpinfo[4], *bp, *bpf20;
    int bp_num, bpf20_rows, bpf20_cols;

    if (x3f_get_camf_matrix(x3f, "KeepImageArea", 4, 0, 0, M_UINT, keep) &&
	x3f_get_camf_matrix_var(x3f, "BadPixels", &bp_num, NULL, NULL,
				M_UINT, (void **)&bp))
      for (i=0; i < bp_num; i++)
	MARK_PIX(bad_pixel_list, bad_pixel_vec,
		 ((bp[i] & 0x000fff00) >> 8) - keep[0],
		 ((bp[i] & 0xfff00000) >> 20) - keep[1],
		 image->columns, image->rows);

    /* NOTE: the numbers of rows and cols in this matrix are
       interchanged due to bug in camera firmware */
    if (x3f_get_camf_matrix_var(x3f, "BadPixelsF20",
				&bpf20_cols, &bpf20_rows, NULL,
				M_UINT, (void **)&bpf20) && bpf20_cols == 3)
      for (row=0; row < bpf20_rows; row++)
	MARK_PIX(bad_pixel_list, bad_pixel_vec,
		 bpf20[3*row + 1], bpf20[3*row + 0],
		 image->columns, image->rows);

    /* NOTE: the numbers of rows and cols in this matrix are
       interchanged due to bug in camera firmware
       TODO: should Jpeg_BadClutersF20 really be used for RAW? It works
       though. */
    if (x3f_get_camf_matrix_var(x3f, "Jpeg_BadClusters",
				&bpf20_cols, &bpf20_rows, NULL,
				M_UINT, (void **)&bpf20) && bpf20_cols == 3)
      for (row=0; row < bpf20_rows; row++)
	MARK_PIX(bad_pixel_list, bad_pixel_vec,
		 bpf20[3*row + 1], bpf20[3*row + 0],
		 image->columns, image->rows);

    /* TODO: should those really be interpolated over, or should they be
       rescaled instead? */
    if (x3f_get_camf_matrix(x3f, "HighlightPixelsInfo", 2, 2, 0, M_UINT,
			    hpinfo))
      for (row = hpinfo[1]; row < image->rows; row += hpinfo[3])
	for (col = hpinfo[0]; col < image->columns; col += hpinfo[2])
	  MARK_PIX(bad_pixel_list, bad_pixel_vec,
		   col, row, image->columns, image->rows);
  } /* colors == 3 */

  if ((colors == 1 && x3f_get_camf_matrix_var(x3f, "BadPixelsLumaF23",
					      &bpf23_len, NULL, NULL,
					      M_UINT, (void **)&bpf23)) ||
      (colors == 3 && x3f_get_camf_matrix_var(x3f, "BadPixelsChromaF23",
					      &bpf23_len, NULL, NULL,
					      M_UINT, (void **)&bpf23)))
    for (i=0, row=-1; i < bpf23_len; i++)
      if (row == -1) row = bpf23[i];
      else if (bpf23[i] == 0) row = -1;
      else {MARK_PIX(bad_pixel_list, bad_pixel_vec,
		     bpf23[i], row,
		     image->columns, image->rows); i++;}

  /* Interpolate over autofocus pixels for sd Quattro and sd Quattro H.
     TODO: The positions shouldn't really be hardcoded. */
  
  if (x3f_get_camf_unsigned(x3f, "CAMERAID", &cameraid)) {
    const grid_t *g = NULL;

    if (cameraid == X3F_CAMERAID_SDQ) {
      static const grid_t sdq_af_luma    = {217, 5641, 16, 1, 464, 3312, 32, 2};
      static const grid_t sdq_af_chroma  = {108, 2820,  8, 1, 232, 1656, 16, 1};
      g = (colors == 1 ? &sdq_af_luma : &sdq_af_chroma);
    } else if (cameraid == X3F_CAMERAID_SDQH) {
      static const grid_t sdqh_af_luma   = {233, 6425, 16, 1, 592, 3888, 32, 2};
      static const grid_t sdqh_af_chroma = {116, 2820,  8, 1, 296, 1944, 16, 1};
      g = (colors == 1 ? &sdqh_af_luma : &sdqh_af_chroma);
    }

    if (g != NULL) {
      int r, c;

      x3f_printf(DEBUG, "Create AF grid for removing bad pixels\n");

      for (row = g->ri; row <= g->rf; row += g->rp)
	for (col = g->ci; col <= g->cf; col += g->cp)
	  for (r = 0; r < g->rs; r++)
	    for (c = 0; c < g->cs; c++)
	      MARK_PIX(bad_pixel_list, bad_pixel_vec, col+c, row+r,
		       image->columns, image->rows);
    }
  }

  /* END - collecting bad pixels */


  /* BEGIN - fixing bad pixels. This part fixes all bad pixels
     collected in the list 'bad_pixel_list', using the mirror data in
     the vector 'bad_pixel_vec'.  This is made in passes. In each pass
     all pixels that can be interpolated are interpolated and also
     removed from the list of bad pixels.  Eventually the list of bad
     pixels is going to be empty. */

  if (bad_pixel_list)
    x3f_printf(DEBUG, "There are bad pixels to fix\n");

  while (bad_pixel_list) {
    bad_pixel_t *p, *pn;
    bad_pixel_t *fixed = NULL;	/* Contains all, in this pass, fixed pixels */
    struct {
      int all_four, two_linear, two_corner, left; /* Statistics */
    } stats = {0,0,0,0};

    /* Iterate over all pixels in the bad pixel list, in this pass */
    for (p=bad_pixel_list; p && (pn=p->next, 1); p=pn) {
      uint16_t *outp =
	&image->data[p->r*image->row_stride + p->c*image->channels];
      uint16_t *inp[4] = {NULL, NULL, NULL, NULL};
      int num = 0;

      /* Collect status of neighbor pixels */
      if (!TEST_PIX(bad_pixel_vec, p->c - 1, p->r, image->columns, image->rows))
	num++, inp[0] =
	  &image->data[p->r*image->row_stride + (p->c - 1)*image->channels];
      if (!TEST_PIX(bad_pixel_vec, p->c + 1, p->r, image->columns, image->rows))
	num++, inp[1] =
	  &image->data[p->r*image->row_stride + (p->c + 1)*image->channels];
      if (!TEST_PIX(bad_pixel_vec, p->c, p->r - 1, image->columns, image->rows))
	num++, inp[2] =
	  &image->data[(p->r - 1)*image->row_stride + p->c*image->channels];
      if (!TEST_PIX(bad_pixel_vec, p->c, p->r + 1, image->columns, image->rows))
	num++, inp[3] =
	  &image->data[(p->r + 1)*image->row_stride + p->c*image->channels];

      /* Test if interpolation is possible ... */
      if (inp[0] && inp[1] && inp[2] && inp[3])
	/* ... all four neighbors are OK */
	stats.all_four++;
      else if (inp[0] && inp[1])
	/* ... left and right are OK */
	inp[2] = inp[3] = NULL, num = 2, stats.two_linear++;
      else if (inp[2] && inp[3])
	/* ... above and under are OK */
	inp[0] = inp[1] = NULL, num = 2, stats.two_linear++;
      else if (fix_corner && num == 2)
	/* ... corner (plus nothing else to do) are OK */
	stats.two_corner++;
      else
	/* ... nope - it was not possible. Look at next without doing
	   interpolation.  */
	{stats.left++; continue;};

      /* Interpolate the actual pixel */
      for (color=0; color < colors; color++) {
	uint32_t sum = 0;
	for (i=0; i<4; i++)
	  if (inp[i]) sum += inp[i][color];
	outp[color] = (sum + num/2)/num;
      }

      /* Remove p from bad_pixel_list */
      if (p->prev) p->prev->next  = p->next;
      else         bad_pixel_list = p->next;
      if (p->next) p->next->prev = p->prev;

      /* Add p to fixed list */
      p->prev = NULL;
      p->next = fixed;
      fixed = p;
    }

    x3f_printf(DEBUG, "Bad pixels pass %d: %d fixed (%d all_four, %d linear, %d corner), %d left\n",
	       stat_pass,
	       stats.all_four + stats.two_linear + stats.two_corner,
	       stats.all_four,
	       stats.two_linear,
	       stats.two_corner,
	       stats.left);

    if (!fixed) {
      /* If nothing else to do, accept corners */
      if (!fix_corner) fix_corner = 1;
      else {
	x3f_printf(WARN, "Failed to interpolate %d bad pixels\n",
		   stats.left);
	fixed = bad_pixel_list;	/* Free remaining list entries */
	bad_pixel_list = NULL;	/* Force termination */
      }
    }

    /* Clear the bad pixel vector and free the list */
    for (p=fixed; p && (pn=p->next, 1); p=pn) {
      CLEAR_PIX(bad_pixel_vec, p->c, p->r, image->columns, image->rows);
      free(p);
    }

    stat_pass++;
  }

  /* END - fixing bad pixels */

  free(bad_pixel_vec);
}

static int preprocess_data(x3f_t *x3f, int fix_bad, char *wb, x3f_image_levels_t *ilevels)
{
  x3f_area16_t image, qtop;
  int row, col, color;
  uint32_t max_raw[3];
  double scale[3], black_level[3], black_dev[3], intermediate_bias;
  int quattro = x3f_image_area_qtop(x3f, &qtop);
  int colors_in = quattro ? 2 : 3;

  if (!x3f_image_area(x3f, &image) || image.channels < 3) return 0;
  if (quattro && (qtop.channels < 1 ||
		  qtop.rows < 2*image.rows || qtop.columns < 2*image.columns))
    return 0;

  if (!get_black_level(x3f, &image, 1, colors_in, black_level, black_dev) ||
      (quattro && !get_black_level(x3f, &qtop, 0, 1,
				   &black_level[2], &black_dev[2]))) {
    x3f_printf(ERR, "Could not get black level\n");
    return 0;
  }
  x3f_printf(DEBUG, "black_level = {%g,%g,%g}, black_dev = {%g,%g,%g}\n",
	     black_level[0], black_level[1], black_level[2],
	     black_dev[0], black_dev[1], black_dev[2]);

  if (!x3f_get_max_raw(x3f, max_raw)) {
    x3f_printf(ERR, "Could not get maximum RAW level\n");
    return 0;
  }
  x3f_printf(DEBUG, "max_raw = {%u,%u,%u}\n",
	     max_raw[0], max_raw[1], max_raw[2]);

  if (!get_intermediate_bias(x3f, wb, black_level, black_dev,
			     &intermediate_bias)) {
    x3f_printf(ERR, "Could not get intermediate bias\n");
    return 0;
  }
  x3f_printf(DEBUG, "intermediate_bias = %g\n", intermediate_bias);
  ilevels->black[0] = ilevels->black[1] = ilevels->black[2] = intermediate_bias;

  if (!get_max_intermediate(x3f, wb, intermediate_bias, ilevels->white)) {
    x3f_printf(ERR, "Could not get maximum intermediate level\n");
    return 0;
  }
  x3f_printf(DEBUG, "max_intermediate = {%u,%u,%u}\n",
	     ilevels->white[0], ilevels->white[1], ilevels->white[2]);

  for (color = 0; color < 3; color++)
    scale[color] = (ilevels->white[color] - ilevels->black[color]) /
      (max_raw[color] - black_level[color]);

  /* Preprocess image data (HUF/TRU->x3rgb16) */
  for (row = 0; row < image.rows; row++)
    for (col = 0; col < image.columns; col++)
      for (color = 0; color < colors_in; color++) {
	uint16_t *valp =
	  &image.data[image.row_stride*row + image.channels*col + color];
	int32_t out =
	  (int32_t)round(scale[color] * (*valp - black_level[color]) +
			 ilevels->black[color]);

	if (out < 0) *valp = 0;
	else if (out > 65535) *valp = 65535;
	else *valp = out;
      }

  if (quattro) {
    /* Preprocess and downsample Quattro top layer (Q->top16) */
    for (row = 0; row < image.rows; row++)
      for (col = 0; col < image.columns; col++) {
	uint16_t *outp =
	  &image.data[image.row_stride*row + image.channels*col + 2];
	uint16_t *row1 =
	  &qtop.data[qtop.row_stride*2*row + qtop.channels*2*col];
	uint16_t *row2 =
	  &qtop.data[qtop.row_stride*(2*row+1) + qtop.channels*2*col];
	uint32_t sum =
	  row1[0] + row1[qtop.channels] + row2[0] + row2[qtop.channels];
	int32_t out = (int32_t)round(scale[2] * (sum/4.0 - black_level[2]) +
				     ilevels->black[2]);

	if (out < 0) *outp = 0;
	else if (out > 65535) *outp = 65535;
	else *outp = out;
      }

    /* Preprocess Quattro top layer (Q->top16) at full resolution */
    for (row = 0; row < qtop.rows; row++)
      for (col = 0; col < qtop.columns; col++) {
	uint16_t *valp = &qtop.data[qtop.row_stride*row + qtop.channels*col];
	int32_t out = (int32_t)round(scale[2] * (*valp - black_level[2]) +
				     ilevels->black[2]);

	if (out < 0) *valp = 0;
	else if (out > 65535) *valp = 65535;
	else *valp = out;
      }
    if (fix_bad) interpolate_bad_pixels(x3f, &qtop, 1);
  }

  if (fix_bad) interpolate_bad_pixels(x3f, &image, 3);

  return 1;
}

static int get_conv(x3f_t *x3f, x3f_color_encoding_t encoding, char *wb,
		    int lutsize, uint16_t max_out, double *lut,
		    double *conv_matrix)
{
  double raw_to_xyz[9];	/* White point for XYZ is assumed to be D65 */
  double xyz_to_rgb[9];
  double raw_to_rgb[9];
  double sensor_iso, capture_iso, iso_scaling;

  if (x3f_get_camf_float(x3f, "SensorISO", &sensor_iso) &&
      x3f_get_camf_float(x3f, "CaptureISO", &capture_iso)) {
    x3f_printf(DEBUG, "SensorISO = %g\n", sensor_iso);
    x3f_printf(DEBUG, "CaptureISO = %g\n", capture_iso);
    iso_scaling = capture_iso/sensor_iso;
  }
  else {
    iso_scaling = 1.0;
    x3f_printf(WARN, "Could not calculate ISO scaling, assuming %g\n",
	       iso_scaling);
  }

  if (!x3f_get_raw_to_xyz(x3f, wb, raw_to_xyz)) {
    x3f_printf(ERR, "Could not get raw_to_xyz for white balance: %s\n", wb);
    return 0;
  }

  switch (encoding) {
  case SRGB:
    x3f_sRGB_LUT(lut, lutsize, max_out);
    x3f_XYZ_to_sRGB(xyz_to_rgb);
    break;
  case ARGB:
    x3f_gamma_LUT(lut, lutsize, max_out, 2.2);
    x3f_XYZ_to_AdobeRGB(xyz_to_rgb);
    break;
  case PPRGB:
    {
      double xyz_to_prophotorgb[9], d65_to_d50[9];

      x3f_gamma_LUT(lut, lutsize, max_out, 1.8);
      x3f_XYZ_to_ProPhotoRGB(xyz_to_prophotorgb);
      /* The standad white point for ProPhoto RGB is D50 */
      x3f_Bradford_D65_to_D50(d65_to_d50);
      x3f_3x3_3x3_mul(xyz_to_prophotorgb, d65_to_d50, xyz_to_rgb);
    }
    break;
  default:
    x3f_printf(ERR, "Unknown color space %d\n", encoding);
    return 0;
  }

  x3f_3x3_3x3_mul(xyz_to_rgb, raw_to_xyz, raw_to_rgb);
  x3f_scalar_3x3_mul(iso_scaling, raw_to_rgb, conv_matrix);

  x3f_printf(DEBUG, "raw_to_rgb\n");
  x3f_3x3_print(DEBUG, raw_to_rgb);
  x3f_printf(DEBUG, "conv_matrix\n");
  x3f_3x3_print(DEBUG, conv_matrix);

  return 1;
}

/* Converts the data in place */

#define LUTSIZE 1024

static int convert_data(x3f_t *x3f,
			x3f_area16_t *image, x3f_image_levels_t *ilevels,
			x3f_color_encoding_t encoding,
			int apply_sgain,
			char *wb)
{
  int row, col, color;
  uint16_t max_out = 65535; /* TODO: should be possible to adjust */

  double conv_matrix[9];
  double lut[LUTSIZE];
  x3f_spatial_gain_corr_t sgain[MAXCORR];
  int sgain_num;

  if (image->channels < 3) return 0;

  if (!get_conv(x3f, encoding, wb, LUTSIZE, max_out, lut, conv_matrix))
    return 0;

  if (apply_sgain) {
    sgain_num = x3f_get_spatial_gain(x3f, wb, sgain);
    if (sgain_num == 0)
      x3f_printf(WARN, "Could not get spatial gain\n");
  } else {
    sgain_num = 0;
  }

  for (row = 0; row < image->rows; row++) {
    for (col = 0; col < image->columns; col++) {
      uint16_t *valp[3];
      double input[3], output[3];

      /* Get the data */
      for (color = 0; color < 3; color++) {
	valp[color] =
	  &image->data[image->row_stride*row + image->channels*col + color];
	input[color] = x3f_calc_spatial_gain(sgain, sgain_num,
					     row, col, color,
					     image->rows, image->columns) *
	  (*valp[color] - ilevels->black[color]) /
	  (ilevels->white[color] - ilevels->black[color]);
      }

      /* Do color conversion */
      x3f_3x3_3x1_mul(conv_matrix, input, output);

      /* Write back the data, doing non linear coding */
      for (color = 0; color < 3; color++)
	*valp[color] = x3f_LUT_lookup(lut, LUTSIZE, output[color]);
    }
  }

  x3f_cleanup_spatial_gain(sgain, sgain_num);

  ilevels->black[0] = ilevels->black[1] = ilevels->black[2] = 0.0;
  ilevels->white[0] = ilevels->white[1] = ilevels->white[2] = max_out;

  return 1;
}

static int run_denoising(x3f_t *x3f)
{
  x3f_area16_t original_image, image;
  x3f_denoise_type_t type = X3F_DENOISE_STD;
  char *sensorid;

  if (!x3f_image_area(x3f, &original_image)) return 0;
  if (!x3f_crop_area_camf(x3f, "ActiveImageArea", &original_image, 1, &image)) {
    image = original_image;
    x3f_printf(WARN, "Could not get active area, denoising entire image\n");
  }

  if (x3f_get_prop_entry(x3f, "SENSORID", &sensorid) &&
      !strcmp(sensorid, "F20"))
    type = X3F_DENOISE_F20;

  x3f_denoise(&image, type);
  return 1;
}

static int expand_quattro(x3f_t *x3f, int denoise, x3f_area16_t *expanded)
{
  x3f_area16_t image, active, qtop, qtop_crop, active_exp;
  uint32_t rect[4];

  if (!x3f_image_area_qtop(x3f, &qtop)) return 0;
  if (!x3f_image_area(x3f, &image)) return 0;
  if (denoise &&
      !x3f_crop_area_camf(x3f, "ActiveImageArea", &image, 1, &active)) {
    active = image;
    x3f_printf(WARN, "Could not get active area, denoising entire image\n");
  }

  rect[0] = 0;
  rect[1] = 0;
  rect[2] = 2*image.columns - 1;
  rect[3] = 2*image.rows - 1;
  if (!x3f_crop_area(rect, &qtop, &qtop_crop)) return 0;

  expanded->columns = qtop_crop.columns;
  expanded->rows = qtop_crop.rows;
  expanded->channels = 3;
  expanded->row_stride = expanded->columns*expanded->channels;
  expanded->data = expanded->buf =
    malloc(expanded->rows*expanded->row_stride*sizeof(uint16_t));

  if (denoise && !x3f_crop_area_camf(x3f, "ActiveImageArea", expanded, 0,
				     &active_exp)) {
    active_exp = *expanded;
    x3f_printf(WARN, "Could not get active area, denoising entire image\n");
  }

  x3f_expand_quattro(&image, denoise ? &active : NULL, &qtop_crop,
		     expanded, denoise ? &active_exp : NULL);

  return 1;
}

/* extern */ int x3f_get_image(x3f_t *x3f,
			       x3f_area16_t *image,
			       x3f_image_levels_t *ilevels,
			       x3f_color_encoding_t encoding,
			       int crop,
			       int fix_bad,
			       int denoise,
			       int apply_sgain,
			       char *wb)
{
  x3f_area16_t original_image, expanded;
  x3f_image_levels_t il;

  if (wb == NULL) wb = x3f_get_wb(x3f);

  if (encoding == QTOP) {
    x3f_area16_t qtop;

    if (!x3f_image_area_qtop(x3f, &qtop)) return 0;
    if (!crop || !x3f_crop_area_camf(x3f, "ActiveImageArea", &qtop, 0, image))
      *image = qtop;

    return ilevels == NULL;
  }

  if (!x3f_image_area(x3f, &original_image)) return 0;
  if (!crop || !x3f_crop_area_camf(x3f, "ActiveImageArea", &original_image, 1,
				   image))
    *image = original_image;

  if (encoding == UNPROCESSED) return ilevels == NULL;

  if (!preprocess_data(x3f, fix_bad, wb, &il)) return 0;

  if (expand_quattro(x3f, denoise, &expanded)) {
    /* NOTE: expand_quattro destroys the data of original_image */
    if (!crop ||
	!x3f_crop_area_camf(x3f, "ActiveImageArea", &expanded, 0, image))
      *image = expanded;
    original_image = expanded;
  }
  else if (denoise && !run_denoising(x3f)) return 0;

  if (encoding != NONE &&
      !convert_data(x3f, &original_image, &il, encoding, apply_sgain, wb)) {
    free(image->buf);
    return 0;
  }

  if (ilevels) *ilevels = il;
  return 1;
}

/* extern */ int x3f_get_preview(x3f_t *x3f,
				 x3f_area16_t *image,
				 x3f_image_levels_t *ilevels,
				 x3f_color_encoding_t encoding,
				 int apply_sgain,
				 char *wb,
				 uint32_t max_width,
				 x3f_area8_t *preview)
{
  int row, col, color;
  uint16_t max_out = 255;

  double conv_matrix[9];
  double lut[LUTSIZE];
  x3f_spatial_gain_corr_t sgain[MAXCORR];
  int sgain_num;

  int reduction, reduction2;

  if (image->channels < 3) return 0;

  if (!get_conv(x3f, encoding, wb, LUTSIZE, max_out, lut, conv_matrix))
    return 0;

  if (apply_sgain) {
    sgain_num = x3f_get_spatial_gain(x3f, wb, sgain);
    if (sgain_num == 0)
      x3f_printf(WARN, "Could not get spatial gain\n");
  } else {
    sgain_num = 0;
  }

  reduction = (image->columns + max_width - 1)/max_width;
  reduction2 = reduction*reduction;
  preview->columns = image->columns/reduction;
  preview->rows = image->rows/reduction;
  preview->channels = 3;
  preview->row_stride = preview->columns*preview->channels;
  preview->data = preview->buf =
    malloc(preview->rows*preview->row_stride*sizeof(uint8_t));

  for (row = 0; row < preview->rows; row++) {
    for (col = 0; col < preview->columns; col++) {
      double input[3], output[3];

      /* Get the data */
      for (color = 0; color < 3; color++) {
	uint32_t acc = 0;
	int r, c;

	for (r=0; r<reduction; r++)
	  for (c=0; c<reduction; c++)
	    acc += image->data[image->row_stride*(row*reduction + r) +
			       image->channels*(col*reduction + c) + color];

	input[color] = x3f_calc_spatial_gain(sgain, sgain_num,
					     row, col, color,
					     preview->rows, preview->columns) *
	  ((double)acc/reduction2 - ilevels->black[color]) /
	  (ilevels->white[color] - ilevels->black[color]);
      }

      /* Do color conversion */
      x3f_3x3_3x1_mul(conv_matrix, input, output);

      /* Write back the data, doing non linear coding */
      for (color = 0; color < 3; color++)
	preview->data[preview->row_stride*row + preview->channels*col + color] =
	  x3f_LUT_lookup(lut, LUTSIZE, output[color]);
    }
  }

  x3f_cleanup_spatial_gain(sgain, sgain_num);

  x3f_crop_area8_camf(x3f, "ActiveImageArea", preview, 1, preview);

  return 1;
}


<|start_filename|>src/x3f_dngtags.h<|end_filename|>
/* X3F_DNGTAGS.H
 *
 * libtiff tag extender for DNG tags.
 *
 * The current version of libtiff (4.0.3) only supports tags for DNG
 * versions <= 1.1.0.0
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#ifndef X3F_DNGTAGS_H
#define X3F_DNGTAGS_H

#include <inttypes.h>
#include <assert.h>

#define PHOTOMETRIC_CFA 32803
#define PHOTOMETRIC_LINEARRAW 34892

#define TIFFTAG_EXTRACAMERAPROFILES 50933
#define TIFFTAG_PROFILENAME 50936
#define TIFFTAG_ASSHOTPROFILENAME 50934
#define TIFFTAG_OPCODELIST1 51008
#define TIFFTAG_OPCODELIST2 51009
#define TIFFTAG_OPCODELIST3 51022
#define TIFFTAG_DEFAULTBLACKRENDER 51110
#define TIFFTAG_FORWARDMATRIX1 50964
#define TIFFTAG_FORWARDMATRIX2 50965

#define DNG_OPCODEFLAG_OPTIONAL 1
#define DNG_OPCODEFLAG_SKIPPREVIEW 2

#define DNG_OPCODE_GAINMAP_ID 9
#define DNG_OPCODE_GAINMAP_VER 0x01030000

typedef struct __attribute__((packed)) {
  uint32_t count;
} dng_opcodelist_header_t;

typedef struct __attribute__((packed)) {
  uint32_t id;
  uint32_t ver;
  uint32_t flags;
  uint32_t parsize;
} dng_opcode_header_t;

typedef struct __attribute__((packed)) {
  dng_opcode_header_t header;
  uint32_t Top;
  uint32_t Left;
  uint32_t Bottom;
  uint32_t Right;
  uint32_t Plane;
  uint32_t Planes;
  uint32_t RowPitch;
  uint32_t ColPitch;
  uint32_t MapPointsV;
  uint32_t MapPointsH;
  double MapSpacingV;
  double MapSpacingH;
  double MapOriginV;
  double MapOriginH;
  uint32_t MapPlanes;
  float MapGain[0];
} dng_opcode_GainMap_t;

typedef uint8_t uint8_alias_t __attribute__((may_alias));
typedef uint16_t uint16_alias_t __attribute__((may_alias));
typedef uint32_t uint32_alias_t __attribute__((may_alias));
typedef uint64_t uint64_alias_t __attribute__((may_alias));

#define PUT_BIG_16(_dst, _src)				\
do {							\
  assert(sizeof(_dst) == 2);				\
  typeof(_dst) _tmp = (_src);	/* For type checking */	\
  uint16_t _srcu = *(uint16_alias_t *)&_tmp;		\
  uint8_alias_t *_dstp = (uint8_alias_t *)&(_dst);	\
							\
  _dstp[0] = (_srcu&0xff00) >> 8;			\
  _dstp[1] = (_srcu&0x00ff) >> 0;			\
 } while (0)

#define PUT_BIG_32(_dst, _src)				\
do {							\
  assert(sizeof(_dst) == 4);				\
  typeof(_dst) _tmp = (_src);	/* For type checking */	\
  uint32_t _srcu = *(uint32_alias_t *)&_tmp;		\
  uint8_alias_t *_dstp = (uint8_alias_t *)&(_dst);	\
							\
  _dstp[0] = (_srcu&0xff000000) >> 24;			\
  _dstp[1] = (_srcu&0x00ff0000) >> 16;			\
  _dstp[2] = (_srcu&0x0000ff00) >> 8;			\
  _dstp[3] = (_srcu&0x000000ff) >> 0;			\
 } while (0)

#define PUT_BIG_64(_dst, _src)				\
do {							\
  assert(sizeof(_dst) == 8);				\
  typeof(_dst) _tmp = (_src);	/* For type checking */	\
  uint64_t _srcu = *(uint64_alias_t *)&_tmp;		\
  uint8_alias_t *_dstp = (uint8_alias_t *)&(_dst);	\
							\
  _dstp[0] = (_srcu&0xff00000000000000) >> 56;		\
  _dstp[1] = (_srcu&0x00ff000000000000) >> 48;		\
  _dstp[2] = (_srcu&0x0000ff0000000000) >> 40;		\
  _dstp[3] = (_srcu&0x000000ff00000000) >> 32;		\
  _dstp[4] = (_srcu&0x00000000ff000000) >> 24;		\
  _dstp[5] = (_srcu&0x0000000000ff0000) >> 16;		\
  _dstp[6] = (_srcu&0x000000000000ff00) >> 8;		\
  _dstp[7] = (_srcu&0x00000000000000ff) >> 0;		\
 } while (0)

#endif	/* X3F_DNGTAGS_H */


<|start_filename|>src/x3f_matrix.c<|end_filename|>
/* X3F_MATRIX.C
 *
 * Library for manipulating matrices.
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#include <math.h>

#include "x3f_matrix.h"
#include "x3f_printf.h"

#define M0(m) *(m+0)
#define M1(m) *(m+1)
#define M2(m) *(m+2)

#define M00(m) *(m+0)
#define M01(m) *(m+1)
#define M02(m) *(m+2)
#define M10(m) *(m+3)
#define M11(m) *(m+4)
#define M12(m) *(m+5)
#define M20(m) *(m+6)
#define M21(m) *(m+7)
#define M22(m) *(m+8)

void x3f_3x1_invert(double *a, double *ainv)
{
  M0(ainv) = 1.0/M0(a); M1(ainv) = 1.0/M1(a); M2(ainv) = 1.0/M2(a);
}

void x3f_3x1_comp_mul(double *a, double *b, double *c)
{
  M0(c) = M0(a)*M0(b); M1(c) = M1(a)*M1(b); M2(c) = M2(a)*M2(b);
}

/* Multiply a scalar with a 3x1 matrix, giving a 3x1 matrix */
void x3f_scalar_3x1_mul(double a, double *b, double *c)
{
  M0(c) = a*M0(b); M1(c) = a*M1(b); M2(c) = a*M2(b);
}

/* Multiply a scalar with a 3x3 matrix, giving a 3x3 matrix */
void x3f_scalar_3x3_mul(double a, double *b, double *c)
{
  M00(c) = a*M00(b); M01(c) = a*M01(b); M02(c) = a*M02(b);
  M10(c) = a*M10(b); M11(c) = a*M11(b); M12(c) = a*M12(b);
  M20(c) = a*M20(b); M21(c) = a*M21(b); M22(c) = a*M22(b);
}

/* Multiply a 3x3 matrix with a 3x1 matrix, giving a 3x1 matrix */
void x3f_3x3_3x1_mul(double *a, double *b, double *c)
{
  M0(c) = M00(a)*M0(b) + M01(a)*M1(b) + M02(a)*M2(b);
  M1(c) = M10(a)*M0(b) + M11(a)*M1(b) + M12(a)*M2(b);
  M2(c) = M20(a)*M0(b) + M21(a)*M1(b) + M22(a)*M2(b);
}

/* Multiply a 3x3 matrix with a 3x3 matrix, giving a 3x3 matrix */
void x3f_3x3_3x3_mul(double *a, double *b, double *c)
{
  M00(c) = M00(a)*M00(b) + M01(a)*M10(b) + M02(a)*M20(b);
  M01(c) = M00(a)*M01(b) + M01(a)*M11(b) + M02(a)*M21(b);
  M02(c) = M00(a)*M02(b) + M01(a)*M12(b) + M02(a)*M22(b);

  M10(c) = M10(a)*M00(b) + M11(a)*M10(b) + M12(a)*M20(b);
  M11(c) = M10(a)*M01(b) + M11(a)*M11(b) + M12(a)*M21(b);
  M12(c) = M10(a)*M02(b) + M11(a)*M12(b) + M12(a)*M22(b);

  M20(c) = M20(a)*M00(b) + M21(a)*M10(b) + M22(a)*M20(b);
  M21(c) = M20(a)*M01(b) + M21(a)*M11(b) + M22(a)*M21(b);
  M22(c) = M20(a)*M02(b) + M21(a)*M12(b) + M22(a)*M22(b);
}

/* Calculate the inverse of a 3x3 matrix
   http://en.wikipedia.org/wiki/Invertible_matrix */
void x3f_3x3_inverse(double *a, double *ainv)
{
  double A, B, C, D, E, F, G, H, I;
  double det;

  A = +(M11(a)*M22(a)-M12(a)*M21(a));
  B = -(M10(a)*M22(a)-M12(a)*M20(a));
  C = +(M10(a)*M21(a)-M11(a)*M20(a));

  D = -(M01(a)*M22(a)-M02(a)*M21(a));
  E = +(M00(a)*M22(a)-M02(a)*M20(a));
  F = -(M00(a)*M21(a)-M01(a)*M20(a));

  G = +(M01(a)*M12(a)-M02(a)*M11(a));
  H = -(M00(a)*M12(a)-M02(a)*M10(a));
  I = +(M00(a)*M11(a)-M01(a)*M10(a));

  det = M00(a)*A + M01(a)*B + M02(a)*C;

  M00(ainv) = A/det; M01(ainv) = D/det; M02(ainv) = G/det;
  M10(ainv) = B/det; M11(ainv) = E/det; M12(ainv) = H/det;
  M20(ainv) = C/det; M21(ainv) = F/det; M22(ainv) = I/det;
}

/* Convert a 3x1 matrix to a 3x3 diagonal matrix */
void x3f_3x3_diag(double *a, double *b)
{
  M00(b) = M0(a); M01(b) = 0.0;   M02(b) = 0.0;
  M10(b) = 0.0;   M11(b) = M1(a); M12(b) = 0.0;
  M20(b) = 0.0;   M21(b) = 0.0;   M22(b) = M2(a);
}

void x3f_3x3_identity(double *a)
{
  M00(a) = 1.0; M01(a) = 0.0; M02(a) = 0.0;
  M10(a) = 0.0; M11(a) = 1.0; M12(a) = 0.0;
  M20(a) = 0.0; M21(a) = 0.0; M22(a) = 1.0;
}

void x3f_3x3_ones(double *a)
{
  M00(a) = 1.0; M01(a) = 1.0; M02(a) = 1.0;
  M10(a) = 1.0; M11(a) = 1.0; M12(a) = 1.0;
  M20(a) = 1.0; M21(a) = 1.0; M22(a) = 1.0;
}

/* Print a 3x1 matrix */
void x3f_3x1_print(x3f_verbosity_t level, double *a)
{
  x3f_printf(level, "%10g\n", M0(a));
  x3f_printf(level, "%10g\n", M1(a));
  x3f_printf(level, "%10g\n", M2(a));
}

/* Print a 3x3 matrix */
void x3f_3x3_print(x3f_verbosity_t level, double *a)
{
  x3f_printf(level, "%10g %10g %10g\n", M00(a), M01(a), M02(a));
  x3f_printf(level, "%10g %10g %10g\n", M10(a), M11(a), M12(a));
  x3f_printf(level, "%10g %10g %10g\n", M20(a), M21(a), M22(a));
}

/* Conversion matrices for color spaces. NOTE - the result is, of
   course, linear and needs some gamma encoding or similar.  */

/* http://www.color.org/ROMMRGB.pdf */
/* http://en.wikipedia.org/wiki/ProPhoto_RGB_color_space */
void x3f_XYZ_to_ProPhotoRGB(double *a)
{
  M00(a) = +1.3460; M01(a) = -0.2556; M02(a) = -0.0511;
  M10(a) = -0.5446; M11(a) = +1.5082; M12(a) = +0.0205;
  M20(a) = +0.0000; M21(a) = +0.0000; M22(a) = +1.2123;
}

void x3f_ProPhotoRGB_to_XYZ(double *a)
{
  M00(a) = 0.7977; M01(a) = 0.1352; M02(a) = 0.0313;
  M10(a) = 0.2880; M11(a) = 0.7119; M12(a) = 0.0001;
  M20(a) = 0.0000; M21(a) = 0.0000; M22(a) = 0.8249;
}

/* http://en.wikipedia.org/wiki/Adobe_RGB_color_space */
/* http://www.adobe.com/digitalimag/pdfs/AdobeRGB1998.pdf */
void x3f_XYZ_to_AdobeRGB(double *a)
{
  M00(a) = +2.04159; M01(a) = -0.56501; M02(a) = -0.34473;
  M10(a) = -0.96924; M11(a) = +1.87597; M12(a) = +0.04156;
  M20(a) = +0.01344; M21(a) = -0.11836; M22(a) = +1.01517;
}

void x3f_AdobeRGB_to_XYZ(double *a)
{
  M00(a) = 0.57667; M01(a) = 0.18556; M02(a) = 0.18823;
  M10(a) = 0.29737; M11(a) = 0.62736; M12(a) = 0.07529;
  M20(a) = 0.02703; M21(a) = 0.07069; M22(a) = 0.99134;
}

/* http://en.wikipedia.org/wiki/SRGB */
void x3f_XYZ_to_sRGB(double *a)
{
  M00(a) = +3.2406; M01(a) = -1.5372; M02(a) = -0.4986;
  M10(a) = -0.9689; M11(a) = +1.8758; M12(a) = +0.0415;
  M20(a) = +0.0557; M21(a) = -0.2040; M22(a) = +1.0570;
}

void x3f_sRGB_to_XYZ(double *a)
{
  M00(a) = 0.4124; M01(a) = 0.3576; M02(a) = 0.1805;
  M10(a) = 0.2126; M11(a) = 0.7152; M12(a) = 0.0722;
  M20(a) = 0.0193; M21(a) = 0.1192; M22(a) = 0.9505;
}

/* http://en.wikipedia.org/wiki/CIE_1931_color_space */
void x3f_CIERGB_to_XYZ(double *a)
{
  M00(a) = 0.49    ; M01(a) = 0.31    ; M02(a) = 0.20    ;
  M10(a) = 0.17697 ; M11(a) = 0.81240 ; M12(a) = 0.01063 ;
  M20(a) = 0.00    ; M21(a) = 0.01    ; M22(a) = 0.99    ;
}

/* http://www.brucelindbloom.com/index.html?Eqn_ChromAdapt.html */
void x3f_Bradford_D50_to_D65(double *a)
{
  M00(a) = +0.9555766 ; M01(a) = -0.0230393; M02(a) = +0.0631636;
  M10(a) = -0.0282895 ; M11(a) = +1.0099416; M12(a) = +0.0210077;
  M20(a) = +0.0122982 ; M21(a) = -0.0204830; M22(a) = +1.3299098;
}

void x3f_Bradford_D65_to_D50(double *a)
{
  M00(a) = +1.0478112 ; M01(a) = +0.0228866; M02(a) = -0.0501270;
  M10(a) = +0.0295424 ; M11(a) = +0.9904844; M12(a) = -0.0170491;
  M20(a) = -0.0092345 ; M21(a) = +0.0150436; M22(a) = +0.7521316;
}

void x3f_sRGB_LUT(double *lut, int size, uint16_t max)
{
  double a = 0.055;
  double thres = 0.0031308;
  int i;

  for (i=0; i<size; i++) {
    double lin = (double)i/(size - 1);
    double srgb;

    if (lin <= thres)
      srgb = 12.92 * lin;
    else
      srgb = (1 + a)*pow(lin, 1/2.4) - a;

    srgb *= max;

    if (srgb < 0)
      lut[i] = 0;
    else if (srgb > max)
      lut[i] = max;
    else
      lut[i] = srgb;
  }
}

void x3f_gamma_LUT(double *lut, int size, uint16_t max, double gamma)
{
  int i;

  for (i=0; i<size; i++) {
    double lin = (double)i/(size - 1);
    double gam;

    gam = pow(lin, 1/gamma);
    gam *= max;

    if (gam < 0)
      lut[i] = 0;
    else if (gam > max)
      lut[i] = max;
    else
      lut[i] = gam;
  }
}

uint16_t x3f_LUT_lookup(double *lut, int size, double val)
{
  double index = val*(size - 1);
  int i = (int)floor(index);
  double frac = index - i;

  if (i<0)
    return (uint16_t)round(lut[0]);
  else if (i>=(size - 1))
    return (uint16_t)round(lut[size-1]);
  else
    return (uint16_t)round(lut[i] + frac*(lut[i+1] - lut[i]));
}


<|start_filename|>src/x3f_spatial_gain.c<|end_filename|>
/* X3F_SPATIAL_GAIN.C
 *
 * Library for adjusting for spatial gain in X3F images.
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#include "x3f_spatial_gain.h"
#include "x3f_meta.h"
#include "x3f_printf.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

static double lens_position(double focal_length, double object_distance)
{
  return 1.0/(1.0/focal_length - 1.0/object_distance);
}

static double get_focal_length(x3f_t *x3f)
{
  char *flength;
  double focal_length;

  if (x3f_get_prop_entry(x3f, "FLENGTH", &flength))
    focal_length = atof(flength);
  else {
    focal_length = 30.0;
    x3f_printf(WARN, "Could not get focal length, assuming %g mm\n",
	       focal_length);
  }

  return focal_length;
}

static double get_object_distance(x3f_t *x3f)
{
  double object_distance;

  if (x3f_get_camf_float(x3f, "ObjectDistance", &object_distance))
    object_distance *= 10.0;	/* Convert cm to mm */
  else {
    object_distance = INFINITY;
    x3f_printf(WARN, "Could not get object distance, assuming %g mm\n",
	       object_distance);
  }

  return object_distance;
}

/* Get the minimum object distance (MOD)  */
static double get_MOD(x3f_t *x3f)
{
  int32_t lens_information;
  double mod;

  if (!x3f_get_camf_signed(x3f, "LensInformation", &lens_information))
    lens_information = -1;

  /* TODO: is there any better way of obtaining MOD? */
  switch (lens_information) {
  case 1003:			/* DP1 Merrill */
    mod = 200.0;
    break;
  case 1004:			/* DP2 Merrill */
    mod = 280.0;
    break;
  case 1005:			/* DP3 Merrill */
    mod = 226.0;
    break;
  default:
    mod = 280.0;
    x3f_printf(WARN, "Could not get MOD, assuming %g mm\n", mod);
  }

  return mod;
}

static int get_merrill_type_gains_table(x3f_t *x3f, char *name, char *chan,
					uint32_t **mgain, int *rows, int *cols,
					double *mingain, double *delta)
{
  char table[32], *val;
  int rows_tmp, cols_tmp;

  sprintf(table, "GainsTable%s", chan);
  if (!x3f_get_camf_property(x3f, name, table, &val) ||
      !x3f_get_camf_matrix_var(x3f, val, &rows_tmp, &cols_tmp, NULL,
			       M_UINT, (void **)mgain) ||
      (*rows != -1 && *rows != rows_tmp) ||
      (*cols != -1 && *cols != cols_tmp))
    return 0;
  *rows = rows_tmp;
  *cols = cols_tmp;

  sprintf(table, "MinGains%s", chan);
  if (!x3f_get_camf_property(x3f, name, table, &val)) return 0;
  *mingain = atof(val);

  sprintf(table, "Delta%s", chan);
  if (!x3f_get_camf_property(x3f, name, table, &val)) return 0;
  *delta = atof(val);

  return 1;
}

typedef struct merrill_spatial_gain_s {
  char *name;
  double x,y;
  struct merrill_spatial_gain_s *next;
} merrill_spatial_gain_t;


int x3f_get_merrill_type_spatial_gain(x3f_t *x3f, int hp_flag,
				      x3f_spatial_gain_corr_t *corr)
{
  char **include_blocks, **include_blocks_val;
  uint32_t include_blocks_num;
  double capture_aperture;
  double *spatial_gain_fstop;
  int num_fstop;
  int corr_num = 3;

  merrill_spatial_gain_t *gains = NULL, *g;
  merrill_spatial_gain_t *q_closest[4] = {NULL, NULL, NULL, NULL};
  double q_closest_dx[4] = {INFINITY, -INFINITY, -INFINITY, INFINITY};
  double q_closest_dy[4] = {INFINITY, INFINITY, -INFINITY, -INFINITY};
  double q_closest_d2[4] = {INFINITY, INFINITY, INFINITY, INFINITY};
  double q_weight_x[4], q_weight_y[4], q_weight[4];
  double x,y;
  int i;

  if (!x3f_get_camf_float(x3f, "CaptureAperture", &capture_aperture)) return 0;

  if (!x3f_get_camf_property_list(x3f, "IncludeBlocks",
				  &include_blocks, &include_blocks_val,
				  &include_blocks_num))
    return 0;

  if (hp_flag) {		/* Quattro HP */
    if (!x3f_get_camf_matrix_var(x3f, "SpatialGainHP_Fstop",
				 &num_fstop, NULL, NULL,
				 M_FLOAT, (void **)&spatial_gain_fstop))
      return 0;
    corr_num = 6; /* R, G, B0, B1, B2, B3 */

    for (i=0; i<include_blocks_num; i++) {
      char **names, **values;
      uint32_t num, aperture_index;
      char dummy;

      if (sscanf(include_blocks[i], "SpatialGainHPProps_%d%c",
		 &aperture_index, &dummy) == 1 &&
	  x3f_get_camf_property_list(x3f, include_blocks[i],
				     &names, &values, &num) &&
	  aperture_index >= 0 && aperture_index < num_fstop) {
	g = alloca(sizeof(merrill_spatial_gain_t));
	g->name = include_blocks[i];
	g->x = 1.0/spatial_gain_fstop[aperture_index];
	g->y = 0.0;		/* unused */
	g->next = gains;
	gains = g;
      }
    }

    x = 1.0/capture_aperture;
    y = 0.0;			/* unused */
  }
  else {			/* Merrill and Quattro */
    if (x3f_get_camf_matrix_var(x3f, "SpatialGain_Fstop",
				&num_fstop, NULL, NULL,
				M_FLOAT, (void **)&spatial_gain_fstop))
      for (i=0; i<include_blocks_num; i++) {
	char **names, **values;
	uint32_t num, aperture_index;
	char focus_distance[4], dummy;

	if (sscanf(include_blocks[i], "SpatialGainsProps_%d_%3s%c",
		   &aperture_index, focus_distance, &dummy) == 2 &&
	    x3f_get_camf_property_list(x3f, include_blocks[i],
				       &names, &values, &num) &&
	    aperture_index >= 0 && aperture_index < num_fstop) {
	  double lenspos;

	  if (!strcmp(focus_distance, "INF"))
	    lenspos = lens_position(get_focal_length(x3f), INFINITY);
	  else if (!strcmp(focus_distance, "MOD"))
	    lenspos = lens_position(get_focal_length(x3f), get_MOD(x3f));
	  else continue;

	  g = alloca(sizeof(merrill_spatial_gain_t));
	  g->name = include_blocks[i];
	  g->x = 1.0/spatial_gain_fstop[aperture_index];
	  g->y = lenspos;
	  g->next = gains;
	  gains = g;
	}
      }
    else
      for (i=0; i<include_blocks_num; i++) {
	char **names, **values;
	uint32_t num;
	double aperture, lenspos;
	char dummy;

	if (sscanf(include_blocks[i], "SpatialGainsProps_%lf_%lf%c",
		   &aperture, &lenspos, &dummy) == 2 &&
	    x3f_get_camf_property_list(x3f, include_blocks[i],
				       &names, &values, &num)) {
	  g = alloca(sizeof(merrill_spatial_gain_t));
	  g->name = include_blocks[i];
	  g->x = 1.0/aperture;
	  g->y = lenspos;
	  g->next = gains;
	  gains = g;
	}
      }

    /* TODO: doing bilinear interpolation with respect to 1/aperture
       and lens position. Is that correct? */
    x = 1.0/capture_aperture;
    y = lens_position(get_focal_length(x3f), get_object_distance(x3f));
  }

  for (g=gains; g; g=g->next) {
    double dx = g->x - x;
    double dy = g->y - y;
    double d2 = dx*dx + dy*dy;
    int q;

    if      (dx > 0.0 && dy > 0.0) q = 0;
    else if (dx > 0.0)             q = 3;
    else if (dy > 0.0)             q = 1;
    else                           q = 2;

    if (d2 < q_closest_d2[q]) {
      q_closest[q] = g;
      q_closest_dx[q] = dx;
      q_closest_dy[q] = dy;
      q_closest_d2[q] = d2;
    }
  }

  /* TODO: bilinear interpolation requires that the points be laid out
     in a more or less rectilinear grid. Is that assumption good
     enough? It appears to be valid for the current cameras. */
  q_weight_x[0] = q_closest_dx[1]/(q_closest_dx[1] - q_closest_dx[0]);
  q_weight_x[1] = q_closest_dx[0]/(q_closest_dx[0] - q_closest_dx[1]);
  q_weight_x[2] = q_closest_dx[3]/(q_closest_dx[3] - q_closest_dx[2]);
  q_weight_x[3] = q_closest_dx[2]/(q_closest_dx[2] - q_closest_dx[3]);

  q_weight_y[0] = q_closest_dy[3]/(q_closest_dy[3] - q_closest_dy[0]);
  q_weight_y[1] = q_closest_dy[2]/(q_closest_dy[2] - q_closest_dy[1]);
  q_weight_y[2] = q_closest_dy[1]/(q_closest_dy[1] - q_closest_dy[2]);
  q_weight_y[3] = q_closest_dy[0]/(q_closest_dy[0] - q_closest_dy[3]);

  x3f_printf(DEBUG, "x = %f y = %f\n", x, y);
  for (i=0; i<4; i++) {
    if (q_closest[i])
      x3f_printf(DEBUG, "q = %d name = %s x = %f y = %f\n",
		 i, q_closest[i]->name, q_closest[i]->x, q_closest[i]->y);
    else
      x3f_printf(DEBUG, "q = %d name = NULL\n", i);
    x3f_printf(DEBUG, "q = %d dx = %f dy = %f d2 = %f wx = %f wy = %f\n",
	       i, q_closest_dx[i], q_closest_dy[i], q_closest_d2[i],
	       q_weight_x[i], q_weight_y[i]);
  }

  for (i=0; i<4; i++) {
    if (isnan(q_weight_x[i])) q_weight_x[i] = 1.0;
    if (isnan(q_weight_y[i])) q_weight_y[i] = 1.0;
    q_weight[i] = q_weight_x[i]*q_weight_y[i];
  }

  for (i=0; i<4; i++)
    x3f_printf(DEBUG, "q = %d name = %s w = %f\n",
	       i, q_closest[i] ? q_closest[i]->name : "NULL", q_weight[i]);

  for (i=0; i<corr_num; i++) {
    x3f_spatial_gain_corr_t *c = &corr[i];
    c->gain = NULL;		/* No interploated gains present */
    c->malloc = 0;

    c->rows = c->cols = -1;
    c->rowoff = c->coloff = 0;
    c->rowpitch = c->colpitch = 1;
    c->chan = i;
    c->channels = 1;

    c->mgain_num = 0;
  }

  if (hp_flag) {
    corr[2].rowoff = 0;
    corr[2].coloff = 0;
    corr[2].rowpitch = corr[2].colpitch = 2;
    corr[2].chan = 2;

    corr[3].rowoff = 0;
    corr[3].coloff = 1;
    corr[3].rowpitch = corr[3].colpitch = 2;
    corr[3].chan = 2;

    corr[4].rowoff = 1;
    corr[4].coloff = 0;
    corr[4].rowpitch = corr[4].colpitch = 2;
    corr[4].chan = 2;

    corr[5].rowoff = 1;
    corr[5].coloff = 1;
    corr[5].rowpitch = corr[5].colpitch = 2;
    corr[5].chan = 2;
  }

  for (i=0; i<4; i++) {
    char *name;
    char *channels_normal[3] = {"R", "G", "B"};
    char *channels_hp[6] = {"R", "G", "B0", "B1", "B2", "B3"};
    char **channels = hp_flag ? channels_hp : channels_normal;
    int j;

    if (!q_closest[i]) continue;
    name = q_closest[i]->name;

    for (j=0; j<corr_num; j++) {
      x3f_spatial_gain_corr_t *c = &corr[j];
      x3f_spatial_gain_corr_merrill_t *m = &c->mgain[c->mgain_num++];

      m->weight = q_weight[i];
      if (!get_merrill_type_gains_table(x3f, name, channels[j],
					&m->gain, &c->rows, &c->cols,
					&m->mingain, &m->delta))
	return 0;
    }
  }

  for (i=0; i<corr_num; i++)
    if (corr[i].mgain_num == 0) return 0;

  return corr_num;
}

int x3f_get_interp_merrill_type_spatial_gain(x3f_t *x3f, int hp_flag,
					     x3f_spatial_gain_corr_t *corr)
{
  int corr_num = x3f_get_merrill_type_spatial_gain(x3f, hp_flag, corr);
  int i;

  for (i=0; i<corr_num; i++) {
    x3f_spatial_gain_corr_t *c = &corr[i];
    int j, num = c->rows*c->cols*c->channels;

    c->gain = malloc(num*sizeof(double));
    c->malloc = 1;

    for (j=0; j<num; j++) {
      int g;

      c->gain[j] = 0.0;
      for (g=0; g < c->mgain_num; g++) {
	x3f_spatial_gain_corr_merrill_t *m = &c->mgain[g];
	c->gain[j] += m->weight*(m->mingain + m->delta*m->gain[j]);
      }
    }
  }

  return corr_num;
}

int x3f_get_classic_spatial_gain(x3f_t *x3f, char *wb,
				 x3f_spatial_gain_corr_t *corr)
{
  char *gain_name;

  if ((x3f_get_camf_property(x3f, "SpatialGainTables", wb, &gain_name) &&
       x3f_get_camf_matrix_var(x3f, gain_name,
			       &corr->rows, &corr->cols, &corr->channels,
			       M_FLOAT, (void **)&corr->gain)) ||
      x3f_get_camf_matrix_var(x3f, "SpatialGain",
			      &corr->rows, &corr->cols, &corr->channels,
			      M_FLOAT, (void **)&corr->gain)) {
    corr->malloc = 0;
    corr->rowoff = corr->coloff = 0;
    corr->rowpitch = corr->colpitch = 1;
    corr->chan = 0;
    corr->mgain_num = 0;

    return 1;
  }

  return 0;
}

int x3f_get_spatial_gain(x3f_t *x3f, char *wb, x3f_spatial_gain_corr_t *corr)
{
  int corr_num = 0;

  corr_num += x3f_get_interp_merrill_type_spatial_gain(x3f, 0, &corr[corr_num]);
  if (corr_num == 0)
    corr_num += x3f_get_classic_spatial_gain(x3f, wb, &corr[corr_num]);

  return corr_num;
}

void x3f_cleanup_spatial_gain(x3f_spatial_gain_corr_t *corr, int corr_num)
{
  int i;

  for (i=0; i<corr_num; i++)
    if (corr[i].malloc) free(corr[i].gain);
}

double x3f_calc_spatial_gain(x3f_spatial_gain_corr_t *corr, int corr_num,
			     int row, int col, int chan, int rows, int cols)
{
  double gain = 1.0;
  double rrel = (double)row/rows;
  double crel = (double)col/cols;
  int i;

  for (i=0; i<corr_num; i++) {
    x3f_spatial_gain_corr_t *c = &corr[i];
    double rc, cc;
    int ri, ci;
    double rf, cf;
    double *r1, *r2;
    int co1, co2;
    double gr1, gr2;
    int ch = chan - c->chan;

    if (ch < 0 || ch >= c->channels) continue;

    if (row%c->rowpitch != c->rowoff) continue;
    if (col%c->colpitch != c->coloff) continue;

    rc = rrel*(c->rows-1);
    ri = (int)floor(rc);
    rf = rc - ri;

    cc = crel*(c->cols-1);
    ci = (int)floor(cc);
    cf = cc - ci;

    if (ri < 0)
      r1 = r2 = &c->gain[0];
    else if (ri >= c->rows-1)
      r1 = r2 = &c->gain[(c->rows-1)*c->cols*c->channels];
    else {
      r1 = &c->gain[ri*c->cols*c->channels];
      r2 = &c->gain[(ri+1)*c->cols*c->channels];
    }

    if (ci < 0)
      co1 = co2 = ch;
    if (ci >= c->cols-1)
      co1 = co2 = (c->cols-1)*c->channels + ch;
    else {
      co1 = ci*c->channels + ch;
      co2 = (ci+1)*c->channels + ch;
    }

    gr1 = r1[co1] + cf*(r1[co2]-r1[co1]);
    gr2 = r2[co1] + cf*(r2[co2]-r2[co1]);

    gain *= gr1 + rf*(gr2-gr1);
  }

  return gain;
}


<|start_filename|>src/x3f_dngtags.c<|end_filename|>
/* X3F_DNGTAGS.C
 *
 * libtiff tag extender for DNG tags.
 *
 * The current version of libtiff (4.0.3) only supports tags for DNG
 * versions <= 1.1.0.0
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#include "x3f_dngtags.h"

#include <stddef.h>
#include <tiffio.h>

static TIFFExtendProc previous_tag_extender = NULL;

static void tag_extender(TIFF *tiff)
{
  static const TIFFFieldInfo tags[] = {
    {TIFFTAG_EXTRACAMERAPROFILES, -1, -1, TIFF_LONG, FIELD_CUSTOM, 1, 1,
     "ExtraCameraProfiles"},
    {TIFFTAG_PROFILENAME, -1, -1, TIFF_ASCII, FIELD_CUSTOM, 1, 0,
     "ProfileName"},
    {TIFFTAG_ASSHOTPROFILENAME, -1, -1, TIFF_ASCII, FIELD_CUSTOM, 1, 0,
     "AsShotProfileName"},
    {TIFFTAG_OPCODELIST1, -1, -1, TIFF_UNDEFINED, FIELD_CUSTOM, 1, 1,
     "OpcodeList1"},
    {TIFFTAG_OPCODELIST2, -1, -1, TIFF_UNDEFINED, FIELD_CUSTOM, 1, 1,
     "OpcodeList2"},
    {TIFFTAG_OPCODELIST3, -1, -1, TIFF_UNDEFINED, FIELD_CUSTOM, 1, 1,
     "OpcodeList3"},
    {TIFFTAG_DEFAULTBLACKRENDER, 1, 1, TIFF_LONG, FIELD_CUSTOM, 1, 0,
     "DefaultBlackRender"},
    {TIFFTAG_FORWARDMATRIX1, -1, -1, TIFF_SRATIONAL, FIELD_CUSTOM, 1, 1,
     "ForwardMatrix1"},
    {TIFFTAG_FORWARDMATRIX2, -1, -1, TIFF_SRATIONAL, FIELD_CUSTOM, 1, 1,
     "ForwardMatrix2"},
  };

  TIFFMergeFieldInfo(tiff, tags, sizeof(tags)/sizeof(tags[0]));

  if (previous_tag_extender)
    previous_tag_extender(tiff);
}

__attribute__((constructor)) static void install_tag_extender(void)
{
  previous_tag_extender = TIFFSetTagExtender(tag_extender);
}


<|start_filename|>src/x3f_dump.h<|end_filename|>
/* X3F_DUMP.H
 *
 * Library for writing unprocessed RAW and JPEG data.
 *
 * Copyright 2015 - <NAME> <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#ifndef X3F_DUMP_H
#define X3F_DUMP_H

#include "x3f_io.h"

extern x3f_return_t x3f_dump_raw_data(x3f_t *x3f, char *outfilename);
extern x3f_return_t x3f_dump_jpeg(x3f_t *x3f, char *outfilename);

#endif


<|start_filename|>src/x3f_printf.c<|end_filename|>
/* X3F_PRINTF.C
 *
 * Library for printout, depending on seriousity level.
 *
 * Copyright 2015 - <NAME> <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#include "x3f_printf.h"

#include <stdio.h>
#include <stdarg.h>

x3f_verbosity_t x3f_printf_level = INFO;

extern void x3f_printf(x3f_verbosity_t level, const char *fmt, ...)
{
  va_list ap;
  FILE *f = level > WARN ? stdout : stderr;

  if (level > x3f_printf_level) return;

  switch(level) {
  case ERR:
    fprintf(f, "ERR: ");
    break;
  case WARN:
    fprintf(f, "WRN: ");
    break;
  case INFO:
    fprintf(f, "   : ");
    break;
  case DEBUG:
    fprintf(f, "dbg: ");
    break;
  }
  va_start(ap, fmt);
  vfprintf(f, fmt, ap);
  va_end(ap);
}


<|start_filename|>src/x3f_version.h<|end_filename|>
/* X3F_VERSION.H
 *
 * Contains the version of X3F Tools
 *
 * Copyright 2016 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#ifndef X3F_VERSION_H
#define X3F_VERSION_H

extern char *version;

#endif


<|start_filename|>src/x3f_printf.h<|end_filename|>
/* X3F_PRINTF.H
 *
 * Library for printout, depending on seriousity level.
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#ifndef X3F_PRINTF_H
#define X3F_PRINTF_H

#ifdef __cplusplus
extern "C" {
#endif

typedef enum {ERR=0, WARN=1, INFO=2, DEBUG=3} x3f_verbosity_t;

extern x3f_verbosity_t x3f_printf_level;

extern void x3f_printf(x3f_verbosity_t level, const char *fmt, ...)
  __attribute__((format(printf, 2, 3)));

#ifdef __cplusplus
}
#endif

#endif


<|start_filename|>src/x3f_print_meta.h<|end_filename|>
/* X3F_PRINT_META.H
 *
 * Library for printing meta data found in X3F files.
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#ifndef X3F_PRINT_META_H
#define X3F_PRINT_META_H

#include "x3f_io.h"

extern uint32_t max_printed_matrix_elements;

extern void x3f_print_meta(x3f_t *x3f);
extern x3f_return_t x3f_dump_meta_data(x3f_t *x3f, char *outfilename);

#endif


<|start_filename|>src/x3f_extract.c<|end_filename|>
/* X3F_EXTRACT.C
 *
 * Extracting images from X3F files.
 * Also converting to histogram and color image.
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#include "x3f_version.h"
#include "x3f_io.h"
#include "x3f_process.h"
#include "x3f_output_dng.h"
#include "x3f_output_tiff.h"
#include "x3f_output_ppm.h"
#include "x3f_histogram.h"
#include "x3f_print_meta.h"
#include "x3f_dump.h"
#include "x3f_denoise.h"
#include "x3f_printf.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

typedef enum
  { META      = 0,
    JPEG      = 1,
    RAW       = 2,
    TIFF      = 3,
    DNG       = 4,
    PPMP3     = 5,
    PPMP6     = 6,
    HISTOGRAM = 7}
  output_file_type_t;

static char *extension[] =
  { ".meta",
    ".jpg",
    ".raw",
    ".tif",
    ".dng",
    ".ppm",
    ".ppm",
    ".csv" };

static void usage(char *progname)
{
  fprintf(stderr,
          "usage: %s <SWITCHES> <file1> ...\n"
          "   -o <DIR>        Use <DIR> as output directory\n"
          "   -v              Verbose output for debugging\n"
          "   -q              Suppress all messages except errors\n"
	  "ONE OFF THE FORMAT SWITCHWES\n"
	  "   -meta           Dump metadata\n"
          "   -jpg            Dump embedded JPEG\n"
          "   -raw            Dump RAW area undecoded\n"
          "   -tiff           Dump RAW/color as 3x16 bit TIFF\n"
          "   -dng            Dump RAW as DNG LinearRaw (default)\n"
          "   -ppm-ascii      Dump RAW/color as 3x16 bit PPM/P3 (ascii)\n"
          "                   NOTE: 16 bit PPM/P3 is not generally supported\n"
          "   -ppm            Dump RAW/color as 3x16 bit PPM/P6 (binary)\n"
          "   -histogram      Dump histogram as csv file\n"
          "   -loghist        Dump histogram as csv file, with log exposure\n"
	  "APPROPRIATE COMBINATIONS OF MODIFIER SWITCHES\n"
	  "   -color <COLOR>  Convert to RGB color space\n"
	  "                   (none, sRGB, AdobeRGB, ProPhotoRGB)\n"
	  "                   'none' means neither scaling, applying gamma\n"
	  "                   nor converting color space.\n"
	  "                   This switch does not affect DNG output\n"
          "   -unprocessed    Dump RAW without any preprocessing\n"
          "   -qtop           Dump Quattro top layer without preprocessing\n"
          "   -no-crop        Do not crop to active area\n"
          "   -no-denoise     Do not denoise RAW data\n"
          "   -no-sgain       Do not apply spatial gain (color compensation)\n"
          "   -no-fix-bad     Do not fix bad pixels\n"
          "   -sgain          Apply spatial gain (default except for Quattro)\n"
          "   -wb <WB>        Select white balance preset\n"
          "   -compress       Enable ZIP compression for DNG and TIFF output\n"
          "   -ocl            Use OpenCL\n"
	  "\n"
	  "STRANGE STUFF\n"
          "   -offset <OFF>   Offset for SD14 and older\n"
          "                   NOTE: If not given, then offset is automatic\n"
          "   -matrixmax <M>  Max num matrix elements in metadata (def=100)\n",
          progname);
  exit(1);
}

#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <errno.h>

static int check_dir(char *Path)
{
  struct stat filestat;
  int ret;

  if ((ret = stat(Path, &filestat)) != 0)
    return ret;

  if (!S_ISDIR(filestat.st_mode)) {
    errno = ENOTDIR;
    return -1;
  }

  return 0;
}

#define MAXPATH 1000
#define EXTMAX 10
#define MAXOUTPATH (MAXPATH+EXTMAX)
#define MAXTMPPATH (MAXOUTPATH+EXTMAX)

static int safecpy(char *dst, const char *src, int dst_size)
{
  if (strnlen(src, dst_size+1) > dst_size) {
    x3f_printf(DEBUG, "safecpy: String too large\n");
    return 1;
  } else {
    strcpy(dst, src);
    return 0;
  }
}

static int safecat(char *dst, const char *src, int dst_size)
{
  if (strnlen(dst, dst_size+1) + strnlen(src, dst_size+1) > dst_size) {
    x3f_printf(DEBUG, "safecat: String too large\n");
    return 1;
  } else {
    strcat(dst, src);
    return 0;
  }
}

#if defined(_WIN32) || defined(_WIN64)
#define PATHSEP "\\"
static const char pathseps[] = PATHSEP "/:";
#else
#define PATHSEP "/"
static const char pathseps[] = PATHSEP;
#endif

static int make_paths(const char *inpath, const char *outdir,
		      const char *ext,
		      char *tmppath, char *outpath)
{
  int err = 0;

  if (outdir && *outdir) {
    const char *ptr = inpath, *sep, *p;

    for (sep=pathseps; *sep; sep++)
      if ((p = strrchr(inpath, *sep)) && p+1 > ptr) ptr = p+1;

    err += safecpy(outpath, outdir, MAXOUTPATH);
    if (!strchr(pathseps, outdir[strlen(outdir)-1]))
      err += safecat(outpath, PATHSEP, MAXOUTPATH);
    err += safecat(outpath, ptr, MAXOUTPATH);
  }
  else err += safecpy(outpath, inpath, MAXOUTPATH);

  err += safecat(outpath, ext, MAXOUTPATH);
  err += safecpy(tmppath, outpath, MAXTMPPATH);
  err += safecat(tmppath, ".tmp", MAXTMPPATH);

  return err;
}

#define Z extract_jpg=0,extract_raw=0,extract_unconverted_raw=0

int main(int argc, char *argv[])
{
  int extract_jpg = 0;
  int extract_meta; /* Always computed */
  int extract_raw = 1;
  int extract_unconverted_raw = 0;
  int crop = 1;
  int fix_bad = 1;
  int denoise = 1;
  int apply_sgain = -1;
  output_file_type_t file_type = DNG;
  x3f_color_encoding_t color_encoding = SRGB;
  int files = 0;
  int errors = 0;
  int log_hist = 0;
  char *wb = NULL;
  int compress = 0;
  int use_opencl = 0;
  char *outdir = NULL;
  x3f_return_t ret;

  int i;

  x3f_printf(INFO, "X3F TOOLS VERSION = %s\n\n", version);

  /* Set stdout and stderr to line buffered mode to avoid scrambling */
  setvbuf(stdout, NULL, _IOLBF, 0);
  setvbuf(stderr, NULL, _IOLBF, 0);

  for (i=1; i<argc; i++)

    /* Only one of those switches is valid, the last one */
    if (!strcmp(argv[i], "-jpg"))
      Z, extract_jpg = 1, file_type = JPEG;
    else if (!strcmp(argv[i], "-meta"))
      Z, file_type = META;
    else if (!strcmp(argv[i], "-raw"))
      Z, extract_unconverted_raw = 1, file_type = RAW;
    else if (!strcmp(argv[i], "-tiff"))
      Z, extract_raw = 1, file_type = TIFF;
    else if (!strcmp(argv[i], "-dng"))
      Z, extract_raw = 1, file_type = DNG;
    else if (!strcmp(argv[i], "-ppm-ascii"))
      Z, extract_raw = 1, file_type = PPMP3;
    else if (!strcmp(argv[i], "-ppm"))
      Z, extract_raw = 1, file_type = PPMP6;
    else if (!strcmp(argv[i], "-histogram"))
      Z, extract_raw = 1, file_type = HISTOGRAM;
    else if (!strcmp(argv[i], "-loghist"))
      Z, extract_raw = 1, file_type = HISTOGRAM, log_hist = 1;

    else if (!strcmp(argv[i], "-color") && (i+1)<argc) {
      char *encoding = argv[++i];
      if (!strcmp(encoding, "none"))
	color_encoding = NONE;
      else if (!strcmp(encoding, "sRGB"))
	color_encoding = SRGB;
      else if (!strcmp(encoding, "AdobeRGB"))
	color_encoding = ARGB;
      else if (!strcmp(encoding, "ProPhotoRGB"))
	color_encoding = PPRGB;
      else {
	fprintf(stderr, "Unknown color encoding: %s\n", encoding);
	usage(argv[0]);
      }
    }
    else if (!strcmp(argv[i], "-o") && (i+1)<argc)
      outdir = argv[++i];
    else if (!strcmp(argv[i], "-v"))
      x3f_printf_level = DEBUG;
    else if (!strcmp(argv[i], "-q"))
      x3f_printf_level = ERR;
    else if (!strcmp(argv[i], "-unprocessed"))
      color_encoding = UNPROCESSED;
    else if (!strcmp(argv[i], "-qtop"))
      color_encoding = QTOP;
    else if (!strcmp(argv[i], "-no-crop"))
      crop = 0;
    else if (!strcmp(argv[i], "-no-fix-bad"))
      fix_bad = 0;
    else if (!strcmp(argv[i], "-no-denoise"))
      denoise = 0;
    else if (!strcmp(argv[i], "-no-sgain"))
      apply_sgain = 0;
    else if (!strcmp(argv[i], "-sgain"))
      apply_sgain = 1;
    else if ((!strcmp(argv[i], "-wb")) && (i+1)<argc)
      wb = argv[++i];
    else if (!strcmp(argv[i], "-compress"))
      compress = 1;
    else if (!strcmp(argv[i], "-ocl"))
      use_opencl = 1;

  /* Strange Stuff */
    else if ((!strcmp(argv[i], "-offset")) && (i+1)<argc)
      legacy_offset = atoi(argv[++i]), auto_legacy_offset = 0;
    else if ((!strcmp(argv[i], "-matrixmax")) && (i+1)<argc)
      max_printed_matrix_elements = atoi(argv[++i]);
    else if (!strncmp(argv[i], "-", 1))
      usage(argv[0]);
    else
      break;			/* Here starts list of files */

  if (outdir != NULL && check_dir(outdir) != 0) {
    x3f_printf(ERR, "Could not find outdir %s\n", outdir);
    usage(argv[0]);
  }

  x3f_set_use_opencl(use_opencl);

  extract_meta =
    file_type == META ||
    file_type == DNG ||
    (extract_raw &&
     (crop || (color_encoding != UNPROCESSED && color_encoding != QTOP)));

  for (; i<argc; i++) {
    char *infile = argv[i];
    FILE *f_in = fopen(infile, "rb");
    x3f_t *x3f = NULL;

    char tmpfile[MAXTMPPATH+1];
    char outfile[MAXOUTPATH+1];
    x3f_return_t ret_dump;
    int sgain;

    files++;

    if (f_in == NULL) {
      x3f_printf(ERR, "Could not open infile %s\n", infile);
      goto found_error;
    }

    x3f_printf(INFO, "READ THE X3F FILE %s\n", infile);
    x3f = x3f_new_from_file(f_in);

    if (x3f == NULL) {
      x3f_printf(ERR, "Could not read infile %s\n", infile);
      goto found_error;
    }

    if (extract_jpg) {
      if (X3F_OK != (ret = x3f_load_data(x3f, x3f_get_thumb_jpeg(x3f)))) {
	x3f_printf(ERR, "Could not load JPEG thumbnail from %s (%s)\n",
		   infile, x3f_err(ret));
	goto found_error;
      }
    }

    if (extract_meta) {
      x3f_directory_entry_t *DE = x3f_get_prop(x3f);

      if (X3F_OK != (ret = x3f_load_data(x3f, x3f_get_camf(x3f)))) {
	x3f_printf(ERR, "Could not load CAMF from %s (%s)\n",
		   infile, x3f_err(ret));
	goto found_error;
      }
      if (DE != NULL)
	/* Not for Quattro */
	if (X3F_OK != (ret = x3f_load_data(x3f, DE))) {
	  x3f_printf(ERR, "Could not load PROP from %s (%s)\n",
		     infile, x3f_err(ret));
	  goto found_error;
	}
      /* We do not load any JPEG meta data */
    }

    if (extract_raw) {
      x3f_directory_entry_t *DE;

      if (NULL == (DE = x3f_get_raw(x3f))) {
	x3f_printf(ERR, "Could not find any matching RAW format\n");
	goto found_error;
      }

      if (X3F_OK != (ret = x3f_load_data(x3f, DE))) {
	x3f_printf(ERR, "Could not load RAW from %s (%s)\n",
		   infile, x3f_err(ret));
	goto found_error;
      }
    }

    if (extract_unconverted_raw) {
      x3f_directory_entry_t *DE;

      if (NULL == (DE = x3f_get_raw(x3f))) {
	x3f_printf(ERR, "Could not find any matching RAW format\n");
	goto found_error;
      }

      if (X3F_OK != (ret = x3f_load_image_block(x3f, DE))) {
	x3f_printf(ERR, "Could not load unconverted RAW from %s (%s)\n",
		   infile, x3f_err(ret));
	goto found_error;
      }
    }

    if (make_paths(infile, outdir, extension[file_type], tmpfile, outfile)) {
      x3f_printf(ERR, "Too large outfile path for infile %s and outdir %s\n",
		 infile, outdir);
      goto found_error;
    }

    unlink(tmpfile);

    /* TODO: Quattro files seem to be already corrected for spatial
       gain. Is that assumption correct? Applying it only worsens the
       result anyhow, so it is disabled by default. */
    sgain =
      apply_sgain == -1 ? x3f->header.version < X3F_VERSION_4_0 : apply_sgain;

    switch (file_type) {
    case META:
      x3f_printf(INFO, "Dump META DATA to %s\n", outfile);
      ret_dump = x3f_dump_meta_data(x3f, tmpfile);
      break;
    case JPEG:
      x3f_printf(INFO, "Dump JPEG to %s\n", outfile);
      ret_dump = x3f_dump_jpeg(x3f, tmpfile);
      break;
    case RAW:
      x3f_printf(INFO, "Dump RAW block to %s\n", outfile);
      ret_dump = x3f_dump_raw_data(x3f, tmpfile);
      break;
    case TIFF:
      x3f_printf(INFO, "Dump RAW as TIFF to %s\n", outfile);
      ret_dump = x3f_dump_raw_data_as_tiff(x3f, tmpfile,
					   color_encoding,
					   crop, fix_bad, denoise, sgain, wb,
					   compress);
      break;
    case DNG:
      x3f_printf(INFO, "Dump RAW as DNG to %s\n", outfile);
      ret_dump = x3f_dump_raw_data_as_dng(x3f, tmpfile,
					  fix_bad, denoise, sgain, wb,
					  compress);
      break;
    case PPMP3:
    case PPMP6:
      x3f_printf(INFO, "Dump RAW as PPM to %s\n", outfile);
      ret_dump = x3f_dump_raw_data_as_ppm(x3f, tmpfile,
					  color_encoding,
					  crop, fix_bad, denoise, sgain, wb,
					  file_type == PPMP6);
      break;
    case HISTOGRAM:
      x3f_printf(INFO, "Dump RAW as CSV histogram to %s\n", outfile);
      ret_dump = x3f_dump_raw_data_as_histogram(x3f, tmpfile,
						color_encoding,
						crop, fix_bad, denoise, sgain, wb,
						log_hist);
      break;
    }

    if (X3F_OK != ret_dump) {
      x3f_printf(ERR, "Could not dump to %s: %s\n", tmpfile, x3f_err(ret_dump));
      errors++;
    } else {
      if (rename(tmpfile, outfile) != 0) {
	x3f_printf(ERR, "Could not rename %s to %s\n", tmpfile, outfile);
	errors++;
      }
    }

    goto clean_up;

  found_error:

    errors++;

  clean_up:

    x3f_delete(x3f);

    if (f_in != NULL)
      fclose(f_in);
  }

  if (files == 0) {
    x3f_printf(ERR, "No files given\n");
    usage(argv[0]);
  }

  x3f_printf(INFO, "Files processed: %d\terrors: %d\n", files, errors);

  return errors > 0;
}


<|start_filename|>src/x3f_image.h<|end_filename|>
/* X3F_IMAGE.H
 *
 * Library for access to image data.
 *
 * Copyright 2015 - Roland and <NAME>
 * BSD-style - see doc/copyright.txt
 *
 */

#ifndef X3F_IMAGE_H
#define X3F_IMAGE_H

#include "x3f_io.h"

typedef enum {COL_SIDE_LEFT = 0,
	      COL_SIDE_RIGHT = 1,
	      COL_SIDE_WRONG = 2} col_side_t;

extern int x3f_image_area(x3f_t *x3f, x3f_area16_t *image);
extern int x3f_image_area_qtop(x3f_t *x3f, x3f_area16_t *image);
extern int x3f_crop_area(uint32_t *coord, x3f_area16_t *image,
			 x3f_area16_t *crop);
extern int x3f_crop_area8(uint32_t *coord, x3f_area8_t *image,
			  x3f_area8_t *crop);
extern int x3f_get_camf_rect(x3f_t *x3f, char *name,
			     x3f_area16_t *image, int rescale,
			     uint32_t *rect);
extern int x3f_crop_area_column(x3f_t *x3f, col_side_t which_side,
				x3f_area16_t *image, int rescale,
				x3f_area16_t *crop);
extern int x3f_crop_area_camf(x3f_t *x3f, char *name,
			      x3f_area16_t *image, int rescale,
			      x3f_area16_t *crop);
extern int x3f_crop_area8_camf(x3f_t *x3f, char *name,
			       x3f_area8_t *image, int rescale,
			       x3f_area8_t *crop);

#endif
