Skip to content

Commit ff5c860

Browse files
daduccidaducci
committed
Changed filenames and bit of logic
1 parent b9821dd commit ff5c860

4 files changed

Lines changed: 46 additions & 40 deletions

File tree

‎commit/core.pyx‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1700,6 +1700,7 @@ cdef class Evaluation :
17001700
xic[ self.DICTIONARY['TRK']['kept']==1 ] *= self.DICTIONARY['TRK']['lenTot'] / self.DICTIONARY['TRK']['len']
17011701

17021702
# call (potential) postprocessing required by the specific model
1703+
#TODO: in this the right place to call such a function?
17031704
if hasattr(self.model, '_postprocess') :
17041705
ret_subinfo = logger.subinfo('Calling model-specific postprocessing', indent_lvl=2, indent_char='-', with_progress=True)
17051706
with ProgressBar(disable=self.verbose<3, hide_on_exit=True, subinfo=ret_subinfo):

‎commit/models.pyx‎

Lines changed: 19 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
#!python
22
#cython: language_level=3, boundscheck=False, wraparound=False, profile=False
33
from os.path import join as pjoin
4+
import os
45
import numpy as np
56
import nibabel
67
from amico.models import BaseModel, StickZeppelinBall as _StickZeppelinBall, CylinderZeppelinBall as _CylinderZeppelinBall
@@ -104,33 +105,36 @@ class ScalarMap( BaseModel ) :
104105
return xic
105106

106107
RESULTS_path = evaluation.get_config('RESULTS_path')
107-
niiISO_img = np.asanyarray( nibabel.load( pjoin(RESULTS_path,'compartment_ISO.nii.gz') ).dataobj ).astype(np.float32)
108-
xLES = niiISO_img[evaluation.DICTIONARY['MASK_ix'], evaluation.DICTIONARY['MASK_iy'], evaluation.DICTIONARY['MASK_iz']]
108+
109+
# lesion contribution are stored in the ISO file
110+
niiISO = nibabel.load( pjoin(RESULTS_path,'compartment_ISO.nii.gz') )
111+
affine = niiISO.affine if nibabel.__version__ >= '2.0.0' else niiISO.get_affine()
112+
xLES = np.asarray( niiISO.dataobj, np.float32 )[evaluation.DICTIONARY['MASK_ix'], evaluation.DICTIONARY['MASK_iy'], evaluation.DICTIONARY['MASK_iz']]
109113
if np.count_nonzero(xLES>0) == 0:
110114
logger.warning('No lesions found')
111115
return xic
112116

113117
# rescale the input scalar map in each voxel according to estimated lesion contributions
114-
niiIC_img = np.asanyarray( nibabel.load( pjoin(RESULTS_path,'compartment_IC.nii.gz') ).dataobj ).astype(np.float32) #IC ma puo essere anche myelin ---> signal?
115-
IC = niiIC_img[evaluation.DICTIONARY['MASK_ix'], evaluation.DICTIONARY['MASK_iy'], evaluation.DICTIONARY['MASK_iz']]
116-
R_scaling_factor = np.zeros_like(xLES, dtype=np.float32)
117-
R_scaling_factor[xLES>0] = (IC[xLES>0] - xLES[xLES>0]) / IC[xLES>0]
118-
niiR_img = np.zeros_like(niiISO_img, dtype=np.float32)
119-
niiR_img[evaluation.DICTIONARY['MASK_ix'], evaluation.DICTIONARY['MASK_iy'], evaluation.DICTIONARY['MASK_iz']] = R_scaling_factor
120-
affine = evaluation.niiDWI.affine if nibabel.__version__ >= '2.0.0' else evaluation.niiDWI.get_affine()
121-
nibabel.save(nibabel.Nifti1Image(niiR_img, affine), pjoin(RESULTS_path,'R_scaling_factor.nii.gz'))
118+
niiIC = nibabel.load( pjoin(RESULTS_path,'compartment_IC.nii.gz') )
119+
IC = np.asarray(niiIC.dataobj, np.float32)[ evaluation.DICTIONARY['MASK_ix'], evaluation.DICTIONARY['MASK_iy'], evaluation.DICTIONARY['MASK_iz'] ]
120+
R = np.zeros_like(xLES, dtype=np.float32)
121+
R[xLES>0] = (IC[xLES>0] - xLES[xLES>0]) / IC[xLES>0]
122+
123+
niiR_img = np.zeros(evaluation.get_config('dim'), dtype=np.float32)
124+
niiR_img[evaluation.DICTIONARY['MASK_ix'], evaluation.DICTIONARY['MASK_iy'], evaluation.DICTIONARY['MASK_iz']] = R
125+
nibabel.Nifti1Image(niiR_img, affine).to_filename( pjoin(RESULTS_path,'R_scaling_factor.nii.gz') )
122126

123127
# save the map of local tissue damage estimated in each voxel
124-
nibabel.save( nibabel.Nifti1Image( niiISO_img, affine ), pjoin(RESULTS_path,'compartment_Lesion.nii.gz') )
128+
os.rename( pjoin(RESULTS_path,'compartment_ISO.nii.gz'), pjoin(RESULTS_path,'compartment_LESION.nii.gz') )
125129

126-
# override ISO map and set it to 0
127-
nibabel.save( nibabel.Nifti1Image( 0*niiISO_img, affine), pjoin(RESULTS_path,'compartment_ISO.nii.gz') )
130+
# set ISO map to 0
131+
nibabel.Nifti1Image( 0*niiR_img, affine).to_filename( pjoin(RESULTS_path,'compartment_ISO.nii.gz') )
128132

129133
# rescale each streamline weight
130134
kept = evaluation.DICTIONARY['TRK']['kept']
131135
cdef double [::1] xic_view = xic[kept==1]
132136
cdef double [::1] xic_scaled_view = xic[kept==1].copy()
133-
cdef float [::1] R_scaled_view = R_scaling_factor
137+
cdef float [::1] R_view = R
134138
cdef unsigned int [::1] idx_v_view = evaluation.DICTIONARY['IC']['v']
135139
cdef unsigned int [::1] idx_f_view = evaluation.DICTIONARY['IC']['fiber']
136140
cdef size_t i, idx_v, idx_f
@@ -139,7 +143,7 @@ class ScalarMap( BaseModel ) :
139143
# Rescaling streamline weights accounting for lesions
140144
for i in range(evaluation.DICTIONARY['IC']['v'].shape[0]):
141145
idx_v = idx_v_view[i]
142-
val = R_scaled_view[idx_v]
146+
val = R_view[idx_v]
143147
if val > 0:
144148
idx_f = idx_f_view[i]
145149
#TODO: allow considering other than the min value

‎commit/trk2dictionary/trk2dictionary.pyx‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -480,7 +480,7 @@ cpdef run( filename_tractogram=None, path_out=None, filename_peaks=None, filenam
480480
cdef float* ptrISO
481481
cdef float [:, :, ::1] niiISO_img
482482
if filename_lesion_mask is not None :
483-
logger.subinfo( 'Restricted ISO map', indent_lvl=2, indent_char='-' )
483+
logger.subinfo( 'Restricting ISO compartments to lesion mask:', indent_lvl=2, indent_char='-' )
484484
niiISO = nibabel.load( filename_lesion_mask )
485485
niiISO_hdr = _get_header( niiISO )
486486
logger.subinfo( f'{niiISO.shape[0]} x {niiISO.shape[1]} x {niiISO.shape[2]}', indent_lvl=3, indent_char='-' )
@@ -491,7 +491,7 @@ cpdef run( filename_tractogram=None, path_out=None, filename_peaks=None, filenam
491491
niiISO_img = np.ascontiguousarray( np.asanyarray( niiISO.dataobj ).astype(np.float32) )
492492
ptrISO = &niiISO_img[0,0,0]
493493
else :
494-
logger.subinfo( 'No ISO map specified, using tdi', indent_lvl=2, indent_char='-' )
494+
logger.subinfo( 'No mask specified to filter ISO compartments', indent_lvl=2, indent_char='-' )
495495
ptrISO = NULL
496496

497497
# write dictionary information info file

‎commit/trk2dictionary/trk2dictionary_c.cpp‎

Lines changed: 24 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -76,7 +76,7 @@ float minSegLen, minFiberLen, maxFiberLen;
7676

7777
// Threads variables
7878
vector<thread> threads;
79-
vector<unsigned long int> totICSegments;
79+
vector<unsigned long int> totICSegments;
8080
vector<unsigned int> totFibers;
8181
unsigned int totECVoxels = 0;
8282
unsigned int totECSegments = 0;
@@ -94,7 +94,7 @@ unsigned int read_fiberTCK( FILE* fp, float fiber[3][MAX_FIB_LEN] , float* toVOX
9494

9595
// ---------- Parallel fuction --------------
9696
int ICSegments( char* str_filename, int isTRK, int n_count, int nReplicas, int n_scalars, int n_properties, float* ptrToVOXMM,
97-
double* ptrTDI , double* ptrBlurRho, double* ptrBlurAngle, double* ptrBlurWeights, bool* ptrBlurApplyTo, short* ptrHashTable, char* path_out,
97+
double* ptrTDI , double* ptrBlurRho, double* ptrBlurAngle, double* ptrBlurWeights, bool* ptrBlurApplyTo, short* ptrHashTable, char* path_out,
9898
unsigned long long int offset, int idx, unsigned int startpos, unsigned int endpos );
9999

100100
int ECSegments(float* ptrPEAKS, int Np, float vf_THR, int ECix, int ECiy, int ECiz,
@@ -148,7 +148,7 @@ int trk2dictionary(
148148
batch_size[i%p_size] += 1;
149149
}
150150

151-
151+
152152
// Compute the starting position
153153
// -----------------------------------------
154154
unsigned int elements = threads_count + 1;
@@ -172,7 +172,7 @@ int trk2dictionary(
172172
else if (strcmp(ext,".tck")==0) // for .tck file
173173
isTRK = 0;
174174
else
175-
return 0;
175+
return 0;
176176

177177

178178
// Open tractogram file and compute the offset for each thread
@@ -185,7 +185,7 @@ int trk2dictionary(
185185

186186
FILE* fpTractogram = fopen(str_filename,"rb");
187187
if (fpTractogram == NULL) return 0;
188-
fseek( fpTractogram, data_offset, SEEK_SET ); // skip the header
188+
fseek( fpTractogram, data_offset, SEEK_SET ); // skip the header
189189

190190
OffsetArr[0] = ftell( fpTractogram );
191191

@@ -200,7 +200,7 @@ int trk2dictionary(
200200
f++;
201201
current = ftell( fpTractogram );
202202
for( int i = 1; i < threads_count; i++ ){
203-
if( f == Pos[i] )
203+
if( f == Pos[i] )
204204
OffsetArr[i] = current;
205205
}
206206
}
@@ -215,7 +215,7 @@ int trk2dictionary(
215215
current = ftell( fpTractogram );
216216

217217
for( int i = 1; i < threads_count; i++ ){
218-
if( f == Pos[i] )
218+
if( f == Pos[i] )
219219
OffsetArr[i] = current;
220220
}
221221

@@ -232,7 +232,7 @@ int trk2dictionary(
232232
// ==========================================
233233
// Parallel IC compartments
234234
// ==========================================
235-
235+
236236
std::cout << " * Exporting IC compartments:" << std::endl;
237237
// unsigned int width = 25;
238238
// PROGRESS = new ProgressBar( (unsigned int) n_count, (unsigned int) width);
@@ -244,7 +244,7 @@ int trk2dictionary(
244244
// ---- Original ------
245245
for( int i = 0; i<threads_count; i++ ){
246246
threads.push_back( thread( ICSegments, str_filename, isTRK, n_count, nReplicas, n_scalars, n_properties, ptrToVOXMM,
247-
ptrTDI[i] , ptrBlurRho, ptrBlurAngle, ptrBlurWeights, ptrBlurApplyTo, ptrHashTable, path_out, OffsetArr[i],
247+
ptrTDI[i] , ptrBlurRho, ptrBlurAngle, ptrBlurWeights, ptrBlurApplyTo, ptrHashTable, path_out, OffsetArr[i],
248248
i, Pos[i], Pos[i+1] ) );
249249
}
250250

@@ -410,10 +410,11 @@ int ISOcompartments(double** ptrTDI, char* path_out, int threads){
410410
OUTPUT_path = OUTPUT_path.substr (0,OUTPUT_path.size()-5);
411411
unsigned int totISOVoxels = 0, v=0;
412412

413-
filename = OUTPUT_path+"/dictionary_ISO_v.dict"; FILE* pDict_ISO_v = fopen( filename.c_str(), "wb" );
413+
filename = OUTPUT_path+"/dictionary_ISO_v.dict";
414+
FILE* pDict_ISO_v = fopen( filename.c_str(), "wb" );
414415

415-
int ix, iy, iz, id, atLeastOne;
416-
int skip = 0;
416+
int ix, iy, iz, id, atLeastOne;
417+
int skip = 0;
417418

418419
for(iz=0; iz<dim.z ;iz++){
419420
for(iy=0; iy<dim.y ;iy++)
@@ -437,8 +438,8 @@ int ISOcompartments(double** ptrTDI, char* path_out, int threads){
437438
}
438439
skip = 0;
439440
v = ix + dim.x * ( iy + dim.y * iz );
440-
fwrite( &v, 4, 1, pDict_ISO_v );
441-
totISOVoxels++;
441+
fwrite( &v, 4, 1, pDict_ISO_v );
442+
totISOVoxels++;
442443
}
443444
}
444445
fclose( pDict_ISO_v );
@@ -453,9 +454,9 @@ int ISOcompartments(double** ptrTDI, char* path_out, int threads){
453454
/* Parallel Function */
454455
/********************************************************************************************************************/
455456

456-
int ICSegments( char* str_filename, int isTRK, int n_count, int nReplicas, int n_scalars, int n_properties, float* ptrToVOXMM, double* ptrTDI, double* ptrBlurRho,
457-
double* ptrBlurAngle, double* ptrBlurWeights, bool* ptrBlurApplyTo, short* ptrHashTable, char* path_out,
458-
unsigned long long int offset, int idx, unsigned int startpos, unsigned int endpos )
457+
int ICSegments( char* str_filename, int isTRK, int n_count, int nReplicas, int n_scalars, int n_properties, float* ptrToVOXMM, double* ptrTDI, double* ptrBlurRho,
458+
double* ptrBlurAngle, double* ptrBlurWeights, bool* ptrBlurApplyTo, short* ptrHashTable, char* path_out,
459+
unsigned long long int offset, int idx, unsigned int startpos, unsigned int endpos )
459460
{
460461

461462
// Variables definition
@@ -474,7 +475,7 @@ unsigned long long int offset, int idx, unsigned int startpos, unsigned int endp
474475
map<segInVoxKey,float> FiberNorm;
475476
map<segInVoxKey,float>::iterator itNorm;
476477

477-
segInVoxKey inVoxKey;
478+
segInVoxKey inVoxKey;
478479

479480
P.resize(nReplicas);
480481

@@ -506,8 +507,8 @@ unsigned long long int offset, int idx, unsigned int startpos, unsigned int endp
506507
int incr_old = 0;
507508
// Iterate over streamlines
508509

509-
for(int f=startpos; f<endpos; f++)
510-
{
510+
for(int f=startpos; f<endpos; f++)
511+
{
511512

512513
if ( isTRK )
513514
N = read_fiberTRK( fpTractogram1, fiber, n_scalars, n_properties );
@@ -528,13 +529,13 @@ unsigned long long int offset, int idx, unsigned int startpos, unsigned int endp
528529
{
529530
// NB: please note inverted ordering for 'v'
530531
v = it->first.x + dim.x * ( it->first.y + dim.y * it->first.z );
531-
o = it->first.o;
532+
o = it->first.o;
532533

533534
fwrite( &sumFibers, 4, 1, pDict_IC_f );
534535
fwrite( &v, 4, 1, pDict_IC_v );
535536
fwrite( &o, 2, 1, pDict_IC_o );
536-
fwrite( &(it->second), 4, 1, pDict_IC_len );
537-
537+
fwrite( &(it->second), 4, 1, pDict_IC_len );
538+
538539
ptrTDI[ it->first.z + dim.z * ( it->first.y + dim.y * it->first.x ) ] += it->second;
539540

540541
inVoxKey.set( it->first.x, it->first.y, it->first.z );

0 commit comments

Comments
 (0)