Image Processing in C++/Go

transform.cc

//===--------------------------------------------
// File: transform.cc
// Usage: display a given image from a specific
//        rgb color space and simulate the image
//        from another color rgb space.
// Author: Xing GUO <[email protected]>
//===-------------------------------------------

#include <iostream>
#include <cmath>
#include <string>
#include "bitmap_image.h"

typedef struct {
  double xr;
  double yr;
  double xg;
  double yg;
  double xb;
  double yb;
  double xw;
  double yw;
} Coefficient;

const Coefficient sRGB = {
  /* xr= */ 0.64,   /* yr= */ 0.33,
  /* xg= */ 0.29,   /* yg= */ 0.60,
  /* xb= */ 0.15,   /* yb= */ 0.06,
  /* xw= */ 0.3127, /* yw= */ 0.3291,
};

const Coefficient screenRGB = {
  /* xr= */ 0.507,  /* yr= */ 0.312,
  /* xg= */ 0.315,  /* yg= */ 0.537,
  /* xb= */ 0.158,  /* yb= */ 0.118,
  /* xw= */ 0.279,  /* yw= */ 0.312,
};

/* helper functions */
void matrix_inv_3x3(
    const double matrix_A[3][3],
    double       inv_matrix_A[3][3]
    ) {
  double det = 0.0;
  for (unsigned i = 0; i < 3; ++ i)
    det = det +
      (matrix_A[0][i] *
        (matrix_A[1][(i+1)%3] * matrix_A[2][(i+2)%3] -
         matrix_A[1][(i+2)%3] * matrix_A[2][(i+1)%3]));
  for (unsigned i = 0; i < 3; ++ i)
    for (unsigned j = 0; j < 3; ++ j)
      inv_matrix_A[i][j] =
        ((matrix_A[(j+1)%3][(i+1)%3] * matrix_A[(j+2)%3][(i+2)%3]) -
         (matrix_A[(j+1)%3][(i+2)%3] * matrix_A[(j+2)%3][(i+1)%3])) / det;
}

void matrix_mul_3x3(
    const double matrix_A[3][3],
    const double matrix_B[3][3],
    double       matrix_C[3][3]
    ) {
  for (unsigned i = 0; i < 3; ++ i)
    for (unsigned j = 0; j < 3; ++ j)
      matrix_C[i][j] =
        matrix_A[i][0] * matrix_B[0][j] +
        matrix_A[i][1] * matrix_B[1][j] +
        matrix_A[i][2] * matrix_B[2][j];
}

void matrix_vec_mul_3(
    const double matrix_A[3][3],
    const double vec_b[3],
    double       vec_c[3]
    ) {
  for (unsigned i = 0; i < 3; ++ i)
    vec_c[i] =
      matrix_A[i][0] * vec_b[0] +
      matrix_A[i][1] * vec_b[1] +
      matrix_A[i][2] * vec_b[2];
}

/* See: http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html */
void get_rgb2xyz_matrix(
    const Coefficient coeff,
    double rgb2xyz_m[3][3]
    ) {
  double XYZ_w[3] = {
    coeff.xw,
    coeff.yw,
    1.0 - coeff.xw - coeff.yw };
  double Xr = coeff.xr / coeff.yr,
         Yr = 1.0,
         Zr = (1.0 - coeff.xr - coeff.yr) / coeff.yr,
         Xg = coeff.xg / coeff.yg,
         Yg = 1.0,
         Zg = (1 - coeff.xg - coeff.yg) / coeff.yg,
         Xb = coeff.xb / coeff.yb,
         Yb = 1.0,
         Zb = (1.0 - coeff.xb - coeff.yb) / coeff.yb;
  double XYZ_rgb_m[3][3] = {
    { Xr, Xg, Xb },
    { Yr, Yg, Yb },
    { Zr, Zg, Zb }
  };
  double inv_XYZ_rgb_m[3][3];
  matrix_inv_3x3(XYZ_rgb_m, inv_XYZ_rgb_m);
  double S_rgb[3];
  matrix_vec_mul_3(inv_XYZ_rgb_m, XYZ_w, S_rgb);
  rgb2xyz_m[0][0] = S_rgb[0] * Xr;
  rgb2xyz_m[0][1] = S_rgb[1] * Xg;
  rgb2xyz_m[0][2] = S_rgb[2] * Xb;
  rgb2xyz_m[1][0] = S_rgb[0] * Yr;
  rgb2xyz_m[1][1] = S_rgb[1] * Yg;
  rgb2xyz_m[1][2] = S_rgb[2] * Yb;
  rgb2xyz_m[2][0] = S_rgb[0] * Zr;
  rgb2xyz_m[2][1] = S_rgb[1] * Zg;
  rgb2xyz_m[2][2] = S_rgb[2] * Zb;
}

bitmap_image transform_im(
    const Coefficient from_coeff,
    const Coefficient to_coeff,
    bitmap_image image_i) {
  double from_coeff_mat[3][3], inv_from_coeff_mat[3][3],
         to_coeff_mat[3][3];
  get_rgb2xyz_matrix(from_coeff, from_coeff_mat);
  get_rgb2xyz_matrix(to_coeff, to_coeff_mat);
  matrix_inv_3x3(from_coeff_mat, inv_from_coeff_mat);
  double transform_matrix[3][3];
  matrix_mul_3x3(inv_from_coeff_mat, to_coeff_mat, transform_matrix);
  unsigned im_width  = image_i.width(),
           im_height = image_i.height();
  for (unsigned i = 0; i < im_width; ++ i) {
    for (unsigned j = 0; j < im_height; ++ j) {
      unsigned char rgb_uc[3];
      double rgb_float_t[3], rgb_float[3];
      image_i.get_pixel(i, j,
          rgb_uc[0], rgb_uc[1], rgb_uc[2]);
      for (unsigned i = 0; i < 3; ++ i)
        rgb_float_t[i] = static_cast<double>(rgb_uc[i]);
      matrix_vec_mul_3(transform_matrix, rgb_float_t, rgb_float);
      for (unsigned i = 0; i < 3; ++ i)
        rgb_uc[i] = static_cast<unsigned char>(rgb_float[i]);
      image_i.set_pixel(i, j, rgb_uc[0], rgb_uc[1], rgb_uc[2]);
    }
  }
  return image_i;
}

int main(int argc, char **argv) {
  if (argc < 2) {
    std::cerr << "ERROR: Too few arguments.\n"
              << "Usage:\n"
              << "  ./transform image.bmp";
    return 1;
  }
  std::string image_name(argv[1]);
  bitmap_image image_i(image_name);
  bitmap_image image_o = transform_im(sRGB, screenRGB, image_i);
  image_o.save_image(image_name + ".cc.o.bmp");
  return 0;
}

transform.go

//===--------------------------------------------
// File: transform.go
// Usage: display a given image from a specific
//        rgb color space and simulate the image
//        from another color rgb space.
// Author: Xing GUO <[email protected]>
//===-------------------------------------------

package main

import (
  "fmt"
  "os"
  "image"
  "image/color"
  "golang.org/x/image/bmp"
)

type Coefficient struct {
  xr float64
  yr float64
  xg float64
  yg float64
  xb float64
  yb float64
  xw float64
  yw float64
}

var sRGB = Coefficient{
  /* xr= */ 0.64,   /* yr= */ 0.33,
  /* xg= */ 0.29,   /* yg= */ 0.60,
  /* xb= */ 0.15,   /* yb= */ 0.06,
  /* xw= */ 0.3127, /* yw= */ 0.3291,
}

var screenRGB = Coefficient{
  /* xr= */ 0.507,  /* yr= */ 0.312,
  /* xg= */ 0.315,  /* yg= */ 0.537,
  /* xb= */ 0.158,  /* yb= */ 0.118,
  /* xw= */ 0.279,  /* yw= */ 0.312,
}

// help functions
func matrix_mul_3x3(matrix_a [3][3]float64,
    matrix_b [3][3]float64) (matrix_c [3][3]float64) {
  for i := 0; i < 3; i ++ {
    for j := 0; j < 3; j ++ {
      matrix_c[i][j] = matrix_a[i][0] * matrix_b[0][j] +
                       matrix_a[i][1] * matrix_b[1][j] +
                       matrix_a[i][2] * matrix_b[2][j]
    }
  }
  return
}

func matrix_vec_mul_3(matrix_a [3][3]float64,
    vec_b [3]float64) (vec_c [3]float64) {
  for i := 0; i < 3; i ++ {
    vec_c[i] = matrix_a[i][0] * vec_b[0] +
               matrix_a[i][1] * vec_b[1] +
               matrix_a[i][2] * vec_b[2]
  }
  return
}

func matrix_inv_3x3(matrix_a [3][3]float64) (inv_matrix_a [3][3]float64) {
  var det float64 = 0
  for i := 0; i < 3; i ++ {
    det = det +
      (matrix_a[0][i] *
        (matrix_a[1][(i+1)%3] * matrix_a[2][(i+2)%3] -
         matrix_a[1][(i+2)%3] * matrix_a[2][(i+1)%3]))
  }
  for i := 0; i < 3; i ++ {
    for j := 0; j < 3; j ++ {
      inv_matrix_a[i][j] =
        ((matrix_a[(j+1)%3][(i+1)%3] * matrix_a[(j+2)%3][(i+2)%3]) -
         (matrix_a[(j+1)%3][(i+2)%3] * matrix_a[(j+2)%3][(i+1)%3])) / det;
    }
  }
  return
}

func get_rgb2xyz_matrix(coeff Coefficient) (rgb2xyz_m [3][3]float64) {
  XYZ_w := [3]float64{
    coeff.xw,
    coeff.yw,
    1.0 - coeff.xw - coeff.yw,
  }
  var Xr float64 = coeff.xr / coeff.yr
  var Yr float64 = 1.0
  var Zr float64 = (1.0 - coeff.xr - coeff.yr) / coeff.yr
  var Xg float64 = coeff.xg / coeff.yg
  var Yg float64 = 1.0
  var Zg float64 = (1 - coeff.xg - coeff.yg) / coeff.yg
  var Xb float64 = coeff.xb / coeff.yb
  var Yb float64 = 1.0
  var Zb float64 = (1.0 - coeff.xb - coeff.yb) / coeff.yb
  XYZ_rgb_m := [3][3]float64{
    { Xr, Xg, Xb },
    { Yr, Yg, Yb },
    { Zr, Zg, Zb },
  }
  inv_XYZ_rgb_m := matrix_inv_3x3(XYZ_rgb_m)
  S_rgb := matrix_vec_mul_3(inv_XYZ_rgb_m, XYZ_w)
  rgb2xyz_m[0][0] = S_rgb[0] * Xr
  rgb2xyz_m[0][1] = S_rgb[1] * Xg
  rgb2xyz_m[0][2] = S_rgb[2] * Xb
  rgb2xyz_m[1][0] = S_rgb[0] * Yr
  rgb2xyz_m[1][1] = S_rgb[1] * Yg
  rgb2xyz_m[1][2] = S_rgb[2] * Yb
  rgb2xyz_m[2][0] = S_rgb[0] * Zr
  rgb2xyz_m[2][1] = S_rgb[1] * Zg
  rgb2xyz_m[2][2] = S_rgb[2] * Zb
  return
}

func transform_im(from_coeff Coefficient, to_coeff Coefficient,
                  im_inp image.Image) image.Image {
  from_coeff_mat := get_rgb2xyz_matrix(from_coeff)
  to_coeff_mat   := get_rgb2xyz_matrix(to_coeff)
  inv_from_coeff_mat := matrix_inv_3x3(from_coeff_mat)
  transform_matrix := matrix_mul_3x3(inv_from_coeff_mat, to_coeff_mat)
  im_bounds := im_inp.Bounds()
  im_outp := image.NewRGBA(
    image.Rect(im_bounds.Min.X, im_bounds.Min.Y, im_bounds.Max.X, im_bounds.Max.Y))
  for i := im_bounds.Min.X; i < im_bounds.Max.X; i ++ {
    for j := im_bounds.Min.Y; j < im_bounds.Max.Y; j ++ {
      r_uc, g_uc, b_uc, _ := im_inp.At(i, j).RGBA()
      rgb_float := [3]float64{
        float64(r_uc >> 8),
        float64(g_uc >> 8),
        float64(b_uc >> 8) }
      rgb_float = matrix_vec_mul_3(transform_matrix, rgb_float)
      rgba := color.RGBA{
        uint8(rgb_float[0]),
        uint8(rgb_float[1]),
        uint8(rgb_float[2]), 255 }
      im_outp.Set(i, j, rgba)
    }
  }
  return im_outp
}

func main() {
  argv := os.Args[1:]
  if len(argv) < 1 {
    fmt.Println("ERROR: Too few arguments.")
    fmt.Println("  Usage: ./transform image.bmp")
    return
  }
  image_file_i, err := os.Open(argv[0])
  if err != nil {
    fmt.Println(err)
    return
  }
  defer image_file_i.Close()
  image_data_i, err := bmp.Decode(image_file_i)
  image_outp := transform_im(sRGB, screenRGB, image_data_i)
  image_file_o, err := os.Create(argv[0]+".go.o.bmp")
  if err != nil {
    fmt.Println(err)
    return
  }
  err = bmp.Encode(image_file_o, image_outp)
  if err != nil {
    fmt.Println(err)
    return
  }
  defer image_file_o.Close()
}