Skip to content
Snippets Groups Projects
Commit 409d65b7 authored by Mark Jenkinson's avatar Mark Jenkinson
Browse files

Working p-thresholding cluster program for FEAT 4 

Must tidy up options and general coding later
parent 114ee4a1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 159 additions and 30 deletions
cluster.cc
+ 159
− 30
View file @ 409d65b7
......@@ -6,17 +6,19 @@
/* CCOPYRIGHT */
#include "fmribmain.h"
#include "newimageall.h"
#include <vector>
#include <algorithm>
#include "fmribmain.h"
#include "newimageall.h"
#include "infer.h"
using namespace NEWIMAGE;
void print_usage(char *argv[])
{
cerr << "Usage: " << argv[0] << " <input file> <output file> "
<< "<threshold> <output_format> [options]" << endl;
cerr << "Usage: " << argv[0] << " <input file> <output file> <cope-image>"
<< " <z-threshold> <p-threshold> <output_format> <DetLambdaHalf>"
<< " <No. Voxels in Mask> [options]" << endl;
cerr << "[-f] : threshold is fractional (of robust range) instead of "
<< "absolute"<<endl;
cerr << "[-6,18,26] : number of connections at each point (default is 26)"
......@@ -27,7 +29,8 @@ void print_usage(char *argv[])
cerr << "1 - size order"<<endl;
cerr << "2 - highest / lowest original value"<<endl;
cerr << "3 - mean original value"<<endl;
// cerr << "4 - unprocessed unique integer label"<<endl;
cerr << "4 - save thresholded map and print out stats"<<endl;
// cerr << "5 - unprocessed unique integer label"<<endl;
}
......@@ -35,52 +38,90 @@ template <class T>
void get_stats(const volume<int>& labelim, const volume<T>& origim,
std::vector<int>& size,
std::vector<T>& maxvals, std::vector<T>& meanvals,
std::vector<int>& maxx, std::vector<int>& maxy,
std::vector<int>& maxz, std::vector<float>& cogx,
std::vector<float>& cogy, std::vector<float>& cogz,
bool invert)
{
int labelnum = labelim.max();
size.resize(labelnum+1,0);
maxvals.resize(labelnum+1, (T) 0);
meanvals.resize(labelnum+1,(T) 0);
maxx.resize(labelnum+1,0);
maxy.resize(labelnum+1,0);
maxz.resize(labelnum+1,0);
cogx.resize(labelnum+1,0.0f);
cogy.resize(labelnum+1,0.0f);
cogz.resize(labelnum+1,0.0f);
std::vector<float> sum(labelnum+1,0.0);
for (int z=labelim.minz(); z<=labelim.maxz(); z++) {
for (int y=labelim.miny(); y<=labelim.maxy(); y++) {
for (int x=labelim.minx(); x<=labelim.maxx(); x++) {
int idx = labelim(x,y,z);
T oxyz = origim(x,y,z);
size[idx]++;
sum[idx]+=(float) origim(x,y,z);
cogx[idx]+=((float) oxyz)*x;
cogy[idx]+=((float) oxyz)*y;
cogz[idx]+=((float) oxyz)*z;
sum[idx]+=(float) oxyz;
if (invert) {
if ((size[idx]==1) || (maxvals[idx]>origim(x,y,z))) {
maxvals[idx] = origim(x,y,z);
if ((size[idx]==1) || (maxvals[idx]>oxyz)) {
maxvals[idx] = oxyz;
maxx[idx] = x;
maxy[idx] = y;
maxz[idx] = z;
}
} else {
if ((size[idx]==1) || (maxvals[idx]<origim(x,y,z))) {
maxvals[idx] = origim(x,y,z);
if ((size[idx]==1) || (maxvals[idx]<oxyz)) {
maxvals[idx] = oxyz;
maxx[idx] = x;
maxy[idx] = y;
maxz[idx] = z;
}
}
}
}
}
for (int n=0; n<=labelnum; n++) {
if (size[n]>0.0) meanvals[n] = (T) (sum[n]/((float) size[n]));
if (size[n]>0.0) {
meanvals[n] = (T) (sum[n]/((float) size[n]));
cogx[n] /= sum[n];
cogy[n] /= sum[n];
cogz[n] /= sum[n];
}
}
}
void get_sizeorder(const std::vector<int>& size, std::vector<int>& sizeorder)
std::vector<int> get_sortindex(const std::vector<int>& vals)
{
// find order of labels so that sizes are in *ascending* order
int length=size.size()-1;
std::vector<pair<int, int> > sortlist(length+1);
for (int n=0; n<=length; n++) {
sortlist[n] = pair<int, int>(size[n],n);
// return the mapping of old indices to new indices in the
// new *ascending* sort of vals
int length=vals.size();
std::vector<pair<int, int> > sortlist(length);
for (int n=0; n<length; n++) {
sortlist[n] = pair<int, int>(vals[n],n);
}
sortlist[0].first = 0; // force the background component to be in 0th posn
sort(sortlist.begin(),sortlist.end()); // O(N.log(N))
std::vector<int> idx(length);
for (int n=0; n<length; n++) {
idx[n] = sortlist[n].second;
}
return idx;
}
void get_sizeorder(const std::vector<int>& size, std::vector<int>& sizeorder)
{
std::vector<int> sizecopy(size), idx;
sizecopy[0] = 0; // force the background component to be in 0th posn
idx = get_sortindex(sizecopy);
// second part of pair is now the prior-index of the sorted values
sizeorder.resize(length+1,0);
for (int n=0; n<=length; n++) {
sizeorder[sortlist[n].second] = n; // maps old index to new
int length = size.size();
sizeorder.resize(length,0);
for (int n=0; n<length; n++) {
sizeorder[idx[n]] = n; // maps old index to new
}
}
......@@ -104,19 +145,23 @@ template <class T>
int fmrib_main(int argc, char *argv[])
{
volume<int> labelim;
int outform = atoi(argv[4]);
std::vector<int> size;
float p_thresh = atof(argv[5]);
int outform = atoi(argv[6]);
float dLh = atof(argv[7]);
int V = atoi(argv[8]);
std::vector<int> size, maxx, maxy, maxz;
std::vector<float> cogx, cogy, cogz;
std::vector<T> maxvals, meanvals;
float th;
bool invert=false, fractionalthresh=false;
int numconnected=26;
{
// read in the volume
volume<T> vin;
read_volume(vin,argv[1]);
// parse the options
bool invert=false, fractionalthresh=false;
int numconnected=26;
for (int n=5; n<argc; n++) {
for (int n=9; n<argc; n++) {
string opt=argv[n];
if (opt=="-f") { fractionalthresh=true; }
else if (opt=="-6") { numconnected=6; }
......@@ -130,8 +175,7 @@ int fmrib_main(int argc, char *argv[])
}
// threshold the volume
float th;
th = atof(argv[3]);
th = atof(argv[4]);
if (th<0.0) {
invert=true;
th=-th;
......@@ -149,7 +193,8 @@ int fmrib_main(int argc, char *argv[])
// process according to the output format
read_volume(vin,argv[1]); // re-read the original volume
get_stats(labelim,vin,size,maxvals,meanvals,invert);
get_stats(labelim,vin,size,maxvals,meanvals,
maxx,maxy,maxz,cogx,cogy,cogz,invert);
}
cout << "CLUSTERS " << size.size()-1 << endl;
......@@ -175,6 +220,90 @@ int fmrib_main(int argc, char *argv[])
relabel_image(labelim,relabeledim,meanvals);
save_volume(relabeledim,argv[2]);
} else if (outform==4) {
Infer infer(dLh, th, V);
std::vector<int> idx;
// size[0] = 0;
// idx = get_sortindex(size);
int length = size.size();
// // re-threshold for p
// std::vector<int> pthreshlabels = size;
// for (int n=1; n<length; n++) {
// int index=idx[n];
// int k = size[index];
// float p = infer(k);
// if (p>p_thresh) { pthreshlabels[n] = 0; }
// }
// {
// volume<int> oldlabels(labelim);
// relabel_image(oldlabels,labelim,pthreshlabels);
// }
// volume<T> vin;
// copyconvert(labelim,vin);
// vin.threshold(0.5); // redundant?
// labelim = connected_components(vin,numconnected);
// read_volume(vin,argv[1]); // re-read the original volume
// get_stats(labelim,vin,size,maxvals,meanvals,
// maxx,maxy,maxz,cogx,cogy,cogz,invert);
// process the COPE image
volume<T> vin;
read_volume(vin,argv[3]);
std::vector<int> copesize, copemaxx, copemaxy, copemaxz;
std::vector<float> copecogx, copecogy, copecogz;
std::vector<T> copemax, copemean;
get_stats(labelim,vin,copesize,copemax,copemean,
copemaxx,copemaxy,copemaxz,copecogx,copecogy,copecogz,invert);
// re-threshold for p
size[0] = 0;
std::vector<int> pthreshsize;
pthreshsize = size;
int nozeroclust=0;
for (int n=1; n<length; n++) {
int k = size[n];
float p = infer(k);
if (p>p_thresh) {
pthreshsize[n] = 0;
nozeroclust++;
}
}
cerr << "No of sub-threshold p-clusters = " << nozeroclust << endl;
idx = get_sortindex(pthreshsize);
// Print results
cout << "Cluster\tNumber\tMaximum\tp\t-log(p)\tMax Pos\tMax Pos\tMax Pos"
<< "\tCOG\tCOG\tCOG\tMaximum\tMax Pos\tMax Pos\tMax Pos\tMean" << endl;
cout << "Index\tVoxels\tValue\t\t\t X \t Y \t Z \t X \t Y \t Z \tCope"
<< "\t X \t Y \t Z \tCope" << endl;
for (int n=length-1; n>=1; n--) {
int index=idx[n];
int k = size[index];
float p = infer(k);
if (pthreshsize[index]>0) {
cout << n - nozeroclust << "\t" << k << "\t" << maxvals[index] << "\t"
<< p << "\t" << -log10(p) << "\t" << maxx[index] << "\t"
<< maxy[index] << "\t" << maxz[index] << "\t"
<< cogx[index] << "\t" << cogy[index] << "\t"
<< cogz[index] << "\t" << copemax[index] << "\t"
<< copecogx[index] << "\t" << copecogy[index] << "\t"
<< copecogz[index] << "\t" << copemean[index] << endl;
}
}
// Threshold the input volume st it is 0 for all non-clusters
// and maintains the same values otherwise
labelim.threshold(1);
volume<T> lcopy;
read_volume(vin,argv[1]);
copyconvert(labelim,lcopy);
vin *= lcopy;
save_volume(vin,argv[2]);
} else if (outform==5) {
save_volume(labelim,argv[2]);
} else {
cerr << "Unrecognised output format" << endl;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment