← Back to branch summary

~ubuntu-branches/ubuntu/raring/voxbo/raring

« back to all changes in this revision

Viewing changes to stats/vbscoregen.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Michael Hanke
  • Date: 2010-06-06 11:33:11 UTC
  • Revision ID: james.westby@ubuntu.com-20100606113311-v3c13imdkkd5n7ae
Tags: upstream-1.8.5~svn1172
Import upstream version 1.8.5~svn1172

Show diffs side-by-side

added added

removed removed

Lines of Context:
stats/vbscoregen.cpp
 
1
 
 
2
// vbscoregen.cpp
 
3
// simulate scores from a mask
 
4
// Copyright (c) 2010 by The VoxBo Development Team
 
5
 
 
6
// VoxBo is free software: you can redistribute it and/or modify it
 
7
// under the terms of the GNU General Public License as published by
 
8
// the Free Software Foundation, either version 3 of the License, or
 
9
// (at your option) any later version.
 
10
// 
 
11
// VoxBo is distributed in the hope that it will be useful, but
 
12
// WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
14
// General Public License for more details.
 
15
// 
 
16
// You should have received a copy of the GNU General Public License
 
17
// along with VoxBo.  If not, see <http://www.gnu.org/licenses/>.
 
18
// 
 
19
// For general information on VoxBo, including the latest complete
 
20
// source code and binary distributions, manual, and associated files,
 
21
// see the VoxBo home page at: http://www.voxbo.org/
 
22
//
 
23
// original version written by Dan Kimberg
 
24
 
 
25
using namespace std;
 
26
 
 
27
#include <fstream>
 
28
#include "vbutil.h"
 
29
#include "vbio.h"
 
30
#include "vbscoregen.hlp.h"
 
31
 
 
32
void vbscoregen_help();
 
33
void vbscoregen_version();
 
34
 
 
35
 
 
36
gsl_rng *rng;
 
37
double gaussian_random(double sigma);
 
38
 
 
39
class VBscoregen {
 
40
public:
 
41
  VBscoregen(string cfg) {parseconfig(cfg);}
 
42
private:
 
43
  Tes lesions;
 
44
  map<string,VBRegion> maskmap;
 
45
  map<string,VBVoxel> voxelmap;
 
46
  double mean1,sd1,mean2,sd2;
 
47
  VB_Vector scores;
 
48
  void parseconfig(string configfile);
 
49
  double genscore(double pct);
 
50
  void do_ifvoxel(string v1);
 
51
  void do_ifvoxelor(string v1,string v2);
 
52
  void do_ifvoxeland(string v1,string v2);
 
53
  void do_region(string r1);
 
54
  void do_regionor(string r1,string r2);
 
55
  void do_regionand(string r1,string r2);
 
56
  void do_regionpct(string r1);
 
57
  void do_regionpctor(string r1,string r2);
 
58
  void do_regionpctand(string r1,string r2);
 
59
};
 
60
 
 
61
int
 
62
main(int argc,char **argv)
 
63
{
 
64
  if (argc!=2) {
 
65
    vbscoregen_help();
 
66
    exit(0);
 
67
  }
 
68
  string arg=argv[1];
 
69
  if (arg=="-v") {
 
70
    vbscoregen_version();
 
71
    exit(0);
 
72
  }
 
73
  VBscoregen sg(arg);
 
74
}
 
75
 
 
76
 
 
77
 
 
78
double
 
79
VBscoregen::genscore(double pct)
 
80
{
 
81
  double mean=mean1+(pct*(mean2-mean1));
 
82
  double sd=sd1+(pct*(sd2-sd1));
 
83
  return mean+gaussian_random(sd);
 
84
}
 
85
 
 
86
void
 
87
VBscoregen::parseconfig(string configfile)
 
88
{
 
89
  ifstream infile;
 
90
  infile.open(configfile.c_str());
 
91
  if (!infile) {
 
92
    // FIXME
 
93
    exit(111);
 
94
  }
 
95
  chdir(xdirname(configfile).c_str());
 
96
 
 
97
  char buf[STRINGLEN];
 
98
  tokenlist args;
 
99
  while (infile.getline(buf,STRINGLEN,'\n')) {
 
100
    args.ParseLine(buf);
 
101
    if (args.size()==0) continue;
 
102
    if (args[0][0]=='#' || args[0][0] == ';')
 
103
      continue;
 
104
    if (args[0]=="intact" && args.size()==3) {
 
105
      mean1=strtod(args[1]);
 
106
      sd1=strtod(args[2]);
 
107
    }
 
108
    else if (args[0]=="damaged" && args.size()==3) {
 
109
      mean2=strtod(args[1]);
 
110
      sd2=strtod(args[2]);
 
111
    }
 
112
    else if (args[0]=="mask" && args.size()==3) {
 
113
      Cube m;
 
114
      if (m.ReadFile(args[2])) {
 
115
        // 
 
116
        exit (202);
 
117
      }
 
118
      VBRegion rr=findregion_mask(m,vb_agt,0.0);
 
119
      maskmap[args[1]]=rr;
 
120
      // cout loaded mask xxxx stats
 
121
    }
 
122
    else if (args[0]=="voxel" && args.size()==5) {
 
123
      VBVoxel vv(strtol(args[2]),strtol(args[3]),strtol(args[4]));
 
124
      voxelmap[args[1]]=vv;
 
125
      // cout loaded voxel xxxx location
 
126
    }
 
127
    else if (args[0]=="lesions" && args.size()==2) {
 
128
      if (lesions.ReadFile(args[1])) {
 
129
        // FIXME
 
130
        exit(191);
 
131
      }
 
132
      scores.resize(lesions.dimt);
 
133
    }
 
134
    else if (args[0]=="ifvoxel" && args.size()==2) {
 
135
      do_ifvoxel(args[1]);
 
136
    }
 
137
    else if (args[0]=="ifvoxelor" && args.size()==3) {
 
138
      do_ifvoxelor(args[1],args[2]);
 
139
    }
 
140
    else if (args[0]=="ifvoxeland" && args.size()==3) {
 
141
      do_ifvoxeland(args[1],args[2]);
 
142
    }
 
143
    else if (args[0]=="region" && args.size()==2) {
 
144
      do_region(args[1]);
 
145
    }
 
146
    else if (args[0]=="regionor" && args.size()==3) {
 
147
      do_regionor(args[1],args[2]);
 
148
    }
 
149
    else if (args[0]=="regionand" && args.size()==3) {
 
150
      do_regionand(args[1],args[2]);
 
151
    }
 
152
    else if (args[0]=="regionpct" && args.size()==2) {
 
153
      do_regionpct(args[1]);
 
154
    }
 
155
    else if (args[0]=="regionpctor" && args.size()==3) {
 
156
      do_regionpctor(args[1],args[2]);
 
157
    }
 
158
    else if (args[0]=="regionpctand" && args.size()==3) {
 
159
      do_regionpctand(args[1],args[2]);
 
160
    }
 
161
    else if (args[0]=="generate" && args.size()==2) {
 
162
      scores.WriteFile(args[1]);
 
163
    }
 
164
    else {
 
165
      infile.close();
 
166
      // 
 
167
      exit(128);
 
168
    }
 
169
  }
 
170
  infile.close();
 
171
  return;
 
172
}
 
173
 
 
174
 
 
175
void
 
176
VBscoregen::do_ifvoxel(string v1)
 
177
{
 
178
  if (!voxelmap.count(v1))
 
179
    exit(999);
 
180
  for (int i=0; i<lesions.dimt; i++) {
 
181
    if (abs(lesions.GetValue(voxelmap[v1],i))>DBL_MIN)
 
182
      scores[i]=genscore(1.0);
 
183
    else
 
184
      scores[i]=genscore(0.0);
 
185
  }
 
186
}
 
187
 
 
188
void
 
189
VBscoregen::do_ifvoxelor(string v1,string v2)
 
190
{
 
191
  if (!voxelmap.count(v1) || !voxelmap.count(v2))
 
192
    exit(999);
 
193
  for (int i=0; i<lesions.dimt; i++) {
 
194
    if (abs(lesions.GetValue(voxelmap[v1],i))>DBL_MIN
 
195
        ||abs(lesions.GetValue(voxelmap[v2],i))>DBL_MIN)
 
196
      scores[i]=genscore(1.0);
 
197
    else
 
198
      scores[i]=genscore(0.0);
 
199
  }
 
200
}
 
201
 
 
202
void
 
203
VBscoregen::do_ifvoxeland(string v1,string v2)
 
204
{
 
205
  if (!voxelmap.count(v1) || !voxelmap.count(v2))
 
206
    exit(999);
 
207
  for (int i=0; i<lesions.dimt; i++) {
 
208
    if (abs(lesions.GetValue(voxelmap[v1],i))>DBL_MIN
 
209
        &&abs(lesions.GetValue(voxelmap[v2],i))>DBL_MIN)
 
210
      scores[i]=genscore(1.0);
 
211
    else
 
212
      scores[i]=genscore(0.0);
 
213
  }
 
214
}
 
215
 
 
216
void
 
217
VBscoregen::do_region(string r1)
 
218
{
 
219
  if (!maskmap.count(r1))
 
220
    exit(999);
 
221
  VBRegion reg=maskmap[r1];
 
222
  for (int i=0; i<lesions.dimt; i++) {
 
223
    bool dflg=0;
 
224
    for (VI v=reg.begin(); v!=reg.end(); v++) {
 
225
      if (abs(lesions.GetValue(v->second,i))>DBL_MIN) {
 
226
        dflg=1;
 
227
        break;
 
228
      }
 
229
    }
 
230
    if (dflg)
 
231
      scores[i]=genscore(1.0);
 
232
    else
 
233
      scores[i]=genscore(0.0);
 
234
  }
 
235
}
 
236
 
 
237
void
 
238
VBscoregen::do_regionor(string r1,string r2)
 
239
{
 
240
  if (!maskmap.count(r1) || !maskmap.count(r2))
 
241
    exit(999);
 
242
  VBRegion reg1=maskmap[r1];
 
243
  VBRegion reg2=maskmap[r2];
 
244
  for (int i=0; i<lesions.dimt; i++) {
 
245
    bool dflg1=0,dflg2=0;
 
246
    for (VI v=reg1.begin(); v!=reg1.end(); v++) {
 
247
      if (abs(lesions.GetValue(v->second,i))>DBL_MIN) {
 
248
        dflg1=1;
 
249
        break;
 
250
      }
 
251
    }
 
252
    if (!dflg1)
 
253
      for (VI v=reg2.begin(); v!=reg2.end(); v++) {
 
254
        if (abs(lesions.GetValue(v->second,i))>DBL_MIN) {
 
255
          dflg2=1;
 
256
          break;
 
257
        }
 
258
      }
 
259
    if (dflg1 || dflg2)
 
260
      scores[i]=genscore(1.0);
 
261
    else
 
262
      scores[i]=genscore(0.0);
 
263
  }
 
264
}
 
265
 
 
266
void
 
267
VBscoregen::do_regionand(string r1,string r2)
 
268
{
 
269
  if (!maskmap.count(r1) || !maskmap.count(r2))
 
270
    exit(999);
 
271
  VBRegion reg1=maskmap[r1];
 
272
  VBRegion reg2=maskmap[r2];
 
273
  for (int i=0; i<lesions.dimt; i++) {
 
274
    bool dflg1=0,dflg2=0;
 
275
    for (VI v=reg1.begin(); v!=reg1.end(); v++) {
 
276
      if (abs(lesions.GetValue(v->second,i))>DBL_MIN) {
 
277
        dflg1=1;
 
278
        break;
 
279
      }
 
280
    }
 
281
    if (dflg1)
 
282
      for (VI v=reg2.begin(); v!=reg2.end(); v++) {
 
283
        if (abs(lesions.GetValue(v->second,i))>DBL_MIN) {
 
284
          dflg2=1;
 
285
          break;
 
286
        }
 
287
      }
 
288
    if (dflg1 && dflg2)
 
289
      scores[i]=genscore(1.0);
 
290
    else
 
291
      scores[i]=genscore(0.0);
 
292
  }
 
293
}
 
294
 
 
295
void
 
296
VBscoregen::do_regionpct(string r1)
 
297
{
 
298
  if (!maskmap.count(r1))
 
299
    exit(999);
 
300
  VBRegion reg1=maskmap[r1];
 
301
  for (int i=0; i<lesions.dimt; i++) {
 
302
    uint32 cnt=0;
 
303
    for (VI v=reg1.begin(); v!=reg1.end(); v++) {
 
304
      if (abs(lesions.GetValue(v->second,i))>DBL_MIN)
 
305
        cnt++;
 
306
    }
 
307
    scores[i]=genscore((double)cnt/(double)reg1.size());
 
308
  }
 
309
}
 
310
 
 
311
void
 
312
VBscoregen::do_regionpctor(string r1,string r2)
 
313
{
 
314
  if (!maskmap.count(r1) || !maskmap.count(r2))
 
315
    exit(999);
 
316
  VBRegion reg1=maskmap[r1];
 
317
  VBRegion reg2=maskmap[r2];
 
318
  for (int i=0; i<lesions.dimt; i++) {
 
319
    double pct1,pct2;
 
320
    uint32 cnt=0;
 
321
    for (VI v=reg1.begin(); v!=reg1.end(); v++) {
 
322
      if (abs(lesions.GetValue(v->second,i))>DBL_MIN)
 
323
        cnt++;
 
324
    }
 
325
    pct1=(double)cnt/(double)reg1.size();
 
326
    for (VI v=reg2.begin(); v!=reg2.end(); v++) {
 
327
      if (abs(lesions.GetValue(v->second,i))>DBL_MIN)
 
328
        cnt++;
 
329
    }
 
330
    pct2=(double)cnt/(double)reg2.size();
 
331
    scores[i]=genscore((pct2>pct1 ? pct2 : pct1));
 
332
  }
 
333
}
 
334
 
 
335
void
 
336
VBscoregen::do_regionpctand(string r1,string r2)
 
337
{
 
338
  if (!maskmap.count(r1) || !maskmap.count(r2))
 
339
    exit(999);
 
340
  VBRegion reg1=maskmap[r1];
 
341
  VBRegion reg2=maskmap[r2];
 
342
  for (int i=0; i<lesions.dimt; i++) {
 
343
    double pct1,pct2;
 
344
    uint32 cnt=0;
 
345
    for (VI v=reg1.begin(); v!=reg1.end(); v++) {
 
346
      if (abs(lesions.GetValue(v->second,i))>DBL_MIN)
 
347
        cnt++;
 
348
    }
 
349
    pct1=(double)cnt/(double)r1.size();
 
350
    for (VI v=reg2.begin(); v!=reg2.end(); v++) {
 
351
      if (abs(lesions.GetValue(v->second,i))>DBL_MIN)
 
352
        cnt++;
 
353
    }
 
354
    pct2=(double)cnt/(double)reg2.size();
 
355
    scores[i]=genscore((pct2>pct1 ? pct1 : pct2));
 
356
  }
 
357
}
 
358
 
 
359
 
 
360
 
 
361
double
 
362
gaussian_random(double sigma)
 
363
{
 
364
  if (!rng) {
 
365
    rng=gsl_rng_alloc(gsl_rng_mt19937);
 
366
    assert(rng);
 
367
    gsl_rng_set(rng,VBRandom());
 
368
  }
 
369
  return gsl_ran_gaussian(rng,sigma);
 
370
}
 
371
 
 
372
 
 
373
 
 
374
 
 
375
void
 
376
vbscoregen_help()
 
377
{
 
378
  cout << boost::format(myhelp) % vbversion;
 
379
}
 
380
 
 
381
void
 
382
vbscoregen_version()
 
383
{
 
384
  printf("VoxBo vbscoregen (v%s)\n",vbversion.c_str());
 
385
}