A computer-implemented method for identifying clouds in a digital image, comprising pixels, of a scene, for at least one digital image of at least one respective scene S, the method comprising: For each scene S: quantifying at least one heterogeneity/homogeneity characteristic, h, in each of a multiplicity of locations within a monochromatic or panchromatic digital image of the scene, by computing, for at least one individual location which defines a vicinity of locations spatially adjacent thereto within said digital image, at least one function which varies monotonically as a function of the individual location's vicinity's homogeneity thereby to define at least one characteristic, indicative of local heterogeneity/homogeneity within said digital image, for the individual location; using a controller for generating an output identifying clouds in the digital image; characterised in that generating an output identifying clouds in the digital image includes performing feature-space and/or image-space transformations for each of plural subsets into which the digital image is divided, including computing a function of a slope of a linear fitted model of a feature cluster in an [h,LNDCI] feature space, thereby to yield transformation outputs, and applying predetermined thresholds for distinguishing between cloud pixels and background pixels to said transformation outputs, wherein said function comprises one of the following: PCOT = h × sin gain − LNDCI × w × cos gain or: PCOT = h × sin gain + LNDCI × w × cos gain or: PCOT = h × sin gain − LNCI × w × cos gain or: PCOT = h × sin gain + LNCI × w × cos gain wherein gain comprises a slope of a linear fitted model of said feature cluster in said [h,LNDCI] feature space and where, for each pixel in each of the plural subsets, LNDCI = log 10 b − h b + h and LNCI = log 10 b − h /b where b denotes brightness of a pixel, w is an empirical model weight, and wherein the cloud pixels are identified for PCOT values below a predetermined threshold. For each scene S: quantifying at least one heterogeneity/homogeneity characteristic, h, in each of a multiplicity of locations within a monochromatic or panchromatic digital image of the scene, by computing, for at least one individual location which defines a vicinity of locations spatially adjacent thereto within said digital image, at least one function which varies monotonically as a function of the individual location's vicinity's homogeneity thereby to define at least one characteristic, indicative of local heterogeneity/homogeneity within said digital image, for the individual location; using a controller for generating an output identifying clouds in the digital image; characterised in that generating an output identifying clouds in the digital image includes quantifying at least one heterogeneity/homogeneity characteristic, h, in each of a multiplicity of locations within a monochromatic or panchromatic digital image of the scene, by computing, for at least one individual location which defines a vicinity of locations spatially adjacent thereto within said digital image, at least one function which varies monotonically as a function of the individual location's vicinity's homogeneity thereby to define at least one characteristic, indicative of local heterogeneity/homogeneity within said digital image, for the individual location; using a controller for generating an output identifying clouds in the digital image; characterised in that generating an output identifying clouds in the digital image includes performing feature-space and/or image-space transformations for each of plural subsets into which the digital image is divided, including computing a function of a slope of a linear fitted model of a feature cluster in an [h,LNDCI] feature space, thereby to yield transformation outputs, and applying predetermined thresholds for distinguishing between cloud pixels and background pixels to said transformation outputs, wherein said function comprises one of the following: PCOT = h × sin gain − LNDCI × w × cos gain or: PCOT = h × sin gain + LNDCI × w × cos gain or: PCOT = h × sin gain − LNCI × w × cos gain or: PCOT = h × sin gain + LNCI × w × cos gain wherein gain comprises a slope of a linear fitted model of said feature cluster in said [h,LNDCI] feature space