← Back to branch summary

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

« back to all changes in this revision

Viewing changes to munge/vbmaskcompare.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:
munge/vbmaskcompare.cpp
 
1
 
 
2
// questions re: fiez et al.
 
3
 
 
4
// fig 3a shows sizes of 16 and 19 and calculates the percent
 
5
// difference as 16%.  the calculation given in the article (diff
 
6
// divided by mean size) gives 17%.  which is right?
 
7
 
 
8
// fig 3b shows surface distances. the algorithm isn't clear, and not
 
9
// spelled out in the refs, so here's my guess.  accumulate all the
 
10
// surface voxels for each volume.  surface-surface matches are 0's.
 
11
// other surface voxels that don't overlap with non-surface voxels are
 
12
// given the distance to the nearest surface voxel (positive for A-B,
 
13
// negative for B-A).  calculate euclidean distances, but round.
 
14
 
 
15
 
 
16
 
 
17
 
 
18
// vbmaskcompare.cpp
 
19
// mask comparison utility for doing inter-rater reliability
 
20
// Copyright (c) 2005-2009 by The VoxBo Development Team
 
21
 
 
22
// VoxBo is free software: you can redistribute it and/or modify it
 
23
// under the terms of the GNU General Public License as published by
 
24
// the Free Software Foundation, either version 3 of the License, or
 
25
// (at your option) any later version.
 
26
// 
 
27
// VoxBo is distributed in the hope that it will be useful, but
 
28
// WITHOUT ANY WARRANTY; without even the implied warranty of
 
29
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
30
// General Public License for more details.
 
31
// 
 
32
// You should have received a copy of the GNU General Public License
 
33
// along with VoxBo.  If not, see <http://www.gnu.org/licenses/>.
 
34
// 
 
35
// For general information on VoxBo, including the latest complete
 
36
// source code and binary distributions, manual, and associated files,
 
37
// see the VoxBo home page at: http://www.voxbo.org/
 
38
//
 
39
// original version written by Dan Kimberg
 
40
// based on Fiez et al.
 
41
 
 
42
#include <math.h>
 
43
#include "vbutil.h"
 
44
#include "vbio.h"
 
45
#include "vbversion.h"
 
46
#include "vbmaskcompare.hlp.h"
 
47
 
 
48
class subrect {
 
49
public:
 
50
  subrect() {xmin=-1;}
 
51
  int32 xmin,xmax,ymin,ymax,zmin,zmax;
 
52
  void include(int x,int y,int z);
 
53
};
 
54
 
 
55
void
 
56
subrect::include(int x,int y,int z)
 
57
{
 
58
  if (xmin==-1) {
 
59
    xmin=xmax=x;
 
60
    ymin=ymax=y;
 
61
    zmin=zmax=z;
 
62
    return;
 
63
  }
 
64
  if (x<xmin) xmin=x; else if (x>xmax) xmax=x;
 
65
  if (y<ymin) ymin=y; else if (y>ymax) ymax=y;
 
66
  if (z<zmin) zmin=z; else if (z>zmax) zmax=z;
 
67
}
 
68
 
 
69
class VBComp {
 
70
private:
 
71
  Cube mask1,mask2;
 
72
  int m_union,m_intersection,m_just1,m_just2,m_neither;
 
73
  int m_discrepant1,m_discrepant2;
 
74
  VBRegion surface1,surface2;   // surfaces for the two volumes
 
75
  VBRegion unique1,unique2;     // unique voxels for the two volumes
 
76
  VBRegion discrepants1,discrepants2,nondiscrepants;
 
77
  VBRegion allvoxels1,allvoxels2;
 
78
  vector<int> distances;
 
79
public:
 
80
  VBComp();
 
81
  int surfaceflag;
 
82
  // double and integer versions of the maximum allowable non-discrepant distance
 
83
  double f_ddist;
 
84
  int32 f_idist;
 
85
  string stem;
 
86
  int load_masks(string m1,string m2);
 
87
  void check_volumes();
 
88
  void check_volumes2();
 
89
  void check_surfaces();
 
90
  void check_discrepants();
 
91
  void write_discrepants();
 
92
  void summarize();
 
93
};
 
94
 
 
95
void vbmaskcompare_help();
 
96
void vbmaskcompare_version();
 
97
 
 
98
int
 
99
main(int argc,char *argv[])
 
100
{
 
101
  if (argc==1) {
 
102
    vbmaskcompare_help();
 
103
    exit(0);
 
104
  }
 
105
  VBComp cmp;
 
106
  cmp.surfaceflag=0;
 
107
  tokenlist args;
 
108
  vector<string> filelist;
 
109
  args.Transfer(argc-1,argv+1);
 
110
  for (int i=0; i<args.size(); i++) {
 
111
    if (args[i]=="-s")
 
112
      cmp.surfaceflag=1;
 
113
    else if (args[i]=="-d")    // OBSOLETED
 
114
      continue;
 
115
    else if (args[i]=="-dd" && i<args.size()-1) {
 
116
      cmp.f_ddist=strtod(args[i+1]);
 
117
      cmp.f_idist=strtol(args[i+1]);
 
118
      i++;
 
119
    }
 
120
    else if (args[i]=="-o" && i<args.size()-1)
 
121
      cmp.stem=args[++i];
 
122
    else if (args[i]=="-h") {
 
123
      vbmaskcompare_help();
 
124
      exit(0);
 
125
    }
 
126
    else if (args[i]=="-v") {
 
127
      vbmaskcompare_version();
 
128
      exit(0);
 
129
    }
 
130
    else
 
131
      filelist.push_back(args[i]);
 
132
  }
 
133
 
 
134
  if (filelist.size()!=2) {
 
135
    vbmaskcompare_help();
 
136
    exit(5);
 
137
  }
 
138
 
 
139
  if (cmp.load_masks(filelist[0],filelist[1])) {
 
140
    printf("[E] vbmaskcompare: couldn't load masks (or discrepant sizes)\n");
 
141
    exit(10);
 
142
  }
 
143
  // printf("[I] vbmaskcompare: checking volumes\n");
 
144
  cmp.check_volumes();
 
145
  // cmp.check_volumes2();  // little extra work to establish bounding box
 
146
  if (cmp.surfaceflag) {
 
147
    printf("[I] vbmaskcompare: checking surfaces\n");
 
148
    cmp.check_surfaces();
 
149
  }
 
150
  cmp.check_discrepants();
 
151
  cmp.summarize();
 
152
  if (cmp.stem.size()) {
 
153
    cmp.write_discrepants();
 
154
  }
 
155
  exit(0);
 
156
}
 
157
 
 
158
VBComp::VBComp()
 
159
{
 
160
  m_intersection=0;
 
161
  m_union=0;
 
162
  m_just1=0;
 
163
  m_just2=0;
 
164
  m_neither=0;
 
165
  f_ddist=2.0;
 
166
  f_idist=2;
 
167
}
 
168
 
 
169
int
 
170
VBComp::load_masks(string m1,string m2)
 
171
{
 
172
  mask1.ReadFile(m1);
 
173
  mask2.ReadFile(m2);
 
174
  if (!(mask1.data))
 
175
    return 101;
 
176
  if (!(mask2.data))
 
177
    return 102;
 
178
  if (mask1.dimx!=mask2.dimx ||
 
179
      mask1.dimy!=mask2.dimy ||
 
180
      mask1.dimz!=mask2.dimz)
 
181
    return 103;
 
182
  return 0;
 
183
}
 
184
 
 
185
void
 
186
VBComp::check_volumes2()
 
187
{
 
188
  // bool found1=0,found2=0;
 
189
  bool v1,v2;
 
190
  subrect r1,r2;
 
191
  int ind=0;
 
192
  for (int k=0; k<mask1.dimz; k++) {
 
193
    for (int j=0; j<mask1.dimy; j++) {
 
194
      for (int i=0; i<mask1.dimx;  i++) {
 
195
        v1=mask1.testValue(i,j,k);
 
196
        v2=mask2.testValue(i,j,k);
 
197
        if (v1) r1.include(i,j,k);
 
198
        if (v2) r2.include(i,j,k);
 
199
        ind++;
 
200
      }
 
201
    }
 
202
  }
 
203
}
 
204
 
 
205
void
 
206
VBComp::check_volumes()
 
207
{
 
208
  bool v1,v2;
 
209
  for (int k=0; k<mask1.dimz; k++) {
 
210
    for (int j=0; j<mask1.dimy; j++) {
 
211
      for (int i=0; i<mask1.dimx;  i++) {
 
212
        v1=(int)mask1.testValue(i,j,k);
 
213
        v2=(int)mask2.testValue(i,j,k);
 
214
        if (v1 && v2) {
 
215
          m_intersection++;
 
216
          nondiscrepants.add(i,j,k,1);
 
217
        }
 
218
        else if (v1) {
 
219
          unique1.add(i,j,k,0.0);
 
220
          m_just1++;
 
221
        }
 
222
        else if (v2) {
 
223
          unique2.add(i,j,k,0.0);
 
224
          m_just2++;
 
225
        }
 
226
        else
 
227
          m_neither++;
 
228
        if (v1) {
 
229
          allvoxels1.add(i,j,k,0.0);
 
230
          if (mask1.is_surface(i,j,k))
 
231
            surface1.add(i,j,k,0.0);
 
232
        }
 
233
        if (v2) {
 
234
          allvoxels2.add(i,j,k,0.0);
 
235
          if (mask2.is_surface(i,j,k))
 
236
            surface2.add(i,j,k,0.0);
 
237
        }
 
238
      }
 
239
    }
 
240
  }
 
241
  m_union=m_intersection+m_just1+m_just2;
 
242
}
 
243
 
 
244
// here's my understanding of the algorithm.  matching surface voxels
 
245
// count as zeros.  for unmatched surface voxels, that don't overlap
 
246
// with the other mask, find the nearest surface voxel on the other
 
247
// surface.
 
248
 
 
249
void
 
250
VBComp::check_surfaces()
 
251
{
 
252
  int mindist;
 
253
  // distances for first surface (including matching surface voxels)
 
254
  // for (int i=0; i<(int)surface1.size(); i++) {
 
255
  for (VI vox1=surface1.begin(); vox1!=surface1.end(); vox1++) {
 
256
    mindist=mask1.dimx+mask1.dimy+mask1.dimz;  // suitably large value
 
257
    for (VI vox2=surface2.begin();vox2!=surface2.end(); vox2++) {
 
258
      int dist=(int)round(voxeldistance(vox1->second,vox2->second));
 
259
      if (dist<mindist)
 
260
        mindist=dist;
 
261
    }
 
262
    if (mindist==0)
 
263
      distances.push_back(0);
 
264
    if (mindist>0 && !(allvoxels2.contains(vox1->second.x,vox1->second.y,vox1->second.z)))
 
265
      distances.push_back(mindist);
 
266
  }
 
267
  // distances for second surface (excluding matching surface voxels)
 
268
  // for (int i=0; i<(int)surface2.size(); i++) {
 
269
  for (VI vox2=surface2.begin(); vox2!=surface2.end(); vox2++) {
 
270
    mindist=mask1.dimx+mask1.dimy+mask1.dimz;  // suitably large value
 
271
    for (VI vox1=surface1.begin(); vox1!=surface1.end(); vox1++) {
 
272
      int dist=(int)round(voxeldistance(vox1->second,vox2->second));
 
273
      if (dist < mindist)
 
274
        mindist=dist;
 
275
    }
 
276
    if (mindist!=0 && !(allvoxels1.contains(vox2->second.x,vox2->second.y,vox2->second.z)))
 
277
      distances.push_back(0-mindist);
 
278
  }
 
279
  double meandist=0;
 
280
  for (size_t i=0; i<distances.size(); i++)
 
281
    meandist+=distances[i];
 
282
  meandist/=distances.size();
 
283
  printf("[I] vbmaskcompare: mean inter-surface distance: %g\n",meandist);
 
284
  
 
285
}
 
286
 
 
287
void
 
288
VBComp::check_discrepants()
 
289
{
 
290
  m_discrepant1=m_discrepant2=0;
 
291
  discrepants1.clear();
 
292
  discrepants2.clear();
 
293
  double dist,mindist;
 
294
  for (VI vox1=unique1.begin(); vox1!=unique1.end(); vox1++) {
 
295
    mindist=999999999;
 
296
    int x1=vox1->second.x-f_idist;
 
297
    int x2=vox1->second.x+f_idist;
 
298
    int y1=vox1->second.y-f_idist;
 
299
    int y2=vox1->second.y+f_idist;
 
300
    int z1=vox1->second.z-f_idist;
 
301
    int z2=vox1->second.z+f_idist;
 
302
    for (int i=x1; i<=x2; i++) {
 
303
      for (int j=y1; j<=y2; j++) {
 
304
        for (int k=z1; k<=z2; k++) {
 
305
          if (!(mask2.testValue(i,j,k))) continue;
 
306
          dist=voxeldistance(i,j,k,vox1->second.x,vox1->second.y,vox1->second.z);
 
307
          if (dist<mindist)
 
308
            mindist=dist;
 
309
          if (dist<=f_ddist)
 
310
            break;
 
311
        }
 
312
      }
 
313
    }
 
314
    if (mindist-f_ddist>DBL_MIN) {
 
315
      m_discrepant1++;
 
316
      discrepants1.add(vox1->second);
 
317
    }
 
318
    else
 
319
      nondiscrepants.add(vox1->second);
 
320
  }
 
321
  for (VI vox2=unique2.begin(); vox2!=unique2.end(); vox2++) {
 
322
    mindist=999999999;
 
323
    int x1=vox2->second.x-f_idist;
 
324
    int x2=vox2->second.x+f_idist;
 
325
    int y1=vox2->second.y-f_idist;
 
326
    int y2=vox2->second.y+f_idist;
 
327
    int z1=vox2->second.z-f_idist;
 
328
    int z2=vox2->second.z+f_idist;
 
329
    for (int i=x1; i<=x2; i++) {
 
330
      for (int j=y1; j<=y2; j++) {
 
331
        for (int k=z1; k<=z2; k++) {
 
332
          if (!(mask1.testValue(i,j,k))) continue;
 
333
          dist=voxeldistance(i,j,k,vox2->second.x,vox2->second.y,vox2->second.z);
 
334
          if (dist<mindist)
 
335
            mindist=dist;
 
336
          if (dist<=f_ddist)
 
337
            break;
 
338
        }
 
339
      }
 
340
    }
 
341
    if (mindist-f_ddist>DBL_MIN) {
 
342
      m_discrepant2++;
 
343
      discrepants2.add(vox2->second);
 
344
    }
 
345
    else
 
346
      nondiscrepants.add(vox2->second);
 
347
  }
 
348
}
 
349
 
 
350
void
 
351
VBComp::write_discrepants()
 
352
{
 
353
  Cube cb;
 
354
  string fname;
 
355
  cb.SetVolume(mask1.dimx,mask1.dimy,mask1.dimz,vb_byte);
 
356
  cb.zero();
 
357
  for (VI vox=discrepants1.begin(); vox!=discrepants1.end(); vox++)
 
358
    cb.setValue<char>(vox->second.x,vox->second.y,vox->second.z,1);
 
359
  if (discrepants1.size()==0)
 
360
    cout << format("[W] vbmaskcompare: no discrepant voxels in first mask\n");
 
361
  fname=stem;
 
362
  replace_string(fname,"XXX","d1");
 
363
  if (cb.WriteFile(fname))
 
364
    cout << format("[E] vbmaskcompare: couldn't write output file %s\n") % fname;
 
365
  else
 
366
    cout << format("[I] vbmaskcompare: wrote output file %s\n") % fname;
 
367
 
 
368
  cb.zero();
 
369
  for (VI vox=discrepants2.begin(); vox!=discrepants2.end(); vox++)
 
370
    cb.setValue<char>(vox->second.x,vox->second.y,vox->second.z,1);
 
371
  if (discrepants2.size()==0)
 
372
    cout << format("[W] vbmaskcompare: no discrepant voxels in second mask\n");
 
373
  fname=stem;
 
374
  replace_string(fname,"XXX","d2");
 
375
  if (cb.WriteFile(fname))
 
376
    cout << format("[E] vbmaskcompare: couldn't write output file %s\n") % fname;
 
377
  else
 
378
    cout << format("[I] vbmaskcompare: wrote output file %s\n") % fname;
 
379
 
 
380
  cb.zero();
 
381
  for (VI vox=nondiscrepants.begin(); vox!=nondiscrepants.end(); vox++)
 
382
    cb.setValue<char>(vox->second.x,vox->second.y,vox->second.z,1);
 
383
  if (nondiscrepants.size()==0)
 
384
    cout << format("[W] vbmaskcompare: no nondiscrepant voxels anywhere\n");
 
385
  fname=stem;
 
386
  replace_string(fname,"XXX","nond");
 
387
  if (cb.WriteFile(fname))
 
388
    cout << format("[E] vbmaskcompare: couldn't write output file %s\n") % fname;
 
389
  else
 
390
    cout << format("[I] vbmaskcompare: wrote output file %s\n") % fname;
 
391
}
 
392
 
 
393
void
 
394
VBComp::summarize()
 
395
{
 
396
  if (m_union==0) {
 
397
    printf("[E] vbmaskcompare: no voxels found in either mask\n");
 
398
    return;
 
399
  }
 
400
  printf("[I] vbmaskcompare: total voxels: %d\n",m_union+m_neither);
 
401
  printf("[I] vbmaskcompare: %s: %d total, %d unique\n",
 
402
         mask1.GetFileName().c_str(),m_intersection+m_just1,m_just1);
 
403
  printf("[I] vbmaskcompare: %s: %d total, %d unique\n",
 
404
         mask2.GetFileName().c_str(),m_intersection+m_just2,m_just2);
 
405
  printf("[I] vbmaskcompare: total mask voxels (mask union): %d\n",m_union);
 
406
  printf("[I] vbmaskcompare: common (mask intersection): %d\n",m_intersection);
 
407
  float dice=(2.0*m_intersection)/(m_union+m_intersection);
 
408
  printf("[I] vbmaskcompare: dice overlap: %.4f\n",dice);
 
409
  float disc=(float)abs(m_union-m_intersection)/(float)m_union;
 
410
  printf("[I] vbmaskcompare: percent non-overlapping: %.2f%%\n",100.0*disc);
 
411
  printf("[I] vbmaskcompare: discrepant voxels in %s: %d\n",
 
412
         mask1.GetFileName().c_str(),m_discrepant1);
 
413
  printf("[I] vbmaskcompare: discrepant voxels in %s: %d\n",
 
414
         mask2.GetFileName().c_str(),m_discrepant2);
 
415
  printf("[I] vbmaskcompare: total discrepant voxels: %d (%.2f%%)\n",m_discrepant1+m_discrepant2,
 
416
         (float)(m_discrepant1+m_discrepant2)/(float)m_union);
 
417
}
 
418
 
 
419
void
 
420
vbmaskcompare_help()
 
421
{
 
422
  cout << boost::format(myhelp) % vbversion;
 
423
}
 
424
 
 
425
void
 
426
vbmaskcompare_version()
 
427
{
 
428
  printf("VoxBo vbmaskcompare (v%s)\n",vbversion.c_str());
 
429
}