Add FDR

package.json

@@ -28,5 +28,8 @@
     "url": "https://github.com/rordenlab/niimath/issues"
   },
   "homepage": "https://github.com/rordenlab/niimath#readme",
-  "devDependencies": {}
-}
+  "dependencies": {
+    "nifti-reader-js": "^0.6.8",
+    "pako": "^2.1.0"
+  }
+}

src/JavaScript.md

@@ -165,3 +165,11 @@ WASM memory management is limited. JavaScript protects memory, improving securit
 | -roi [xmin] [xsize] [ymin] [ysize] [zmin] [zsize] [tmin] [tsize]  | NO  |  zero outside roi (using voxel coordinates). Inputting -1 for a size will set it to the full image extent for that dimension |
 | -bptf  [hp_sigma] [lp_sigma]  | NO  |  (-t in ip.c) Bandpass temporal filtering; nonlinear highpass and Gaussian linear lowpass (with sigmas in volumes, not seconds); set either sigma[0 to skip that filter |
 | -roc [AROC-thresh] [outfile] [4Dnoiseonly] [truth]  | NO  |  take (normally binary) truth and test c |
+
+### Alternative
+
+One can use emscriptens built in file system and memory support if you compile as follows:
+
+```
+emcc -s USE_ZLIB=1 -O1 -DHAVE_ZLIB -DFSLSTYLE -DREJECT_COMPLEX -DNII2MESH -std=gnu99 -o niimathwasm.js niimath.c MarchingCubes.c meshify.c quadric.c base64.c radixsort.c fdr.c bwlabel.c bw.c core.c tensor.c core32.c core64.c niftilib/nifti2_io.c znzlib/znzlib.c -I./niftilib -I./znzlib -lm -lz
+```

src/Makefile

@@ -3,7 +3,7 @@ CNAME=gcc
 #linker flags
 LFLAGS=-lm -lz
 #c flags
-CFLAGS=-O3 -DHAVE_ZLIB -DFSLSTYLE -DPIGZ -DREJECT_COMPLEX -DNII2MESH
+CFLAGS=-O1 -DHAVE_ZLIB -DFSLSTYLE -DPIGZ -DREJECT_COMPLEX -DNII2MESH
 # Westmere required for fast gz https://github.com/tikv/tikv/issues/4843
 #CFLAGS=-O3 -DHAVE_ZLIB -march=westmere
 
@@ -22,5 +22,5 @@ ifeq "$(OMP)" "1"
 endif
 
 all:
-	$(CNAME) $(CFLAGS) -o niimath niimath.c MarchingCubes.c meshify.c quadric.c base64.c radixsort.c bwlabel.c bw.c core.c tensor.c core32.c core64.c niftilib/nifti2_io.c znzlib/znzlib.c -I./niftilib -I./znzlib $(LFLAGS)
+	$(CNAME) $(CFLAGS) -o niimath niimath.c MarchingCubes.c meshify.c quadric.c base64.c radixsort.c fdr.c bwlabel.c bw.c core.c tensor.c core32.c core64.c niftilib/nifti2_io.c znzlib/znzlib.c -I./niftilib -I./znzlib $(LFLAGS)
 

src/bw.c

@@ -52,6 +52,11 @@
   #include <immintrin.h>
 #endif
 
+#ifndef _mm_malloc
+  #define _mm_malloc(size, alignment) malloc(size)
+  #define _mm_free(ptr) free(ptr)
+#endif
+
 #ifndef M_PI //not defined by older gcc unless compiled with `-std=gnu99`:
   #define M_PI 3.14159265358979323846
 #endif

src/core.c

@@ -8,7 +8,9 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <ctype.h>
-#define USE_MESH
+#ifndef USING_WASM
+	#define USE_MESH
+#endif
 #ifdef USE_MESH
 	#include "meshify.h"
 	#include "quadric.h"
@@ -18,6 +20,10 @@
 #else
 	#include <immintrin.h>
 #endif
+#ifndef _mm_malloc
+  #define _mm_malloc(size, alignment) malloc(size)
+  #define _mm_free(ptr) free(ptr)
+#endif
 
 #ifdef USING_WASM
 #define WASM_EXPORT(name) \
@@ -503,7 +509,7 @@ int nifti_image_change_datatype(nifti_image *nim, int dt, in_hdr *ihdr) {
 		double *o64 = (double *)dat;
 		if (idt == DT_FLOAT32) {
 			for (size_t i = 0; i < nim->nvox; i++)
-				o64[i] = (u32[i] * scl) + inter;
+				o64[i] = (f32[i] * scl) + inter; //<<<<<<<<>>>>
 			ok = 0;
 		}
 		if (idt == DT_UINT64) {
@@ -1193,19 +1199,12 @@ static double Mitchell_filter(double t) {
 	return (0.0);
 }
 
-static double filter(double t) {
-	/* f(t) = 2|t|^3 - 3|t|^2 + 1, -1 <= t <= 1 */
-	if (t < 0.0)
-		t = -t;
-	if (t < 1.0)
-		return ((2.0 * t - 3.0) * t * t + 1.0);
-	return (0.0);
-}
-
 CLIST *createFilter(int srcXsize, int dstXsize, int filterMethod) {
-	//CLIST	* createFilter(int srcXsize, int dstXsize, double (*filterf)(), double fwidth) {
-	double (*filterf)() = filter;
-	double fwidth;
+	//Schumacher's resampler in Graphics Gems 3, for improvements see
+	//see https://github.com/richgel999/imageresampler/blob/master/resampler.cpp
+	// method 1 is the default: linear
+	double (*filterf)(double t) = triangle_filter;
+	double fwidth = triangle_support;
 	if (filterMethod == 0) {
 		filterf = box_filter;
 		fwidth = box_support;
@@ -1218,9 +1217,6 @@ CLIST *createFilter(int srcXsize, int dstXsize, int filterMethod) {
 	} else if (filterMethod == 4) {
 		filterf = Mitchell_filter;
 		fwidth = Mitchell_support;
-	} else { //method 1 is the default: linear
-		filterf = triangle_filter;
-		fwidth = triangle_support;
 	}
 	CLIST *contrib = (CLIST *)calloc(dstXsize, sizeof(CLIST));
 	double xscale = (double)dstXsize / (double)srcXsize;

src/coreFLT.c

@@ -48,6 +48,19 @@
 	#define tensor_decomp //tensor_decomp support is optional
 #endif
 
+#ifndef USING_WASM
+	#ifdef __x86_64__
+		#include <immintrin.h>
+	#else
+		#include "arm_malloc.h"
+	#endif
+#endif
+#ifndef _mm_malloc
+	#define _mm_malloc(size, alignment) malloc(size)
+	#define _mm_free(ptr) free(ptr)
+	#undef SIMD
+#endif
+
 #ifdef SIMD //explicitly vectorize (SSE,AVX,Neon)
 	#ifdef __x86_64__
 		#ifdef DT32
@@ -64,13 +77,6 @@
 		#endif
 	#endif
 #endif
-#ifndef USING_WASM
-	#ifdef __x86_64__
-		#include <immintrin.h>
-	#else
-		#include "arm_malloc.h"
-	#endif
-#endif
 
 #ifdef bandpass
 	#include "bw.h"
@@ -81,7 +87,7 @@
 #ifdef tensor_decomp
 	#include "tensor.h"
 #endif
-
+#include "fdr.h"
 //#define TFCE //formerly we used Christian Gaser's tfce, new bespoke code handles connectivity
 //#ifdef TFCE //we now use in-built tfce function
 //	#include "tfce_pthread.h"
@@ -551,7 +557,7 @@ staticx void blurS(flt *img, int nx, int ny, flt xmm, flt Sigmamm, flt kernelWid
 	flt *kWeight = (flt *)_mm_malloc(nx * sizeof(flt), 64); //ensure sum of kernel = 1.0
 	for (int i = 0; i < nx; i++) {
 		kStart[i] = MAX(-cutoffvox, -i); //do not read below 0
-		kEnd[i] = MIN(cutoffvox, nx - i - 1); //do not read beyond final columnn
+		kEnd[i] = MIN(cutoffvox, nx - i - 1); //do not read beyond final column
 		if ((i > 0) && (kStart[i] == (kStart[i - 1])) && (kEnd[i] == (kEnd[i - 1]))) { //reuse weight
 			kWeight[i] = kWeight[i - 1];
 			continue;
@@ -608,7 +614,7 @@ staticx void blurP(flt *img, int nx, int ny, flt xmm, flt FWHMmm, flt kernelWid)
 	flt *kWeight = (flt *)_mm_malloc(nx * sizeof(flt), 64); //ensure sum of kernel = 1.0
 	for (int i = 0; i < nx; i++) {
 		kStart[i] = MAX(-cutoffvox, -i); //do not read below 0
-		kEnd[i] = MIN(cutoffvox, nx - i - 1); //do not read beyond final columnn
+		kEnd[i] = MIN(cutoffvox, nx - i - 1); //do not read beyond final column
 		if ((i > 0) && (kStart[i] == (kStart[i - 1])) && (kEnd[i] == (kEnd[i - 1]))) { //reuse weight
 			kWeight[i] = kWeight[i - 1];
 			continue;
@@ -4596,6 +4602,29 @@ staticx int nifti_binary(nifti_image *nim, char *fin, enum eOp op) {
 	return 0;
 } // nifti_binary()
 
+//bork
+staticx int nifti_fdr(nifti_image *nim, double qval) {
+	if (nim->datatype != DT_FLOAT32) {
+		printfx("nifti_fdr: Unsupported datatype %d\n", nim->datatype);
+		exit(1);
+	}
+	if (nim->intent_code != NIFTI_INTENT_PVAL) {
+		printfx("nifti_fdr intent_code should be 22 (NIFTI_INTENT_PVAL) not %d\n", nim->intent_code);
+	}
+	#ifdef DT32
+	size_t nvox = nim->nx * nim->ny * nim->nz;
+	flt *img = (flt *)nim->data;
+	float thresh = fdr(img, qval, nvox);
+	printfx("nifti_fdr %g %g\n", qval, thresh);
+	#else
+	printfx("'-dt double' does not support fdr\n" );
+	exit(1);
+	#endif
+
+
+	exit(1);
+}
+
 staticx void nifti_compare(nifti_image *nim, char *fin) {
 	if (nim->nvox < 1)
 		exit(1);
@@ -5205,6 +5234,10 @@ int mainWASM(nifti_image *nim, int argc, char *argv[]) {
 		} else if (!strcmp(argv[ac], "--compare")) { //--function terminates without saving image
 			ac++;
 			nifti_compare(nim, argv[ac]); //always terminates
+		} else if (!strcmp(argv[ac], "-fdr")) { //--function terminates without saving image
+				ac++;
+				double qval = strtod(argv[ac], &end);
+			nifti_fdr(nim, qval); //always terminates
 		}
 		#endif //WASM does not (yet) resize images
 		else if (!strcmp(argv[ac], "-edt"))
@@ -5449,7 +5482,7 @@ static char* splitArgv(char **str, char **word){
 		(*str)++;
 		inquotes = true;
 	}
-	// Set phrase begining
+	// Set phrase beginning
 	*word = *str;
 	// Skip all chars if in quotes
 	if( inquotes ){

src/fdr.c

@@ -0,0 +1,65 @@
+#include <stdlib.h>
+#include <limits.h>
+#include <stdbool.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include "radixsort.h"
+#ifdef _MSC_VER
+
+#else 
+ #include <unistd.h>
+#endif
+#include <stdint.h>
+
+#define printfd(...) fprintf(stderr, __VA_ARGS__)
+
+#ifndef MAX
+#define MAX(A,B) ((A) > (B) ? (A) : (B))
+#endif
+
+#ifndef MIN
+#define MIN(A,B) ((A) > (B) ? (B) : (A))
+#endif
+
+float fdr(float *pvals, double qval, size_t nvox) {
+	double cV = 1.0;
+	if ((qval < 0.0) || (qval > 1.0)) {
+		printfd("qval out of range (0-1).\n");
+		return NAN;
+	}
+	//========[PART 1: FDR THRESHOLD]========================================
+	//Sort p-values
+	float* dx_out = (float*)malloc(nvox*sizeof(float));
+	uint32_t* idx_in = (uint32_t*)malloc(nvox*sizeof(uint32_t));
+	uint32_t* idx_out = (uint32_t*)malloc(nvox*sizeof(uint32_t));
+	float mn = pvals[0];
+	float mx = pvals[0];
+	for (int i = 0; i < nvox; i++) {
+		idx_in[i] = i;
+		mn = MIN(mn, pvals[i]);
+		mx = MAX(mx, pvals[i]);
+	}
+	if ((mn < 0.0) || (mx > 1.0)) {
+		printfd("Values out of range (0-1).\n");
+		return NAN;
+	}
+	radix11sort_f32(pvals, dx_out, idx_in, idx_out, nvox);
+	free(idx_in);
+	free(idx_out);
+
+	// Order (indices), in the same size as the pvalues
+	int v = 0;
+	do {
+		int v2 = v;
+		do {
+			v2++;
+		} while ( v2 < nvox );
+		v++;
+	} while ( v < nvox );
+
+	free(dx_out);
+	float thresh = 66.0;
+	return thresh;
+}
+

src/fdr.h

@@ -0,0 +1,16 @@
+#ifndef EX_FDR_H
+#define EX_FDR_H
+
+#ifdef  __cplusplus
+//extern "C" {
+#endif
+
+#include <stdbool.h>
+
+float fdr(float *pvals, double qval, size_t nvox);
+
+#ifdef  __cplusplus
+//}
+#endif
+
+#endif // EX_FDR_H