Esto es una idea de Ester Esteban, asi que toda la culpa es suya.

Siguiendo Zeestraten 2017 Neurology queremos determinar la altura del pico de los histogramas normalizados de MD.

Se supone que ya tenemos la imagen MD por el preprocesamiento previo. Ahora hemos de sacar el histograma de esta imagen,

[osotolongo@brick03 nph]$ fslstats 0005_dti_MD.nii.gz -l 0 -u 0.004 -h 1000 > hist.txt; sed -i "1ibins" hist.txt

Esto produce un perfil de los valores de MD bastante claro, con un pico muy determinado.

Obviamente, si cambiamos el numero de bins, la altura cambia,

Asi que vamos a normalizar por el area bajo la curva (A.K.A. densidad)

La altura de este pico es lo que vamos a tomar como el valor de MD-NPH,

> max(h$bins, na.rm = TRUE)/sum(h$bins, na.rm=TRUE)
[1] 0.08052718

Pero con una salvedad, no voy a normalizar por el area la curva del histograma individual sino por el valor medio de esta area entre todos los histogramas, para que todos los valores de densidad tengan un mismo nivel de referencia.

Voy a,

1. reutilizar los componentes del pipeline
2. seguir el paper en cuanto a los valores de threshold y demas
3. usar List::Util para manipular los arrays

asi que todo deberia ser sencillo.

Ok. Pero no taaaan sencillo.

Antes que nada, no voy a sacar el valor medio de toda la imagen MD sino de la porcion correspondiente a la WM. Es decir, que tengo que empezar por hacer una mascara de la WM. Es to no deberia ser dificil si reutilizo el aseg de FS,

get_wm_mask.sh

#!/bin/sh
subject=$1
shift
 
tmp_dir=$1
shift
 
if [ ! -f ${tmp_dir}/register.dat ]; then
        if [  ! -d ${tmp_dir} ]; then mkdir -p ${tmp_dir}; fi;
        tkregister2 --mov $SUBJECTS_DIR/${subject}/mri/rawavg.mgz --noedit --s ${subject} --regheader --reg ${tmp_dir}/register.dat;
fi
if [ ! -f ${tmp_dir}/all_aseg.nii.gz ]; then
        mri_label2vol --seg $SUBJECTS_DIR/${subject}/mri/aseg.mgz --temp $SUBJECTS_DIR/${subject}/mri/rawavg.mgz --o ${tmp_dir}/all_aseg.nii.gz --reg ${tmp_dir}/register.dat;
fi
${FSLDIR}/bin/fslmaths ${tmp_dir}/all_aseg.nii.gz -uthr 2  -thr 2 -div 2 ${tmp_dir}/lhwm
${FSLDIR}/bin/fslmaths ${tmp_dir}/all_aseg.nii.gz -uthr 41  -thr 41 -div 41 ${tmp_dir}/rhwm
${FSLDIR}/bin/fslmaths ${tmp_dir}/lhwm -add ${tmp_dir}/rhwm -bin ${tmp_dir}/wm_mask.nii.gz
~                                                                                                                      

Pero ojo, que esa mascara esta en espacio nativo T1, asi que tendria que hacer algo como,

antsApplyTransforms -d 3 -i ${tmp_dir}/wm_mask.nii.gz -r ${tmp_dir}/hifi_b0.nii.gz -t ${tmp_dir}/$subject_epi_reg_ANTS.mat -o ${tmp_dir}/wm_tmp.nii.gz

y ahora a esto hay que hacerle un treshold,

fslmaths ${tmp_dir}/wm_tmp.nii.gz -thr 0.5 ${tmp_dir}/wm_mask.nii.gz

La mascara resultante hay que aplicarla a la imagen MD y a esta ultima imagen es a la que hay que sacarle el histograma.

Primero leo la DB del proyecto y por cada unos de las imagenes MD que he llegado a procesar ejecuto el procedimiento anterior,

foreach my $subject (@dtis){
        my $dti_md = $w_dir.'/'.$subject.'_dti_MD.nii.gz';
        if($subject and -e $dti_md){
                my $order = $ENV{'FSLDIR'}.'/bin/fslstats '.$dti_md.' -l 0 -u '.$up_threshold.' -h '.$nbins;
                my @shist = qx/$order/;
                chomp @shist;
                push @allsums, sum(@shist);
                $nph{$subject} = max(@shist);
        }
}

ahora, el valor para normalizar es el medio de todas las areas,

my $norm = sum(@allsums)/scalar(@allsums);

y saco todo por STDOUT,

print "Subject,NPH\n";
foreach my $subject (sort keys %nph){
        my $nphv = $nph{$subject}/$norm;
        print "$subject,$nphv\n";
}

y el codigo completo online en https://github.com/asqwerty666/acenip/blob/main/tools/nph.pl

Esto es bastante facil de usar,

[osotolongo@brick03 nph]$ ./nph.pl bioface > nph.csv
[osotolongo@brick03 nph]$ sed 's/;/,/;1iSubject,Subject_ID' ../bioface_mri.csv > bioface_codes.csv
[osotolongo@brick03 nph]$ join -t, -j 1 bioface_codes.csv nph.csv | awk -F"," '{print $2","$3}' | csvsort -c 1 > normalised_peak_height.csv
[osotolongo@brick03 nph]$ head normalised_peak_height.csv
Subject_ID,NPH
B002,0.0731249323687369
B003,0.0732767676009448
B004,0.0745456763272544
B005,0.0746812434988687
B006,0.0697899799470258
B007,0.0788675577583173
B008,0.0932810594443464
B009,0.0787428359604322
B010,0.0681414831401963

Nota: Para simplificar he usado csvsort de csvkit pero tambien se puede dar una vuelta y usar gnu sort.