diff --git a/README.md b/README.md index f2ba929..daed196 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@

SpectralUnmixing

-A general, fast, flexible, and including spectral unmixing package. Oriented towards VSWIR imaging spectroscopy data but applicable for different sensor types. Includes options for different treatments of endmember library assemblages, including MESMA and bootstrapping (aka monte carlo) strategies. +A general, fast, flexible, and including spectral unmixing package. Oriented towards VSWIR imaging spectroscopy data but applicable for different sensor types. Includes options for different treatments of endmember library assemblages, including MESMA and bootstrapping (aka Monte Carlo) strategies. ## Installation This package has not been registered yet, but will be soon. In the interim, after cloning and navigating into the repository, it can be installed from the Julia REPL, by running @@ -12,17 +12,19 @@ export JULIA_PROJECT=${PWD} ``` ## Using the script +Currently the package supports reading and writing ENVI raster data. + Basic: ``` -julia unmix.jl REFLECTANCE_IMAGE ENDMEMBER_LIBRARY ENDMEMBER_COLUMN OUTPUT_BASE --mode sma +julia unmix.jl REFLECTANCE_IMAGE ENDMEMBER_LIBRARY ENDMEMBER_COLUMN OUTPUT_BASE --mode sma ``` Parallel implementation (with 10 cores): ``` -julia -p 10 unmix.jl REFLECTANCE_IMAGE ENDMEMBER_LIBRARY ENDMEMBER_COLUMN OUTPUT_BASE --mode sma +julia -p 10 unmix.jl REFLECTANCE_IMAGE ENDMEMBER_LIBRARY ENDMEMBER_COLUMN OUTPUT_BASE --mode sma ``` Bootstrapping uncertainty: @@ -42,4 +44,3 @@ Preset maximum number of endmembers used for unmixing: ``` julia -p 10 unmix.jl REFLECTANCE_IMAGE ENDMEMBER_LIBRARY ENDMEMBER_COLUMN OUTPUT_BASE --mode sma --n_mc 50 --normalization brightness --num_endmembers 10 ``` - diff --git a/examples/simulation_demo.ipynb b/examples/simulation_demo.ipynb index e3f865c..3400180 100644 --- a/examples/simulation_demo.ipynb +++ b/examples/simulation_demo.ipynb @@ -28,7 +28,7 @@ "source": [ "# Spectral Library Preparation\n", "\n", - "Start by loading up a spectral library, filtering it down to the desired classes, interpolating the spectra to the relevant wavelengths, and then removing undesired wavelength regions. By default the ignored regions correspond to atmospheric absorption regions. Here we split the library into the highest level classes: soil, photosynthetic vegetation (pv) and non-photosynthetic vegetation (npv). As the plots below show, all the endmember spectra within each class share general attributes while having small variations between them." + "Start by loading up a spectral library (Ochoa et al. [10.21232/QijXnpFt](https://doi.org/10.21232/QijXnpFt)), filtering it down to the desired classes, interpolating the spectra to the relevant wavelengths, and then removing undesired wavelength regions. By default the ignored regions correspond to atmospheric absorption regions. Here we split the library into the highest level classes: soil, photosynthetic vegetation (pv) and non-photosynthetic vegetation (npv). As the plots below show, all the endmember spectra within each class share general attributes while having small variations between them." ] }, { diff --git a/src/Datasets.jl b/src/Datasets.jl index ff6f2a2..f0aaff1 100644 --- a/src/Datasets.jl +++ b/src/Datasets.jl @@ -17,7 +17,17 @@ using ArchGDAL using GDAL +""" + set_band_names(filename::String, band_names::Vector) +Set the names of the bands in a raster dataset specified by the given filename. + +# Arguments +- `filename::String`: Path to the raster file. The file must be in a format supported by +GDAL. The data file is updated in place. +- `band_names::Vector`: A vector of strings containing the new names for each band in the +raster dataset. The length of this vector must match the number of bands in the raster. +""" function set_band_names(filename::String, band_names::Vector) GC.gc() local_ods = GDAL.gdalopen(filename, GDAL.GA_Update) @@ -27,45 +37,104 @@ function set_band_names(filename::String, band_names::Vector) local_ods = nothing end -function initiate_output_datasets(output_files::Vector{String}, x_len::Int64, y_len::Int64, output_bands::Vector, reference_dataset::ArchGDAL.IDataset) +""" + initiate_output_datasets(output_files::Vector{String}, x_len::Int64, y_len::Int64, + output_bands::Vector, reference_dataset::ArchGDAL.IDataset) + +Initialize multiple output raster datasets in ENVI format with specified dimensions and +band information using a reference dataset for projection and geotransformation. + +# Arguments +- `output_files::Vector{String}`: File paths of the output datasets. Each entry corresponds +to a separate output raster. +- `x_len::Int64`: Width of the output datasets. +- `y_len::Int64`: Height of the output datasets. +- `output_bands::Vector`: Number of bands for each output dataset. The length of this +vector must match the length of `output_files`. +- `reference_dataset::ArchGDAL.IDataset`: Reference dataset object from which projection +and geotransformation information will be copied. + +# Notes +- All bands in the output datasets are initialized with the no-data value of `-9999`. +""" +function initiate_output_datasets(output_files::Vector{String}, x_len::Int64, y_len::Int64, + output_bands::Vector, reference_dataset::ArchGDAL.IDataset) GC.gc() for _o in 1:length(output_files) - outDataset = ArchGDAL.create(output_files[_o], driver=ArchGDAL.getdriver("ENVI"), width=x_len, - height=y_len, nbands=output_bands[_o], dtype=Float64) + outDataset = ArchGDAL.create( + output_files[_o], driver=ArchGDAL.getdriver("ENVI"), width=x_len, + height=y_len, nbands=output_bands[_o], dtype=Float64 + ) ArchGDAL.setproj!(outDataset, ArchGDAL.getproj(reference_dataset)) - @info "Output Image Size (x,y,b): $x_len, $y_len, $output_bands. Creating output fractional cover dataset." + @info "Output Image Size (x,y,b): $x_len, $y_len, $output_bands. + Creating output fractional cover dataset." try - ArchGDAL.setgeotransform!(outDataset, ArchGDAL.getgeotransform(reference_dataset)) + ArchGDAL.setgeotransform!( + outDataset, ArchGDAL.getgeotransform(reference_dataset) + ) catch e - println("No geotransorm avaialble, proceeding without") + println("No geotransorm available, proceeding without") end for _b in 1:length(output_bands) - ArchGDAL.setnodatavalue!(ArchGDAL.getband(outDataset,_b), -9999) + ArchGDAL.setnodatavalue!(ArchGDAL.getband(outDataset, _b), -9999) end output = zeros(x_len, y_len, output_bands[_o]) .- 9999 - ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2]) + ArchGDAL.write!( + outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], + size(output)[2] + ) outDataset = nothing - end end +""" + write_results(output_files::Vector{String}, output_bands::Vector, x_len::Int64, + y_len::Int64, results, args) + +Write processed results to specified output ENVI raster files, including primary outputs, +uncertainty estimates, and complete fractions, based on the provided parameters. + +# Arguments +- `output_files::Vector{String}`: File paths of the output datasets to be created or +updated. Each entry corresponds to a separate output raster. +- `output_bands::Vector`: Number of bands for each output dataset. The length of this +vector must match the number of output files. +- `x_len::Int64`: Width of the output datasets. +- `y_len::Int64`: Height of the output datasets. +- `results`: A collection of results, where each result is expected to be a tuple or +similar structure. Specific elements contain the data to be written to the output datasets. +- `args`: An object or structure containing additional parameters that dictate the +writing behavior. Available options are: + - `args.n_mc`: If greater than 1, an uncertainty output is also written to the second + output file. + - `args.write_complete_fractions`: If set to 1, write complete fractions to the + subsequent output file. + +# Notes +- The primary output is initialized with a no-data value of `-9999` and written to the + first output file specified in `output_files`. +- Each output array is permuted to match the desired dimensions before writing. +""" +function write_results(output_files::Vector{String}, output_bands::Vector, x_len::Int64, + y_len::Int64, results, args) -function write_results(output_files::Vector{String}, output_bands::Vector, x_len::Int64, y_len::Int64, results, args) GC.gc() @info "Write primary output" output = zeros(y_len, x_len, output_bands[1]) .- 9999 for res in results if isnothing(res[2]) == false - output[res[1],res[3], :] = res[2] + output[res[1], res[3], :] = res[2] end end - output = permutedims( output, (2,1,3)) - outDataset = ArchGDAL.read(output_files[1], flags=1, alloweddrivers =["ENVI"]) - ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2]) + output = permutedims(output, (2, 1, 3)) + outDataset = ArchGDAL.read(output_files[1], flags=1, alloweddrivers=["ENVI"]) + ArchGDAL.write!( + outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2] + ) outDataset = nothing ods_idx = 2 @@ -74,13 +143,16 @@ function write_results(output_files::Vector{String}, output_bands::Vector, x_len output = zeros(y_len, x_len, output_bands[ods_idx]) .- 9999 for res in results if isnothing(res[4]) == false - output[res[1],res[3], :] = res[4] + output[res[1], res[3], :] = res[4] end end - output = permutedims( output, (2,1,3)) - outDataset = ArchGDAL.read(output_files[ods_idx], flags=1, alloweddrivers =["ENVI"]) - ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2]) + output = permutedims(output, (2, 1, 3)) + outDataset = ArchGDAL.read(output_files[ods_idx], flags=1, alloweddrivers=["ENVI"]) + ArchGDAL.write!( + outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], + size(output)[2] + ) outDataset = nothing ods_idx += 1 end @@ -90,29 +162,53 @@ function write_results(output_files::Vector{String}, output_bands::Vector, x_len output = zeros(y_len, x_len, output_bands[ods_idx]) .- 9999 for res in results if isnothing(res[5]) == false - output[res[1],res[3], :] = res[5] + output[res[1], res[3], :] = res[5] end end - output = permutedims( output, (2,1,3)) - outDataset = ArchGDAL.read(output_files[ods_idx], flags=1, alloweddrivers =["ENVI"]) - ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2]) + output = permutedims(output, (2, 1, 3)) + outDataset = ArchGDAL.read(output_files[ods_idx], flags=1, alloweddrivers=["ENVI"]) + ArchGDAL.write!( + outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], + size(output)[2] + ) outDataset = nothing ods_idx += 1 end - end -function write_line_results(output_files::Vector{String}, results, n_mc::Int64, write_complete_fractions::Bool) +""" + write_line_results(output_files::Vector{String}, results, n_mc::Int64, + write_complete_fractions::Bool) + +Write line-based results to specified output ENVI raster files, including primary results, +uncertainty estimates, and complete fractions based on provided parameters. + +# Arguments +- `output_files::Vector{String}`: File paths of the output datasets to be created or +updated. Each entry corresponds to a separate output raster. +- `results`: A collection of results, where each entry is expected to be a tuple or +similar structure. Specific elements contain the data to be written to the output datasets. +- `n_mc::Int64`: If greater than one, uncertainty outputs will be written to the +subsequent output file. +- `write_complete_fractions::Bool`: Indicates whether to write complete fractions +to the subsequent output file. +""" +function write_line_results(output_files::Vector{String}, results, n_mc::Int64, + write_complete_fractions::Bool) + GC.gc() - + if isnothing(results[2]) == false - + output = zeros(size(results[3])[1], 1, size(results[2])[2]) .- 9999 output[results[3], 1, :] = results[2] - outDataset = ArchGDAL.read(output_files[1], flags=1, alloweddrivers =["ENVI"]) - ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, results[1]-1, size(output)[1], 1) + outDataset = ArchGDAL.read(output_files[1], flags=1, alloweddrivers=["ENVI"]) + ArchGDAL.write!( + outDataset, output, [1:size(output)[end];], 0, results[1] - 1, + size(output)[1], 1 + ) outDataset = nothing end ods_idx = 2 @@ -121,9 +217,14 @@ function write_line_results(output_files::Vector{String}, results, n_mc::Int64, if isnothing(results[4]) == false output = zeros(size(results[3])[1], 1, size(results[4])[2]) .- 9999 output[results[3], 1, :] = results[4] - - outDataset = ArchGDAL.read(output_files[ods_idx], flags=1, alloweddrivers =["ENVI"]) - ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, results[1]-1, size(output)[1], 1) + + outDataset = ArchGDAL.read( + output_files[ods_idx], flags=1, alloweddrivers=["ENVI"] + ) + ArchGDAL.write!( + outDataset, output, [1:size(output)[end];], 0, results[1] - 1, + size(output)[1], 1 + ) outDataset = nothing end ods_idx += 1 @@ -133,36 +234,72 @@ function write_line_results(output_files::Vector{String}, results, n_mc::Int64, if isnothing(results[5]) == false output = zeros(size(results[3])[1], 1, size(results[5])[2]) .- 9999 output[results[3], 1, :] = results[5] - - outDataset = ArchGDAL.read(output_files[ods_idx], flags=1, alloweddrivers =["ENVI"]) - ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, results[1]-1, size(output)[1], 1) + + outDataset = ArchGDAL.read( + output_files[ods_idx], flags=1, alloweddrivers=["ENVI"] + ) + ArchGDAL.write!( + outDataset, output, [1:size(output)[end];], 0, results[1] - 1, + size(output)[1], 1 + ) outDataset = nothing end ods_idx += 1 end - end -function load_line(reflectance_file::String, reflectance_uncertainty_file::String, line::Int64, - good_bands::Array{Bool}, refl_nodata::Float64) - - reflectance_dataset = ArchGDAL.read(reflectance_file, alloweddrivers =["ENVI"]) - img_dat = convert(Array{Float64},ArchGDAL.readraster(reflectance_file, alloweddrivers =["ENVI"])[:,line,:]) - img_dat = img_dat[:, good_bands] - good_data = .!all(img_dat .== refl_nodata, dims=2)[:,1] - img_dat = img_dat[good_data,:] - - if sum(good_data) >= 1 - if reflectance_uncertainty_file != "" - unc_dat = convert(Array{Float64},ArchGDAL.readraster(reflectance_uncertainty_file, alloweddrivers =["ENVI"])[:,line,:]) - unc_dat = unc_dat[:, good_bands] - unc_dat = unc_dat[good_data,:] - else - unc_dat = nothing - end +""" + load_line(reflectance_file::String, reflectance_uncertainty_file::String, + line::Int64, good_bands::Array{Bool}, refl_nodata::Float64) + +Load a specific line (row) of reflectance data and its associated uncertainty from raster +files, filtering based on specified good bands and no-data values. + +# Arguments +- `reflectance_file::String`: Path to the ENVI reflectance raster file. +- `reflectance_uncertainty_file::String`: Path to the uncertainty raster file corresponding +to the reflectance data. An empty string indicates no uncertainty data. +- `line::Int64`: Line (row) of the raster data to be loaded. +- `good_bands::Array{Bool}`: Boolean array indicating which bands of the data to keep and +process. +- `refl_nodata::Float64`: The no-data value used in the reflectance data. + +# Returns +- A tuple containing: + - `img_dat::Array{Float64}`: The filtered reflectance data for the specified line and + good bands. If no valid data exists, this will be `nothing`. + - `unc_dat::Array{Float64}`: The uncertainty data, filtered similarly. If no uncertainty + data is loaded or valid, this will be `nothing`. + - `good_data`: A boolean array indicating which pixels contain valid reflectance data. +""" +function load_line(reflectance_file::String, reflectance_uncertainty_file::String, + line::Int64, good_bands::Array{Bool}, refl_nodata::Float64) + + reflectance_dataset = ArchGDAL.read(reflectance_file, alloweddrivers=["ENVI"]) + img_dat = convert( + Array{Float64}, + ArchGDAL.readraster(reflectance_file, alloweddrivers=["ENVI"])[:, line, :] + ) + img_dat = img_dat[:, good_bands] + good_data = .!all(img_dat .== refl_nodata, dims=2)[:, 1] + img_dat = img_dat[good_data, :] + + if sum(good_data) >= 1 + if reflectance_uncertainty_file != "" + unc_dat = convert( + Array{Float64}, + ArchGDAL.readraster( + reflectance_uncertainty_file, alloweddrivers=["ENVI"] + )[:, line, :] + ) + unc_dat = unc_dat[:, good_bands] + unc_dat = unc_dat[good_data, :] else - return nothing, nothing, good_data + unc_dat = nothing end + else + return nothing, nothing, good_data + end - return img_dat, unc_dat, good_data + return img_dat, unc_dat, good_data end diff --git a/src/EndmemberLibrary.jl b/src/EndmemberLibrary.jl index b6fbfce..a61bdf3 100644 --- a/src/EndmemberLibrary.jl +++ b/src/EndmemberLibrary.jl @@ -25,24 +25,49 @@ using NMF using MultivariateStats using Clustering +""" + nanargmin(input::Array) + +Return the index of the minimum value in the input array, ignoring NaN values. +""" function nanargmax(input::Array) x = copy(input) - x[isnan.(x)] .= -Inf; - return argmax(x); + x[isnan.(x)] .= -Inf + return argmax(x) end +""" + nanargmin(input::Array) + +Return the index of the maximum value in the input array, ignoring NaN values. +""" function nanargmin(input::Array) x = copy(input) - x[isnan.(x)] .= Inf; - return argmin(x); + x[isnan.(x)] .= Inf + return argmin(x) end +""" + read_envi_wavelengths(filename::String, nm::Bool=true) + +Read wavelength data from an ENVI header file with option to convert from microns to +nanometers. + +- If no wavelengths are found, returns `nothing`. + +# Notes +- The function constructs the expected header filename by replacing the extension of the +given filename with `.hdr`. +- This function logs an error if no wavelength data is found in the header file. +""" function read_envi_wavelengths(filename::String, nm::Bool=true) header_name = splitext(filename)[1] * ".hdr" header = readlines(header_name) found = false for line in header - if occursin("wavelength = {", line) || occursin("wavelength= {", line) || occursin("wavelength={", line) + if occursin("wavelength = {", line) || + occursin("wavelength= {", line) || + occursin("wavelength={", line) header = line found = true break @@ -54,32 +79,46 @@ function read_envi_wavelengths(filename::String, nm::Bool=true) return nothing end - wavelengths = [parse(Float64, strip(x)) for x in split(split(split(header, "{")[2], "}")[1],",")] + wavelengths = [ + parse(Float64, strip(x)) for x in split(split(split(header, "{")[2], "}")[1], ",") + ] if nm && all(wavelengths .< 25) wavelengths = wavelengths .* 1000 @info "Converting wavelengths read from $filename to nm from microns." end - + return wavelengths end +""" + get_good_bands_mask(wavelengths::Array{Float64}, wavelength_pairs) + +Return a boolean mask for `wavelengths` *outside* the specified wavelength ranges given +by `wavelength_pairs`. +""" function get_good_bands_mask(wavelengths::Array{Float64}, wavelength_pairs) good_bands = ones(Bool, length(wavelengths)) - + for wvp in wavelength_pairs wavelength_diff = wavelengths .- wvp[1] - wavelength_diff[wavelength_diff .< 0] .= maximum(filter(!isnan, wavelength_diff)) + wavelength_diff[wavelength_diff.<0] .= maximum(filter(!isnan, wavelength_diff)) lower_index = nanargmin(wavelength_diff) wavelength_diff = wvp[2] .- wavelengths - wavelength_diff[wavelength_diff .< 0] .= maximum(filter(!isnan, wavelength_diff)) + wavelength_diff[wavelength_diff.<0] .= maximum(filter(!isnan, wavelength_diff)) upper_index = nanargmin(wavelength_diff) good_bands[lower_index:upper_index] .= false end return good_bands end +""" + mutable struct SpectralLibrary + +A structure for managing spectral data from a file, supporting class labels, scaling, and +wavelength filtering. +""" mutable struct SpectralLibrary file_name::String class_header_name::String @@ -88,49 +127,97 @@ mutable struct SpectralLibrary class_valid_keys scale_factor::Float64 wavelength_regions_ignore - SpectralLibrary(file_name::String, class_header_name::String, spectral_starting_column::Int64 = 2, truncate_end_columns::Int64 = 0, class_valid_keys = nothing, - scale_factor = 1.0, wavelength_regions_ignore= [0,440,1310,1490,1770,2050,2440,2880]) = - new(file_name, class_header_name, spectral_starting_column, truncate_end_columns, class_valid_keys, scale_factor, wavelength_regions_ignore) - + @doc """ + SpectralLibrary(file_name::String, + class_header_name::String, + spectral_starting_column::Int64=2, + truncate_end_columns::Int64=0, + class_valid_keys=nothing, + scale_factor=1.0, + wavelength_regions_ignore=[0, 440, 1310, 1490, 1770, 2050, 2440, 2880]) + + Constructor for the `SpectralLibrary` struct. + + # Arguments + - `file_name::String`: Endmember Library data file name, CSV-like. + - `class_header_name::String`: Column for class labels in data file. + - `spectral_starting_column::Int64=2`: Column for start of spectral data + - `truncate_end_columns::Int64=0`: Number of columns to ignore from end of data. + - `class_valid_keys=nothing`: List of valid class labels. + - `scale_factor::Float64=1.0`: Spectral scaling factor. + - `wavelength_regions_ignore=[0,440,1310,1490,1770,2050,2440,2880]`: List of wavelength + regions to ignore as beginning/end pairs in list. + """ + SpectralLibrary(file_name::String, + class_header_name::String, + spectral_starting_column::Int64=2, + truncate_end_columns::Int64=0, + class_valid_keys=nothing, + scale_factor=1.0, + wavelength_regions_ignore=[0, 440, 1310, 1490, 1770, 2050, 2440, 2880]) = + new(file_name, class_header_name, spectral_starting_column, + truncate_end_columns, class_valid_keys, scale_factor, wavelength_regions_ignore) spectra classes good_bands wavelengths end +""" + load_data!(library::SpectralLibrary) + +Load and preprocess spectral data associated with the `SpectralLibrary` object. + +- Load spectral data from a CSV-like file `library.file_name`. +- Process wavelength ignore regions into paired sets `library.wavelength_regions_ignore`. +- Load `library.spectra`, `library.classes`, and `library.class_valid_keys`. +- Sort and load `library.wavelengths` as nanometers and rearrange spectral data accordingly. +- Generate mask `library.good_bands` based on `library.wavelength_regions_ignore`. + +# Returns +- The updated `SpectralLibrary` instance with loaded and processed data. +""" function load_data!(library::SpectralLibrary) df = DataFrame(CSV.File(library.file_name)) paired_sets = [] for i in 1:2:size(library.wavelength_regions_ignore)[1] - push!(paired_sets, [library.wavelength_regions_ignore[i],library.wavelength_regions_ignore[i+1]]) + push!(paired_sets, [library.wavelength_regions_ignore[i], + library.wavelength_regions_ignore[i+1]]) end library.wavelength_regions_ignore = paired_sets @info "Ignoring wavelength regions: " * string(library.wavelength_regions_ignore) try - library.spectra = Matrix(df[:,library.spectral_starting_column:end-library.truncate_end_columns]) + library.spectra = Matrix( + df[:, library.spectral_starting_column:end-library.truncate_end_columns]) catch e - throw(ArgumentError("Could not read spectra from endmember library. Try adjusting spectral_starting_column or truncate_end_columns, or reworking the library.")) + throw(ArgumentError("Could not read spectra from endmember library. Try adjusting + spectral_starting_column or truncate_end_columns, or reworking the library.")) end try - library.classes = convert(Array{AbstractString}, df[!,library.class_header_name]) + library.classes = convert(Array{AbstractString}, df[!, library.class_header_name]) catch e - throw(ArgumentError("Could not read classes from endmember library using class key "*library.class_header_name * - ". Try adjusting class_header_name, or reworking the library.")) + throw(ArgumentError("Could not read classes from endmember library using class key" + * library.class_header_name * + ". Try adjusting class_header_name, or reworking the library." + )) end try - library.wavelengths = parse.(Float64, names(df[:,library.spectral_starting_column:end-library.truncate_end_columns])) + library.wavelengths = parse.(Float64, names( + df[:, library.spectral_starting_column:end-library.truncate_end_columns]) + ) catch e - throw(ArgumentError("Could not read wavelengths from endmember library. Try adjusting class_header_name, or reworking the library.")) + throw(ArgumentError("Could not read wavelengths from endmember library. Try + adjusting class_header_name, or reworking the library.")) end wl_order = sortperm(library.wavelengths) - library.spectra = library.spectra[:,wl_order] + library.spectra = library.spectra[:, wl_order] library.wavelengths = library.wavelengths[wl_order] if isnothing(library.class_valid_keys) @@ -149,24 +236,47 @@ function load_data!(library::SpectralLibrary) if ignore_regions_nm && all(library.wavelengths .< 25) library.wavelengths = library.wavelengths .* 1000 - @info "Unit mismatch between library wavelengths and wavelength regions to ignore detected. Converting library wavelengths to nm from microns" + @info "Unit mismatch between library wavelengths and wavelength regions to ignore + detected. Converting library wavelengths to nm from microns." end return library end +""" + save_data(library::SpectralLibrary, output_filename::String; + class_label_name::String="Label") -function save_data(library::SpectralLibrary, output_filename::String; class_label_name::String = "Label") +Wrte the spectra and classes from a `SpectralLibrary` to a CSV file. + +- The output CSV will have the class labels as the first column under `class_label_name`, +followed by columns for each wavelength with the spectra values filling the corresponding +rows. +""" +function save_data(library::SpectralLibrary, output_filename::String; + class_label_name::String="Label") df = DataFrame() insertcols!(df, 1, class_label_name => library.classes) for (_wv, wv) in enumerate(library.wavelengths) - insertcols!(df, 1 + _wv, string(wv) => library.spectra[:,_wv]) + insertcols!(df, 1 + _wv, string(wv) => library.spectra[:, _wv]) end CSV.write(output_filename, df) - end +""" + filter_by_class!(library::SpectralLibrary) + +Filter spectra in the `SpectralLibrary` based on `library.class_valid_keys`. + +- Update `library.spectra` and `library.classes` to include only those matching `library. +class_valid_keys`. + +# Notes +- This function modifies the `library` in place. +- If `library.class_valid_keys` is `nothing`, logs a message indicating that no filtering +will occur, and exits without making any changes. +""" function filter_by_class!(library::SpectralLibrary) if isnothing(library.class_valid_keys) @info "No class valid keys provided, no filtering occuring" @@ -175,16 +285,26 @@ function filter_by_class!(library::SpectralLibrary) valid_classes = zeros(Bool, size(library.spectra)[1]) for cla in library.class_valid_keys - valid_classes[library.classes .== cla] .= true + valid_classes[library.classes.==cla] .= true end - library.spectra = library.spectra[valid_classes,:] + library.spectra = library.spectra[valid_classes, :] library.classes = library.classes[valid_classes] end -function remove_wavelength_region_inplace!(library::SpectralLibrary, set_as_nans::Bool=false) +""" + remove_wavelength_region_inplace!(library::SpectralLibrary, set_as_nans::Bool=false) + +Remove wavelength regions from `library` outside `library.good_bands` either by setting +them to NaN or by filtering them out. + +# Notes +- This function modifies the `library` in place. +""" +function remove_wavelength_region_inplace!(library::SpectralLibrary, + set_as_nans::Bool=false) if set_as_nans - library.spectra[:,.!library.good_bands] .= NaN + library.spectra[:, .!library.good_bands] .= NaN library.wavelengths[.!library.good_bands] .= NaN else library.spectra = library.spectra[:, library.good_bands] @@ -193,13 +313,29 @@ function remove_wavelength_region_inplace!(library::SpectralLibrary, set_as_nans end end -function interpolate_library_to_new_wavelengths!(library::SpectralLibrary, new_wavelengths::Array{Float64}) +""" + interpolate_library_to_new_wavelengths!(library::SpectralLibrary, new_wavelengths::Array + {Float64}) + +Interpolate `library` spectra to new wavelengths. + +- Perform linear interpolation, applying a flat extrapolation boundary condition. +- Update `library.wavelengths`, `library.good_bands`, and `library.spectra` matrix to the +resampled spectra. + +# Notes +- This function modifies the `library` in place +""" +function interpolate_library_to_new_wavelengths!(library::SpectralLibrary, + new_wavelengths::Array{Float64}) old_spectra = copy(library.spectra) library.spectra = zeros((size(library.spectra)[1], length(new_wavelengths))) for _s in 1:size(old_spectra)[1] - fit = LinearInterpolation(library.wavelengths, old_spectra[_s,:], extrapolation_bc=Flat()); - library.spectra[_s,:] = fit(new_wavelengths) + fit = LinearInterpolation( + library.wavelengths, old_spectra[_s, :], extrapolation_bc=Flat() + ) + library.spectra[_s, :] = fit(new_wavelengths) end library.wavelengths = new_wavelengths @@ -207,19 +343,54 @@ function interpolate_library_to_new_wavelengths!(library::SpectralLibrary, new_w library.good_bands = good_bands end +""" + scale_library!(library::SpectralLibrary, scaling_factor=nothing) + +Scale the spectral data in the `SpectralLibrary` by the library's defined scale factor +(default) or a specified scaling factor. +""" function scale_library!(library::SpectralLibrary, scaling_factor=nothing) if isnothing(scaling_factor) - library.spectra = library.spectra ./ library.scale_factor + library.spectra = library.spectra ./ library.scale_factor else library.spectra /= scaling_factor end end +""" + brightness_normalize!(library::SpectralLibrary) + +Normalize `library.spectra` based on the RMS brightness of `library.good_bands`. + +# Notes +- This function modifies the `library` in place. +- Only considers `library.good_bands` in the brightness calculation. +""" function brightness_normalize!(library::SpectralLibrary) - library.spectra = library.spectra ./ sqrt.(mean(library.spectra[:,library.good_bands].^2, dims=2)) + library.spectra = library.spectra ./ + sqrt.(mean(library.spectra[:, library.good_bands] .^ 2, dims=2)) end +""" + split_library(library::SpectralLibrary, split_fraction::Float64) + +Split a `SpectralLibrary` into two new libraries based on a specified fraction of the +total spectra. + +# Returns +- `Tuple{SpectralLibrary, SpectralLibrary}` with each output library a random, +mutually exclusive subset of the original library containing `split_fraction` and `1 - +split_fraction` of the spectra, respectively. + +# Notes +- The split is random; consecutive calls with the same `library` may yield different +results. +""" function split_library(library::SpectralLibrary, split_fraction::Float64) + + if !(0 < split_fraction < 1) + throw(ArgumentError("split_fraction must be between 0 and 1 (exclusive).")) + end perm = randperm(size(library.spectra)[1]) split_1 = perm[1:Int(round(split_fraction * length(perm)))] @@ -228,8 +399,8 @@ function split_library(library::SpectralLibrary, split_fraction::Float64) output_library_1 = deepcopy(library) output_library_2 = deepcopy(library) - output_library_1.spectra = output_library_1.spectra[split_1,:] - output_library_2.spectra = output_library_2.spectra[split_2,:] + output_library_1.spectra = output_library_1.spectra[split_1, :] + output_library_2.spectra = output_library_2.spectra[split_2, :] output_library_1.classes = output_library_1.classes[split_1] output_library_2.classes = output_library_2.classes[split_2] @@ -238,62 +409,109 @@ function split_library(library::SpectralLibrary, split_fraction::Float64) end +""" + reduce_endmembers_nmf!(library::SpectralLibrary, max_endmembers_per_class::Int64) + +Reduce the number of endmembers in the spectral library using Non-negative Matrix +Factorization (NMF). + +- For each class in `library.class_valid_keys`, apply NMF to the spectra subset to identify +the specified maximum number of endmembers. +- Update `library.classes` and `library.spectra` to contain the reduced set of endmembers. + +# Notes +- This function modifies the `library` in place. +- Only `library.good_bands` will be considered in the NMF computation. +- The NMF uses a maximum of 500 iterations with a tolerance of 1.0e-2 for convergence. +""" function reduce_endmembers_nmf!(library::SpectralLibrary, max_endmembers_per_class::Int64) reduced_library = [] - new_classes = copy(library.classes)[1:max_endmembers_per_class*size(library.class_valid_keys)[1]] + new_classes = + copy(library.classes)[1:max_endmembers_per_class*size(library.class_valid_keys)[1]] for (_c, cla) in enumerate(library.class_valid_keys) - library_subset = library.spectra[library.classes .== cla, library.good_bands] + library_subset = library.spectra[library.classes.==cla, library.good_bands] out_spectra = zeros(max_endmembers_per_class, size(library.spectra)[2]) - out_spectra[:,:] .= NaN + out_spectra[:, :] .= NaN r = nnmf(library_subset, max_endmembers_per_class; maxiter=500, tol=1.0e-2) - out_spectra[:,library.good_bands] = r.H + out_spectra[:, library.good_bands] = r.H - push!(reduced_library,out_spectra) - new_classes[(_c-1) * max_endmembers_per_class + 1:_c * max_endmembers_per_class] .= cla + push!(reduced_library, out_spectra) + new_classes[(_c-1)*max_endmembers_per_class+1:_c*max_endmembers_per_class] .= cla end - library.spectra = cat(reduced_library...,dims=1) + library.spectra = cat(reduced_library..., dims=1) library.classes = new_classes end +""" + reduce_endmembers_pca!(library::SpectralLibrary, max_endmembers_per_class::Int64) + +Reduce the number of endmembers in the spectral library using Principal Component Analysis +(PCA). + +- For each class in `library.class_valid_keys`, apply PCA to the spectra subset to identify +the specified maximum number of endmembers. +- Update `library.classes` and `library.spectra` to contain the reduced set of endmembers. + +# Notes +- This function modifies the `library` in place. +- Only `library.good_bands` will be considered in the PCA computation. +""" function reduce_endmembers_pca!(library::SpectralLibrary, max_endmembers_per_class::Int64) reduced_library = [] - new_classes = copy(library.classes)[1:max_endmembers_per_class*size(library.class_valid_keys)[1]] + new_classes = + copy(library.classes)[1:max_endmembers_per_class*size(library.class_valid_keys)[1]] for (_c, cla) in enumerate(library.class_valid_keys) - library_subset = library.spectra[library.classes .== cla, library.good_bands] + library_subset = library.spectra[library.classes.==cla, library.good_bands] out_spectra = zeros(max_endmembers_per_class, size(library.spectra)[2]) - out_spectra[:,:] .= NaN + out_spectra[:, :] .= NaN model = fit(PCA, library_subset', maxoutdim=max_endmembers_per_class, pratio=1) - out_spectra[:,library.good_bands] = projection(model)' + out_spectra[:, library.good_bands] = projection(model)' - push!(reduced_library,out_spectra) - new_classes[(_c-1) * max_endmembers_per_class + 1:_c * max_endmembers_per_class] .= cla + push!(reduced_library, out_spectra) + new_classes[(_c-1)*max_endmembers_per_class+1:_c*max_endmembers_per_class] .= cla end - library.spectra = cat(reduced_library...,dims=1) + library.spectra = cat(reduced_library..., dims=1) library.classes = new_classes end -function reduce_endmembers_kmeans!(library::SpectralLibrary, max_endmembers_per_class::Int64) +""" + reduce_endmembers_kmeans!(library::SpectralLibrary, max_endmembers_per_class::Int64) + +Reduce the number of endmembers in the spectral library using K-means clustering. + +- For each class in `library.class_valid_keys`, apply K-means clustering to the spectra +subset to identify the specified maximum number of endmembers from cluster centers. +- Update `library.classes` and `library.spectra` to contain the reduced set of endmembers. + +# Notes +- This function modifies the `library` in place. +- Only `library.good_bands` will be considered in the K-means computation. +- The K-means algorithm is run with a maximum of 1000 iterations. +""" +function reduce_endmembers_kmeans!(library::SpectralLibrary, + max_endmembers_per_class::Int64) reduced_library = [] - new_classes = copy(library.classes)[1:max_endmembers_per_class*size(library.class_valid_keys)[1]] + new_classes = + copy(library.classes)[1:max_endmembers_per_class*size(library.class_valid_keys)[1]] for (_c, cla) in enumerate(library.class_valid_keys) - library_subset = library.spectra[library.classes .== cla, library.good_bands] + library_subset = library.spectra[library.classes.==cla, library.good_bands] out_spectra = zeros(max_endmembers_per_class, size(library.spectra)[2]) - out_spectra[:,:] .= NaN + out_spectra[:, :] .= NaN model = kmeans(library_subset', max_endmembers_per_class, maxiter=1000) - out_spectra[:,library.good_bands] = model.centers' + out_spectra[:, library.good_bands] = model.centers' - push!(reduced_library,out_spectra) - new_classes[(_c-1) * max_endmembers_per_class + 1:_c * max_endmembers_per_class] .= cla + push!(reduced_library, out_spectra) + new_classes[(_c-1)*max_endmembers_per_class+1:_c*max_endmembers_per_class] .= cla end - library.spectra = cat(reduced_library...,dims=1) + library.spectra = cat(reduced_library..., dims=1) library.classes = new_classes end diff --git a/src/Plotting.jl b/src/Plotting.jl index 7c16382..448e6da 100644 --- a/src/Plotting.jl +++ b/src/Plotting.jl @@ -1,13 +1,35 @@ +# Copyright 2022 California Institute of Technology +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# Author: Philip G. Brodrick, philip.g.brodrick@jpl.nasa.gov using Plots using Statistics +""" + plot_mean_endmembers(endmember_library::SpectralLibrary; output_name::String="") -function plot_mean_endmembers(endmember_library; output_name = "") +Plot the mean spectra of endmembers of each class from a spectral library. + +- Saves file to `output_name` if provided, otherwise defaults to no saving. +""" +function plot_mean_endmembers(endmember_library::SpectralLibrary; output_name::String="") plots = [] for (_u, u) in enumerate(endmember_library.class_valid_keys) - mean_spectra = mean(endmember_library.spectra[endmember_library.classes .== u,:],dims=1)[:] - plot!(endmember_library.wavelengths, mean_spectra, label=u, xlim=[300,3200]) + mean_spectra = + mean(endmember_library.spectra[endmember_library.classes.==u, :], dims=1)[:] + plot!(endmember_library.wavelengths, mean_spectra, label=u, xlim=[300, 3200]) end xlabel!("Wavelength [nm]") ylabel!("Reflectance") @@ -19,20 +41,38 @@ function plot_mean_endmembers(endmember_library; output_name = "") plot!() end -function plot_endmembers(endmember_library::SpectralLibrary; output_name::String = "") +""" + plot_endmembers(endmember_library::SpectralLibrary; output_name::String="") + +Plot all endmember spectra from a spectral library. +- Endmember spectra are colored and labeled by class. +- Saves file to `output_name` if provided, otherwise defaults to no saving. +""" +function plot_endmembers(endmember_library::SpectralLibrary; output_name::String="") for (_u, u) in enumerate(endmember_library.class_valid_keys) if _u == 1 - plot(endmember_library.wavelengths, endmember_library.spectra[endmember_library.classes .== u,:]', lab="", xlim=[300,3200], color=palette(:tab10)[Mod(_u)],dpi=200) + plot(endmember_library.wavelengths, + endmember_library.spectra[endmember_library.classes.==u, :]', + lab="", + xlim=[300, 3200], + color=palette(:tab10)[Mod(_u)], + dpi=200 + ) else - plot!(endmember_library.wavelengths, endmember_library.spectra[endmember_library.classes .== u,:]', lab="",xlim=[300,3200], color=palette(:tab10)[Mod(_u)]) + plot!(endmember_library.wavelengths, + endmember_library.spectra[endmember_library.classes.==u, :]', + lab="", + xlim=[300, 3200], + color=palette(:tab10)[Mod(_u)] + ) end end xlabel!("Wavelenth [nm]") ylabel!("Reflectance") xticks!([500, 1000, 1500, 2000, 2500, 3000]) for (_u, u) in enumerate(endmember_library.class_valid_keys) - plot!([1:2],[0,0.3], color=palette(:tab10)[Mod(_u)], labels=u) + plot!([1:2], [0, 0.3], color=palette(:tab10)[Mod(_u)], labels=u) end if output_name != "" savefig(output_name) @@ -40,23 +80,48 @@ function plot_endmembers(endmember_library::SpectralLibrary; output_name::String plot!() end -function plot_endmembers_individually(endmember_library::SpectralLibrary; output_name::String = "", legend_on::Bool = false) +""" + plot_endmembers_individually(endmember_library::SpectralLibrary; + output_name::String="", legend_on::Bool=false) + +Plot all endmember spectra individually from a spectral library, grouped by class. + +- Create seperate plots for each class, legends are disabled by default. +- Saves file to `output_name` if provided, otherwise defaults to no saving. +""" +function plot_endmembers_individually(endmember_library::SpectralLibrary; + output_name::String="", legend_on::Bool=false) + plots = [] spectra = endmember_library.spectra classes = endmember_library.classes for (_u, u) in enumerate(endmember_library.class_valid_keys) - sp = spectra[classes .== u,:] - sp[broadcast(isnan,sp)] .= 0 - toplot = spectra[classes .== u,:] + sp = spectra[classes.==u, :] + sp[broadcast(isnan, sp)] .= 0 + toplot = spectra[classes.==u, :] if size(toplot)[2] == sum(endmember_library.good_bands) - push!(plots, plot(endmember_library.wavelengths[endmember_library.good_bands], toplot', title=u, xlabel="Wavelength [nm]", ylabel="Reflectance")) - else - push!(plots, plot(endmember_library.wavelengths, toplot', title=u, xlabel="Wavelength [nm]", ylabel="Reflectance")) + push!(plots, + plot(endmember_library.wavelengths[endmember_library.good_bands], + toplot', + title=u, + xlabel="Wavelength [nm]", + ylabel="Reflectance" + ) + ) + else + push!(plots, + plot(endmember_library.wavelengths, + toplot', + title=u, + xlabel="Wavelength [nm]", + ylabel="Reflectance" + ) + ) end xticks!([500, 1000, 1500, 2000, 2500]) end - plot(plots...,size=(1000,600),dpi=200) + plot(plots..., size=(1000, 600), dpi=200) if legend_on plot!(legend=:outertopright) else diff --git a/src/Solvers.jl b/src/Solvers.jl index 9f7afd4..f676ace 100644 --- a/src/Solvers.jl +++ b/src/Solvers.jl @@ -18,45 +18,119 @@ using JuMP using NLopt using LinearAlgebra -function opt_solve(A::Matrix{Float64}, b::Vector{Float64}, x0::Vector{Float64}, lb::Vector{Float64}, ub::Vector{Float64}) +""" + opt_solve(A::Matrix{Float64}, b::Vector{Float64}, x0::Vector{Float64}, + lb::Vector{Float64}, ub::Vector{Float64}) + +Solve a nonlinear least squares problem, finding the vector `x` subject to lower (`lb`) +and upper (`ub`) bounds, that minimizes the residual ||Ax - b||² using the NLopt library's +implementation of the SLSQP algorithm. + +# Arguments +- `A::Matrix{Float64}`: Coefficient matrix of size (m, n). +- `b::Vector{Float64}`: Target vector of size (m,). +- `x0::Vector{Float64}`: Initial guess vector of size (n,) for the optimization variables. +Values will be clipped to the range [0, 1]. +- `lb::Vector{Float64}`: Vector of size (n,) specifying the lower bounds for the +optimization variables. +- `ub::Vector{Float64}`: Vector of size (n,) specifying the upper bounds for the +optimization variables. + +# Returns +- A tuple containing: + - `x`: The optimized values of the variables `x` (vector of size (n,)). + - `mse_opt`: The exponential of the objective function value at the optimum. +""" +function opt_solve(A::Matrix{Float64}, b::Vector{Float64}, x0::Vector{Float64}, + lb::Vector{Float64}, ub::Vector{Float64}) #mle = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)) mle = Model(NLopt.Optimizer) set_optimizer_attribute(mle, "algorithm", :LD_SLSQP) #LD_SLSQP - x0[x0 .< 0] .= 0 - x0[x0 .> 1] .= 1 + x0[x0.<0] .= 0 + x0[x0.>1] .= 1 @variable(mle, lb[i] <= x[i=1:length(x0)] <= ub[i]) #@constraint(mle, sum(x) == 1) - for n in 1:length(x0) set_start_value(x[n], x0[n]) end + for n in 1:length(x0) + set_start_value(x[n], x0[n]) + end #@NLexpression(mle, mse , sum( (b .- sum(x[i] .* A[:,i] for i in 1:length(x) )).^2 ) ) - @expression(mle, mse , sum( (b .- sum(x[i] .* A[:,i] for i in 1:length(x) )).^2 ) ) - + @expression(mle, mse, sum((b .- sum(x[i] .* A[:, i] for i in 1:length(x))) .^ 2)) + @NLobjective(mle, Min, mse) JuMP.optimize!(mle) - return value.(mle[:x]) , exp(objective_value(mle)) + return value.(mle[:x]), exp(objective_value(mle)) end +""" + dolsq(A, b; method::String="default") + +Solve a least squares problem, finding the vector `x` that minimizes the residual +||Ax - b||². +# Arguments +- `A`: Coefficient matrix of size (m, n). +- `b`: Target vector of size (m,). +- `method`: An optional keyword argument specifying the solver method. +Available options are: + - `"default"`: Solve the system using the backslash operator (`\`). + - `"pinv"`: Use the pseudoinverse of `A` to compute the solution. + - `"qr"`: Use QR decomposition to solve the least squares problem. + +# Returns +- `x`: Vector of size (n,) that minimizes the least squares residual. +""" function dolsq(A, b; method::String="default") if method == "default" x = A \ b elseif method == "pinv" - x = pinv(A)*b + x = pinv(A) * b elseif method == "qr" #Q,R = qr(A) #x = inv(R)*Q'*b qrA = qr(A) - x = qrA\b + x = qrA \ b end return x end +""" + bvls(A, b, x_lsq, lb, ub, tol::Float64, max_iter::Int64, verbose::Int64, + inverse_method::String) + +Solve a bounded value least squares problem, finding the vector `x` subject to lower (`lb`) +and upper (`ub`) bounds, that minimizes the residual ||Ax - b||². + +- The function employs an iterative approach, adjusting the solution based on the specified + bounds. +- Convergence is determined by whether KKT optimality condition (see +[`compute_kkt_optimality`](@ref)) is within the specified tolerance. + +# Arguments +- `A`: Coefficient matrix of size (m, n). +- `b`: Target vector of size (m,). +- `x_lsq`: Initial guess vector of size (n,) for the optimization variables. +- `lb`: Vector of size (n,) specifying the lower bounds for the +optimization variables. +- `ub`: Vector of size (n,) specifying the upper bounds for the +optimization variables. +- `tol::Float64`: Tolerance for convergence of the optimization. +- `max_iter::Int64`: Maximum number of iterations allowed. If not specified, defaults to n. +- `verbose::Int64`: Verbosity of output (0, 1, or 2). +- `inverse_method::String`: Least squares solver method, see [`dolsq`](@ref) for options. + +# Returns +- A tuple containing: + - `x`: Vector of size (n,) that minimizes the least squares residual. + - `cost`: The final cost function value of the minimized sum of squared residuals. +""" +function bvls(A, b, x_lsq, lb, ub, tol::Float64, max_iter::Int64, verbose::Int64, + inverse_method::String) -function bvls(A, b, x_lsq, lb, ub, tol::Float64, max_iter::Int64, verbose::Int64, inverse_method::String) n_iter = 0 m, n = size(A) @@ -73,10 +147,10 @@ function bvls(A, b, x_lsq, lb, ub, tol::Float64, max_iter::Int64, verbose::Int64 free_set = on_bound .== 0 active_set = .!free_set - free_set = (1:size(free_set)[1])[free_set .!= 0] + free_set = (1:size(free_set)[1])[free_set.!=0] r = A * x - b - cost = 0.5 * dot(r , r) + cost = 0.5 * dot(r, r) initial_cost = cost g = A' * r @@ -118,12 +192,12 @@ function bvls(A, b, x_lsq, lb, ub, tol::Float64, max_iter::Int64, verbose::Int64 ind = free_set[.!v] x[ind] = z[.!v] - r = A * x -b - cost_new = 0.5 * dot(r , r) + r = A * x - b + cost_new = 0.5 * dot(r, r) cost_change = cost - cost_new cost = cost_new g = A' * r - step_norm = sum((x[free_set] .- x_free_old).^2) + step_norm = sum((x[free_set] .- x_free_old) .^ 2) if any(v) free_set = free_set[.!v] @@ -159,7 +233,7 @@ function bvls(A, b, x_lsq, lb, ub, tol::Float64, max_iter::Int64, verbose::Int64 free_set = on_bound .== 0 sum(free_set) active_set = .!free_set - free_set = (1:size(free_set)[1])[free_set .!= 0] + free_set = (1:size(free_set)[1])[free_set.!=0] x_free = x[free_set] x_free_old = copy(x_free) @@ -170,12 +244,16 @@ function bvls(A, b, x_lsq, lb, ub, tol::Float64, max_iter::Int64, verbose::Int64 b_free = b - A * (x .* active_set) z = dolsq(A_free, b_free, method=inverse_method) - lbv = (1:size(free_set)[1])[ z .< lb_free] - ubv = (1:size(free_set)[1])[ z .> ub_free] + lbv = (1:size(free_set)[1])[z.ub_free] v = cat(lbv, ubv, dims=1) if size(v)[1] > 0 - alphas = cat(lb_free[lbv] - x_free[lbv], ub_free[ubv] - x_free[ubv],dims=1) ./ (z[v] - x_free[v]) + alphas = cat( + lb_free[lbv] - x_free[lbv], + ub_free[ubv] - x_free[ubv], + dims=1 + ) ./ (z[v] - x_free[v]) i = argmin(alphas) i_free = v[i] @@ -198,10 +276,10 @@ function bvls(A, b, x_lsq, lb, ub, tol::Float64, max_iter::Int64, verbose::Int64 end end #while @label start - step_norm = sum((x_free .- x_free_old).^2) + step_norm = sum((x_free .- x_free_old) .^ 2) r = A * x - b - cost_new = 0.5 * dot(r , r) + cost_new = 0.5 * dot(r, r) cost_change = cost - cost_new combo = tol * cost @@ -218,16 +296,30 @@ function bvls(A, b, x_lsq, lb, ub, tol::Float64, max_iter::Int64, verbose::Int64 termination_status = 0 end - x[x .< 1e-5] .= 0 + x[x.<1e-5] .= 0 return x, cost end -function compute_kkt_optimality(g::Vector{Float64}, on_bound::Vector) - g_kkt = g .* on_bound - free_set = on_bound .== 0 - g_kkt[free_set] = broadcast(abs, g[free_set]) +""" + compute_kkt_optimality(g::Vector{Float64}, on_bound::Vector) - return maximum(g_kkt) +Computes the Karush-Kuhn-Tucker (KKT) optimality condition value for a given gradient +vector and a vector indicating which variables are on bounds. -end +- Returns a value representing the maximum KKT condition across all variables +# Arguments +- `g`: A vector of size (n,) representing the gradient of the objective function + at the current point. +- `on_bound`: A vector of size (n,) indicating the status of each variable: + - `-1` if the variable is at its lower bound, + - `1` if the variable is at its upper bound, + - `0` if the variable is free (not constrained). +""" +function compute_kkt_optimality(g::Vector{Float64}, on_bound::Vector) + g_kkt = g .* on_bound + free_set = on_bound .== 0 + g_kkt[free_set] = broadcast(abs, g[free_set]) + + return maximum(g_kkt) +end diff --git a/src/SpectralUnmixing.jl b/src/SpectralUnmixing.jl index ebb3447..a539772 100644 --- a/src/SpectralUnmixing.jl +++ b/src/SpectralUnmixing.jl @@ -45,42 +45,92 @@ export initiate_output_datasets, set_band_names, write_results # Unmixing and simulation functions export unmix_line, unmix_pixel, simulate_pixel, unmix_and_write_line +""" + wl_index(wavelengths::Vector{Float64}, target::Float64) + +Return the index of the wavelength in a vector that is closest to the specified target +wavelength. +""" function wl_index(wavelengths::Vector{Float64}, target) argmin(abs.(wavelengths .- target)) end +""" + scale_data(refl::Matrix{Float64}, wavelengths::Vector{Float64}, criteria::String) + +Scale the reflectance data based on specified criteria. + +# Arguments +- `refl::Matrix{Float64}`: Reflectance data either size (n,) or (m,n), where n is the +number of bands and m the number of pixels. +- `wavelengths::Vector{Float64}`: A vector of wavelengths of size (n,). +- `criteria::String`: Normalization options: + - `"none"`: No scaling will be applied, and the original reflectance data will be + returned. + - `"brightness"`: The data will be normalized based on brightness, excluding specified + bad wavelength regions [1300, 1500] and [1800, 2000]. + - A specific wavelength (as a string): The data will be normalized using the reflectance + value at this wavelength. + +# Returns +- A matrix of scaled reflectance data with same dimensions as `refl`. +""" function scale_data(refl::Matrix{Float64}, wavelengths::Vector{Float64}, criteria::String) if criteria == "none" return refl elseif criteria == "brightness" - bad_regions_wl = [[1300,1500],[1800,2000]] + bad_regions_wl = [[1300, 1500], [1800, 2000]] good_bands = convert(Vector{Bool}, ones(length(wavelengths))) for br in bad_regions_wl good_bands[wl_index(wavelengths, br[1]):wl_index(wavelengths, br[2])] .= false end if length(size(refl)) == 2 - norm = sqrt.(mean(refl[:,good_bands].^2, dims=2)) + norm = sqrt.(mean(refl[:, good_bands] .^ 2, dims=2)) else - norm = sqrt.(mean(refl[good_bands].^2)) + norm = sqrt.(mean(refl[good_bands] .^ 2)) end else try - target_wl = parse(Float64,criteria) + target_wl = parse(Float64, criteria) if length(size(refl)) == 2 - norm = refl[:,wl_index(wavelengths, target_wl)] ./ 0.5 + norm = refl[:, wl_index(wavelengths, target_wl)] ./ 0.5 else norm = refl[wl_index(wavelengths, target_wl)] ./ 0.5 end catch e - throw(ArgumentError(string("normalization must be [none, brightness, or a specific wavelength]. Provided:", criteria))) + throw(ArgumentError(string("normalization must be [none, brightness, or a + specific wavelength]. Provided:", criteria))) end end return refl ./ norm end -function get_sma_permutation(class_idx, num_endmembers::Vector{Int64}, combination_type::String, library_length::Int64) +""" + get_sma_permutation(class_idx, num_endmembers::Vector{Int64}, + combination_type::String, library_length::Int64) + +Generate a permutation of endmember indices for SMA. + +# Arguments +- `class_idx`: A collection (such as a vector of vectors) of size n where each element +is a set of indices corresponding to endmembers belonging to one class. +- `num_endmembers::Vector{Int64}`: Desired number of endmembers to select. The first +element specifies the number of endmembers to permute; if it is `-1`, all endmembers in +library will be included. +- `combination_type::String`: Permutation type options: + - `"class-even"`: Randomly select one endmember from each class until the required number + of endmembers is selected. + - Any other types: Randomly selects endmembers from the entire library. +- `library_length::Int64`: Size of the entire endmember library. + +# Returns +- A vector of integers representing the permuted endmember indices. +""" +function get_sma_permutation(class_idx, num_endmembers::Vector{Int64}, + combination_type::String, library_length::Int64) + if num_endmembers[1] != -1 if combination_type == "class-even" @@ -94,10 +144,10 @@ function get_sma_permutation(class_idx, num_endmembers::Vector{Int64}, combinati while selector <= num_endmembers[1] _p = mod(selector, length(perm_class_idx)) + 1 push!(perm, perm_class_idx[_p][1]) - deleteat!(perm_class_idx[_p],1) + deleteat!(perm_class_idx[_p], 1) if length(perm_class_idx[_p]) == 0 - deleteat!(perm_class_idx,_p) + deleteat!(perm_class_idx, _p) end selector += 1 end @@ -106,33 +156,100 @@ function get_sma_permutation(class_idx, num_endmembers::Vector{Int64}, combinati perm = randperm(library_length)[1:num_endmembers[1]] end else - perm = convert(Vector{Int64},1:library_length) + perm = convert(Vector{Int64}, 1:library_length) end return perm end +""" + results_from_mc(results::Matrix{Float64}, cost::Vector{Float64}, mode::String) + +Process Monte Carlo simulations and return variance and result based on `mode`. + +# Arguments +- `results::Matrix{Float64}`: Outcomes of MC simulations, where each row represents a +simulation and each column corresponds to a different variable. +- `cost::Vector{Float64}`: A vector of Float64 values representing the cost associated +with each simulation, used to determine the best result when applicable. +- `mode::String`: How to process the results. It can be: + - Contain `"best"`: Returns the result corresponding to the simulation with the lowest + cost. + - Any other value: Computes and returns the mean of the results across all simulations. + +# Returns +- A tuple containing: + - `output::Vector{Float64}`: MC simulation results + - `output_var::Vector{Float64}`: Standard deviation of the results for each variable, or + `nothing` if there is only one simulation. +""" function results_from_mc(results::Matrix{Float64}, cost::Vector{Float64}, mode::String) if size(results)[1] == 1 output_var = nothing else - output_var = std(results, dims=1)[1,:] + output_var = std(results, dims=1)[1, :] end if occursin("best", mode) best_idx = argmin(cost) - output = results[best_idx,:] + output = results[best_idx, :] else - output = mean(results, dims=1)[1,:] + output = mean(results, dims=1)[1, :] end - - return output, output_var + return output, output_var end -function unmix_pixel(library::SpectralLibrary, img_dat_input::Array{Float64}, unc_dat, class_idx, options, mode::String, n_mc::Int64, num_endmembers::Vector{Int64}, normalization::String, optimization::String, max_combinations::Int64, combination_type::String) - +""" + unmix_pixel(library::SpectralLibrary, img_dat_input::Array{Float64}, unc_dat, + class_idx, options, mode::String, n_mc::Int64, + num_endmembers::Vector{Int64}, normalization::String, optimization::String, + max_combinations::Int64, combination_type::String) + +Unmix a pixel's spectral data using a given spectral library with options for SMA or MESMA +and various Monte Carlo and optimization approaches. + +# Arguments +- `library::SpectralLibrary`: Endmember library for unmixing. +- `img_dat_input::Array{Float64}`: Spectral data of the pixel. +- `unc_dat`: An optional array of uncertainty in the pixel spectral data. +- `class_idx`: A collection of indices indicating class memberships for different +endmembers in `library.spectra`. See also [`prepare_combinations`](@ref), +[`prepare_options`](@ref) +- `options`: A collection of potential endmember combinations for the unmixing process. See +also [`prepare_options`](@ref). +- `mode::String`: Determines the unmixing approach: + - `"sma"`: Select endmembers randomly for MC unmixing, returns mean fractions. + - `"sma-best"`: SMA and output best (lowest cost) MC fraction. + - `"mesma"`: Evaluate different combinations of endmembers representing each class. + - `"mesma-best"`: MESMA and output best (lowest cost) MC fraction. +- `n_mc::Int64`: Number of MC iterations. +- `num_endmembers::Vector{Int64}`: The number of endmembers to consider. +- `normalization::String`: The normalization method to apply to the spectral data. See +[`scale_data`](@ref). +- `optimization::String`: The optimization approach for unmixing: (e.g., `"bvls"`, +`"ldsqp"`, or `"inverse"`). Can optionally contain `"pinv"` or `"qr"` to specify the +inverse method, if applicable. +- `max_combinations::Int64`: Maximum number of combinations to consider. +- `combination_type::String`: The type of combination (e.g., `"class-even"`). See also +[`prepare_combinations`](@ref), [`prepare_options`](@ref), [`get_sma_permutation`](@ref). + +# Returns +- A tuple containing: + - `output_mixture::Vector{Float64}`: The estimated fraction of each class in the + `library`, appended with the brightness. + - `output_mixture_var::Vector{Float64}`: The variance of each class in the `library`, + appended with the brightness variance. + - `output_comp_frac::Vector{Float64}`: The estimated fraction of each endmember in the + `library`, appended with the brightness. + - `output_comp_frac_var::Vector{Float64}`: The variance of each endmember in the + `library`, appended with the brightness variance. +""" +function unmix_pixel(library::SpectralLibrary, img_dat_input::Array{Float64}, unc_dat, + class_idx, options, mode::String, n_mc::Int64, num_endmembers::Vector{Int64}, + normalization::String, optimization::String, max_combinations::Int64, + combination_type::String) if length(size(img_dat_input)) == 1 img_dat = reshape(img_dat_input, 1, length(img_dat_input)) @@ -140,7 +257,7 @@ function unmix_pixel(library::SpectralLibrary, img_dat_input::Array{Float64}, un img_dat = img_dat_input end - mc_comp_frac = zeros(n_mc, size(library.spectra)[1]+1) + mc_comp_frac = zeros(n_mc, size(library.spectra)[1] + 1) scores = zeros(n_mc) for mc in 1:n_mc #monte carlo loop Random.seed!(mc) @@ -159,7 +276,9 @@ function unmix_pixel(library::SpectralLibrary, img_dat_input::Array{Float64}, un end if mode == "sma" || mode == "sma-best" - perm = get_sma_permutation(class_idx, num_endmembers, combination_type, size(library.spectra)[1]) + perm = get_sma_permutation( + class_idx, num_endmembers, combination_type, size(library.spectra)[1] + ) G = library.spectra[perm, library.good_bands] G = scale_data(G, library.wavelengths[library.good_bands], normalization)' @@ -169,35 +288,44 @@ function unmix_pixel(library::SpectralLibrary, img_dat_input::Array{Float64}, un x0 = x0[:] res = nothing if occursin("bvls", optimization) - res, cost = bvls(G, d[:], x0, zeros(size(x0)), ones(size(x0)), 1e-3, 100, 1, inverse_method) + res, cost = bvls( + G, d[:], x0, zeros(size(x0)), ones(size(x0)), 1e-3, 100, 1, + inverse_method + ) elseif occursin("ldsqp", optimization) - res, cost = opt_solve(G, d[:], x0, zeros(length(x0)), ones(length(x0)) ) + res, cost = opt_solve(G, d[:], x0, zeros(length(x0)), ones(length(x0))) elseif occursin("inverse", optimization) res = x0 r = G * x0 - d[:] - cost = dot(r,r) + cost = dot(r, r) end mc_comp_frac[mc, perm] = res scores[mc] = cost elseif mode == "mesma" || mode == "mesma-best" solutions = [] - costs = zeros(size(options)[1]).+1e12 + costs = zeros(size(options)[1]) .+ 1e12 if max_combinations != -1 && length(options) > max_combinations perm = randperm(length(options))[1:max_combinations] else - perm = convert(Vector{Int64},1:length(options)) + perm = convert(Vector{Int64}, 1:length(options)) end - + for (_comb, comb) in enumerate(options[perm]) comb = [c for c in comb] #G = hcat(library.spectra[comb,:], ones(size(library.spectra[comb,:])[1],1)) - G = scale_data(library.spectra[comb, library.good_bands], library.wavelengths[library.good_bands], normalization)' + G = scale_data( + library.spectra[comb, library.good_bands], + library.wavelengths[library.good_bands], normalization + )' x0 = dolsq(G, d') x0 = x0[:] ls = nothing if optimization == "bvls" - ls, lc = bvls(G, d[:], x0, zeros(size(x0)), ones(size(x0)), 1e-3, 10, 1, inverse_method) + ls, lc = bvls( + G, d[:], x0, zeros(size(x0)), ones(size(x0)), 1e-3, 10, 1, + inverse_method + ) costs[_comb] = lc elseif optimization == "ldsqp" ls, lc = opt_solve(G, d[:], x0, 0, 1) @@ -205,9 +333,9 @@ function unmix_pixel(library::SpectralLibrary, img_dat_input::Array{Float64}, un elseif optimization == "inverse" ls = x0 r = G * x0 - d[:] - costs[_comb] = dot(r,r) + costs[_comb] = dot(r, r) end - push!(solutions,ls) + push!(solutions, ls) end best = argmin(costs) @@ -221,37 +349,55 @@ function unmix_pixel(library::SpectralLibrary, img_dat_input::Array{Float64}, un end # Calculate the sum of values (inverse of shade), and then normalize - mc_comp_frac[mc_comp_frac .< 0] .= 0 - mc_comp_frac[:,end] = sum(mc_comp_frac,dims=2) - mc_comp_frac[:,1:end-1] = mc_comp_frac[:,1:end-1] ./ mc_comp_frac[:,end] + mc_comp_frac[mc_comp_frac.<0] .= 0 + mc_comp_frac[:, end] = sum(mc_comp_frac, dims=2) + mc_comp_frac[:, 1:end-1] = mc_comp_frac[:, 1:end-1] ./ mc_comp_frac[:, end] # Aggregate results from per-library to per-unique-class mixture_results = zeros(size(mc_comp_frac)[1], length(library.class_valid_keys) + 1) for _i in 1:size(mc_comp_frac)[1] for (_class, cl) in enumerate(library.class_valid_keys) - mixture_results[_i, _class] = sum(mc_comp_frac[_i,1:end-1][cl .== library.classes]) + mixture_results[_i, _class] = + sum(mc_comp_frac[_i, 1:end-1][cl.==library.classes]) end - mixture_results[_i, end] = mc_comp_frac[_i,end] + mixture_results[_i, end] = mc_comp_frac[_i, end] end output_mixture, output_mixture_var = results_from_mc(mixture_results, scores, mode) output_comp_frac, output_comp_frac_var = results_from_mc(mc_comp_frac, scores, mode) return output_mixture, output_mixture_var, output_comp_frac, output_comp_frac_var - end +""" + prepare_combinations(library::SpectralLibrary, combination_type::String) + +Prepare the set of class indices for the specified combination type. + +- Currently only supports `"class-even"`: returns the list of indices grouped by class. +""" function prepare_combinations(library::SpectralLibrary, combination_type::String) class_idx = [] if combination_type == "class-even" for uc in library.class_valid_keys - push!(class_idx, (1:size(library.classes)[1])[library.classes .== uc]) + push!(class_idx, (1:size(library.classes)[1])[library.classes.==uc]) end end return class_idx end -function prepare_options(library::SpectralLibrary, combination_type::String, num_endmembers::Vector{Int64}, class_idx) +""" + prepare_options(library::SpectralLibrary, combination_type::String, + num_endmembers::Vector{Int64}, class_idx) + +Prepare combinations of endmembers based on the specified combination type. + +- Combination type options: + - `"class-even"``: Generate all combinations where one endmember is selected from each class. + - `"all"`: Generate all possible combinations of `num_endmembers` spectra. +""" +function prepare_options(library::SpectralLibrary, combination_type::String, + num_endmembers::Vector{Int64}, class_idx) # Prepare combinations if relevant options = [] @@ -260,7 +406,7 @@ function prepare_options(library::SpectralLibrary, combination_type::String, num elseif combination_type == "all" for num in num_endmembers combo = [c for c in combinations(1:length(library.classes), num)] - push!(options,combo...) + push!(options, combo...) end else error("Invalid combiation string") @@ -269,27 +415,75 @@ function prepare_options(library::SpectralLibrary, combination_type::String, num return options end - - -function unmix_line(line::Int64, reflectance_file::String, mode::String, refl_nodata::Float64, - refl_scale::Float64, normalization::String, library::SpectralLibrary, - reflectance_uncertainty_file::String = "", n_mc::Int64 = 1, - combination_type::String = "all", num_endmembers::Vector{Int64} = [2,3], - max_combinations::Int64 = -1, optimization="bvls") +""" + unmix_line(line::Int64, reflectance_file::String, mode::String, + refl_nodata::Float64, refl_scale::Float64, + normalization::String, library::SpectralLibrary, + reflectance_uncertainty_file::String="", n_mc::Int64=1, + combination_type::String="all", + num_endmembers::Vector{Int64}=[2, 3], + max_combinations::Int64=-1, optimization="bvls") + +Unmix a specific line (row of pixels) of reflectance data. See also [`load_line`](@ref), +[`unmix_pixel`](@ref). + +# Arguments +- `line::Int64`: Index of the line to unmix. +- `reflectance_file::String`: Path to the reflectance data. +- `mode::String`: The mode of unmixing to be used (e.g., "sma", "mesma-best"). See +[`unmix_pixel`](@ref)). +- `refl_nodata::Float64`: The no data value in the reflectance file. +- `refl_scale::Float64`: Scaling factor for the reflectance data. +- `normalization::String`: The normalization method to apply to the reflectance data. See +[`scale_data`](@ref). +- `library::SpectralLibrary`: Spectral library object containing endmembers for unmixing. +- `reflectance_uncertainty_file::String`: Optional path to the reflectance uncertainty file. +- `n_mc::Int64=1`: Number of Monte Carlo iterations to perform. +- `combination_type::String="all"`: The type of endmember combinations to prepare for +unmixing. See also [`prepare_combinations`](@ref), [`prepare_options`](@ref). +- `num_endmembers::Vector{Int64}=[2, 3]`: The number of endmembers to consider in +combinations. +- `max_combinations::Int64=-1`: Maximum number of combinations to consider, defaults to no +limit. +- `optimization::String`: The optimization method to use (default is "bvls"). See also +[`unmix_pixel`](@ref). + +# Returns +- A tuple containing (see also [`unmix_pixel`](@ref)): + - The line index. + - Estimated class fractions for each pixel in the line. + - A boolean mask of pixels with valid data. + - Standard deviations of the mixture results (if applicable). + - Complete fractions for each endmember in the unmixing of each pixel. + +# Notes +- Initializes a random seed for reproducibility. +- Logs the execution time for unmixing a line. +- Returns `nothing` if input data is invalid or missing. +""" +function unmix_line(line::Int64, reflectance_file::String, mode::String, + refl_nodata::Float64, refl_scale::Float64, normalization::String, + library::SpectralLibrary, reflectance_uncertainty_file::String="", n_mc::Int64=1, + combination_type::String="all", num_endmembers::Vector{Int64}=[2, 3], + max_combinations::Int64=-1, optimization="bvls") Random.seed!(13) - - img_dat, unc_dat, good_data = load_line(reflectance_file, reflectance_uncertainty_file, line, library.good_bands, refl_nodata) + + img_dat, unc_dat, good_data = load_line( + reflectance_file, reflectance_uncertainty_file, line, library.good_bands, + refl_nodata + ) if isnothing(img_dat) return line, nothing, good_data, nothing, nothing end mixture_results = fill(-9999.0, sum(good_data), size(library.class_valid_keys)[1] + 1) complete_fractions = zeros(size(img_dat)[1], size(library.spectra)[1] + 1) - - + + if n_mc > 1 - mixture_results_std = fill(-9999.0, sum(good_data), size(library.class_valid_keys)[1] + 1) + mixture_results_std = + fill(-9999.0, sum(good_data), size(library.class_valid_keys)[1] + 1) complete_fractions_std = zeros(size(img_dat)[1], size(library.spectra)[1] + 1) else mixture_results_std = nothing @@ -311,49 +505,78 @@ function unmix_line(line::Int64, reflectance_file::String, mode::String, refl_no # Solve for each pixel for _i in 1:size(img_dat)[1] # Pixel loop - lid = img_dat[_i:_i,:] - if isnothing(unc_dat) - lud = nothing - else - lud = unc_dat[_i:_i,:] + lid = img_dat[_i:_i, :] + if isnothing(unc_dat) + lud = nothing + else + lud = unc_dat[_i:_i, :] end - loc_mixture_res, loc_mixture_var, loc_cf_res, loc_cf_var = unmix_pixel(library, lid, - lud, class_idx, options, mode, n_mc, num_endmembers, normalization, optimization, - max_combinations, combination_type) + loc_mixture_res, loc_mixture_var, loc_cf_res, loc_cf_var = + unmix_pixel(library, lid, + lud, class_idx, options, mode, n_mc, num_endmembers, normalization, + optimization, max_combinations, combination_type + ) - complete_fractions[_i,:] = loc_cf_res - mixture_results[_i,:] = loc_mixture_res + complete_fractions[_i, :] = loc_cf_res + mixture_results[_i, :] = loc_mixture_res if n_mc > 1 - complete_fractions_std[_i,:] = loc_cf_var - mixture_results_std[_i,:] = loc_mixture_var + complete_fractions_std[_i, :] = loc_cf_var + mixture_results_std[_i, :] = loc_mixture_var end end - elapsed_time = time() - start_time - if line != 1 - @info string("Line " , line , " run in " , round(elapsed_time, digits=4) , " seconds") + elapsed_time = time() - start_time + if line != 1 + @info string("Line ", line, " run in ", round(elapsed_time, digits=4), " seconds") end return line, mixture_results, good_data, mixture_results_std, complete_fractions - end -function unmix_and_write_line(line::Int64, reflectance_file::String, mode::String, refl_nodata::Float64, - refl_scale::Float64, normalization::String, library::SpectralLibrary, output_files::Vector{String}, - write_complete_fractions::Bool, - reflectance_uncertainty_file::String = "", n_mc::Int64 = 1, - combination_type::String = "all", num_endmembers::Vector{Int64} = [2,3], - max_combinations::Int64 = -1, optimization="bvls") - - line_results = unmix_line(line, reflectance_file, mode, refl_nodata, refl_scale, normalization, library, +""" + unmix_and_write_line(line::Int64, reflectance_file::String, mode::String, + refl_nodata::Float64, refl_scale::Float64, + normalization::String, library::SpectralLibrary, + output_files::Vector{String}, write_complete_fractions::Bool, + reflectance_uncertainty_file::String="", n_mc::Int64=1, + combination_type::String="all", + num_endmembers::Vector{Int64}=[2, 3], + max_combinations::Int64=-1, optimization="bvls") + +Unmix the specified `line` of `reflectance_file` and write the results to `output_files`. + +- Calls [`unmix_line`(@ref)] and [`write_line_results`](@ref), see for arguments. +""" +function unmix_and_write_line(line::Int64, reflectance_file::String, mode::String, + refl_nodata::Float64, refl_scale::Float64, normalization::String, + library::SpectralLibrary, output_files::Vector{String}, write_complete_fractions::Bool, + reflectance_uncertainty_file::String="", n_mc::Int64=1, combination_type::String="all", + num_endmembers::Vector{Int64}=[2, 3], max_combinations::Int64=-1, optimization="bvls") + + line_results = unmix_line( + line, reflectance_file, mode, refl_nodata, refl_scale, normalization, library, reflectance_uncertainty_file, n_mc, combination_type, num_endmembers, - max_combinations, optimization) + max_combinations, optimization + ) - write_line_results(output_files, line_results, n_mc, write_complete_fractions) - + write_line_results(output_files, line_results, n_mc, write_complete_fractions) end +""" + class_assign_fractions(complete_fractions, library::SpectralLibrary) + +Aggregate endmember fractions into per-unique-class fractions. + +# Arguments +- `complete_fractions`: A matrix or vector containing fractions for each pixel, +where each column represents an endmember. +- `library::SpectralLibrary`: The spectral library of the endemembers used in unmixing. + +# Returns +- A matrix or vector of class fractions, where each row corresponds to a pixel and each +column corresponds to a unique class from the library. +""" function class_assign_fractions(complete_fractions, library::SpectralLibrary) # Aggregate results from per-library to per-unique-class if length(size(complete_fractions)) == 1 @@ -365,30 +588,53 @@ function class_assign_fractions(complete_fractions, library::SpectralLibrary) for _i in 1:size(cf)[1] for (_class, cl) in enumerate(library.class_valid_keys) - mixture_results[_i, _class] = sum(complete_fractions[_i,1:end-1][cl .== library.classes]) + mixture_results[_i, _class] = + sum(complete_fractions[_i, 1:end-1][cl.==library.classes]) end end - + if length(size(complete_fractions)) == 1 - return mixture_results[1,:] + return mixture_results[1, :] else return mixture_results end - end - -function simulate_pixel(library::SpectralLibrary, max_components::Int64, combination_type::String, seed::Int64) +""" + simulate_pixel(library::SpectralLibrary, max_components::Int64, + combination_type::String, seed::Int64) + +Simulate reflectance data for a pixel from a random normalized distribution of +endmembers in `library`. + +# Arguments +- `library::SpectralLibrary`: Spectral library containing endmembers and class +information. +- `max_components::Int64`: Maximum number of endmembers to use in the simulation. +- `combination_type::String`: The type of combinations to consider when selecting +endmembers (e.g., "all" or "class-even"). +- `seed::Int64`: Seed for the random number generator. + +# Returns +- A tuple containing: + - `simulated_rfl`: Simulated reflectance spectrum of the pixel. + - `output_distribution`: A vector of endmember contributions for the pixel. + - `output_distribution_classes`: A vector of contributions of each unique class. +""" +function simulate_pixel(library::SpectralLibrary, max_components::Int64, + combination_type::String, seed::Int64) Random.seed!(seed) - + output_mixture = zeros(size(library.spectra)[2]) output_mixture[:] .= NaN class_idx = prepare_combinations(library, combination_type) - perm = get_sma_permutation(class_idx, [max_components], combination_type, size(library.spectra)[1]) + perm = get_sma_permutation( + class_idx, [max_components], combination_type, size(library.spectra)[1] + ) - G = library.spectra[perm,:] + G = library.spectra[perm, :] distribution = rand(max_components) distribution = distribution ./ sum(distribution) @@ -398,13 +644,11 @@ function simulate_pixel(library::SpectralLibrary, max_components::Int64, combina output_distribution_classes = zeros(size(library.class_valid_keys)) for (_class, cl) in enumerate(library.class_valid_keys) - output_distribution_classes[_class] = sum(output_distribution[cl .== library.classes]) + output_distribution_classes[_class] = sum(output_distribution[cl.==library.classes]) end simulated_rfl = G' * distribution return simulated_rfl, output_distribution, output_distribution_classes - end - end