I was very confused, because I would make a suggested change, such as using FLoops or rewriting it for matrices, and I was getting inconsistent results. A run consisted of processing 6 images, and sometimes the processing times for the images were all over the place.

Someone suggested that I had a global variable problem, but I didn’t have any global variables in my file. Nonetheless, while using JET to find type inconsistencies, it found a global variable in another file (that I include()-ed in my main file) that had a const global variable that was being redefined. I cleaned up the types, and used @isdefined to prevent the constants from being redefined, and the result is lightening-fast code that went from 38 seconds down to 0.58 seconds to process 6 images! It’s probably faster than that, as I still have some diagnostic code left in that reports timings.

Thanks everyone! The cleaned-up coded is below:

```
function processHyperSpectral( images::Vector{Matrix{Gray{N0f16}}}, CMF::Vector{XYZ{Float64}}, Illuminant::Vector{Float64})
numWL = length(CMF)
cmfArray = zeros( (3,numWL))
for i in eachindex(CMF)
cmfArray[1, i] = CMF[i].x
cmfArray[2, i] = CMF[i].y
cmfArray[3,i] = CMF[i].z
end
dimensions = size(images[1])
horiz = dimensions[1]
vert = dimensions[2]
XYZlayers =Vector{Matrix{Float64}}()
for j in 1:3
push!( XYZlayers, fill(0.0, horiz, vert) )
end
#----------------
# Precalculate Illuminant x CMF
IllumiCMF = zeros( (3,numWL))
for i in eachindex(CMF)
for j in 1:3
IllumiCMF[j, i] = cmfArray[j, i] .* Illuminant[i]
end
end
#------------------
# 5) Multiply hyperspectral by cmf
for i in eachindex( images)
for j in 1:3
XYZlayers[j] .+= images[i] .* IllumiCMF[j,i] # CHANGE 1B
end
end
#---------------
# We'll norm these by XYZ, which is the average area under the 3 curves
XYZ = sum(cmfArray)
XYZ /= 3
for j in 1:3
XYZlayers[j] ./= XYZ # CHANGE 1c
end
#-----------------------------
# Next, convert XYZ to RGB
#colorLayers = fill(0.5, 3, horiz, vert)
colorLayers = Vector{Matrix{Float64}}()
for j in 1:3
push!( colorLayers, fill(0.0, horiz, vert) )
end
for xy in CartesianIndices(colorLayers[1])
A = XYZlayers[1][ xy]
B = XYZlayers[2][ xy]
C = XYZlayers[3][ xy]
colorLayers[1][ xy] = muladd( 3.2404542,
A,
muladd(-1.5371385,
B,
-0.4985314 * C)
)
colorLayers[2][ xy] = muladd( -0.969266,
A,
muladd(1.8760108,
XYZlayers[2][ xy],
0.041556 * C)
)
colorLayers[3][ xy] = muladd( 0.0556434,
A,
muladd(-0.2040259,
B,
1.0572252 * C)
)
end
#------------------
# 6) perform gamma correction and composit images into the final image
for j in 1:3
clamp!(colorLayers[j], 0, 1)
end
newImage = fill(RGB(0.0, .0, 0.0), horiz, vert)
for xy in CartesianIndices(colorLayers[1])
redTemp, greenTemp, blueTemp = sRGBGamma(colorLayers[1][ xy], colorLayers[2][ xy], colorLayers[3][ xy])
newImage[xy] = RGB( redTemp, greenTemp, blueTemp ) # collect changes tuple to an array
end
return newImage
end
```