Skip to content

ipca

IteratedPCA

Bases: MetaAlgo

Iterated PCA Change Classifier.

The number of unchanged pixels usually outnumber changed pixels and since they're highly correlated over time, they should lie along the first principal axis while the changed pixels lie along the second axis. However, since the principal components are calculated from the covariance matrix computed for all pixels, the no-change axis might be poorly defined. Iterated PCA solves this by calculating the principal components iteratively and weighting each pixel by its probability to be no change pixels.

Accepted flags
  • niter = Number of iterations IPCA should be run
References
  • Wiemker, R. (1997). An iterative spectral-spatial Bayesian labeling approach for unsupervised robust change detection on remotely sensed multispectral imagery. In Proceedings of the 7th International Conference on Computer Analysis of Images and Patterns, volume LCNS 1296, pages 263–370.
  • Canty, M.J. (2019). Image Analysis, Classification, and Change Detection in Remote Sensing: With Algorithms for Python (4th ed.). CRC Press. https://doi.org/10.1201/9780429464348
Source code in changedet/algos/ipca.py 
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 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
106
107
108
@AlgoCatalog.register("ipca")
class IteratedPCA(MetaAlgo):
    """Iterated PCA Change Classifier.

    The number of unchanged pixels usually outnumber changed pixels and since they're highly
    correlated over time, they should lie along the first principal axis while the changed pixels
    lie along the second axis. However, since the principal components are calculated from the
    covariance matrix computed for **all** pixels, the no-change axis might be poorly defined.
    Iterated PCA solves this by calculating the principal components iteratively and weighting
    each pixel by its probability to be no change pixels.

    Accepted flags
    --------------
    - niter = Number of iterations IPCA should be run

    References
    ----------
    - Wiemker, R. (1997). An iterative spectral-spatial Bayesian labeling approach for unsupervised
    robust change detection on remotely sensed multispectral imagery. In Proceedings of the 7th
    International Conference on Computer Analysis of Images and Patterns, volume LCNS 1296, pages
    263–370.
    - Canty, M.J. (2019). Image Analysis, Classification, and Change Detection in Remote Sensing:
    With Algorithms for Python (4th ed.). CRC Press. https://doi.org/10.1201/9780429464348
    """

    @classmethod
    def run(cls, im1: np.ndarray, im2: np.ndarray, **flags: Any) -> np.ndarray:
        """Run IPCA algorithm.

        Args:
            im1 (np.ndarray): Image 1 array
            im2 (np.ndarray): Image 2 array
            **flags (dict): Flags for the algorithm

        Run `changedet --algo ipca algo --help` for information on flags
        """
        niter = flags.get("niter", 5)
        band_idx = flags.get("band")
        logger = flags["logger"]

        band_str = "all bands" if band_idx == -1 else "band " + str(band_idx)

        logger.info("Running IPCA algorithm for %d iteration(s) on %s", niter, band_str)

        bands, rows, cols = im1.shape

        cmaps = []
        n_clusters = 2
        for band in range(bands):
            # TODO: Allow band sublist
            if band_idx == -1:
                logger.info(f"Processing band {band + 1}")
            else:
                logger.info(f"Processing band {band_idx}")

            fim1 = im1[band].flatten()
            fim2 = im2[band].flatten()

            X = np.zeros((rows * cols, 2))
            X[:, 0] = fim1 - np.mean(fim1)
            X[:, 1] = fim2 - np.mean(fim2)

            # Initial PCA
            ows = OnlineWeightStats(2)
            ows.update(X)
            eivals, eigvecs = np.linalg.eigh(ows.cov)
            eivals = eivals[::-1]
            eigvecs = eigvecs[:, ::-1]

            pcs = X @ eigvecs

            gmm = GMM(n_clusters, 60)

            for _ in range(niter):
                # Calculate responsibility matrix
                U = np.random.rand(X.shape[0], n_clusters)
                for _ in range(n_clusters):
                    sigma = np.sqrt(eivals[1])
                    unfrozen = np.where(np.abs(pcs[:, 1]) >= sigma)[0]
                    frozen = np.where(np.abs(pcs[:, 1]) < sigma)[0]
                    U[frozen, 0] = 1.0
                    U[frozen, 1] = 0.0
                    for j in range(2):
                        U[:, j] = U[:, j] / np.sum(U, 1)

                r, _, _, _ = gmm.fit(X, U, unfrozen)

                tmp = r[:, 0] + 1 - 1
                ows.update(X, tmp)
                eivals, eigvecs = np.linalg.eigh(ows.cov)
                eivals = eivals[::-1]
                eigvecs = eigvecs[:, ::-1]

                pcs = X @ eigvecs

            cmap = np.reshape(r[:, 1], (rows, cols))
            cmaps.append(cmap)

        return np.array(cmaps)

run(im1, im2, **flags) classmethod

Run IPCA algorithm.

Parameters:

Name Type Description Default
im1 ndarray

Image 1 array

 required
im2 ndarray

Image 2 array

 required
**flags dict

Flags for the algorithm

 {}

Run changedet --algo ipca algo --help for information on flags

Source code in changedet/algos/ipca.py 
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 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
106
107
108
@classmethod
def run(cls, im1: np.ndarray, im2: np.ndarray, **flags: Any) -> np.ndarray:
    """Run IPCA algorithm.

    Args:
        im1 (np.ndarray): Image 1 array
        im2 (np.ndarray): Image 2 array
        **flags (dict): Flags for the algorithm

    Run `changedet --algo ipca algo --help` for information on flags
    """
    niter = flags.get("niter", 5)
    band_idx = flags.get("band")
    logger = flags["logger"]

    band_str = "all bands" if band_idx == -1 else "band " + str(band_idx)

    logger.info("Running IPCA algorithm for %d iteration(s) on %s", niter, band_str)

    bands, rows, cols = im1.shape

    cmaps = []
    n_clusters = 2
    for band in range(bands):
        # TODO: Allow band sublist
        if band_idx == -1:
            logger.info(f"Processing band {band + 1}")
        else:
            logger.info(f"Processing band {band_idx}")

        fim1 = im1[band].flatten()
        fim2 = im2[band].flatten()

        X = np.zeros((rows * cols, 2))
        X[:, 0] = fim1 - np.mean(fim1)
        X[:, 1] = fim2 - np.mean(fim2)

        # Initial PCA
        ows = OnlineWeightStats(2)
        ows.update(X)
        eivals, eigvecs = np.linalg.eigh(ows.cov)
        eivals = eivals[::-1]
        eigvecs = eigvecs[:, ::-1]

        pcs = X @ eigvecs

        gmm = GMM(n_clusters, 60)

        for _ in range(niter):
            # Calculate responsibility matrix
            U = np.random.rand(X.shape[0], n_clusters)
            for _ in range(n_clusters):
                sigma = np.sqrt(eivals[1])
                unfrozen = np.where(np.abs(pcs[:, 1]) >= sigma)[0]
                frozen = np.where(np.abs(pcs[:, 1]) < sigma)[0]
                U[frozen, 0] = 1.0
                U[frozen, 1] = 0.0
                for j in range(2):
                    U[:, j] = U[:, j] / np.sum(U, 1)

            r, _, _, _ = gmm.fit(X, U, unfrozen)

            tmp = r[:, 0] + 1 - 1
            ows.update(X, tmp)
            eivals, eigvecs = np.linalg.eigh(ows.cov)
            eivals = eivals[::-1]
            eigvecs = eigvecs[:, ::-1]

            pcs = X @ eigvecs

        cmap = np.reshape(r[:, 1], (rows, cols))
        cmaps.append(cmap)

    return np.array(cmaps)