WASM extended niimathx() calculates robust range

src/bwlabel.c

@@ -39,7 +39,7 @@
  ***************************************************************/
 
 void fill_tratab(uint32_t  *tt,     /* Translation table */
-                 uint32_t  ttn,     /* Size of translation table */
+                 /*uint32_t  ttn,*/     /* Size of translation table */
                  uint32_t  *nabo,   /* Set of neighbours */
                  uint32_t  nr_set)  /* Number of neighbours in nabo */
 {
@@ -83,8 +83,8 @@ uint32_t check_previous_slice(uint32_t  *il,     /* Initial labelling map */
                                   uint32_t  sl,      /* slice */
                                   size_t        dim[3],  /* dimensions of il */
                                   uint32_t  conn,    /* Connectivity criterion */
-                                  uint32_t  *tt,     /* Translation table */
-                                  uint32_t  ttn)     /* Size of translation table */
+                                  uint32_t  *tt)     /* Translation table */
+//                                  uint32_t  ttn)     /* Size of translation table */
 {
    uint32_t l=0;
    uint32_t nabo[9];
@@ -113,7 +113,7 @@ uint32_t check_previous_slice(uint32_t  *il,     /* Initial labelling map */
 
    if (nr_set) 
    {
-      fill_tratab(tt,ttn,nabo,nr_set);
+      fill_tratab(tt,/*ttn,*/nabo,nr_set);
       return(nabo[0]);
    }
    else {return(0);}
@@ -164,7 +164,7 @@ uint32_t do_initial_labelling(uint8_t        *bw,   /* Binary map */
             nr_set = 0;
             if (bw[idx(r,c,sl,dim)])
             {
-               nabo[0] = check_previous_slice(il,r,c,sl,dim,conn,*tt,ttn);
+               nabo[0] = check_previous_slice(il,r,c,sl,dim,conn,*tt /*,ttn*/);
                if (nabo[0]) {nr_set += 1;}
                /*
                   For six(surface)-connectivity
@@ -198,7 +198,7 @@ uint32_t do_initial_labelling(uint8_t        *bw,   /* Binary map */
                if (nr_set)
                {
                   il[idx(r,c,sl,dim)] = nabo[0];
-                  fill_tratab(*tt,ttn,nabo,nr_set);
+                  fill_tratab(*tt,/*ttn,*/nabo,nr_set);
                }
                else
                {
@@ -298,6 +298,10 @@ int bwlabel(nifti_image *nim, int conn) {
   uint32_t  *tt = NULL;
   uint32_t ttn = do_initial_labelling(bw,dim,conn,il,&tt);
   double nl = translate_labels(il,dim,tt,ttn,l);
+  nim->cal_min = 0;
+  nim->cal_max = nl;
+  nim->scl_inter = 0.0;
+  nim->scl_slope = 1.0;
   _mm_free(il);
   _mm_free(tt);
   if (nim->datatype == DT_FLOAT32) {

src/coreFLT.c

@@ -942,7 +942,7 @@ staticx int nifti_c2h(nifti_image *nim) {
 	flt kUninterestingDarkUnits = 900.0; //e.g. -1000..-100
 	flt kInterestingMidUnits = 200.0; //e.g. unenhanced CT: -100..+100
 	flt kScaleRatio = 10;
-	flt kMax = kInterestingMidUnits * (kScaleRatio+1);
+	//flt kMax = kInterestingMidUnits * (kScaleRatio+1);
 	if (nim->nvox < 1)
 		return 1;
 	flt mn = darkest_voxel(nim);
@@ -1008,7 +1008,6 @@ staticx int nifti_otsu(nifti_image *nim, int mode, int makeBinary) { //binarize
 	//Create histogram of intensity frequency
 	// hist[0..kOtsuBins-1]: each bin is number of pixels with this intensity
 	flt mn, mx;
-
 	if (nifti_robust_range(nim, &mn, &mx, 0) != 0)
 		return 1;
 	if (mn >= mx)
@@ -1144,31 +1143,29 @@ staticx int nifti_crop(nifti_image *nim, int tmin, int tsize) {
 	return 0;
 }
 
-staticx void nifti_add2(flt *v, size_t n, flt intercept1) {
+/*staticx void nifti_add2(flt *v, size_t n, flt intercept1) {
 	//#pragma omp parallel for
 	for (size_t i = 0; i < n; i++)
 		v[i] += intercept1;
-} //nifti_add()
+} //nifti_add()*/
 
 staticx int nifti_rescale(nifti_image *nim, double scale, double intercept) {
 	//linear transform of data
 	if (nim->nvox < 1)
 		return 1;
-	flt scl = scale;
-	flt inter = intercept;
 	flt *f32 = (flt *)nim->data;
 	if (intercept == 0.0) {
 		if (scale == 1.0)
 			return 0; //nothing to do
-		nifti_mul(f32, nim->nvox, scl);
+		nifti_mul(f32, nim->nvox, scale);
 		return 0;
 	} else if (scale == 1.0) {
 		nifti_add(f32, nim->nvox, intercept);
 		return 0;
 	}
-	nifti_fma(f32, nim->nvox, scl, inter);
+	nifti_fma(f32, nim->nvox, scale, intercept);
 	//for (size_t i = 0; i < nim->nvox; i++ )
-	//	f32[i] = (f32[i] * scl) + inter;
+	//	f32[i] = (f32[i] * scl) + intercept;
 	return 0;
 }
 
@@ -2394,7 +2391,7 @@ staticx int nifti_tensor_decomp(nifti_image *nim, int isUpperTriangle) {
 			in6[v] = in32[iv];
 			iv += nvox3D;
 		}
-		EIG_tsfunc(0.0, 0.0, 0, in6, 0.0, 0.0, NULL, 0, out14, isUpperTriangle);
+		EIG_tsfunc(0, in6, out14, isUpperTriangle);
 		size_t ov = i;
 		for (int v = 0; v < 14; v++) {
 			out32[ov] = out14[v];
@@ -3108,7 +3105,7 @@ staticx int nifti_roi(nifti_image *nim, int xmin, int xsize, int ymin, int ysize
 	return 0;
 }
 
-staticx int nifti_sobel(nifti_image *nim, int offc, int isBinary) {
+staticx int nifti_sobel(nifti_image *nim, int isBinary) {
 	//sobel is simply one kernel pass per dimension.
 	// this could be achieved with successive passes of "-kernel"
 	// here it is done in a single pass for cache efficiency
@@ -5124,9 +5121,9 @@ int mainWASM(nifti_image *nim, int argc, char *argv[]) {
 		else if (!strcmp(argv[ac], "-c2h"))
 			ok = nifti_c2h(nim);
 		else if (!strcmp(argv[ac], "-sobel_binary"))
-			ok = nifti_sobel(nim, 1, 1);
+			ok = nifti_sobel(nim, 1);
 		else if (!strcmp(argv[ac], "-sobel"))
-			ok = nifti_sobel(nim, 1, 0);
+			ok = nifti_sobel(nim, 0);
 		else if (!strcmp(argv[ac], "-demean"))
 			ok = nifti_demean(nim);
 		else if (!strcmp(argv[ac], "-detrend"))
@@ -5433,7 +5430,52 @@ float clampf(float d, float min, float max) {
 	return t > max ? max : t;
 }
 
-__attribute__((used)) int niimath (void *img, int datatype, int nx, int ny, int nz, int nt, float dx, float dy, float dz, float dt, char * cmdstr){
+staticx int nifti_unscale(nifti_image *nim, double scale, double intercept) {
+	//invert linear transform of data
+	if ((nim->nvox < 1) || (scale == 0.0))
+		return 1;
+	flt *f32 = (flt *)nim->data;
+	if (intercept == 0.0) {
+		if (scale == 1.0)
+			return 0; //nothing to do
+		nifti_mul(f32, nim->nvox, 1.0 / scale);
+		return 0;
+	} else if (scale == 1.0) {
+		nifti_add(f32, nim->nvox, - intercept);
+		return 0;
+	}
+	flt scl = 1.0 / scale;
+	for (size_t i = 0; i < nim->nvox; i++ )
+		f32[i] = (f32[i] - intercept ) * scl;
+	return 0;
+}
+
+int mainWASM2(nifti_image *nim, int argc, char *argv[], float *vars) {
+	nim->scl_slope = vars[0];
+	nim->scl_inter = vars[1];
+	nim->cal_min = vars[2];
+	nim->cal_max = vars[3];
+	if (nim->scl_slope != 0.0)
+		nifti_rescale(nim, nim->scl_slope, nim->scl_inter);
+	int ret = mainWASM(nim, argc, argv);
+	if (nim->scl_slope == 0.0) //bogus vars*, do not rescale, do not waste time calculating robust range
+		return ret;
+	nifti_unscale(nim, nim->scl_slope, nim->scl_inter);
+	flt mn, mx;
+	nifti_robust_range(nim, &mn, &mx, 0);
+	vars[0] = nim->scl_slope;
+	vars[1] = nim->scl_inter;
+	vars[2] = mn;
+	vars[3] = mx;
+	return ret;
+}
+
+//niimathx extends 'niimath' call to have float vars that may be changed by computations
+//vars[0] = scl_slope
+//vars[1] = scl_inter
+//vars[2] = cal_min
+//vars[3] = cal_max
+__attribute__((used)) int niimathx (void *img, float *vars, int datatype, int nx, int ny, int nz, int nt, float dx, float dy, float dz, float dt, char * cmdstr){
 	int nvox = nx * ny * nz * MAX(nt, 1);
 	if (nvox < 1) return 101;
 	#define argc_MAX 128
@@ -5459,14 +5501,14 @@ __attribute__((used)) int niimath (void *img, int datatype, int nx, int ny, int
 	nim.cal_min = 0.0;
 	nim.datatype = DT_FLOAT;
 	if (datatype == DT_FLOAT) 
-		return mainWASM(&nim, argc, argv);//niimath_core(&nim, cmdstr);
+		return mainWASM2(&nim, argc, argv, vars);//niimath_core(&nim, cmdstr);
 	if (datatype == DT_SIGNED_SHORT) {
 		int16_t *img16 = (int16_t *)img;
 		float *img32 = (float *) _mm_malloc(nvox * sizeof(float), 64);
 		for (int i = 0; i < nvox; ++i)
 			img32[i] = img16[i];
 		nim.data = img32;
-		int ret = mainWASM(&nim, argc, argv);
+		int ret = mainWASM2(&nim, argc, argv, vars);
 		for (int i = 0; i < nvox; ++i)
 			img16[i] = clampf(img32[i], -32768.0, 32767.0);//img32[i];
 		_mm_free(img32);
@@ -5478,12 +5520,17 @@ __attribute__((used)) int niimath (void *img, int datatype, int nx, int ny, int
 		for (int i = 0; i < nvox; ++i)
 			img32[i] = img8[i];
 		nim.data = img32;
-		int ret = mainWASM(&nim, argc, argv);
+		int ret = mainWASM2(&nim, argc, argv, vars);
 		for (int i = 0; i < nvox; ++i)
 			img8[i] = clampf(img32[i], 0.0, 255.0); //img32[i]; 
 		_mm_free(img32);
 		return ret;
 	}
 	return 88;
 }
+
+__attribute__((used)) int niimath (void *img, int datatype, int nx, int ny, int nz, int nt, float dx, float dy, float dz, float dt, char * cmdstr){
+	float vars[] = {0.0, 0.0, 0.0, 0.0};
+	return niimathx (img, &vars[0], datatype, nx, ny, nz, nt, dx, dy, dz, dt, cmdstr);
+}
 #endif //WASM code

src/niimath.js

@@ -4,8 +4,6 @@
 // emcc -O2 -s ALLOW_MEMORY_GROWTH -s MAXIMUM_MEMORY=4GB -s WASM=1 -DUSING_WASM -I. core32.c nifti2_wasm.c core.c walloc.c -o funcx.js
 //Test on image
 // node niimath.js T1.nii -sqr sT1.nii
-//Note: since we renew the ArrayBuffer we do not need to pre-allocate worst case TOTAL_MEMORY:
-// emcc -O2 -s ALLOW_MEMORY_GROWTH -s MAXIMUM_MEMORY=4GB -s TOTAL_MEMORY=268435456 -s WASM=1 -DUSING_WASM -I. core32.c nifti2_wasm.c core.c walloc.c -o funcx.js
 
 const fs = require('fs')
 const nifti = require('nifti-reader-js')
@@ -94,7 +92,7 @@ function niimath_shim(instance, memory, hdr, img8, cmd) {
       console.log('Only dataypes 2,4,16 supported')
       return
     }
-    let {niimath, walloc, wfree} = instance.exports
+    let {niimathx, niimath, walloc, wfree} = instance.exports
     let nx = hdr.dims[1] //number of columns
     let ny = hdr.dims[2] //number of rows
     let nz = hdr.dims[3] //number of slices
@@ -103,7 +101,7 @@ function niimath_shim(instance, memory, hdr, img8, cmd) {
     let dy = hdr.pixDims[2] //space between rows
     let dz = hdr.pixDims[3] //space between slices
     let dt = hdr.pixDims[4] //time between volumes
-    let bpv = Math.floor(hdr.numBitsPerVoxel/8) 
+    let bpv = Math.floor(hdr.numBitsPerVoxel/8)
     console.log(`dims: ${nx}x${ny}x${nz}x${nt} nbyper: ${bpv}`)
     //allocate WASM null-terminated command string
     let cptr = walloc(cmd.length + 1)
@@ -113,27 +111,43 @@ function niimath_shim(instance, memory, hdr, img8, cmd) {
         cmdstr[i] = cmd.charCodeAt(i)
     let cstr = new Uint8Array(instance.exports.memory.buffer, cptr, cmd.length + 1)
     cstr.set(cmdstr)
+    //allocate WASM float variables:
+    let fptr = walloc(4 * 4) //4 = sizeof(float)
+    memory.record_malloc(fptr, 4 * 4)
+    fvars = new Float32Array(instance.exports.memory.buffer, fptr, 4)
+    fvars[0] = hdr.scl_slope
+    fvars[1] = hdr.scl_inter
+    fvars[2] = hdr.cal_min
+    fvars[3] = hdr.cal_max
     //allocate WASM image data
-    let nvox = nx * ny * nz * nt; 
+    let nvox = nx * ny * nz * nt
     let ptr = walloc(nvox * bpv)
     memory.record_malloc(ptr, nvox * bpv)
     cimg = new Uint8Array(instance.exports.memory.buffer, ptr, nvox * bpv)
     cimg.set(img8)
     //run WASM
     startTime = new Date()
-    let ok = niimath(ptr, datatype, nx, ny, nz, nt, dx, dy, dz, dt, cptr)
+    //let ok = niimath(ptr, datatype, nx, ny, nz, nt, dx, dy, dz, dt, cptr)
+    let ok = niimathx(ptr, fptr, datatype, nx, ny, nz, nt, dx, dy, dz, dt, cptr)
     if (ok != 0) {
         console.log(" -> '", cmd, " generated a fatal error: ", ok)
         return
     }
     console.log("'", cmd, "' required ", Math.round((new Date()) - startTime), "ms")
+    //renew fvars
+    fvars = new Float32Array(instance.exports.memory.buffer, fptr, 4)
+    vars = new Float32Array(4)
+    vars.set(fvars)
+    console.log("slope", vars[0], "intercept", vars[1], "cal_min", vars[2], "cal_max", vars[3])
     //Problem: JavaScript will generate "detached ArrayBuffer" error if malloc requires memory growth
     //Unintuitive Solution: renew cimg pointer
     // https://depth-first.com/articles/2019/10/16/compiling-c-to-webassembly-and-running-it-without-emscripten/
     cimg = new Uint8Array(instance.exports.memory.buffer, ptr, nvox * bpv)
     //copy WASM image data to JavaScript image data
     img8.set(cimg)
     //free WASM memory
+    memory.record_free(fptr)
+    wfree(fptr)
     memory.record_free(cptr)
     wfree(cptr)
     memory.record_free(ptr)

src/tensor.c

@@ -1692,10 +1692,9 @@ void symeig_double( int n , double *a , double *e ) {
 /* Author: Daniel Glen, 15 Nov 2004 */
 //https://github.com/afni/afni/blob/b6a9f7a21c1f3231ff09efbd861f8975ad48e525/src/3dDTeig.c
 //  Components developed at the US National Institutes of Health (after 15 Jan 2001) are not copyrighted.
-void EIG_tsfunc( double tzero, double tdelta ,
+void EIG_tsfunc(
                           int npts, float ts[],
-                          double ts_mean, double ts_slope,
-                          void * ud, int nbriks, float * val, int isUpperTriangle )
+                          float * val, int isUpperTriangle )
 {
   #define SMALLNUMBER 1E-4
   int i,j;

src/tensor.h

@@ -5,10 +5,9 @@
 extern "C" {
 #endif
 
-void EIG_tsfunc( double tzero, double tdelta ,
+void EIG_tsfunc(
                           int npts, float ts[],
-                          double ts_mean, double ts_slope,
-                          void * ud, int nbriks, float * val, int isUpperTriangle );
+                          float * val, int isUpperTriangle );
 
 #ifdef  __cplusplus
 }