Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
F
film
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Container Registry
Model registry
Operate
Environments
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
FSL
film
Commits
480a0d9b
Commit
480a0d9b
authored
17 years ago
by
Matthew Webster
Browse files
Options
Downloads
Patches
Plain Diff
In progress Newimage conversion, DO NOT TAG STABLE
parent
d4b72baa
No related branches found
No related tags found
No related merge requests found
Changes
4
Hide whitespace changes
Inline
Side-by-side
Showing
4 changed files
AutoCorrEstimator.cc
+37
-78
37 additions, 78 deletions
AutoCorrEstimator.cc
AutoCorrEstimator.h
+17
-9
17 additions, 9 deletions
AutoCorrEstimator.h
film_gls.cc
+40
-50
40 additions, 50 deletions
film_gls.cc
film_ols.cc
+1
-1
1 addition, 1 deletion
film_ols.cc
with
95 additions
and
138 deletions
AutoCorrEstimator.cc
+
37
−
78
View file @
480a0d9b
...
...
@@ -29,7 +29,6 @@ namespace FILM {
void
AutoCorrEstimator
::
setDesignMatrix
(
const
Matrix
&
dm
)
{
Tracer
tr
(
"AutoCorrEstimator::setDesignMatrix"
);
int
sizeTS
=
xdata
.
getNumVolumes
();
int
numPars
=
dm
.
Ncols
();
dminFFTReal
.
ReSize
(
zeropad
,
numPars
);
...
...
@@ -62,17 +61,15 @@ namespace FILM {
void
AutoCorrEstimator
::
preWhiten
(
ColumnVector
&
in
,
ColumnVector
&
ret
,
int
i
,
Matrix
&
dmret
,
bool
highfreqremovalonly
)
{
Tracer
tr
(
"AutoCorrEstimator::preWhiten"
);
int
sizeTS
=
xdata
.
getNumVolumes
();
int
numPars
=
dminFFTReal
.
getNumSeries
();
ret
.
ReSize
(
sizeTS
);
dmret
.
ReSize
(
sizeTS
,
numPars
);
// FFT auto corr estimate
dummy
=
0
;
vrow
=
0
;
vrow
.
Rows
(
1
,
sizeTS
/
2
)
=
acEst
.
getSeries
(
i
).
Rows
(
1
,
sizeTS
/
2
);
vrow
.
Rows
(
zeropad
-
sizeTS
/
2
+
2
,
zeropad
)
=
acEst
.
getSeries
(
i
).
Rows
(
2
,
sizeTS
/
2
).
Reverse
();
vrow
.
Rows
(
1
,
sizeTS
/
2
)
=
acEst
.
Column
(
i
).
Rows
(
1
,
sizeTS
/
2
);
vrow
.
Rows
(
zeropad
-
sizeTS
/
2
+
2
,
zeropad
)
=
acEst
.
Column
(
i
).
Rows
(
2
,
sizeTS
/
2
).
Reverse
();
FFT
(
vrow
,
dummy
,
ac_fft_real
,
ac_fft_im
);
...
...
@@ -151,9 +148,6 @@ namespace FILM {
cerr
<<
"Prewhitening... "
;
int
sizeTS
=
xdata
.
getNumVolumes
();
int
numTS
=
xdata
.
getNumSeries
();
ret
.
ReSize
(
sizeTS
,
numTS
);
// make sure p_vrow is cyclic (even function)
...
...
@@ -172,8 +166,8 @@ namespace FILM {
// FFT auto corr estimate
dummy
=
0
;
vrow
=
0
;
vrow
.
Rows
(
1
,
sizeTS
/
2
)
=
acEst
.
getSeries
(
i
).
Rows
(
1
,
sizeTS
/
2
);
vrow
.
Rows
(
zeropad
-
sizeTS
/
2
+
2
,
zeropad
)
=
acEst
.
getSeries
(
i
).
Rows
(
2
,
sizeTS
/
2
).
Reverse
();
vrow
.
Rows
(
1
,
sizeTS
/
2
)
=
acEst
.
Column
(
i
).
Rows
(
1
,
sizeTS
/
2
);
vrow
.
Rows
(
zeropad
-
sizeTS
/
2
+
2
,
zeropad
)
=
acEst
.
Column
(
i
).
Rows
(
2
,
sizeTS
/
2
).
Reverse
();
FFT
(
vrow
,
dummy
,
ac_fft_real
,
ac_fft_im
);
// FFT x data
...
...
@@ -213,7 +207,7 @@ namespace FILM {
cerr
<<
" Completed"
<<
endl
;
}
void
AutoCorrEstimator
::
fitAutoRegressiveModel
()
Matrix
AutoCorrEstimator
::
fitAutoRegressiveModel
()
{
Tracer
trace
(
"AutoCorrEstimator::fitAutoRegressiveModel"
);
...
...
@@ -221,9 +215,6 @@ namespace FILM {
const
int
maxorder
=
15
;
const
int
minorder
=
1
;
int
sizeTS
=
xdata
.
getNumVolumes
();
int
numTS
=
xdata
.
getNumSeries
();
// setup temp variables
ColumnVector
x
(
sizeTS
);
...
...
@@ -236,7 +227,7 @@ namespace FILM {
for
(
int
i
=
1
;
i
<=
numTS
;
i
++
)
{
x
=
xdata
.
getSeries
(
i
).
As
Column
();
x
=
xdata
.
Column
(
i
);
ColumnVector
betastmp
;
order
(
i
)
=
pacf
(
x
,
minorder
,
maxorder
,
betastmp
);
...
...
@@ -269,8 +260,6 @@ namespace FILM {
else
Kinv
.
SubMatrix
(
j
,
j
,
1
,
j
)
=
Krow
.
Rows
(
sizeTS
-
j
+
1
,
sizeTS
).
t
();
}
//MISCMATHS::write_ascii_matrix(Kinv,"Kinv");
// Kinv now becomes V:
if
(
order
(
i
)
!=
1
)
...
...
@@ -292,21 +281,15 @@ namespace FILM {
// output betas:
write_ascii_matrix
(
LogSingleton
::
getInstance
().
appendDir
(
"betas"
),
betas
);
VolumeInfo
vinfo
=
xdata
.
getInfo
();
vinfo
.
v
=
maxorder
;
betas
.
unthresholdSeries
(
vinfo
,
xdata
.
getPreThresholdPositions
());
betas
.
writeAsFloat
(
LogSingleton
::
getInstance
().
getDir
()
+
"/betas"
);
countLargeE
=
0
;
cerr
<<
" Completed"
<<
endl
;
return
(
betas
);
}
int
AutoCorrEstimator
::
pacf
(
const
ColumnVector
&
x
,
int
minorder
,
int
maxorder
,
ColumnVector
&
betas
)
{
Tracer
ts
(
"pacf"
);
int
order
=
-
1
;
int
sizeTS
=
x
.
Nrows
();
// Set c
Matrix
c
(
1
,
1
);
...
...
@@ -342,11 +325,10 @@ namespace FILM {
return
order
;
}
int
AutoCorrEstimator
::
establishUsanThresh
(
const
V
olum
e
&
epivol
)
int
AutoCorrEstimator
::
establishUsanThresh
(
const
C
olum
nVector
&
epivol
)
{
int
usanthresh
=
100
;
int
num
=
epivol
.
getVolumeSize
();
int
num
=
epivol
.
Nrows
();
Histogram
hist
(
epivol
,
num
/
200
);
hist
.
generate
();
float
mode
=
hist
.
mode
();
...
...
@@ -358,7 +340,7 @@ namespace FILM {
// Work out standard deviation from mode for values greater than mode:
for
(
int
i
=
1
;
i
<=
num
;
i
++
)
{
if
(
epivol
(
i
)
>
mode
)
{
sum
=
sum
+
(
epivol
(
i
)
-
mode
)
*
(
epivol
(
i
)
-
mode
);
sum
+
=
(
epivol
(
i
)
-
mode
)
*
(
epivol
(
i
)
-
mode
);
count
++
;
}
}
...
...
@@ -371,54 +353,39 @@ namespace FILM {
return
usanthresh
;
}
void
AutoCorrEstimator
::
spatiallySmooth
(
const
string
&
usanfname
,
const
V
olum
e
&
epivol
,
int
masksize
,
const
string
&
epifname
,
const
string
&
susanpath
,
int
usan_thresh
,
int
lag
)
{
void
AutoCorrEstimator
::
spatiallySmooth
(
const
string
&
usanfname
,
const
C
olum
nVector
&
epivol
,
int
masksize
,
const
string
&
epifname
,
int
usan_thresh
,
int
lag
)
{
Tracer
trace
(
"AutoCorrEstimator::spatiallySmooth"
);
if
(
xdata
.
getNumSeries
()
<=
1
)
if
(
numTS
<=
1
)
{
cerr
<<
"Warning: Number of voxels = "
<<
xdata
.
getNumSeries
()
<<
". Spatial smoothing of autocorrelation estimates is not carried out"
<<
endl
;
cerr
<<
"Warning: Number of voxels = "
<<
numTS
<<
". Spatial smoothing of autocorrelation estimates is not carried out"
<<
endl
;
}
else
{
if
(
lag
==
0
)
lag
=
MISCMATHS
::
Min
(
40
,
int
(
xdata
.
getNumVolumes
()
/
4
));
if
(
usan_thresh
==
0
)
{
// Establish epi thresh to use:
usan_thresh
=
establishUsanThresh
(
epivol
);
}
lag
=
MISCMATHS
::
Min
(
40
,
int
(
sizeTS
/
4
));
volume
<
float
>
susan_vol
(
xdata
.
getInfo
().
x
,
xdata
.
getInfo
().
y
,
xdata
.
getInfo
().
z
);
volume
<
float
>
usan_area
(
xdata
.
getInfo
().
x
,
xdata
.
getInfo
().
y
,
xdata
.
getInfo
().
z
);
if
(
usan_thresh
==
0
)
usan_thresh
=
establishUsanThresh
(
epivol
);
// Establish epi thresh to use:
volume4D
<
float
>
susan_vol
(
mask
.
xsize
(),
mask
.
ysize
(),
mask
.
zsize
(),
1
);
volume
<
float
>
usan_area
(
mask
.
xsize
(),
mask
.
ysize
(),
mask
.
zsize
());
volume
<
float
>
usan_vol
;
read_volume
(
usan_vol
,
usanfname
);
volume
<
float
>
kernel
;
kernel
=
gaussian_kernel3D
(
masksize
,
xdata
.
getInfo
().
vx
,
xdata
.
getInfo
().
vy
,
xdata
.
getInfo
().
vz
,
2.0
);
// Setup volume for reading and writing volumes:
Volume
vol
(
acEst
.
getNumSeries
(),
xdata
.
getInfo
(),
xdata
.
getPreThresholdPositions
());
int
i
=
2
;
kernel
=
gaussian_kernel3D
(
masksize
,
mask
.
xdim
(),
mask
.
ydim
(),
mask
.
zdim
(),
2.0
);
int
factor
=
10000
;
cerr
<<
"Spatially smoothing auto corr estimates"
<<
endl
;
for
(;
i
<=
lag
;
i
++
)
for
(
int
i
=
2
;
i
<=
lag
;
i
++
)
{
// setup susan input
vol
=
acEst
.
getVolume
(
i
).
AsColumn
()
*
factor
;
vol
.
unthreshold
();
susan_vol
.
insert_vec
(
vol
);
// call susan
susan_vol
=
susan_convolve
(
susan_vol
,
kernel
,
1
,
0
,
1
,
&
usan_area
,
usan_vol
,
usan_thresh
*
usan_thresh
);
susan_vol
.
setmatrix
(
acEst
.
Row
(
i
),
mask
);
susan_vol
*=
factor
;
susan_vol
[
0
]
=
susan_convolve
(
susan_vol
[
0
],
kernel
,
1
,
0
,
1
,
&
usan_area
,
usan_vol
,
usan_thresh
*
usan_thresh
);
// insert output back into acEst
vol
=
susan_vol
.
vec
();
vol
.
threshold
();
acEst
.
setVolume
(
static_cast
<
RowVector
>
((
vol
/
factor
).
AsRow
()),
i
);
susan_vol
/=
factor
;
acEst
.
Row
(
i
)
=
susan_vol
.
matrix
(
mask
);
cerr
<<
"."
;
}
...
...
@@ -453,14 +420,13 @@ namespace FILM {
ColumnVector
fft_im
;
ColumnVector
dummy
(
zeropad
);
ColumnVector
realifft
(
zeropad
);
int
sizeTS
=
xdata
.
getNumVolumes
();
for
(
int
i
=
1
;
i
<=
xdata
.
getNumSeries
()
;
i
++
)
for
(
int
i
=
1
;
i
<=
numTS
;
i
++
)
{
dummy
=
0
;
vrow
=
0
;
vrow
.
Rows
(
1
,
sizeTS
/
2
)
=
acEst
.
getSeries
(
i
).
Rows
(
1
,
sizeTS
/
2
);
vrow
.
Rows
(
zeropad
-
sizeTS
/
2
+
2
,
zeropad
)
=
acEst
.
getSeries
(
i
).
Rows
(
2
,
sizeTS
/
2
).
Reverse
();
vrow
.
Rows
(
1
,
sizeTS
/
2
)
=
acEst
.
Column
(
i
).
Rows
(
1
,
sizeTS
/
2
);
vrow
.
Rows
(
zeropad
-
sizeTS
/
2
+
2
,
zeropad
)
=
acEst
.
Column
(
i
).
Rows
(
2
,
sizeTS
/
2
).
Reverse
();
FFT
(
vrow
,
dummy
,
fft_real
,
fft_im
);
...
...
@@ -478,8 +444,6 @@ namespace FILM {
Tracer
tr
(
"AutoCorrEstimator::multitaper"
);
cerr
<<
"Multitapering... "
;
int
sizeTS
=
xdata
.
getNumVolumes
();
int
numTS
=
xdata
.
getNumSeries
();
Matrix
slepians
;
getSlepians
(
M
,
sizeTS
,
slepians
);
...
...
@@ -504,10 +468,8 @@ namespace FILM {
// Compute FFT for each slepian taper
for
(
int
k
=
1
;
k
<=
slepians
.
Ncols
();
k
++
)
{
x
.
Rows
(
1
,
sizeTS
)
=
SP
(
slepians
.
Column
(
k
),
xdata
.
getSeries
(
i
));
//LogSingleton::getInstance().out("x", xdata.getSeries(i), false);
x
.
Rows
(
1
,
sizeTS
)
=
SP
(
slepians
.
Column
(
k
),
xdata
.
Column
(
i
));
FFT
(
x
,
dummy
,
fft_real
,
fft_im
);
for
(
int
j
=
1
;
j
<=
zeropad
;
j
++
)
{
...
...
@@ -532,7 +494,7 @@ namespace FILM {
//LogSingleton::getInstance().out("fftreal", fft_real);
float
varx
=
MISCMATHS
::
var
(
ColumnVector
(
x
.
Rows
(
1
,
sizeTS
))).
AsScalar
();
acEst
.
setSeries
(
realifft
.
Rows
(
1
,
sizeTS
)
/
varx
,
i
)
;
acEst
.
Column
(
i
)
=
realifft
.
Rows
(
1
,
sizeTS
)
/
varx
;
}
countLargeE
=
0
;
cerr
<<
"Completed"
<<
endl
;
...
...
@@ -571,7 +533,6 @@ namespace FILM {
cerr
<<
"Tukey M = "
<<
M
<<
endl
;
cerr
<<
"Tukey estimates... "
;
int
sizeTS
=
xdata
.
getNumVolumes
();
ColumnVector
window
(
M
);
...
...
@@ -580,7 +541,7 @@ namespace FILM {
window
(
j
)
=
0.5
*
(
1
+
cos
(
M_PI
*
j
/
(
float
(
M
))));
}
for
(
int
i
=
1
;
i
<=
xdata
.
getNumSerie
s
();
i
++
)
{
for
(
int
i
=
1
;
i
<=
xdata
.
Ncol
s
();
i
++
)
{
acEst
.
SubMatrix
(
1
,
M
,
i
,
i
)
=
SP
(
acEst
.
SubMatrix
(
1
,
M
,
i
,
i
),
window
);
acEst
.
SubMatrix
(
M
+
1
,
sizeTS
,
i
,
i
)
=
0
;
...
...
@@ -596,8 +557,7 @@ namespace FILM {
cerr
<<
"Using New PAVA on AutoCorr estimates... "
;
for
(
int
i
=
1
;
i
<=
xdata
.
getNumSeries
();
i
++
)
{
int
sizeTS
=
xdata
.
getNumVolumes
();
for
(
int
i
=
1
;
i
<=
numTS
;
i
++
)
{
int
stopat
=
(
int
)
sizeTS
/
2
;
// 5% point of distribution of autocorr about zero
...
...
@@ -683,9 +643,8 @@ namespace FILM {
cerr
<<
"Applying constraints to AutoCorr estimates... "
;
for
(
int
i
=
1
;
i
<=
xdata
.
getNumSeries
()
;
i
++
)
for
(
int
i
=
1
;
i
<=
numTS
;
i
++
)
{
int
sizeTS
=
xdata
.
getNumVolumes
();
int
j
=
3
;
int
stopat
=
(
int
)
sizeTS
/
4
;
...
...
@@ -740,11 +699,11 @@ namespace FILM {
{
Tracer
tr
(
"AutoCorrEstimator::getMeanEstimate"
);
ret
.
ReSize
(
acEst
.
getNumVolume
s
());
ret
.
ReSize
(
acEst
.
Nrow
s
());
// Calc global Vrow:
for
(
int
i
=
1
;
i
<=
acEst
.
getNumVolume
s
();
i
++
)
for
(
int
i
=
1
;
i
<=
acEst
.
Nrow
s
();
i
++
)
{
ret
(
i
)
=
MISCMATHS
::
mean
(
ColumnVector
(
acEst
.
getVolume
(
i
).
AsColumn
())).
AsScalar
();
ret
(
i
)
=
MISCMATHS
::
mean
(
ColumnVector
(
acEst
.
Row
(
i
).
AsColumn
())).
AsScalar
();
}
}
}
...
...
This diff is collapsed.
Click to expand it.
AutoCorrEstimator.h
+
17
−
9
View file @
480a0d9b
...
...
@@ -14,6 +14,7 @@
#define WANT_STREAM
#define WANT_MATH
#include
"newimage/newimageall.h"
#include
"miscmaths/volumeseries.h"
#include
"miscmaths/volume.h"
#include
"newmatap.h"
...
...
@@ -22,15 +23,18 @@
using
namespace
NEWMAT
;
using
namespace
MISCMATHS
;
using
namespace
NEWIMAGE
;
namespace
FILM
{
class
AutoCorrEstimator
{
public:
AutoCorrEstimator
(
const
VolumeSeries
&
pxdata
)
:
AutoCorrEstimator
(
const
Matrix
&
pxdata
)
:
sizeTS
(
pxdata
.
Nrows
()),
numTS
(
pxdata
.
Ncols
()),
xdata
(
pxdata
),
acEst
(
pxdata
.
getNumVolumes
(),
pxdata
.
getNumSeries
(),
pxdata
.
getInfo
(),
pxdata
.
getPreThresholdPosition
s
()),
acEst
(
pxdata
.
Nrows
(),
pxdata
.
Ncol
s
()),
vrow
(),
xrow
(),
dummy
(),
...
...
@@ -38,7 +42,7 @@ namespace FILM {
dm_mn
(),
zeropad
(
0
)
{
zeropad
=
MISCMATHS
::
nextpow2
(
pxdata
.
getNumVolume
s
());
zeropad
=
MISCMATHS
::
nextpow2
(
pxdata
.
Nrow
s
());
vrow
.
ReSize
(
zeropad
);
xrow
.
ReSize
(
zeropad
);
dummy
.
ReSize
(
zeropad
);
...
...
@@ -50,19 +54,19 @@ namespace FILM {
}
void
calcRaw
(
int
lag
=
0
);
void
spatiallySmooth
(
const
string
&
usanfname
,
const
V
olum
e
&
epivol
,
int
masksize
,
const
string
&
epifname
,
const
string
&
susanpath
,
int
usanthresh
,
int
lag
=
0
);
void
spatiallySmooth
(
const
string
&
usanfname
,
const
C
olum
nVector
&
epivol
,
int
masksize
,
const
string
&
epifname
,
int
usanthresh
,
int
lag
=
0
);
void
applyConstraints
();
void
filter
(
const
ColumnVector
&
filterFFT
);
void
fitAutoRegressiveModel
();
Matrix
fitAutoRegressiveModel
();
void
pava
();
void
preWhiten
(
VolumeSeries
&
in
,
VolumeSeries
&
ret
);
void
preWhiten
(
ColumnVector
&
in
,
ColumnVector
&
ret
,
int
i
,
Matrix
&
dmret
,
bool
highfreqremovalonly
=
false
);
void
setDesignMatrix
(
const
Matrix
&
dm
);
int
establishUsanThresh
(
const
V
olum
e
&
epivol
);
int
establishUsanThresh
(
const
C
olum
nVector
&
epivol
);
void
getMeanEstimate
(
ColumnVector
&
ret
);
VolumeSeries
&
getEstimates
()
{
return
acEst
;
}
Matrix
&
getEstimates
()
{
return
acEst
;
}
VolumeSeries
&
getE
()
{
return
E
;
}
ColumnVector
&
getCountLargeE
(){
return
countLargeE
;
}
...
...
@@ -70,15 +74,19 @@ namespace FILM {
void
tukey
(
int
M
);
void
multitaper
(
int
M
);
int
pacf
(
const
ColumnVector
&
x
,
int
minorder
,
int
maxorder
,
ColumnVector
&
betas
);
volume
<
float
>
mask
;
private
:
const
int
sizeTS
;
const
int
numTS
;
AutoCorrEstimator
();
const
AutoCorrEstimator
&
operator
=
(
AutoCorrEstimator
&
);
AutoCorrEstimator
(
AutoCorrEstimator
&
);
void
getSlepians
(
int
M
,
int
sizeTS
,
Matrix
&
slepians
);
const
VolumeSeries
&
xdata
;
VolumeSeries
acEst
;
const
Matrix
&
xdata
;
Matrix
acEst
;
VolumeSeries
E
;
ColumnVector
countLargeE
;
...
...
This diff is collapsed.
Click to expand it.
film_gls.cc
+
40
−
50
View file @
480a0d9b
...
...
@@ -44,16 +44,24 @@ int main(int argc, char *argv[])
globalopts
.
parse_command_line
(
argc
,
argv
,
logger
);
// load data
volume4D
<
double
>
input_data
;
read_volume4D
(
input_data
,
globalopts
.
inputfname
);
volume4D
<
float
>
input_data
;
volumeinfo
vinfo
;
read_volume4D
(
input_data
,
globalopts
.
inputfname
,
vinfo
);
int
sizeTS
=
input_data
.
tsize
();
Matrix
datam
;
VolumeSeries
x
,
y
;
y
.
read
(
globalopts
.
inputfname
);
x
.
setInfo
(
y
.
getInfo
());
x
=
input_data
.
matrix
();
x
.
thresholdSeries
(
globalopts
.
thresh
,
true
);
//save_volume_dtype(input_data[int(sizeTS/2)-1],logger.getDir() + "/" + globalopts.epifname,DT_SIGNED_INT );
save_volume
(
input_data
[
int
(
sizeTS
/
2
)
-
1
],
logger
.
getDir
()
+
"/"
+
globalopts
.
epifname
);
volume
<
double
>
mask
(
input_data
[
0
]);
volume
<
float
>
mask
(
input_data
[
0
]);
volume4D
<
float
>
epivol2
(
mask
.
xsize
(),
mask
.
ysize
(),
mask
.
zsize
(),
1
);
epivol2
[
0
]
=
input_data
[
int
(
sizeTS
/
2
)
-
1
];
mask
=
0
;
for
(
int
t
=
0
;
t
<
input_data
.
tsize
();
t
++
)
mask
+=
input_data
[
t
];
mask
/=
input_data
.
tsize
();
...
...
@@ -64,32 +72,9 @@ int main(int argc, char *argv[])
for
(
int
t
=
0
;
t
<
input_data
.
tsize
();
t
++
)
input_data
(
x
,
y
,
z
,
t
)
-=
mask
(
x
,
y
,
z
);
mask
.
binarise
(
globalopts
.
thresh
,
mask
.
max
()
+
1
,
exclusive
);
input_data
*=
mask
;
save_volume4D
(
input_data
,
"foo2"
);
datam
=
input_data
.
matrix
(
mask
);
int
sizeTS
=
input_data
.
tsize
();
int
numTS
=
mask
.
sum
();
cerr
<<
datam
.
Nrows
()
<<
" "
<<
datam
.
Ncols
()
<<
endl
;
// if needed for SUSAN spatial smoothing, output the halfway volume for use later
Volume
epivol
;
if
(
globalopts
.
smoothACEst
)
{
epivol
=
x
.
getVolume
(
int
(
sizeTS
/
2
)).
AsColumn
();
epivol
.
setInfo
(
x
.
getInfo
());
epivol
.
writeAsInt
(
logger
.
getDir
()
+
"/"
+
globalopts
.
epifname
);
}
// This also removes the mean from each of the time series:
x
.
thresholdSeries
(
globalopts
.
thresh
,
true
);
x
.
writeThresholdedSeriesAsFloat
(
x
.
getInfo
(),
x
.
getPreThresholdPositions
(),
"foo3"
);
// if needed for SUSAN spatial smoothing, also threshold the epi volume
if
(
globalopts
.
smoothACEst
)
{
epivol
.
setPreThresholdPositions
(
x
.
getPreThresholdPositions
());
epivol
.
threshold
();
}
// Load paradigm:
Paradigm
parad
;
...
...
@@ -122,16 +107,18 @@ int main(int argc, char *argv[])
// Residuals container:
VolumeSeries
residuals
(
sizeTS
,
numTS
,
x
.
getInfo
(),
x
.
getPreThresholdPositions
());
Matrix
residualsm
=
datam
;
// Setup autocorrelation estimator:
AutoCorrEstimator
acEst
(
residuals
);
acEst
.
mask
=
mask
;
if
(
!
globalopts
.
noest
)
{
cout
<<
"Calculating residuals..."
<<
endl
;
for
(
int
i
=
1
;
i
<=
numTS
;
i
++
)
{
//glimGls.setData(x.getSeries(i), parad.getDesignMatrix(), i);
glimGls
.
setData
(
datam
.
Column
(
i
),
parad
.
getDesignMatrix
(),
i
);
residuals
.
setSeries
(
glimGls
.
getResiduals
(),
i
);
}
...
...
@@ -141,7 +128,11 @@ int main(int argc, char *argv[])
if
(
globalopts
.
fitAutoRegressiveModel
)
{
acEst
.
fitAutoRegressiveModel
();
volume4D
<
float
>
beta
;
beta
.
setmatrix
(
acEst
.
fitAutoRegressiveModel
(),
mask
);
copybasicproperties
(
input_data
,
beta
);
FslSetCalMinMax
(
&
vinfo
,
beta
.
min
(),
beta
.
max
());
save_volume4D
(
beta
,
LogSingleton
::
getInstance
().
getDir
()
+
"/betas"
,
vinfo
);
}
else
if
(
globalopts
.
tukey
)
{
...
...
@@ -153,7 +144,8 @@ int main(int argc, char *argv[])
// SUSAN smooth raw estimates:
if
(
globalopts
.
smoothACEst
)
{
acEst
.
spatiallySmooth
(
logger
.
getDir
()
+
"/"
+
globalopts
.
epifname
,
epivol
,
globalopts
.
ms
,
globalopts
.
epifname
,
globalopts
.
susanpath
,
globalopts
.
epith
,
globalopts
.
tukeysize
);
ColumnVector
epivol
=
epivol2
.
matrix
(
mask
).
t
();
acEst
.
spatiallySmooth
(
logger
.
getDir
()
+
"/"
+
globalopts
.
epifname
,
epivol
,
globalopts
.
ms
,
globalopts
.
epifname
,
globalopts
.
epith
,
globalopts
.
tukeysize
);
}
acEst
.
tukey
(
globalopts
.
tukeysize
);
...
...
@@ -170,7 +162,8 @@ int main(int argc, char *argv[])
// Smooth raw estimates:
if
(
globalopts
.
smoothACEst
)
{
acEst
.
spatiallySmooth
(
logger
.
getDir
()
+
"/"
+
globalopts
.
epifname
,
epivol
,
globalopts
.
ms
,
globalopts
.
epifname
,
globalopts
.
susanpath
,
globalopts
.
epith
);
ColumnVector
epivol
=
epivol2
.
matrix
(
mask
).
t
();
acEst
.
spatiallySmooth
(
logger
.
getDir
()
+
"/"
+
globalopts
.
epifname
,
epivol
,
globalopts
.
ms
,
globalopts
.
epifname
,
globalopts
.
epith
);
}
acEst
.
pava
();
...
...
@@ -196,7 +189,7 @@ int main(int argc, char *argv[])
acEst
.
setDesignMatrix
(
parad
.
getDesignMatrix
());
// Use autocorr estimate to prewhiten data:
xprew
=
x
.
getSeries
(
i
);
xprew
=
datam
.
Column
(
i
);
Matrix
designmattw
;
acEst
.
preWhiten
(
xprew
,
xw
,
i
,
designmattw
);
...
...
@@ -207,9 +200,8 @@ int main(int argc, char *argv[])
co
++
;
}
x
.
setSeries
(
xw
,
i
);
glimGls
.
setData
(
x
.
getSeries
(
i
),
designmattw
,
i
);
datam
.
Column
(
i
)
=
xw
;
glimGls
.
setData
(
datam
.
Column
(
i
),
designmattw
,
i
);
if
(
globalopts
.
output_pwdata
||
globalopts
.
verbose
)
{
...
...
@@ -224,8 +216,7 @@ int main(int argc, char *argv[])
cout
.
flush
();
co
++
;
}
glimGls
.
setData
(
x
.
getSeries
(
i
),
parad
.
getDesignMatrix
(),
i
);
glimGls
.
setData
(
datam
.
Column
(
i
),
parad
.
getDesignMatrix
(),
i
);
residuals
.
setSeries
(
glimGls
.
getResiduals
(),
i
);
if
(
globalopts
.
output_pwdata
||
globalopts
.
verbose
)
...
...
@@ -245,16 +236,15 @@ int main(int argc, char *argv[])
cerr
<<
"Saving results... "
<<
endl
;
// write residuals
//residuals.unthresholdSeries(x.getInfo(),x.getPreThresholdPositions());
//residuals.writeAsFloat(logger.getDir() + "/res4d");
residuals
.
writeThresholdedSeriesAsFloat
(
x
.
getInfo
(),
x
.
getPreThresholdPositions
(),
logger
.
getDir
()
+
"/res4d"
);
residuals
.
CleanUp
();
if
(
globalopts
.
output_pwdata
||
globalopts
.
verbose
)
{
// Write out whitened data
x
.
writeThresholdedSeriesAsFloat
(
x
.
getInfo
(),
x
.
getPreThresholdPositions
(),
logger
.
getDir
()
+
"/prewhitened_data"
);
input_data
.
setmatrix
(
datam
,
mask
);
FslSetCalMinMax
(
&
vinfo
,
input_data
.
min
(),
input_data
.
max
());
save_volume4D
(
input_data
,
logger
.
getDir
()
+
"/prewhitened_data"
,
vinfo
);
// Write out whitened design matrix
write_vest
(
logger
.
appendDir
(
"mean_prewhitened_dm.mat"
),
mean_prewhitened_dm
);
...
...
@@ -263,19 +253,19 @@ int main(int argc, char *argv[])
x
.
CleanUp
();
// Write out threshac:
VolumeSeries
&
threshac
=
acEst
.
getEstimates
();
Matrix
&
threshac
m
=
acEst
.
getEstimates
();
int
cutoff
=
sizeTS
/
2
;
if
(
globalopts
.
tukey
)
cutoff
=
globalopts
.
tukeysize
;
cutoff
=
MISCMATHS
::
Max
(
1
,
cutoff
);
threshac
=
threshac
.
Rows
(
1
,
cutoff
);
V
olume
Info
volinfo
=
x
.
getInfo
()
;
volinfo
.
v
=
cutoff
;
volinfo
.
intent_code
=
NIFTI_INTENT_ESTIMATE
;
threshac
.
unthresholdSeries
(
volinfo
,
x
.
getPreThresholdPositions
()
);
threshac
.
writeAsFloat
(
logger
.
getDir
()
+
"/
threshac
1"
);
threshac
.
thresholdSeries
(
);
threshac
.
CleanUp
();
threshac
m
=
threshac
m
.
Rows
(
1
,
cutoff
);
v
olume
4D
<
float
>
threshac
;
threshac
.
setmatrix
(
threshacm
,
mask
)
;
copybasicproperties
(
input_data
,
threshac
)
;
threshac
.
set_intent
(
NIFTI_INTENT_ESTIMATE
,
0
,
0
,
0
);
FslSetCalMinMax
(
&
vinfo
,
threshac
.
min
(),
threshac
.
max
()
);
save_volume4D
(
threshac
,
logger
.
getDir
()
+
"/threshac1"
,
vinfo
);
threshac
m
.
CleanUp
();
// save gls results:
glimGls
.
Save
(
x
.
getInfo
(),
x
.
getPreThresholdPositions
());
...
...
This diff is collapsed.
Click to expand it.
film_ols.cc
+
1
−
1
View file @
480a0d9b
...
...
@@ -132,7 +132,7 @@ int main(int argc, char *argv[])
// Smooth raw estimates:
if
(
globalopts
.
smoothACEst
)
{
acEst
.
spatiallySmooth
(
logger
.
getDir
()
+
"/"
+
globalopts
.
epifname
,
epivol
,
globalopts
.
ms
,
globalopts
.
epifname
,
globalopts
.
susanpath
,
0
);
acEst
.
spatiallySmooth
(
logger
.
getDir
()
+
"/"
+
globalopts
.
epifname
,
epivol
,
globalopts
.
ms
,
globalopts
.
epifname
,
0
);
}
// Apply constraints to estimate autocorr:
...
...
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment