Rcpp 経由で miniCBF 画像を読み込む

Rcpp と CBFlib の練習も兼ねてやってみた。エラー処理は無視!

cbflib-wrapper.cpp:

#include <Rcpp.h>
#include <cbf/cbf.h>

using namespace Rcpp;
using namespace std;

// [[Rcpp::export]]
NumericMatrix readCBF(string filename) {
  FILE *fh = fopen(filename.c_str(), "rb");
  cbf_handle minicbf;
  
  cbf_make_handle(&minicbf);
  cbf_read_widefile(minicbf, fh, MSG_DIGEST);
  
  cbf_find_category (minicbf, "array_data");
  //cbf_find_column   (minicbf, "array_id");
  //cbf_find_row      (minicbf, "null");
  cbf_find_column   (minicbf, "data");
  
  unsigned int compression;
  int binary_id = 0, elsigned = 0, elunsigned = 0, minelement = 0, maxelement = 0;
  size_t elsize = 0, elements = 0, dim1 = 0, dim2 = 0, dim3 = 0, padding = 0;
  const char *byteorder = NULL;
  cbf_get_integerarrayparameters_wdims(minicbf, &compression, &binary_id,
                                      &elsize, &elsigned, &elunsigned, &elements, &minelement, &maxelement,(const char **) &byteorder,
                                      &dim1, &dim2, &dim3, &padding);
  // printf("wdin1 = %ld wdim2 = %ld wdim3 = %ld elsize = %ld elsigned = %d padding = %ld\n", dim1, dim2, dim3, elsize, elsigned, padding);
  
  int *buf = (int*)calloc(elements, sizeof(int)); 
  size_t elements_read;
  cbf_get_integerarray(minicbf, &binary_id, buf, elsize, elsigned, elements, &elements_read);
  cbf_free_handle(minicbf);

  NumericMatrix mat(dim1, dim2); // dim1 is fast scan, which is continuous in memory
  for (int i = 0; i < dim2; i++) {
    for (int j = 0; j < dim1; j++) {
      mat(j, i) = *buf++; // R matrix is column major!
    }
  }
  
  return mat;
}

Rコード:

library(Rcpp)

scale_img <- function(img, min, max) {
  img <- (img - min) / (max - min)
  img <- ifelse(img < 0, 0, img)
  ifelse(img > 1, 1, img)
}

Sys.setenv("PKG_LIBS"="-lcbf")
Rcpp::sourceCpp("cbflib-wrapper.cpp")
minicbf <- readCBF("test.cbf")
image(scale_img(minicbf, 0, 25), useRaster=TRUE, col=grey(seq(0, 1, length.out=256)))

CBFlib へのリンクを、PKG_LIBS 環境変数によって指定するのがポイント。

CIF ファイルは、まず datablock があり (data_XXXX)、その中にカテゴリ (_array_data.XXXX) があって、その中に key-value ペアとしてデータが入っている構成。loop を使っている場合 (mmCIF とか) は、まず row を選択してから、column を選択してデータにアクセスする。loop がない場合は、row は 1 つしかないので、わざわざ選択しなくてもよい。