Esto es una idea de Ester Esteban, asi que todo es culpa 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.

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