

In [95]:

    
%matplotlib inline
import pandas as pd
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
import mia









    



Warning: Cannot change to a different GUI toolkit: qt. Using osx instead.

Loading and Preprocessing

Loading the hologic and synthetic datasets.



In [56]:

    
hologic = pd.DataFrame.from_csv("hologic_blobs.csv")
phantom = pd.DataFrame.from_csv("synthetic_blobs.csv")

Loading the meta data for the real and synthetic datasets.



In [57]:

    
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 [58]:

    
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 [59]:

    
hologic_blob_features = mia.analysis.features_from_blobs(hologic)
phantom_blob_features = mia.analysis.features_from_blobs(phantom)

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



In [60]:

    
hologic_blob_features['patient_id'] = hologic_meta.drop_duplicates()['patient_id']
hologic_blob_features_subset = mia.analysis.create_random_subset(hologic_blob_features, 
                                                                 'patient_id')

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



In [61]:

    
syn_feature_meta = mia.analysis.remove_duplicate_index(phantom_meta)
phantom_blob_features['phantom_name'] = syn_feature_meta.phantom_name.tolist()
phantom_blob_features_subset = mia.analysis.create_random_subset(phantom_blob_features, 'phantom_name')

Combine the features from both datasets.



In [62]:

    
features = pd.concat([hologic_blob_features_subset, phantom_blob_features_subset])
assert features.shape[0] == 96
features.head()









    Out[62]:






  
    
      
      blob_count
      avg_radius
      std_radius
      min_radius
      max_radius
      small_radius_count
      med_radius_count
      large_radius_count
      density
      upper_dist_count
      25%
      50%
      75%
    
  
  
    
      p214-010-60001-cr.png
        78
       19.054538
       17.506086
       8
        90.509668
        68
        4
       6
       40.749811
       22
       8
       11.313708
       22.627417
    
    
      p214-010-60005-cl.png
        97
       20.132590
       23.255605
       8
       181.019336
        94
        2
       1
       41.456308
       27
       8
       11.313708
       22.627417
    
    
      p214-010-60008-ml.png
       310
       18.777405
       24.388840
       8
       181.019336
       299
        5
       6
       37.667995
       68
       8
       11.313708
       16.000000
    
    
      p214-010-60012-cr.png
       185
       13.947419
        7.884251
       8
        45.254834
       152
       28
       5
       47.292289
       63
       8
       11.313708
       16.000000
    
    
      p214-010-60013-cr.png
       141
       21.821567
       30.666648
       8
       181.019336
       134
        1
       6
       44.586708
       38
       8
       11.313708
       22.627417

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



In [63]:

    
selected_features = features.drop(['min_radius'], axis=1)

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 [64]:

    
ks_stats = [list(stats.ks_2samp(hologic_blob_features[col], 
                                phantom_blob_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/blob_features_ks.tex")
ks_test









    Out[64]:






  
    
      
      KS
      p-value
    
  
  
    
      blob_count
       0.341667
       2.774753e-07
    
    
      avg_radius
       0.847222
       1.360953e-42
    
    
      std_radius
       0.711111
       3.929143e-30
    
    
      max_radius
       0.363889
       3.327680e-08
    
    
      small_radius_count
       0.319444
       2.024353e-06
    
    
      med_radius_count
       0.338889
       3.583275e-07
    
    
      large_radius_count
       0.733333
       5.112303e-32
    
    
      density
       0.169444
       4.114705e-02
    
    
      upper_dist_count
       0.345833
       1.883393e-07
    
    
      25%
       0.358333
       5.726005e-08
    
    
      50%
       0.743056
       7.334764e-33
    
    
      75%
       0.777778
       5.796838e-36

Dimensionality Reduction

t-SNE

Running t-SNE to obtain a two dimensional representation.



In [65]:

    
real_index = hologic_blob_features_subset.index
phantom_index = phantom_blob_features_subset.index



In [66]:

    
kwargs = {
    'learning_rate': 200,
    'perplexity': 30,
    'verbose': 1
}



In [67]:

    
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.641093
[t-SNE] Error after 72 iterations with early exaggeration: 12.697723
[t-SNE] Error after 144 iterations: 0.495631



In [68]:

    
mia.plotting.plot_mapping_2d(SNE_mapping_2d, real_index, phantom_index, labels)
plt.savefig('figures/mappings/blob_SNE_mapping_2d.png', dpi=300)

Running t-SNE to obtain a 3 dimensional mapping



In [98]:

    
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.641093
[t-SNE] Error after 100 iterations with early exaggeration: 14.272421
[t-SNE] Error after 300 iterations: 2.104165



In [99]:

    
mia.plotting.plot_mapping_3d(SNE_mapping_3d, real_index, phantom_index, labels)









    Out[99]:





<matplotlib.axes._subplots.Axes3DSubplot at 0x10cca73d0>

Isomap

Running Isomap to obtain a 2 dimensional mapping



In [71]:

    
iso_kwargs = {
    'n_neighbors': 4,
}



In [72]:

    
iso_mapping_2d = mia.analysis.isomap(selected_features, n_components=2, **iso_kwargs)



In [73]:

    
mia.plotting.plot_mapping_2d(iso_mapping_2d, real_index, phantom_index, labels)
plt.savefig('figures/mappings/blob_iso_mapping_2d.png', dpi=300)



In [74]:

    
iso_mapping_3d = mia.analysis.isomap(selected_features, n_components=3, **iso_kwargs)



In [100]:

    
mia.plotting.plot_mapping_3d(iso_mapping_3d, real_index, phantom_index, labels)









    Out[100]:





<matplotlib.axes._subplots.Axes3DSubplot at 0x10c306810>

Locally Linear Embedding

Running locally linear embedding to obtain 2d mapping



In [76]:

    
lle_kwargs = {
    'n_neighbors': 4,
}



In [77]:

    
lle_mapping_2d = mia.analysis.lle(selected_features, n_components=2, **lle_kwargs)



In [78]:

    
mia.plotting.plot_mapping_2d(lle_mapping_2d, real_index, phantom_index, labels)
plt.savefig('figures/mappings/blob_lle_mapping_2d.png', dpi=300)



In [79]:

    
lle_mapping_3d = mia.analysis.lle(selected_features, n_components=3, **lle_kwargs)



In [101]:

    
mia.plotting.plot_mapping_3d(lle_mapping_3d,real_index, phantom_index, labels)









    Out[101]:





<matplotlib.axes._subplots.Axes3DSubplot at 0x10d9d7810>

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 [81]:

    
max_k = 50



In [82]:

    
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 [83]:

    
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 [84]:

    
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/blob_trustworthiness_2d.png', dpi=300)



In [85]:

    
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 [86]:

    
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/blob_continuity_2d.png', dpi=300)



In [87]:

    
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 [88]:

    
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/blob_lcmc_2d.png', dpi=300)

3D Mappings



In [89]:

    
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 [90]:

    
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/blob_trustworthiness_3d.png', dpi=300)



In [91]:

    
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 [92]:

    
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/blob_continuity_3d.png', dpi=300)



In [93]:

    
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 [94]:

    
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/blob_lcmc_3d.png', dpi=300)

	blob_count	avg_radius	std_radius	min_radius	max_radius	small_radius_count	med_radius_count	large_radius_count	density	upper_dist_count	25%	50%	75%
p214-010-60001-cr.png	78	19.054538	17.506086	8	90.509668	68	4	6	40.749811	22	8	11.313708	22.627417
p214-010-60005-cl.png	97	20.132590	23.255605	8	181.019336	94	2	1	41.456308	27	8	11.313708	22.627417
p214-010-60008-ml.png	310	18.777405	24.388840	8	181.019336	299	5	6	37.667995	68	8	11.313708	16.000000
p214-010-60012-cr.png	185	13.947419	7.884251	8	45.254834	152	28	5	47.292289	63	8	11.313708	16.000000
p214-010-60013-cr.png	141	21.821567	30.666648	8	181.019336	134	1	6	44.586708	38	8	11.313708	22.627417

	KS	p-value
blob_count	0.341667	2.774753e-07
avg_radius	0.847222	1.360953e-42
std_radius	0.711111	3.929143e-30
max_radius	0.363889	3.327680e-08
small_radius_count	0.319444	2.024353e-06
med_radius_count	0.338889	3.583275e-07
large_radius_count	0.733333	5.112303e-32
density	0.169444	4.114705e-02
upper_dist_count	0.345833	1.883393e-07
25%	0.358333	5.726005e-08
50%	0.743056	7.334764e-33
75%	0.777778	5.796838e-36