diff --git a/Makefile b/Makefile index 5fada13..fd7e5d2 100755 --- a/Makefile +++ b/Makefile @@ -284,6 +284,9 @@ LN2_REGRESS_OUT: LN3_LAYERS: $(CC) $(CFLAGS) -o LN3_LAYERS src/LN3_LAYERS.cpp $(LIBRARIES) $(LFLAGS) +LN3_NOLAD: + $(CC) $(CFLAGS) -o LN3_NOLAD src/LN3_NOLAD.cpp $(LIBRARIES) $(LFLAGS) + # ============================================================================= clean: diff --git a/dep/laynii_lib.cpp b/dep/laynii_lib.cpp index 3dec54f..99c8d8d 100644 --- a/dep/laynii_lib.cpp +++ b/dep/laynii_lib.cpp @@ -1046,3 +1046,373 @@ nifti_image* iterative_smoothing(nifti_image* nii_in, int iter_smooth, } return nii_smooth; } + +// ============================================================================ +// WIP NOLAD... +// ============================================================================ + +float ln_gaussian(float distance, float sigma) { + return 1. / (sigma * std::sqrt(2. * 3.141592)) + * std::exp(-0.5 * distance * distance / (sigma * sigma)); +} + +void ln_normalize_to_zero_one(float* data, int data_size) { + + // Find minimum and maximum + float temp_max = std::numeric_limits::min(); + float temp_min = std::numeric_limits::max(); + for (int i = 0; i != data_size; ++i) { + if (*(data + i) > temp_max) { + temp_max = *(data + i); + } + if (*(data + i) < temp_min) { + temp_min = *(data + i); + } + } + + // Translate minimum to zero + for (int i = 0; i != data_size; ++i) { + *(data + i) -= temp_min; + } + + // Scale maximum to one + temp_max -= temp_min; + for (int i = 0; i != data_size; ++i) { + *(data + i) /= temp_max; + } +} + + +void ln_smooth_gaussian_iterative_3D(float* data_in, + const int nx, const int ny, const int nz, const int nt, + const float dx, const float dy, const float dz, + const float fwhm, const int nr_iterations, const bool log) { + // NOTE: Overwrites the input with the smoothed image at the end + + int data_size = nx * ny * nz * nt; + + float* data_temp = (float*)malloc(data_size * sizeof(float)); + + // Compute Gaussian weights + float w_0 = ln_gaussian(0, fwhm); + float w_dX = ln_gaussian(dx, fwhm); + float w_dY = ln_gaussian(dy, fwhm); + float w_dZ = ln_gaussian(dz, fwhm); + + // Loop over every data point + int ix, iy, iz, it, j; + for (int n = 0; n != nr_iterations; ++n) { + + if (log) std::printf(" Iteration: %i/%i\n", n+1, nr_iterations); + + for (int i = 0; i != data_size; ++i) { + + std::tie(ix, iy, iz, it) = ind2sub_4D(i, nx, ny, nz); + float new_val = 0, total_weight = 0; + + // Start with the voxel itself + new_val += *(data_in + i) * w_0; + total_weight += w_0; + + // ---------------------------------------------------------------- + // 1-jump neighbours + // ---------------------------------------------------------------- + if (ix > 0) { + j = sub2ind_4D(ix-1, iy, iz, it, nx, ny, nz); + new_val += *(data_in + j) * w_dX; + total_weight += w_dX; + } + if (ix < nx-1) { + j = sub2ind_4D(ix+1, iy, iz, it, nx, ny, nz); + new_val += *(data_in + j) * w_dX; + total_weight += w_dX; + } + if (iy > 0) { + j = sub2ind_4D(ix, iy-1, iz, it, nx, ny, nz); + new_val += *(data_in + j) * w_dY; + total_weight += w_dY; + } + if (iy < ny-1) { + j = sub2ind_4D(ix, iy+1, iz, it, nx, ny, nz); + new_val += *(data_in + j) * w_dY; + total_weight += w_dY; + } + if (iz > 0) { + j = sub2ind_4D(ix, iy, iz-1, it, nx, ny, nz); + new_val += *(data_in + j) * w_dZ; + total_weight += w_dZ; + } + if (iz < nz-1) { + j = sub2ind_4D(ix, iy, iz+1, it, nx, ny, nz); + new_val += *(data_in + j) * w_dZ; + total_weight += w_dZ; + } + + // ---------------------------------------------------------------- + *(data_temp + i) = new_val / total_weight; + } + + // Swap image data for the next iteration + for (int i = 0; i != data_size; ++i) { + *(data_in + i) = *(data_temp + i); + } + } + free(data_temp); +} + + +void ln_compute_gradients_3D(const float* data_in, float* data_grad_x, float* data_grad_y, float* data_grad_z, + const int nx, const int ny, const int nz, const int nt) { + + int data_size = nx * ny * nz * nt; + + // Loop over every data point + int ix, iy, iz, it, j, k; + for (int i = 0; i != data_size; ++i) { + std::tie(ix, iy, iz, it) = ind2sub_4D(i, nx, ny, nz); + + // -------------------------------------------------------------------- + // 1-jump neighbours + // -------------------------------------------------------------------- + if (ix > 0 && ix < nx-1) { + j = sub2ind_4D(ix-1, iy, iz, it, nx, ny, nz); + k = sub2ind_4D(ix+1, iy, iz, it, nx, ny, nz); + *(data_grad_x + i) = *(data_in + j) - *(data_in + k); + } + if (iy > 0 && iy < ny-1) { + j = sub2ind_4D(ix, iy-1, iz, it, nx, ny, nz); + k = sub2ind_4D(ix, iy+1, iz, it, nx, ny, nz); + *(data_grad_y + i) = *(data_in + j) - *(data_in + k); + } + if (iz > 0 && iz < nz-1) { + j = sub2ind_4D(ix, iy, iz-1, it, nx, ny, nz); + k = sub2ind_4D(ix, iy, iz+1, it, nx, ny, nz); + *(data_grad_z + i) = *(data_in + j) - *(data_in + k); + } + } +} + + +void ln_compute_gradients_3D_over_x(const float* data_in, float* data_out, + const int nx, const int ny, const int nz, const int nt) { + + int data_size = nx * ny * nz * nt; + + // Loop over every data point + int ix, iy, iz, it, j, k; + for (int i = 0; i != data_size; ++i) { + std::tie(ix, iy, iz, it) = ind2sub_4D(i, nx, ny, nz); + if (ix > 0 && ix < nx-1) { + j = sub2ind_4D(ix-1, iy, iz, it, nx, ny, nz); + k = sub2ind_4D(ix+1, iy, iz, it, nx, ny, nz); + *(data_out + i) = *(data_in + j) - *(data_in + k); + } + } +} + + +void ln_compute_gradients_3D_over_y(const float* data_in, float* data_out, + const int nx, const int ny, const int nz, const int nt) { + + int data_size = nx * ny * nz * nt; + + // Loop over every data point + int ix, iy, iz, it, j, k; + for (int i = 0; i != data_size; ++i) { + std::tie(ix, iy, iz, it) = ind2sub_4D(i, nx, ny, nz); + if (iy > 0 && iy < ny-1) { + j = sub2ind_4D(ix, iy-1, iz, it, nx, ny, nz); + k = sub2ind_4D(ix, iy+1, iz, it, nx, ny, nz); + *(data_out + i) = *(data_in + j) - *(data_in + k); + } + } +} + + +void ln_compute_gradients_3D_over_z(const float* data_in, float* data_out, + const int nx, const int ny, const int nz, const int nt) { + + int data_size = nx * ny * nz * nt; + + // Loop over every data point + int ix, iy, iz, it, j, k; + for (int i = 0; i != data_size; ++i) { + std::tie(ix, iy, iz, it) = ind2sub_4D(i, nx, ny, nz); + if (iz > 0 && iz < nz-1) { + j = sub2ind_4D(ix, iy, iz-1, it, nx, ny, nz); + k = sub2ind_4D(ix, iy, iz+1, it, nx, ny, nz); + *(data_out + i) = *(data_in + j) - *(data_in + k); + } + } +} + + +void ln_compute_hessian_3D(const float* data_in, float* data_shorthessian, + const int nx, const int ny, const int nz, const int nt, + const int dx, const int dy, const int dz, const int vscale) { + // NOTE: Hessian data should be 6 times larger than the input + // NOTE: Hessian values are saved consecutively for each voxel + // *(data_shorthessian + i*6 + 0) = 2nd derivative xx + // *(data_shorthessian + i*6 + 1) = 2nd derivative xy yx + // *(data_shorthessian + i*6 + 2) = 2nd derivative xz zx + // *(data_shorthessian + i*6 + 3) = 2nd derivative yy + // *(data_shorthessian + i*6 + 4) = 2nd derivative yz zy + // *(data_shorthessian + i*6 + 5) = 2nd derivative zz + + int data_size = nx * ny * nz * nt; + float FWHM = 1.0; + + // Allocate memory (NOTE: I have prioritized RAM optimization) + float* data_grad_1st = (float*)malloc(data_size * sizeof(float)); + float* data_grad_2nd = (float*)malloc(data_size * sizeof(float)); + + // x + ln_compute_gradients_3D_over_x(data_in, data_grad_1st, nx, ny, nz, nt); + if (vscale > 0) { + std::printf("\n Smoothing 1st gradient (iterative 3D Gaussian [FWHM = %f, iterations = %i])...\n", FWHM, vscale); + ln_smooth_gaussian_iterative_3D(data_grad_1st, nx, ny, nz, nt, dx, dy, dz, FWHM, vscale); + } + + // xx + ln_compute_gradients_3D_over_x(data_grad_1st, data_grad_2nd, nx, ny, nz, nt); + for (uint32_t i = 0; i != data_size; ++i) { + *(data_shorthessian + i*6 + 0) = *(data_grad_2nd + i); + } + + // xy + ln_compute_gradients_3D_over_y(data_grad_1st, data_grad_2nd, nx, ny, nz, nt); + for (uint32_t i = 0; i != data_size; ++i) { + *(data_shorthessian + i*6 + 1) = *(data_grad_2nd + i); + } + + // xz + ln_compute_gradients_3D_over_z(data_grad_1st, data_grad_2nd, nx, ny, nz, nt); + for (uint32_t i = 0; i != data_size; ++i) { + *(data_shorthessian + i*6 + 2) = *(data_grad_2nd + i); + } + + // y + ln_compute_gradients_3D_over_y(data_in, data_grad_1st, nx, ny, nz, nt); + if (vscale > 0) { + std::printf("\n Smoothing 2nd gradient (iterative 3D Gaussian [FWHM = %f, iterations = %i])...\n", FWHM, vscale); + ln_smooth_gaussian_iterative_3D(data_grad_1st, nx, ny, nz, nt, dx, dy, dz, FWHM, vscale); + } + + // yy + ln_compute_gradients_3D_over_y(data_grad_1st, data_grad_2nd, nx, ny, nz, nt); + for (uint32_t i = 0; i != data_size; ++i) { + *(data_shorthessian + i*6 + 3) = *(data_grad_2nd + i); + } + + // yz + ln_compute_gradients_3D_over_z(data_grad_1st, data_grad_2nd, nx, ny, nz, nt); + for (uint32_t i = 0; i != data_size; ++i) { + *(data_shorthessian + i*6 + 4) = *(data_grad_2nd + i); + } + + // z + ln_compute_gradients_3D_over_z(data_in, data_grad_1st, nx, ny, nz, nt); + if (vscale > 0) { + std::printf("\n Smoothing 3rd gradient (iterative 3D Gaussian [FWHM = %f, iterations = %i])...\n", FWHM, vscale); + ln_smooth_gaussian_iterative_3D(data_grad_1st, nx, ny, nz, nt, dx, dy, dz, FWHM, vscale); + } + + // zz + ln_compute_gradients_3D_over_z(data_grad_1st, data_grad_2nd, nx, ny, nz, nt); + for (uint32_t i = 0; i != data_size; ++i) { + *(data_shorthessian + i*6 + 5) = *(data_grad_2nd + i); + } + + free(data_grad_1st); + free(data_grad_2nd); +} + +void ln_compute_eigen_values_3D(const float* data_hessian, float* data_eigval1, float* data_eigval2, float* data_eigval3, + const int nx, const int ny, const int nz, const int nt) { + // NOTE: Implementing Delledalle et al. 2017, Hal. + // NOTE: I simplified complex conjugates as I do not have complex values. + + int data_size = nx * ny * nz * nt; + + for (uint32_t i = 0; i != data_size; ++i) { + float a = *(data_hessian + i*6 + 0); // xx + float b = *(data_hessian + i*6 + 3); // yy + float c = *(data_hessian + i*6 + 5); // zz + float d = *(data_hessian + i*6 + 1); // xy, yx + float e = *(data_hessian + i*6 + 4); // yz, zy + float f = *(data_hessian + i*6 + 2); // xz, zx + + // -------------------------------------------------------------------- + // Compute eigen values 3D + // -------------------------------------------------------------------- + float x1 = a*a + b*b + c*c - a*b - a*c - b*c + 3 * (d*d + f*f + e*e); + + float t1 = 2*a - b - c; + float t2 = 2*b - a - c; + float t3 = 2*c - a - b; + float x2 = - t1 * t2 * t3 + 9*( t3*(d*d) + t2*(f*f) + t1*(e*e) ) - 54*( d*e*f ); + + float phi; + if (x2 > 0) { + phi = std::atan2( 0, std::sqrt( 4.*( (x1*x1)*x1 ) - x2*x2 ) / x2 ); + } else if (x2 < 0) { + phi = std::atan2( 0, std::sqrt( 4.*( (x1*x1)*x1 ) - x2*x2 ) / x2 ) + M_PI; + } else { + phi = M_PI / 2.; + } + + float lambda1 = ( a + b + c - 2.*std::sqrt( x1 ) * std::cos( phi/3. ) ) / 3.; + float lambda2 = ( a + b + c + 2.*std::sqrt( x1 ) * std::cos( (phi - M_PI)/3. ) ) / 3.; + float lambda3 = ( a + b + c + 2.*std::sqrt( x1 ) * std::cos( (phi + M_PI)/3. ) ) / 3.; + + if (std::isnan(lambda1)) { + lambda1 = 0.; + } + if (std::isnan(lambda2)) { + lambda2 = 0.; + } + if (std::isnan(lambda3)) { + lambda3 = 0.; + } + + // -------------------------------------------------------------------- + // Sort by magnitude while preserving sign of the eigenvalues + // -------------------------------------------------------------------- + bool lambda1_sign = std::signbit(lambda1); + bool lambda2_sign = std::signbit(lambda2); + bool lambda3_sign = std::signbit(lambda3); + + float lambda1_abs = std::abs(lambda1); + float lambda2_abs = std::abs(lambda2); + float lambda3_abs = std::abs(lambda3); + + *(data_eigval1 + i) = std::min( std::min( lambda1_abs, lambda2_abs ), lambda3_abs ); + *(data_eigval3 + i) = std::max( std::max( lambda1_abs, lambda2_abs ), lambda3_abs ); + *(data_eigval2 + i) = lambda1_abs + lambda2_abs + lambda3_abs - ( *(data_eigval1 + i) + *(data_eigval3 + i) ); + + // Put back the signbit + if (*(data_eigval1 + i) == lambda1_abs) { + *(data_eigval1 + i) *= static_cast(lambda1_sign) * 2 - 1; + } else if (*(data_eigval1 + i) == lambda2_abs) { + *(data_eigval1 + i) *= static_cast(lambda2_sign) * 2 - 1; + } else if (*(data_eigval1 + i) == lambda3_abs) { + *(data_eigval1 + i) *= static_cast(lambda3_sign) * 2 - 1; + } + + if (*(data_eigval2 + i) == lambda1_abs) { + *(data_eigval2 + i) *= static_cast(lambda1_sign) * 2 - 1; + } else if (*(data_eigval2 + i) == lambda2_abs) { + *(data_eigval2 + i) *= static_cast(lambda2_sign) * 2 - 1; + } else if (*(data_eigval2 + i) == lambda3_abs) { + *(data_eigval2 + i) *= static_cast(lambda3_sign) * 2 - 1; + } + + if (*(data_eigval3 + i) == lambda1_abs) { + *(data_eigval3 + i) *= static_cast(lambda1_sign) * 2 - 1; + } else if (*(data_eigval3 + i) == lambda2_abs) { + *(data_eigval3 + i) *= static_cast(lambda2_sign) * 2 - 1; + } else if (*(data_eigval3 + i) == lambda3_abs) { + *(data_eigval3 + i) *= static_cast(lambda3_sign) * 2 - 1; + } + } +} diff --git a/dep/laynii_lib.h b/dep/laynii_lib.h index 0698540..4c4db5b 100644 --- a/dep/laynii_lib.h +++ b/dep/laynii_lib.h @@ -78,3 +78,34 @@ nifti_image* iterative_smoothing(nifti_image* nii_in, int iter_smooth, // Preprocessor macros. // ============================================================================ #define PI 3.14159265; + +// ============================================================================ +// WIP NOLAD... +// ============================================================================ + +float ln_gaussian(float distance, float sigma); + +void ln_normalize_to_zero_one(float* data, int data_size); + +void ln_smooth_gaussian_iterative_3D(float* data_in, + const int nx, const int ny, const int nz, const int nt, + const float dx, const float dy, const float dz, + const float FWHM_val, const int nr_iterations, const bool log = true); + +void ln_compute_gradients_3D(const float* data, float* data_grad_x, float* data_grad_y, float* data_grad_z, + const int nx, const int ny, const int nz, const int nt); + + +void ln_compute_gradients_3D_over_x(const float* data, float* data_out, + const int nx, const int ny, const int nz, const int nt); +void ln_compute_gradients_3D_over_y(const float* data, float* data_out, + const int nx, const int ny, const int nz, const int nt); +void ln_compute_gradients_3D_over_z(const float* data, float* data_out, + const int nx, const int ny, const int nz, const int nt); + +void ln_compute_hessian_3D(const float* data, float* data_shorthessian, + const int nx, const int ny, const int nz, const int nt, + const int dx, const int dy, const int dz, const int nr_smth_iterations); + +void ln_compute_eigen_values_3D(const float* data_shorthessian, float* data_eigval1, float* data_eigval2, float* data_eigval3, + const int nx, const int ny, const int nz, const int nt); diff --git a/src/LN3_NOLAD.cpp b/src/LN3_NOLAD.cpp new file mode 100644 index 0000000..5f2bcd3 --- /dev/null +++ b/src/LN3_NOLAD.cpp @@ -0,0 +1,232 @@ +#include "../dep/laynii_lib.h" +#include + +int show_help(void) { + printf( + "LN3_NOLAD: Nonlinear anisotropic diffusion filter.\n" + "\n" + "Usage:\n" + " LN3_NOLAD -input input.nii\n" + " ../LN3_NOLAD -input input.nii\n" + "\n" + "Options:\n" + " -help : Show this help.\n" + " -input : Nifti image that will be used to compute gradients.\n" + " This can be a 4D nifti. in 4D case, 3D gradients\n" + " will be computed for each volume.\n" + " -nscale : (Optional) Noise scale. Number of Gaussian smoothing iterations applied \n" + " to scalar image. No smoothing ('0') by default.\n" + " -fscale : (Optional) Feature scale. Number of Gaussian smoothing iterations applied \n" + " to first order gradients (vector field). No smoothing ('0') by default.\n" + " -output : (Optional) Output basename for all outputs.\n" + " -debug : (Optional) Save extra intermediate outputs.\n" + "\n" + "\n"); + return 0; +} + +// NOTE[Faruk]: References are: +// - Weickert, J. (1998). Anisotropic diffusion in image processing. Image Rochester NY, 256(3), 170. +// - Mirebeau, J.-M., Fehrenbach, J., Risser, L., & Tobji, S. (2015). Anisotropic Diffusion in ITK, 1-9. + +int main(int argc, char* argv[]) { + nifti_image *nii1 = NULL; + char *fin1 = NULL, *fout = NULL; + int ac; + bool mode_debug = false; + int NSCALE = 0, FSCALE = 0; + float LAMBDA=0.001, ALPHA=0.001, M=4; + + + // Process user options + if (argc < 2) return show_help(); + for (ac = 1; ac < argc; ac++) { + if (!strncmp(argv[ac], "-h", 2)) { + return show_help(); + } else if (!strcmp(argv[ac], "-input")) { + if (++ac >= argc) { + fprintf(stderr, "** missing argument for -input\n"); + return 1; + } + fin1 = argv[ac]; + fout = argv[ac]; + } else if (!strcmp(argv[ac], "-nscale")) { + if (++ac >= argc) { + fprintf(stderr, "** missing argument for -nscale\n"); + } else { + NSCALE = atof(argv[ac]); + } + } else if (!strcmp(argv[ac], "-fscale")) { + if (++ac >= argc) { + fprintf(stderr, "** missing argument for -fscale\n"); + } else { + FSCALE = atof(argv[ac]); + } + } else if (!strcmp(argv[ac], "-debug")) { + mode_debug = true; + } else if (!strcmp(argv[ac], "-output")) { + if (++ac >= argc) { + fprintf(stderr, "** missing argument for -output\n"); + return 1; + } + fout = argv[ac]; + } else { + fprintf(stderr, "** invalid option, '%s'\n", argv[ac]); + return 1; + } + } + + if (!fin1) { + fprintf(stderr, "** missing option '-input'\n"); + return 1; + } + + // Read input dataset, including data + nii1 = nifti_image_read(fin1, 1); + if (!nii1) { + fprintf(stderr, "** failed to read NIfTI from '%s'\n", fin1); + return 2; + } + + log_welcome("LN3_NOLAD"); + log_nifti_descriptives(nii1); + + // Get dimensions of input + const uint32_t nx = nii1->nx; + const uint32_t ny = nii1->ny; + const uint32_t nz = nii1->nz; + const uint32_t nt = nii1->nt; + + const uint32_t end_x = nx - 1; + const uint32_t end_y = ny - 1; + const uint32_t end_z = nz - 1; + + const uint32_t nr_voxels = nz * ny * nx; + const uint32_t data_size = nz * ny * nx * nt; + + const float dx = nii1->pixdim[1]; + const float dy = nii1->pixdim[2]; + const float dz = nii1->pixdim[3]; + + // ======================================================================== + // Fix input datatype issues + // ======================================================================== + nifti_image* nii_input = copy_nifti_as_float32(nii1); + float* data_input = static_cast(nii_input->data); + + // Prepare generic output niftis + nifti_image* nii_out_float32 = copy_nifti_as_float32(nii1); + float* nii_out_float32_data = static_cast(nii_out_float32->data); + + // ======================================================================== + // Normalize by maximum + // ======================================================================== + std::printf("\n Normalizing (minimum to 0 and maximum to 1)...\n"); + ln_normalize_to_zero_one(data_input, data_size); + + if (mode_debug) { + std::printf(" DEBUG: Saving output...\n"); + save_output_nifti(fout, "DEBUG1-normalize_to_zero_one", nii_input, true); + } + + // ======================================================================== + // Noise scale smoothing + // ======================================================================== + float FWHM = 1.0; // Default + if (NSCALE > 0) { + std::printf("\n Smoothing (iterative 3D Gaussian [FWHM = %f, iterations = %i])...\n", FWHM, NSCALE); + ln_smooth_gaussian_iterative_3D(data_input, nx, ny, nz, nt, dx, dy, dz, FWHM, NSCALE); + + if (mode_debug) { + std::printf(" DEBUG: Saving output...\n"); + save_output_nifti(fout, "DEBUG2-smooth_gaussian", nii_input, true); + } + } + + // ======================================================================== + // Compute Hessian + // ======================================================================== + std::printf("\n Computing Hessian matrices...\n"); + + float* data_hessian = (float*)malloc(data_size * 6 * sizeof(float)); + ln_compute_hessian_3D(data_input, data_hessian, nx, ny, nz, nt, dx, dy, dz, FSCALE); + + // ------------------------------------------------------------------------- + if (mode_debug) { + // Compute trace (for testing) + float* data_trace = (float*)malloc(data_size * sizeof(float)); + for (uint32_t i = 0; i != data_size; ++i) { + *(data_trace + i) = *(data_hessian + i*6 + 0) + *(data_hessian + i*6 + 3) + *(data_hessian + i*6 + 5); + } + std::printf(" DEBUG: Saving output...\n"); + for (uint32_t i = 0; i != data_size; ++i) *(nii_out_float32_data + i) = *(data_trace + i); + free(data_trace); + save_output_nifti(fout, "DEBUG3-hessian_trace", nii_out_float32, true); + + // Compute off-trace (for testing) + float* data_offtrace = (float*)malloc(data_size * sizeof(float)); + for (uint32_t i = 0; i != data_size; ++i) { + *(data_offtrace + i) = *(data_hessian + i*6 + 1) + *(data_hessian + i*6 + 2) + *(data_hessian + i*6 + 4); + } + std::printf(" DEBUG: Saving output...\n"); + for (uint32_t i = 0; i != data_size; ++i) *(nii_out_float32_data + i) = *(data_offtrace + i); + free(data_offtrace); + save_output_nifti(fout, "DEBUG3-hessian_offtrace", nii_out_float32, true); + } + + // ======================================================================== + // Eigendecomposition + // ======================================================================== + std::printf("\n Computing Eigen values...\n"); + + float* data_eigval1 = (float*)malloc(data_size * sizeof(float)); + float* data_eigval2 = (float*)malloc(data_size * sizeof(float)); + float* data_eigval3 = (float*)malloc(data_size * sizeof(float)); + ln_compute_eigen_values_3D(data_hessian, data_eigval1, data_eigval2, data_eigval3, nx, ny, nz, nt); + + // ------------------------------------------------------------------------- + if (mode_debug) { + std::printf(" DEBUG: Saving output...\n"); + for (uint32_t i = 0; i != data_size; ++i) *(nii_out_float32_data + i) = *(data_eigval1 + i); + save_output_nifti(fout, "DEBUG4-eigen_value_1", nii_out_float32, true); + for (uint32_t i = 0; i != data_size; ++i) *(nii_out_float32_data + i) = *(data_eigval2 + i); + save_output_nifti(fout, "DEBUG4-eigen_value_2", nii_out_float32, true); + for (uint32_t i = 0; i != data_size; ++i) *(nii_out_float32_data + i) = *(data_eigval3 + i); + save_output_nifti(fout, "DEBUG4-eigen_value_3", nii_out_float32, true); + } + + // ======================================================================== + // Compute diffusion weights + // ======================================================================== + float* data_diffw1 = (float*)malloc(data_size * sizeof(float)); + float* data_diffw2 = (float*)malloc(data_size * sizeof(float)); + float* data_diffw3 = (float*)malloc(data_size * sizeof(float)); + + for (uint32_t i = 0; i != data_size; ++i) { + // Apply compositional closure + float eigvalsum = std::abs(*(data_eigval1 + i)) + std::abs(*(data_eigval2 + i)) + std::abs(*(data_eigval3 + i)); + *(data_diffw1 + i) = 1 - (std::abs(*(data_eigval1 + i)) / eigvalsum); + *(data_diffw2 + i) = 1 - (std::abs(*(data_eigval2 + i)) / eigvalsum); + *(data_diffw3 + i) = 1 - (std::abs(*(data_eigval3 + i)) / eigvalsum); + + // Reclose for balance + eigvalsum = *(data_diffw1 + i) + *(data_diffw2 + i) + *(data_diffw3 + i); + *(data_diffw1 + i) /= eigvalsum; + *(data_diffw2 + i) /= eigvalsum; + *(data_diffw3 + i) /= eigvalsum; + } + + // ------------------------------------------------------------------------- + if (mode_debug) { + std::printf(" DEBUG: Saving output...\n"); + for (uint32_t i = 0; i != data_size; ++i) *(nii_out_float32_data + i) = *(data_diffw1 + i); + save_output_nifti(fout, "DEBUG5-diffweight_1", nii_out_float32, true); + for (uint32_t i = 0; i != data_size; ++i) *(nii_out_float32_data + i) = *(data_diffw2 + i); + save_output_nifti(fout, "DEBUG5-diffweight_2", nii_out_float32, true); + for (uint32_t i = 0; i != data_size; ++i) *(nii_out_float32_data + i) = *(data_diffw3 + i); + save_output_nifti(fout, "DEBUG5-diffweight_3", nii_out_float32, true); + } + + cout << "\n Finished." << endl; + return 0; +}