In [86]:
%matplotlib inline
import pandas as pd
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
import mia

Loading and Preprocessing

Loading the hologic and synthetic datasets.


In [42]:
hologic = pd.DataFrame.from_csv("real-lines.csv")
phantom = pd.DataFrame.from_csv("phantom-lines.csv")

Loading the meta data for the real and synthetic datasets.


In [43]:
hologic_meta = mia.analysis.create_hologic_meta_data(hologic, "meta_data/real_meta.csv")
phantom_meta = mia.analysis.create_synthetic_meta_data(phantom, 
                                                       "meta_data/synthetic_meta.csv")
phantom_meta.index.name = 'img_name'

Prepare the BI-RADS/VBD labels for both datasets.


In [44]:
hologic_labels = hologic_meta.drop_duplicates().BIRADS
phantom_labels = phantom_meta['VBD.1']

class_labels = pd.concat([hologic_labels, phantom_labels])
class_labels.index.name = "img_name"
labels = mia.analysis.remove_duplicate_index(class_labels)[0]

Creating Features

Create blob features from distribution of blobs


In [45]:
hologic_line_features = mia.analysis.features_from_lines(hologic)
phantom_line_features = mia.analysis.features_from_lines(phantom)

Take a random subset of the real mammograms. This is important so that each patient is not over represented.


In [46]:
hologic_line_features['patient_id'] = hologic_meta.drop_duplicates()['patient_id']
hologic_line_features_subset = mia.analysis.create_random_subset(hologic_line_features, 
                                                                 'patient_id')

Take a random subset of the phantom mammograms. This is important so that each case is not over represented.


In [47]:
syn_feature_meta = mia.analysis.remove_duplicate_index(phantom_meta)
phantom_line_features['phantom_name'] = syn_feature_meta.phantom_name.tolist()
phantom_line_features_subset = mia.analysis.create_random_subset(phantom_line_features, 
                                                                 'phantom_name')

Combine the features from both datasets.


In [48]:
features = pd.concat([hologic_line_features_subset, phantom_line_features_subset])
assert features.shape[0] == 96
features.head()


Out[48]:
count mean std min 25% 50% 75% max skew kurtosis upper_dist_count
p214-010-60001-cl.png 72 161.791667 245.194659 1 61.5 94.5 137.25 1744 4.786025 26.793023 16
p214-010-60005-ml.png 124 177.153226 252.405502 1 57.0 91.5 174.75 1725 3.323887 13.873128 31
p214-010-60008-cl.png 105 99.695238 77.001889 1 57.0 76.0 110.00 454 2.537117 7.788977 33
p214-010-60012-mr.png 213 163.037559 315.235968 1 55.0 79.0 158.00 3677 7.384801 74.429175 52
p214-010-60013-cr.png 225 155.368889 180.493203 1 68.0 95.0 162.00 1285 3.584626 15.784599 60

Filter some features, such as the min, to remove noise.


In [49]:
selected_features = features.drop(['min'], axis=1)
selected_features.fillna(0, inplace=True)

Compare Real and Synthetic Features

Compare the distributions of features detected from the real mammograms and the phantoms using the Kolmogorov-Smirnov two sample test.


In [50]:
ks_stats = [list(stats.ks_2samp(hologic_line_features[col], 
                                phantom_line_features[col]))
                                for col in selected_features.columns]

ks_test = pd.DataFrame(ks_stats, columns=['KS', 'p-value'], index=selected_features.columns)
ks_test.to_latex("tables/line_features_ks.tex")
ks_test


Out[50]:
KS p-value
count 0.933333 1.338265e-51
mean 0.593056 4.356206e-21
std 0.654167 1.450514e-25
25% 0.143056 1.255126e-01
50% 0.304167 7.344875e-06
75% 0.508333 1.320817e-15
max 0.737500 2.231475e-32
skew 0.480556 5.426886e-14
kurtosis 0.506944 1.598384e-15
upper_dist_count 0.913889 1.724254e-49

Dimensionality Reduction

t-SNE

Running t-SNE to obtain a two dimensional representation.


In [51]:
real_index = hologic_line_features_subset.index
phantom_index = phantom_line_features_subset.index

In [52]:
kwargs = {
    'learning_rate': 200,
    'perplexity': 20,
    'verbose': 1
}

In [53]:
SNE_mapping_2d = mia.analysis.tSNE(selected_features, n_components=2, **kwargs)


[t-SNE] Computing pairwise distances...
[t-SNE] Computed conditional probabilities for sample 96 / 96
[t-SNE] Mean sigma: 1.347012
[t-SNE] Error after 65 iterations with early exaggeration: 12.391945
[t-SNE] Error after 132 iterations: 0.662618

In [54]:
mia.plotting.plot_mapping_2d(SNE_mapping_2d, real_index, phantom_index, labels)
plt.savefig('figures/mappings/line_SNE_mapping_2d.png', dpi=300)


Running t-SNE to obtain a 3 dimensional mapping


In [55]:
SNE_mapping_3d = mia.analysis.tSNE(selected_features, n_components=3, **kwargs)


[t-SNE] Computing pairwise distances...
[t-SNE] Computed conditional probabilities for sample 96 / 96
[t-SNE] Mean sigma: 1.347012
[t-SNE] Error after 100 iterations with early exaggeration: 16.628509
[t-SNE] Error after 302 iterations: 2.665706

In [83]:
mia.plotting.plot_mapping_3d(SNE_mapping_3d, real_index, phantom_index, labels)


Out[83]:
<matplotlib.axes._subplots.Axes3DSubplot at 0x110b11910>

Isomap

Running Isomap to obtain a 2 dimensional mapping


In [57]:
iso_kwargs = {
    'n_neighbors': 4,
}

In [58]:
iso_mapping_2d = mia.analysis.isomap(selected_features, n_components=2, **iso_kwargs)

In [59]:
mia.plotting.plot_mapping_2d(iso_mapping_2d, real_index, phantom_index, labels)
plt.savefig('figures/mappings/line_iso_mapping_2d.png', dpi=300)



In [60]:
iso_mapping_3d = mia.analysis.isomap(selected_features, n_components=3, **iso_kwargs)

In [87]:
mia.plotting.plot_mapping_3d(iso_mapping_3d, real_index, phantom_index, labels)


Out[87]:
<matplotlib.axes._subplots.Axes3DSubplot at 0x1128091d0>
<matplotlib.figure.Figure at 0x111347bd0>

Locally Linear Embedding

Running locally linear embedding to obtain 2d mapping


In [62]:
lle_kwargs = {
    'n_neighbors': 4,
}

In [63]:
lle_mapping_2d = mia.analysis.lle(selected_features, n_components=2, **lle_kwargs)

In [64]:
mia.plotting.plot_mapping_2d(lle_mapping_2d, real_index, phantom_index, labels)
plt.savefig('figures/mappings/line_lle_mapping_2d.png', dpi=300)



In [65]:
lle_mapping_3d = mia.analysis.lle(selected_features, n_components=3, **lle_kwargs)

In [85]:
mia.plotting.plot_mapping_3d(lle_mapping_3d, real_index, phantom_index, labels)


Out[85]:
<matplotlib.axes._subplots.Axes3DSubplot at 0x10902cf90>

Quality Assessment of Dimensionality Reduction

Assess the quality of the DR against measurements from the co-ranking matrices. First create co-ranking matrices for each of the dimensionality reduction mappings


In [67]:
max_k = 50

In [68]:
SNE_mapping_2d_cm = mia.coranking.coranking_matrix(selected_features, 
                                                   SNE_mapping_2d)
iso_mapping_2d_cm = mia.coranking.coranking_matrix(selected_features, 
                                                   iso_mapping_2d)
lle_mapping_2d_cm = mia.coranking.coranking_matrix(selected_features, 
                                                   lle_mapping_2d)

SNE_mapping_3d_cm = mia.coranking.coranking_matrix(selected_features, 
                                                   SNE_mapping_3d)
iso_mapping_3d_cm = mia.coranking.coranking_matrix(selected_features, 
                                                   iso_mapping_3d)
lle_mapping_3d_cm = mia.coranking.coranking_matrix(selected_features, 
                                                   lle_mapping_3d)

2D Mappings


In [69]:
SNE_trustworthiness_2d = [mia.coranking.trustworthiness(SNE_mapping_2d_cm, k) 
                          for k in range(1, max_k)]
iso_trustworthiness_2d = [mia.coranking.trustworthiness(iso_mapping_2d_cm, k) 
                          for k in range(1, max_k)]
lle_trustworthiness_2d = [mia.coranking.trustworthiness(lle_mapping_2d_cm, k) 
                          for k in range(1, max_k)]

In [70]:
trustworthiness_df = pd.DataFrame([SNE_trustworthiness_2d,
                                   iso_trustworthiness_2d,
                                   lle_trustworthiness_2d], 
                                   index=['SNE', 'Isomap', 'LLE']).T
trustworthiness_df.plot()
plt.savefig('figures/quality_measures/line_trustworthiness_2d.png', dpi=300)



In [71]:
SNE_continuity_2d = [mia.coranking.continuity(SNE_mapping_2d_cm, k) 
                     for k in range(1, max_k)]
iso_continuity_2d = [mia.coranking.continuity(iso_mapping_2d_cm, k) 
                     for k in range(1, max_k)]
lle_continuity_2d = [mia.coranking.continuity(lle_mapping_2d_cm, k) 
                     for k in range(1, max_k)]

In [72]:
continuity_df = pd.DataFrame([SNE_continuity_2d,
                              iso_continuity_2d,
                              lle_continuity_2d], 
                              index=['SNE', 'Isomap', 'LLE']).T
continuity_df.plot()
plt.savefig('figures/quality_measures/line_continuity_2d.png', dpi=300)



In [73]:
SNE_lcmc_2d = [mia.coranking.LCMC(SNE_mapping_2d_cm, k) 
               for k in range(2, max_k)]
iso_lcmc_2d = [mia.coranking.LCMC(iso_mapping_2d_cm, k) 
               for k in range(2, max_k)]
lle_lcmc_2d = [mia.coranking.LCMC(lle_mapping_2d_cm, k) 
               for k in range(2, max_k)]

In [74]:
lcmc_df = pd.DataFrame([SNE_lcmc_2d,
                        iso_lcmc_2d,
                        lle_lcmc_2d], 
                        index=['SNE', 'Isomap', 'LLE']).T
lcmc_df.plot()
plt.savefig('figures/quality_measures/line_lcmc_2d.png', dpi=300)


3D Mappings


In [75]:
SNE_trustworthiness_3d = [mia.coranking.trustworthiness(SNE_mapping_3d_cm, k) 
                          for k in range(1, max_k)]
iso_trustworthiness_3d = [mia.coranking.trustworthiness(iso_mapping_3d_cm, k) 
                          for k in range(1, max_k)]
lle_trustworthiness_3d = [mia.coranking.trustworthiness(lle_mapping_3d_cm, k) 
                          for k in range(1, max_k)]

In [76]:
trustworthiness3d_df = pd.DataFrame([SNE_trustworthiness_3d,
                                   iso_trustworthiness_3d,
                                   lle_trustworthiness_3d], 
                                   index=['SNE', 'Isomap', 'LLE']).T
trustworthiness3d_df.plot()
plt.savefig('figures/quality_measures/line_trustworthiness_3d.png', dpi=300)



In [77]:
SNE_continuity_3d = [mia.coranking.continuity(SNE_mapping_3d_cm, k) 
                     for k in range(1, max_k)]
iso_continuity_3d = [mia.coranking.continuity(iso_mapping_3d_cm, k) 
                     for k in range(1, max_k)]
lle_continuity_3d = [mia.coranking.continuity(lle_mapping_3d_cm, k) 
                     for k in range(1, max_k)]

In [78]:
continuity3d_df = pd.DataFrame([SNE_continuity_3d,
                              iso_continuity_3d,
                              lle_continuity_3d], 
                              index=['SNE', 'Isomap', 'LLE']).T
continuity3d_df.plot()
plt.savefig('figures/quality_measures/line_continuity_3d.png', dpi=300)



In [79]:
SNE_lcmc_3d = [mia.coranking.LCMC(SNE_mapping_3d_cm, k) 
               for k in range(2, max_k)]
iso_lcmc_3d = [mia.coranking.LCMC(iso_mapping_3d_cm, k) 
               for k in range(2, max_k)]
lle_lcmc_3d = [mia.coranking.LCMC(lle_mapping_3d_cm, k) 
               for k in range(2, max_k)]

In [80]:
lcmc3d_df = pd.DataFrame([SNE_lcmc_3d,
                        iso_lcmc_3d,
                        lle_lcmc_3d], 
                        index=['SNE', 'Isomap', 'LLE']).T
lcmc3d_df.plot()
plt.savefig('figures/quality_measures/line_lcmc_3d.png', dpi=300)