Disagreement between miget_real_value and miget_real_value_hyperslab

[cfhammill]

Typically in RMINC we read data in using miget_real_value_hyperslab. However one of our functions, mincGetVoxel uses miget_real_value. I would have assumed that these two functions would agree on any particular voxel's real value, but that doesn't seem to be the case at least for double precision volumes. Internally miget_real_value uses mirw_hyperslab_raw which does no conversion on the value, and then post-hoc calls miconvert_voxel_to_real. Whereas miget_real_value_hyperslab calls mirw_hyperslab_icv which does its own, apparently more complicated scaling internally. So I suspect that miconvert_voxel_to_real is missing something relevant for the rescaling.

Here's a sample C program illustrating the problem

test-miget-real-voxel.c

#include <minc2.h>
#include <stdio.h>
int main(int argc, char** argv){
  char* file=argv[1];
  misize_t location[3];
  misize_t count[3];
  mihandle_t hvol;
  int result;
  int i;
  double single;
  double slab;
  for(i = 0; i < 3; i++){
    location[i] = 50;
    count[i] = 1;
  }
  result = miopen_volume(file,
                         MI2_OPEN_READ, &hvol);
  if (result != MI_NOERROR) {
    printf("Error opening input file: %s.\n", file);
    return 1;
  }
  result = miget_real_value(hvol, location, 3, &single);
  printf("Getting a single value returns: %f.\n", single);
  result = miget_real_value_hyperslab(hvol, MI_TYPE_DOUBLE, location, count, &slab);
  printf("Getting a size 1 slab returns: %f.\n", slab);
  return 0;
}

An easy way to test is to convert the icbm average to unsigned double

mincresample -2 -double mni_icbm152_t1_tal_nlin_sym_09a.mnc mni_icbm152_t1_tal_nlin_sym_09a_signed_double.mnc

Then compile and run the above program to check the difference, on my machine I get

Getting a single value returns: 15.561608.
Getting a size 1 slab returns: 17.435906.

This was noticed by one of our users getting the expected answer when using mincAnova which uses miget_real_value_hyperslab internally, but a very different answer when plucking individual voxels with mincGetVoxel which uses miget_real_value

[vfonov]

miget_real_value calls miget_voxel_value which calls miget_voxel_value_hyperslab
And then applies slice scaling. Perhaps the minc volume in question has wrong slice scaling set?

[cfhammill]

Does miget_hyperslab_icv not implement slice scaling? Because if it does it's not through miconvert_voxel_to_real which I thought probably accounts for the difference. I hope the mincresample call above doesn't screw up slice scaling.

[vfonov]

mincresample uses minc1 API, miget_real_value_hyperslab and miget_real_value is from minc2 API

[cfhammill]

So do you think mincresample is screwing up the slice scaling causing the disagreement?

[vfonov]

what's the file in question?

[cfhammill]

mni_icbm152_t1_tal_nlin_sym_09a.mnc

[vfonov]

so, looks like the problem is in miget_real_value
miget_real_value_hyperslab and miget_hyperslab_with_icv return expected result
$ mincextract -start 50,50,50 -count 1,1,1 mni_icbm152_t1_tal_nlin_sym_09a_signed_double.mnc
17.435906448874884944
$ ./test mni_icbm152_t1_tal_nlin_sym_09a_signed_double.mnc
miget_real_value: 15.561608.
miget_real_value_hyperslab: 17.435906.
miget_hyperslab_with_icv: 17.435906.

[vfonov]

result = miget_hyperslab_with_icv(hvol, MI_TYPE_DOUBLE, location, count, &slab_icv); printf("miget_hyperslab_with_icv: %f.\n", slab_icv);

[vfonov]

so, there is piece of code in https://github.com/BIC-MNI/libminc/blob/develop/libsrc2/hyper.c#L579 (used in miget_hyperslab_with_icv) that is missing in https://github.com/BIC-MNI/libminc/blob/develop/libsrc2/convert.c#L89 , used from miget_real_value

[cfhammill]

Ok, thanks for looking in to this Vlad!

So that code says that scaling is disabled if you have float volume? I guess that kind of makes sense.

[vfonov]

yes, it should be disabled for float and double , but miget_real_value doesn't do that

[cfhammill]

Want me to issue a PR or should I leave this with you?

[vfonov]

it will be better if you make the PR
I don't know if I will have much time working on this anytime soon

[cfhammill]

No problem, I'll open one later this week.

[gdevenyi]

Hi, we're planning a new libminc/minc-toolkit-v2 release, do you have time to put together a PR, or should we do it?