hausdorff_distance.cxx 11.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
/* -----------------------------------------------------------------------
   See COPYRIGHT.TXT and LICENSE.TXT for copyright and license information
   ----------------------------------------------------------------------- */
#include "plm_config.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "itkImage.h"
10
#include "wirth.h"
11

Gregory C. Sharp's avatar
Gregory C. Sharp committed
12
#include "distance_map.h"
13
#include "hausdorff_distance.h"
14
#include "image_boundary.h"
15
16
#include "itk_bbox.h"
#include "itk_crop.h"
17
#include "itk_image_header_compare.h"
18
#include "itk_image_load.h"
19
#include "itk_image_save.h"
Gregory C. Sharp's avatar
Gregory C. Sharp committed
20
#include "itk_resample.h"
21
#include "itk_union.h"
22
#include "logfile.h"
Gregory C. Sharp's avatar
Gregory C. Sharp committed
23
24
25
#include "plm_image.h"
#include "plm_image_header.h"
#include "volume.h"
26
27
28
29

class Hausdorff_distance_private {
public:
    Hausdorff_distance_private () {
30
        pct_hausdorff_distance_fraction = 0.95;
31
        dmap_alg = "";
32
        maximum_distance = FLT_MAX;
33
        vbb = ADAPTIVE_PADDING;
34
35
36
37
        this->clear_statistics ();
    }
public:
    void clear_statistics () {
38
        hausdorff_distance = 0.f;
Gregory C. Sharp's avatar
Gregory C. Sharp committed
39
        min_min_hausdorff_distance = FLT_MAX;
40
41
        avg_avg_hausdorff_distance = 0.f;
        max_avg_hausdorff_distance = 0.f;
42
        pct_hausdorff_distance = 0.f;
43
        boundary_hausdorff_distance = 0.f;
Gregory C. Sharp's avatar
Gregory C. Sharp committed
44
        min_min_boundary_hausdorff_distance = FLT_MAX;
45
46
        avg_avg_boundary_hausdorff_distance = 0.f;
        max_avg_boundary_hausdorff_distance = 0.f;
47
        pct_boundary_hausdorff_distance = 0.f;
48
49
50
    }
public:
    float hausdorff_distance;
Gregory C. Sharp's avatar
Gregory C. Sharp committed
51
    float min_min_hausdorff_distance;
52
53
    float avg_avg_hausdorff_distance;
    float max_avg_hausdorff_distance;
54
    float pct_hausdorff_distance;
55
    float boundary_hausdorff_distance;
Gregory C. Sharp's avatar
Gregory C. Sharp committed
56
    float min_min_boundary_hausdorff_distance;
57
58
    float avg_avg_boundary_hausdorff_distance;
    float max_avg_boundary_hausdorff_distance;
59
    float pct_boundary_hausdorff_distance;
60
61
    float pct_hausdorff_distance_fraction;

62
    std::string dmap_alg;
63
    float maximum_distance;
64
    Volume_boundary_behavior vbb;
65

66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
    UCharImageType::Pointer ref_image;
    UCharImageType::Pointer cmp_image;
};

Hausdorff_distance::Hausdorff_distance ()
{
    d_ptr = new Hausdorff_distance_private;
}

Hausdorff_distance::~Hausdorff_distance ()
{
    delete d_ptr;
}

void 
Hausdorff_distance::set_reference_image (const char* image_fn)
{
    d_ptr->ref_image = itk_image_load_uchar (image_fn, 0);
}

void 
Hausdorff_distance::set_reference_image (
    const UCharImageType::Pointer image)
{
    d_ptr->ref_image = image;
}

void 
Hausdorff_distance::set_compare_image (const char* image_fn)
{
    d_ptr->cmp_image = itk_image_load_uchar (image_fn, 0);
}

void 
Hausdorff_distance::set_compare_image (
    const UCharImageType::Pointer image)
{
    d_ptr->cmp_image = image;
}

106
107
108
109
110
111
112
void 
Hausdorff_distance::set_hausdorff_distance_fraction (
    float hausdorff_distance_fraction)
{
    d_ptr->pct_hausdorff_distance_fraction = hausdorff_distance_fraction;
}

113
114
115
116
117
118
void 
Hausdorff_distance::set_distance_map_algorithm (const std::string& dmap_alg)
{
    d_ptr->dmap_alg = dmap_alg;
}

119
120
121
122
123
124
void 
Hausdorff_distance::set_maximum_distance (float maximum_distance)
{
    d_ptr->maximum_distance = maximum_distance;
}

125
126
127
128
129
130
void
Hausdorff_distance::set_volume_boundary_behavior (Volume_boundary_behavior vbb)
{
    d_ptr->vbb = vbb;
}

Gregory C. Sharp's avatar
Gregory C. Sharp committed
131
132
void 
Hausdorff_distance::run_internal (
133
134
    UCharImageType::Pointer image1,
    UCharImageType::Pointer image2
Gregory C. Sharp's avatar
Gregory C. Sharp committed
135
136
)
{
137
138
    /* Compute distance map */
    Distance_map dmap_filter;
139
    dmap_filter.set_algorithm (d_ptr->dmap_alg);
140
141
142
    dmap_filter.set_input_image (image2);
    dmap_filter.set_inside_is_positive (false);
    dmap_filter.set_use_squared_distance (false);
143
    dmap_filter.set_maximum_distance (d_ptr->maximum_distance);
144
    dmap_filter.set_volume_boundary_behavior (d_ptr->vbb);
145
146
147
    dmap_filter.run ();
    FloatImageType::Pointer dmap = dmap_filter.get_output_image ();

Gregory C. Sharp's avatar
Gregory C. Sharp committed
148
    /* Convert to Plm_image type */
149
    Plm_image pli_uchar (image1);
150
    Volume::Pointer vol_uchar = pli_uchar.get_volume_uchar ();
Gregory C. Sharp's avatar
Gregory C. Sharp committed
151
152
    unsigned char *img_uchar = (unsigned char*) vol_uchar->img;
    Plm_image pli_dmap (dmap);
153
    Volume::Pointer vol_dmap = pli_dmap.get_volume_float ();
Gregory C. Sharp's avatar
Gregory C. Sharp committed
154
155
    float *img_dmap = (float*) vol_dmap->img;

156
157
    /* Find boundary pixels */
    Image_boundary ib;
158
    ib.set_volume_boundary_behavior (d_ptr->vbb);
159
160
161
    ib.set_input_image (image1);
    ib.run ();
    UCharImageType::Pointer itk_ib = ib.get_output_image ();
162
    
163
164
    /* Convert to plm_image */
    Plm_image pli_ib (itk_ib);
165
    Volume::Pointer vol_ib = pli_ib.get_volume_uchar ();
166
167
    unsigned char *img_ib = (unsigned char*) vol_ib->img;

168
    /* Make an array to store the distances */
169
170
    float *h_distance_array = new float[vol_uchar->npix];
    float *bh_distance_array = new float[vol_uchar->npix];
171

Gregory C. Sharp's avatar
Gregory C. Sharp committed
172
    /* Loop through voxels, find distances */
Gregory C. Sharp's avatar
Gregory C. Sharp committed
173
174
    float min_h_distance = FLT_MAX;
    float min_bh_distance = FLT_MAX;
175
176
    float max_h_distance = 0;
    float max_bh_distance = 0;
177
178
179
180
    double sum_h_distance = 0;
    double sum_bh_distance = 0;
    plm_long num_h_vox = 0;
    plm_long num_bh_vox = 0;
Gregory C. Sharp's avatar
Gregory C. Sharp committed
181
    for (plm_long i = 0; i < vol_uchar->npix; i++) {
182
183
184
185
186
187
188
189
190
191
192
193
194
        if (!img_uchar[i]) {
            continue;
        }

        /* Get distance map value for this voxel */
        float h_dist = img_dmap[i];   /* dist for set hausdorff */
        float bh_dist = img_dmap[i];  /* dist for boundary hausdorff */
        if (img_dmap[i] < 0) {
            h_dist = 0;
            bh_dist = - bh_dist;
        }

        /* Update statistics for hausdorff */
195
196
        if (h_dist > max_h_distance) {
            max_h_distance = h_dist;
197
        }
Gregory C. Sharp's avatar
Gregory C. Sharp committed
198
199
200
        if (h_dist < min_h_distance) {
            min_h_distance = h_dist;
        }
201
202
203
204
205
        sum_h_distance += h_dist;
        h_distance_array[num_h_vox] = h_dist;
        num_h_vox ++;
        
        /* Update statistics for boundary hausdorff */
206
        if (img_ib[i]) {
207
208
209
            if (bh_dist > max_bh_distance) {
                max_bh_distance = bh_dist;
            }
Gregory C. Sharp's avatar
Gregory C. Sharp committed
210
211
212
            if (bh_dist < min_bh_distance) {
                min_bh_distance = bh_dist;
            }
213
214
215
            sum_bh_distance += bh_dist;
            bh_distance_array[num_bh_vox] = bh_dist;
            num_bh_vox ++;
Gregory C. Sharp's avatar
Gregory C. Sharp committed
216
217
218
        }
    }

219
    /* Figure out HD95 stuff */
220
221
222
223
224
225
226
227
228
229
    float h_pct = 0, bh_pct = 0;
    if (num_h_vox > 0) {
        int ordinal = (int) floor (
            d_ptr->pct_hausdorff_distance_fraction * num_h_vox-1);
        if (ordinal > num_h_vox - 1) {
            ordinal = num_h_vox - 1;
        }
        h_pct = kth_smallest (h_distance_array, num_h_vox, ordinal);
    }
    if (num_bh_vox > 0) {
230
        int ordinal = (int) floor (
231
232
233
            d_ptr->pct_hausdorff_distance_fraction * num_bh_vox-1);
        if (ordinal > num_bh_vox - 1) {
            ordinal = num_bh_vox - 1;
234
        }
235
        bh_pct = kth_smallest (bh_distance_array, num_bh_vox, ordinal);
236
237
238
    }

    /* Record results */
239
240
241
242
243
    if (max_h_distance > d_ptr->hausdorff_distance) {
        d_ptr->hausdorff_distance = max_h_distance;
    }
    if (max_bh_distance > d_ptr->boundary_hausdorff_distance) {
        d_ptr->boundary_hausdorff_distance = max_bh_distance;
Gregory C. Sharp's avatar
Gregory C. Sharp committed
244
    }
Gregory C. Sharp's avatar
Gregory C. Sharp committed
245
246
247
248
249
250
    if (min_h_distance < d_ptr->min_min_hausdorff_distance) {
        d_ptr->min_min_hausdorff_distance = min_h_distance;
    }
    if (min_bh_distance < d_ptr->min_min_boundary_hausdorff_distance) {
        d_ptr->min_min_boundary_hausdorff_distance = min_bh_distance;
    }
251
    if (num_h_vox > 0) {
252
253
254
255
        float ahd = sum_h_distance / num_h_vox;
        d_ptr->avg_avg_hausdorff_distance += 0.5 * ahd;
        d_ptr->max_avg_hausdorff_distance = std::max (
            d_ptr->max_avg_hausdorff_distance, ahd);
256
257
258
        d_ptr->pct_hausdorff_distance += 0.5 * h_pct;
    }
    if (num_bh_vox > 0) {
259
260
261
262
        float abhd = sum_bh_distance / num_bh_vox;
        d_ptr->avg_avg_boundary_hausdorff_distance += 0.5 * abhd;
        d_ptr->max_avg_boundary_hausdorff_distance = std::max (
            d_ptr->max_avg_boundary_hausdorff_distance, abhd);
263
        d_ptr->pct_boundary_hausdorff_distance += 0.5 * bh_pct;
Gregory C. Sharp's avatar
Gregory C. Sharp committed
264
    }
Gregory C. Sharp's avatar
Gregory C. Sharp committed
265

266
267
    delete[] h_distance_array;
    delete[] bh_distance_array;
Gregory C. Sharp's avatar
Gregory C. Sharp committed
268
269
}

270
271
void 
Hausdorff_distance::run ()
Gregory C. Sharp's avatar
Gregory C. Sharp committed
272
{
273
    /* Resample and/or expand images based on geometry of reference */
Gregory C. Sharp's avatar
Gregory C. Sharp committed
274
    if (!itk_image_header_compare (d_ptr->ref_image, d_ptr->cmp_image)) {
275
        Plm_image_header pih;
276
        
277
278
279
280
281
        pih.set_geometry_to_contain (
            Plm_image_header (d_ptr->cmp_image),
            Plm_image_header (d_ptr->ref_image));
        d_ptr->cmp_image = resample_image (d_ptr->cmp_image, pih, 0, 0);
        d_ptr->ref_image = resample_image (d_ptr->ref_image, pih, 0, 0);
Gregory C. Sharp's avatar
Gregory C. Sharp committed
282
283
    }

284
285
286
287
288
289
290
291
292
293
    /* Crop images to union bounding box containing both structures */
    UCharImageType::Pointer union_image
        = itk_union (d_ptr->cmp_image, d_ptr->ref_image);
    float bbox_coordinates[6];
    int bbox_indices[6];
    itk_bbox (union_image, bbox_coordinates, bbox_indices);
    d_ptr->ref_image = itk_crop_by_index (d_ptr->ref_image, bbox_indices);
    d_ptr->cmp_image = itk_crop_by_index (d_ptr->cmp_image, bbox_indices);
    
    /* Compute distance maps and score the results */
294
295
296
    d_ptr->clear_statistics ();
    this->run_internal (d_ptr->ref_image, d_ptr->cmp_image);
    this->run_internal (d_ptr->cmp_image, d_ptr->ref_image);
Gregory C. Sharp's avatar
Gregory C. Sharp committed
297
298
}

299
300
float 
Hausdorff_distance::get_hausdorff ()
301
302
303
304
{
    return d_ptr->hausdorff_distance;
}

Gregory C. Sharp's avatar
Gregory C. Sharp committed
305
306
307
308
309
310
float 
Hausdorff_distance::get_min_min_hausdorff ()
{
    return d_ptr->min_min_hausdorff_distance;
}

311
float 
312
Hausdorff_distance::get_avg_average_hausdorff ()
313
{
314
315
316
317
318
319
320
    return d_ptr->avg_avg_hausdorff_distance;
}

float 
Hausdorff_distance::get_max_average_hausdorff ()
{
    return d_ptr->max_avg_hausdorff_distance;
321
322
}

323
324
325
326
327
328
float 
Hausdorff_distance::get_percent_hausdorff ()
{
    return d_ptr->pct_hausdorff_distance;
}

329
330
331
332
333
334
float 
Hausdorff_distance::get_boundary_hausdorff ()
{
    return d_ptr->boundary_hausdorff_distance;
}

Gregory C. Sharp's avatar
Gregory C. Sharp committed
335
336
337
338
339
340
float 
Hausdorff_distance::get_min_min_boundary_hausdorff ()
{
    return d_ptr->min_min_boundary_hausdorff_distance;
}

341
float 
342
Hausdorff_distance::get_avg_average_boundary_hausdorff ()
343
{
344
345
346
347
348
349
350
    return d_ptr->avg_avg_boundary_hausdorff_distance;
}

float 
Hausdorff_distance::get_max_average_boundary_hausdorff ()
{
    return d_ptr->max_avg_boundary_hausdorff_distance;
351
352
353
354
355
356
357
358
}

float 
Hausdorff_distance::get_percent_boundary_hausdorff ()
{
    return d_ptr->pct_boundary_hausdorff_distance;
}

359
360
361
362
363
void 
Hausdorff_distance::debug ()
{
    lprintf (
	"Hausdorff distance = %f\n"
364
365
	"Avg average Hausdorff distance = %f\n"
	"Max average Hausdorff distance = %f\n"
366
	"Percent (%.2f) Hausdorff distance = %f\n"
367
	"Hausdorff distance (boundary) = %f\n"
368
369
	"Avg average Hausdorff distance (boundary) = %f\n"
	"Max average Hausdorff distance (boundary) = %f\n"
370
	"Percent (%.2f) Hausdorff distance (boundary) = %f\n",
371
	this->get_hausdorff (),
372
373
	this->get_avg_average_hausdorff (),
	this->get_max_average_hausdorff (),
374
375
        d_ptr->pct_hausdorff_distance_fraction,
	this->get_percent_hausdorff (),
376
	this->get_boundary_hausdorff (),
377
378
	this->get_avg_average_boundary_hausdorff (),
	this->get_max_average_boundary_hausdorff (),
379
380
381
        d_ptr->pct_hausdorff_distance_fraction,
	this->get_percent_boundary_hausdorff ()
    );
382
383
}

384
385
386
387
388
389
390
391
392
void 
do_hausdorff (
    UCharImageType::Pointer image_1, 
    UCharImageType::Pointer image_2
)
{
    Hausdorff_distance hd;
    hd.set_reference_image (image_1);
    hd.set_compare_image (image_2);
393
394
    hd.run ();
    hd.debug ();
395
}