gmca test

Author: jstarck

pycs/sparsity/sparse2d/mr_gmca.py

@@ -19,7 +19,6 @@
 from pycs.misc.cosmostat_init import *
 from pycs.misc.cosmostat_init import writefits
 
-
 ##
 #  Function that calls mr_gmca to perform blind source separation on the
 #  input data.
@@ -75,7 +74,9 @@ def mr_gmca(data, opt=None, path="./", remove_files=True, verbose=False, FileOut
 
     result = readfits(file_out_source)
     est_mixmat = readfits(file_out_mat)
+    est_mixmat = est_mixmat.T
     est_invmixmat = readfits(file_out_invmat)
+    est_invmixmat = est_invmixmat.T
 
     # Return the mr_transform results (and the output file names).
     if remove_files:
@@ -87,3 +88,92 @@ def mr_gmca(data, opt=None, path="./", remove_files=True, verbose=False, FileOut
         return result, est_mixmat, est_invmixmat
     else:
         return result, est_mixmat, est_invmixmat
+
+
+# Main
+if __name__ == "__main__":
+    # to run the test, need to install scikit-image and import the two following pakages
+    from skimage import data, color
+    from skimage.transform import resize
+    import matplotlib.pyplot as plt
+    from pycs.sparsity.sparse2d.bss_eval import *
+
+    sources = load_source_images() 
+    # mixed_images, A = mix_sources_images(sources)
+    #  np.random.seed(0)
+    mixed_images, A = mix_sources_images_noise(sources, noise_level=0.1)
+    info(sources)
+    
+    print("Mixing Matrix (A):\n", A)
+    print("Sources shape:", sources.shape)          # (2, H, W)
+    print("Mixed images shape:", mixed_images.shape)  # (3, H, W)
+
+    optF = '-S2 -K3 -t14 -n5'  # with bi-orthogonal WT abd final denoising at 3sigma
+    optStarlet = '-S2 -K3 -t2 -n5'  # with starlet  and final denoising at 3sigma
+    optCurvelet = '-E1 -S2 -K3  -t28 -n5 '  # with curvelet abd final denoising at 3sigma
+
+    corrPerm=False
+    verbose=False
+    SRec, Emat, Eimat = mr_gmca(mixed_images, opt=optF, remove_files=False, verbose=verbose)
+    SRec =  reorder_and_fix_sign(sources, SRec)
+    # print(" ==> Bi-Orth Wavelet Source err = ", compute_sdr(sources, SRec))
+    # print("A_true shape:", A.shape)
+    # print("A_est  shape:", Emat.shape)
+    # error = amari_error(A, Emat)
+    # print(" ==> Mixing matrix err = ", amari_error(A, Emat))
+    # print("Mixing Matrix (A):\n", Emat)
+    CA, NMSE = evaluate(A, sources, Emat, SRec, corrPerm=corrPerm)
+    print(' ==> Bi-Orth Wavelet Source err: CA = %.4f | NMSE= %.4f' % (CA, NMSE))
+    
+    # Visualization
+    fig, axs = plt.subplots(1, 5, figsize=(15, 5))
+    axs[0].imshow(sources[0], cmap='gray')
+    axs[0].set_title("Source Image 1")
+    axs[1].imshow(sources[1], cmap='gray')
+    axs[1].set_title("Source Image 2")
+    for i in range(3):
+        axs[i+2].imshow(mixed_images[i], cmap='gray')
+        axs[i+2].set_title(f"Mixed Image {i+1}")
+    fig, axs = plt.subplots(1, 2, figsize=(15, 5))
+    axs[0].imshow(SRec[0], cmap='gray')
+    axs[0].set_title("7/9 WT GMCA Image 1")
+    axs[1].imshow(SRec[1], cmap='gray')
+    axs[1].set_title("7/9 WT GMCA Image 2")
+    for ax in axs:
+        ax.axis('off')
+    plt.tight_layout()
+    plt.show()
+    
+    # ----- STARLET ------ 
+    SRec, Emat, Eimat = mr_gmca(mixed_images, opt=optStarlet, remove_files=False, verbose=verbose)
+    SRec =  reorder_and_fix_sign(sources, SRec)
+    CA, NMSE = evaluate(A, sources, Emat, SRec, corrPerm=corrPerm)
+    print('==> Starlet Source err: CA = %.4f | NMSE =  %.4f' % (CA, NMSE))
+    
+    # Visualization
+    fig, axs = plt.subplots(1, 2, figsize=(15, 5))
+    axs[0].imshow(SRec[0], cmap='gray')
+    axs[0].set_title("Starlet GMCA Image 1")
+    axs[1].imshow(SRec[1], cmap='gray')
+    axs[1].set_title("Starlet GMCA Image 2")
+    for ax in axs:
+        ax.axis('off')
+    plt.tight_layout()
+    plt.show()
+
+    # ----- CURVELET ------ 
+    SRec, Emat, Eimat = mr_gmca(mixed_images, opt=optCurvelet, remove_files=False, verbose=verbose)
+    SRec =  reorder_and_fix_sign(sources, SRec)
+    CA, NMSE = evaluate(A, sources, Emat, SRec, corrPerm=corrPerm)
+    print('==> Curvelet Source err: CA = %.4f | NMSE =  %.4f' % (CA, NMSE))
+    
+    # Visualization
+    fig, axs = plt.subplots(1, 2, figsize=(15, 5))
+    axs[0].imshow(SRec[0], cmap='gray')
+    axs[0].set_title("Curvelet GMCA Image 1")
+    axs[1].imshow(SRec[1], cmap='gray')
+    axs[1].set_title("Curvelet GMCA Image 2")
+    for ax in axs:
+        ax.axis('off')
+    plt.tight_layout()
+    plt.show()