Fix sampling of FG and BG signals in RawDataViewer

src/Viewer/RawDataViewer.jl

@@ -244,8 +244,21 @@ end
       # TODO: Ensure that the measurements fit together (num samples / patches)
       # otherwise -> error
 
-      numFGFrames = minimum(acqNumFGFrames.(fs))
-      numBGFrames = minimum(acqNumBGFrames.(fs))    
+      numFGFrames = unique(acqNumFGFrames.(fs))
+      if length(numFGFrames) > 1
+        numFGFrames = minimum(numFGFrames)
+        @warn "Different number of foreground frames in data, limit frames to $numFGFrames"
+      else 
+        numFGFrames = first(numFGFrames)
+      end
+
+      numBGFrames = unique(acqNumBGFrames.(fs))
+      if length(numBGFrames) > 1
+        numBGFrames = minimum(numBGFrames)
+        @warn "Different number of background frames in data, limit frames to $numBGFrames"
+      else 
+        numBGFrames = first(numBGFrames)
+      end
 
       dataFGVec = Any[]
       dataBGVec = Any[]
@@ -271,7 +284,7 @@ end
                     tfCorrection=get_gtk_property(m["cbCorrTF"], :active, Bool))
         push!(dataFGVec, data)
 
-        if acqNumBGFrames(f) > 0
+        if numBGFrames > 0
           dataBG = getMeasurements(f, false, frames=measBGFrameIdx(f),
                 bgCorrection=false, spectralLeakageCorrection = get_gtk_property(m["cbSLCorr"], :active, Bool),
                 tfCorrection=get_gtk_property(m["cbCorrTF"], :active, Bool))
@@ -381,42 +394,30 @@ end
     numFreq = floor(Int, size(data,1) ./ 2 .+ 1)
 
     maxPoints = 5000
-    sp = length(minTP:maxTP) > maxPoints ? round(Int,length(minTP:maxTP) / maxPoints)  : 1
-
-    steps = minTP:sp:maxTP
-    dataCompressed = zeros(length(steps), size(data,2))
-    if sp > 1
-      for j=1:size(data,2)
-        for l=1:length(steps)
-          st = steps[l]
-          en = min(st+sp,steps[end])
-          med_ = median(data[st:en,j])
-          max_ = maximum(data[st:en,j])
-          min_ = minimum(data[st:en,j])
-          dataCompressed[l,j] = rand(Bool) ? max_ : min_ #abs(max_ - med_) > abs(med_ - min_) ? max_ : min_
-        end
-      end
-    else
-      dataCompressed = data[steps,:]
-    end
+    # Reduce each channel to either min or max value if intervals of more than maxpoint points
+    reduceOp = (vals -> map(slice -> rand(Bool) ? maximum(slice) : minimum(slice), eachslice(vals, dims=2)))
+    timeCompressed, dataCompressed = prepareTimeDataForPlotting(data, timePoints, indexInterval=(minTP,maxTP),
+                        maxpoints=maxPoints, reduceOp = reduceOp)
 
-    fTD, axTD, lTD1 = CairoMakie.lines(timePoints[steps], dataCompressed[:,1], 
+    fTD, axTD, lTD1 = CairoMakie.lines(timeCompressed, dataCompressed[:,1], 
                             figure = (; figure_padding=4, resolution = (1000, 800), fontsize = 11),
                             axis = (; title = "Time Domain"),
                             color = CairoMakie.RGBf(colors[1]...),
                             label = labels_[1])
     for j=2:size(data,2)
-      CairoMakie.lines!(axTD, timePoints[steps],dataCompressed[:,j], 
+      CairoMakie.lines!(axTD, timeCompressed, dataCompressed[:,j], 
                         color = CairoMakie.RGBf(colors[j]...),  #linewidth=3)
                         label = labels_[j])
     end
     if length(m.dataBG) > 0 && get_gtk_property(m["cbShowBG"], :active, Bool)
-      CairoMakie.lines!(axTD, timePoints[minTP:sp:maxTP],dataBG[minTP:sp:maxTP,1], color=:black,
+      _, dataBGCompressed = prepareTimeDataForPlotting(dataBG, timePoints, indexInterval=(minTP,maxTP),
+                        maxpoints=maxPoints, reduceOp = reduceOp)
+      CairoMakie.lines!(axTD, timeCompressed, dataBGCompressed[:, 1], color=:black,
              label="BG", linestyle = :dash)
     end    
     CairoMakie.autolimits!(axTD)
-    if timePoints[steps[end]] > timePoints[steps[1]]
-      CairoMakie.xlims!(axTD, timePoints[steps[1]], timePoints[steps[end]])
+    if timeCompressed[end] > timeCompressed[1]
+      CairoMakie.xlims!(axTD, timeCompressed[1], timeCompressed[end])
     end
     if !autoRangingTD && maxValTD > minValTD 
       CairoMakie.ylims!(axTD, minValTD, maxValTD)
@@ -436,33 +437,26 @@ end
       spFr = length(minFr:maxFr) > maxPoints ? round(Int,length(minFr:maxFr) / maxPoints)  : 1
 
       stepsFr = minFr:spFr:maxFr
-      freqDataCompressed = zeros(length(stepsFr), size(freqdata,2))
-      if spFr > 1
-        for j=1:size(freqdata,2)
-          for l=1:length(stepsFr)
-            st = stepsFr[l]
-            en = min(st+spFr,stepsFr[end])
-            freqDataCompressed[l,j] = maximum(freqdata[st:en,j])
-          end
-        end
-      else
-        freqDataCompressed = freqdata[stepsFr,:]
-      end
+      freqsCompressed, freqDataCompressed = prepareTimeDataForPlotting(freqdata, freq, indexInterval=(minFr,maxFr),
+                        maxpoints=maxPoints, reduceOp = vals -> map(maximum, eachslice(vals, dims=2)))
+
 
-      fFD, axFD, lFD1 = CairoMakie.lines(freq[stepsFr],freqDataCompressed[:,1], 
+      fFD, axFD, lFD1 = CairoMakie.lines(freqsCompressed,freqDataCompressed[:,1], 
                               figure = (; figure_padding=4, resolution = (1000, 800), fontsize = 11),
                               axis = (; title = "Frequency Domain", yscale=log10),
                               color = CairoMakie.RGBf(colors[1]...),
                               label=labels_[1])
       for j=2:size(data,2)
-        CairoMakie.lines!(axFD, freq[stepsFr], freqDataCompressed[:,j], 
+        CairoMakie.lines!(axFD, freqsCompressed, freqDataCompressed[:,j], 
                      color = CairoMakie.RGBf(colors[j]...), label=labels_[j])
       end
       if length(m.dataBG) > 0 && get_gtk_property(m["cbShowBG"], :active, Bool)
-        CairoMakie.lines!(axTD, timePoints[minTP:sp:maxTP],dataBG[minTP:sp:maxTP,1], color=:black,
-              label="BG", linestyle = :dash)
+        #CairoMakie.lines!(axTD, timePoints[minTP:sp:maxTP],dataBG[minTP:sp:maxTP,1], color=:black,
+        #      label="BG", linestyle = :dash) # isn't this duplicate?
         if showFD
-          CairoMakie.lines!(axFD, freq[minFr:spFr:maxFr],abs.(rfft(dataBG,1)[minFr:spFr:maxFr,1]) / size(dataBG,1),
+          _, freqDataBGCompressed = prepareTimeDataForPlotting(abs.(rfft(dataBG,1)), freq, indexInterval=(minFr,maxFr),
+                        maxpoints=maxPoints, reduceOp = vals -> map(maximum, eachslice(vals, dims=2)))
+          CairoMakie.lines!(axFD, freqsCompressed, freqDataBGCompressed[:, 1] / size(dataBG,1),
                  color=:black, label="BG", linestyle = :dash)
         end
       end    

src/Viewer/Viewer.jl

@@ -1,3 +1,63 @@
+"""
+    prepareTimeDataForPlotting(data::AbstractMatrix, times = 1:length(data); timeInterval = nothing, maxpoints = length(data), reduceOp = nothing)
+
+Reduces time-series `data` to at most `maxpoints` by either resampling or reducing the data with `reduceOp`. The time-series is defined by `times` and `timeInterval`.
+If `timeInterval` is not given, the whole time-series is used. If `reduceOp` is not given, the data is resampled. 
+The function returns the reduced time-series `t` and the reduced data `d`.
+"""
+function prepareTimeDataForPlotting(data::AbstractMatrix, times = 1:length(data); indexInterval = nothing, timeInterval = nothing, maxpoints = length(data), reduceOp = nothing)
+  idx1, idx2 = timeSlice(times, indexInterval, timeInterval)
+  if length(idx1:idx2) > maxpoints
+    t, d = reducePoints(reduceOp, data, times, idx1, idx2, maxpoints)
+  else
+    t = times[idx1:idx2]
+    d = data[idx1:idx2, :]
+  end
+  # return [[Point2f(t[l], d[l, k]) for l in 1:size(d,1)] for k in 1:size(d, 2)]
+  return t, d
+end
+# Main variant is defined for matrices, so we reshape into Nx1 matrix
+function prepareTimeDataForPlotting(data::AbstractVector, args...; kwargs...)
+  # return first(prepareTimeDataForPlotting(reshape(data, :, 1), args...; kwargs...)) 
+  t, d = prepareTimeDataForPlotting(reshape(data, :, 1), args...; kwargs...)
+  return t, vec(d)
+end
+
+
+timeSlice(times, indexInterval, timeInterval) = error("Cannot specify both indexInterval and timeInterval")
+timeSlice(times, indexInterval, ::Nothing) = first(indexInterval), last(indexInterval)
+function timeSlice(times, ::Nothing, timeInterval)
+  t1, t2 = timeInterval
+  idx1 = findfirst(t -> (t>t1), times)
+  idx2 = findlast(t -> (t<t2), times)
+  idx1 = (idx1 == nothing) ? 1 : idx1 
+  idx2 = (idx2 == nothing) ? length(times) : idx2
+  idx1 = min(max(1,idx1),length(times))
+  idx2 = min(max(1,idx2),length(times))
+  return idx1, idx2
+end
+timeSlice(times, ::Nothing) = 1, length(times)
+
+function reducePoints(reduceOp, data, times, idx1, idx2, maxpoints)
+  stepsize = round(Int, (length(idx1:idx2)) / maxpoints, RoundUp)
+  steps = idx1:stepsize:idx2
+  # Sample time
+  t = range(times[idx1], times[idx2], length=length(steps))
+  # Sample data with reduceOp
+  d = zeros(eltype(data), length(t), size(data, 2))
+  # Zip generates start and end indices for each step
+  for (i, (st, en)) in enumerate(zip(steps, map(i -> min(i + stepsize -1, idx2), steps)))
+    d[i, :] = reduceOp(@view data[st:en, :])
+    #@info "Reducing data from $st to $en with $(d[i, :])"
+  end
+  return t, d
+end
+function reducePoints(::Nothing, data, times, idx1, idx2, maxpoints)
+  d = DSP.resample(data[idx1:idx2], maxpoints / length(idx1:idx2))
+  t = range(times[idx1], times[idx2], length=length(d))
+  return t, d
+end
+
 include("BaseViewer.jl")
 include("SimpleDataViewer.jl")
 include("DataViewer/DataViewer.jl")