flow-plots 0.3.1

Package for drawing and interacting with plots in flow cytometry data
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
use crate::contour::ContourData;
use crate::PlotBytes;
use crate::create_axis_specs;
use crate::density_calc::RawPixelData;
use crate::options::DensityPlotOptions;
use crate::render::{ProgressInfo, RenderConfig};
use flow_fcs::{TransformType, Transformable};
use plotters::prelude::*;

/// Format a value using the transform type
///
/// This replicates the Formattable::format logic since the trait is not exported.
fn format_transform_value(transform: &TransformType, value: &f32) -> String {
    match transform {
        TransformType::Linear => format!("{:.1e}", value),
        TransformType::Arcsinh { cofactor: _ } => {
            // Convert from transformed space back to original space
            let original_value = transform.inverse_transform(value);
            // Make nice rounded labels in original space
            format!("{:.1e}", original_value)
        }
        TransformType::Biexponential { .. } => {
            // Convert from transformed space back to original space using inverse transform
            let original_value = transform.inverse_transform(value);
            // Make nice rounded labels in original space
            format!("{:.1e}", original_value)
        }
    }
}
use anyhow::Result;
use image::RgbImage;
use plotters::{
    backend::BitMapBackend, chart::ChartBuilder, prelude::IntoDrawingArea, style::WHITE,
};

/// Render pixels to a JPEG image using the Plotters backend
///
/// This function handles the complete rendering pipeline:
/// 1. Sets up Plotters chart with axes and mesh
/// 2. Writes pixels directly to the buffer
/// 3. Encodes to JPEG format
///
/// Progress reporting is handled via the RenderConfig if provided.
pub fn render_pixels(
    pixels: Vec<RawPixelData>,
    options: &DensityPlotOptions,
    render_config: &mut RenderConfig,
) -> Result<PlotBytes> {
    use crate::options::PlotOptions;

    let base = options.base();
    let width = base.width;
    let height = base.height;
    let margin = base.margin;
    let x_label_area_size = base.x_label_area_size;
    let y_label_area_size = base.y_label_area_size;

    let setup_start = std::time::Instant::now();
    // Use RGB buffer (3 bytes per pixel) since we'll encode to JPEG which doesn't support alpha
    let mut pixel_buffer = vec![255; (width * height * 3) as usize];

    let (plot_x_range, plot_y_range, x_spec, y_spec) = {
        let backend = BitMapBackend::with_buffer(&mut pixel_buffer, (width, height));
        let root = backend.into_drawing_area();
        root.fill(&WHITE)
            .map_err(|e| anyhow::anyhow!("failed to fill plot background: {e}"))?;

        // Create appropriate ranges based on transform types
        let (x_spec, y_spec) = create_axis_specs(
            &options.x_axis.range,
            &options.y_axis.range,
            &options.x_axis.transform,
            &options.y_axis.transform,
        )?;

        let mut chart = ChartBuilder::on(&root)
            .margin(margin)
            .x_label_area_size(x_label_area_size)
            .y_label_area_size(y_label_area_size)
            .build_cartesian_2d(x_spec.start..x_spec.end, y_spec.start..y_spec.end)?;

        // Clone transforms to avoid lifetime issues with closures
        let x_transform_clone = options.x_axis.transform.clone();
        let y_transform_clone = options.y_axis.transform.clone();

        // Create owned closures for formatters
        let x_formatter =
            move |x: &f32| -> String { format_transform_value(&x_transform_clone, x) };
        let y_formatter =
            move |y: &f32| -> String { format_transform_value(&y_transform_clone, y) };

        let mut mesh = chart.configure_mesh();
        mesh.x_max_light_lines(4)
            .y_max_light_lines(4)
            .x_labels(10)
            .y_labels(10)
            .x_label_formatter(&x_formatter)
            .y_label_formatter(&y_formatter);

        // Add axis labels if provided
        if let Some(ref x_label) = options.x_axis.label {
            mesh.x_desc(x_label);
        }
        if let Some(ref y_label) = options.y_axis.label {
            mesh.y_desc(y_label);
        }

        let mesh_start = std::time::Instant::now();
        mesh.draw()
            .map_err(|e| anyhow::anyhow!("failed to draw plot mesh: {e}"))?;
        eprintln!("    ├─ Mesh drawing: {:?}", mesh_start.elapsed());

        // Get the plotting area bounds (we'll use these after Plotters releases the buffer)
        let plotting_area = chart.plotting_area();
        let (plot_x_range, plot_y_range) = plotting_area.get_pixel_range();

        root.present()
            .map_err(|e| anyhow::anyhow!("failed to present plotters buffer: {e}"))?;

        (plot_x_range, plot_y_range, x_spec, y_spec)
    }; // End Plotters scope - pixel_buffer is now released and we can write to it

    // DIRECT PIXEL BUFFER WRITING - 10-50x faster than Plotters series rendering
    // Now that Plotters has released pixel_buffer, we can write directly
    let series_start = std::time::Instant::now();

    let plot_x_start = plot_x_range.start as f32;
    let plot_y_start = plot_y_range.start as f32;
    let plot_width = (plot_x_range.end - plot_x_range.start) as f32;
    let plot_height = (plot_y_range.end - plot_y_range.start) as f32;

    // Calculate scale factors from data coordinates to screen pixels
    let data_width = x_spec.end - x_spec.start;
    let data_height = y_spec.end - y_spec.start;

    // Stream pixel chunks during rendering using configurable chunk size
    let mut pixel_count = 0;
    let total_pixels = pixels.len();
    let chunk_size = 1000; // Default chunk size for progress reporting

    // Write each pixel directly to the buffer
    for pixel in &pixels {
        let data_x = pixel.x;
        let data_y = pixel.y;

        // Transform data coordinates to screen pixel coordinates
        let rel_x = (data_x - x_spec.start) / data_width;
        let rel_y = (y_spec.end - data_y) / data_height; // Flip Y (screen coords go down)

        let screen_x = (plot_x_start + rel_x * plot_width) as i32;
        let screen_y = (plot_y_start + rel_y * plot_height) as i32;

        // Bounds check
        if screen_x >= plot_x_range.start
            && screen_x < plot_x_range.end
            && screen_y >= plot_y_range.start
            && screen_y < plot_y_range.end
        {
            let px = screen_x as u32;
            let py = screen_y as u32;

            // Write to pixel buffer (RGB format - 3 bytes per pixel)
            let idx = ((py * width + px) * 3) as usize;

            if idx + 2 < pixel_buffer.len() {
                pixel_buffer[idx] = pixel.r;
                pixel_buffer[idx + 1] = pixel.g;
                pixel_buffer[idx + 2] = pixel.b;
            }
        }

        pixel_count += 1;

        // Emit progress every chunk_size pixels
        if pixel_count % chunk_size == 0 || pixel_count == total_pixels {
            let percent = (pixel_count as f32 / total_pixels as f32) * 100.0;

            // Create a small sample of pixels for this chunk (for visualization)
            let chunk_start = (pixel_count - chunk_size.min(pixel_count)).max(0);
            let chunk_end = pixel_count;
            let chunk_pixels: Vec<RawPixelData> = pixels
                .iter()
                .skip(chunk_start)
                .take(chunk_end - chunk_start)
                .map(|p| RawPixelData {
                    x: p.x,
                    y: p.y,
                    r: p.r,
                    g: p.g,
                    b: p.b,
                })
                .collect();

            render_config.report_progress(ProgressInfo {
                pixels: chunk_pixels,
                percent,
            });
        }
    }

    eprintln!(
        "    ├─ Direct pixel writing: {:?} ({} pixels)",
        series_start.elapsed(),
        pixels.len()
    );
    eprintln!("    ├─ Total plotting: {:?}", setup_start.elapsed());

    let img_start = std::time::Instant::now();
    let img: RgbImage = image::ImageBuffer::from_vec(width, height, pixel_buffer)
        .ok_or_else(|| anyhow::anyhow!("plot image buffer had unexpected size"))?;
    eprintln!("    ├─ Image buffer conversion: {:?}", img_start.elapsed());

    let encode_start = std::time::Instant::now();

    // Pre-allocate Vec with estimated JPEG size
    // RGB buffer is (width * height * 3) bytes
    // JPEG at quality 85 typically compresses to ~10-15% of raw size for density plots
    let raw_size = (width * height * 3) as usize;
    let estimated_jpeg_size = raw_size / 8; // Conservative estimate (~12.5% of raw)
    let mut encoded_data = Vec::with_capacity(estimated_jpeg_size);

    // JPEG encoding is faster and produces smaller files for density plots
    // Quality 85 provides good visual quality with ~2x smaller file size vs PNG
    let mut encoder = image::codecs::jpeg::JpegEncoder::new_with_quality(&mut encoded_data, 85);
    encoder
        .encode(img.as_raw(), width, height, image::ExtendedColorType::Rgb8)
        .map_err(|e| anyhow::anyhow!("failed to JPEG encode plot: {e}"))?;
    eprintln!("    └─ JPEG encoding: {:?}", encode_start.elapsed());

    // Return the JPEG-encoded bytes directly
    Ok(encoded_data)
}

/// Render contour plot to JPEG using Plotters LineSeries
///
/// Draws contour lines from KDE density estimation plus optional outlier scatter points.
pub fn render_contour(
    contour_data: ContourData,
    options: &DensityPlotOptions,
    _render_config: &mut RenderConfig,
) -> Result<PlotBytes> {
    use crate::options::PlotOptions;

    let base = options.base();
    let width = base.width;
    let height = base.height;
    let margin = base.margin;
    let x_label_area_size = base.x_label_area_size;
    let y_label_area_size = base.y_label_area_size;

    let (x_spec, y_spec) = create_axis_specs(
        &options.x_axis.range,
        &options.y_axis.range,
        &options.x_axis.transform,
        &options.y_axis.transform,
    )?;

    let mut pixel_buffer = vec![255; (width * height * 3) as usize];

    {
        let backend = BitMapBackend::with_buffer(&mut pixel_buffer, (width, height));
        let root = backend.into_drawing_area();
        root.fill(&WHITE)
            .map_err(|e| anyhow::anyhow!("failed to fill plot background: {e}"))?;

        let x_transform_clone = options.x_axis.transform.clone();
        let y_transform_clone = options.y_axis.transform.clone();
        let x_formatter = move |x: &f64| -> String {
            format_transform_value(&x_transform_clone, &(*x as f32))
        };
        let y_formatter = move |y: &f64| -> String {
            format_transform_value(&y_transform_clone, &(*y as f32))
        };

        let mut chart = ChartBuilder::on(&root)
            .margin(margin)
            .x_label_area_size(x_label_area_size)
            .y_label_area_size(y_label_area_size)
            .build_cartesian_2d(
                x_spec.start as f64..x_spec.end as f64,
                y_spec.start as f64..y_spec.end as f64,
            )?;

        let mut mesh = chart.configure_mesh();
        mesh.x_max_light_lines(4)
            .y_max_light_lines(4)
            .x_labels(10)
            .y_labels(10)
            .x_label_formatter(&x_formatter)
            .y_label_formatter(&y_formatter);
        if let Some(ref x_label) = options.x_axis.label {
            mesh.x_desc(x_label);
        }
        if let Some(ref y_label) = options.y_axis.label {
            mesh.y_desc(y_label);
        }
        mesh.draw()
            .map_err(|e| anyhow::anyhow!("failed to draw plot mesh: {e}"))?;

        let stroke_width = options.contour_line_thickness.max(0.5).min(5.0) as u32;
        let contour_color = RGBColor(60, 60, 60);

        // Draw contour lines
        for path in &contour_data.contours {
            if path.len() < 2 {
                continue;
            }
            let points: Vec<(f64, f64)> = path.iter().copied().collect();
            chart
                .draw_series(LineSeries::new(
                    points,
                    contour_color.stroke_width(stroke_width),
                ))
                .map_err(|e| anyhow::anyhow!("failed to draw contour: {e}"))?;
        }

        // Draw outlier points if present
        if !contour_data.outliers.is_empty() {
            let outlier_color = RGBColor(150, 150, 150);
            chart
                .draw_series(
                    contour_data
                        .outliers
                        .iter()
                        .map(|&(x, y)| Circle::new((x, y), 2, outlier_color.filled())),
                )
                .map_err(|e| anyhow::anyhow!("failed to draw outliers: {e}"))?;
        }

        root.present()
            .map_err(|e| anyhow::anyhow!("failed to present plotters buffer: {e}"))?;
    }

    let img: RgbImage =
        image::ImageBuffer::from_vec(width, height, pixel_buffer)
            .ok_or_else(|| anyhow::anyhow!("plot image buffer had unexpected size"))?;

    let mut encoded_data = Vec::new();
    let mut encoder =
        image::codecs::jpeg::JpegEncoder::new_with_quality(&mut encoded_data, 85);
    encoder
        .encode(img.as_raw(), width, height, image::ExtendedColorType::Rgb8)
        .map_err(|e| anyhow::anyhow!("failed to JPEG encode plot: {e}"))?;

    Ok(encoded_data)
}

/// Render spectral signature plot to JPEG
///
/// Creates a line plot showing normalized spectral signatures (0.0 to 1.0) across detector channels.
pub fn render_spectral_signature(
    data: (Vec<(usize, f64)>, Vec<String>),
    options: &crate::options::spectral::SpectralSignaturePlotOptions,
    _render_config: &mut RenderConfig,
) -> Result<PlotBytes> {
    use crate::options::PlotOptions;
    use plotters::prelude::*;

    let (spectrum_data, channel_names) = data;
    let base = options.base();
    let width = base.width;
    let height = base.height;
    let margin = base.margin;
    let x_label_area_size = base.x_label_area_size;
    let y_label_area_size = base.y_label_area_size;

    // Create RGB buffer
    let mut pixel_buffer = vec![255; (width * height * 3) as usize];

    // Determine x and y ranges (use f32 to match plotters expectations)
    let x_min = 0.0f32;
    let x_max = spectrum_data
        .iter()
        .map(|(idx, _)| *idx as f32)
        .fold(0.0f32, f32::max)
        .max(1.0);
    let y_min = 0.0f32;
    let y_max = 1.0f32;

    // Clone channel_names for the closure
    let channel_names_clone = channel_names.clone();

    {
        let backend = BitMapBackend::with_buffer(&mut pixel_buffer, (width, height));
        let root = backend.into_drawing_area();
        root.fill(&WHITE)
            .map_err(|e| anyhow::anyhow!("failed to fill plot background: {e}"))?;

        let mut chart = ChartBuilder::on(&root)
            .margin(margin)
            .x_label_area_size(x_label_area_size)
            .y_label_area_size(y_label_area_size)
            .build_cartesian_2d(x_min..x_max, y_min..y_max)
            .map_err(|e| anyhow::anyhow!("failed to build chart: {e}"))?;

        // Create formatter for x-axis labels if channel names are provided
        let x_formatter: Option<Box<dyn Fn(&f32) -> String>> = if !channel_names_clone.is_empty()
            && channel_names_clone.len() == spectrum_data.len()
        {
            let channel_names_for_formatter = channel_names_clone.clone();
            Some(Box::new(move |x: &f32| -> String {
                // Find closest channel index
                let idx = x.round() as usize;
                if idx < channel_names_for_formatter.len() {
                    channel_names_for_formatter[idx].clone()
                } else {
                    format!("{:.0}", x)
                }
            }))
        } else {
            None
        };

        // Configure mesh
        let mut mesh = chart.configure_mesh();
        if options.show_grid {
            mesh.x_max_light_lines(4).y_max_light_lines(4);
        }

        // Set axis labels
        if let Some(ref x_axis) = options.x_axis {
            if let Some(ref label) = x_axis.label {
                mesh.x_desc(label);
            }
        } else {
            mesh.x_desc("Channel");
        }

        if let Some(ref y_axis) = options.y_axis {
            if let Some(ref label) = y_axis.label {
                mesh.y_desc(label);
            }
        } else {
            mesh.y_desc("Normalized Intensity");
        }

        // Apply x-axis formatter if provided
        if let Some(ref formatter) = x_formatter {
            mesh.x_label_formatter(formatter);
        }

        // Set number of x labels to match number of channels for better readability
        let x_label_count = if !channel_names_clone.is_empty() {
            channel_names_clone.len().min(20) // Show all channels if <= 20, otherwise show 20
        } else {
            10
        };

        mesh.x_labels(x_label_count)
            .y_labels(10)
            .draw()
            .map_err(|e| anyhow::anyhow!("failed to draw mesh: {e}"))?;

        // Draw the spectral signature line
        if !spectrum_data.is_empty() {
            // Parse hex color (e.g., "#1f77b4" or "1f77b4")
            let line_color = if options.line_color.starts_with('#') {
                let hex = &options.line_color[1..];
                if hex.len() == 6 {
                    let r = u8::from_str_radix(&hex[0..2], 16).unwrap_or(31);
                    let g = u8::from_str_radix(&hex[2..4], 16).unwrap_or(119);
                    let b = u8::from_str_radix(&hex[4..6], 16).unwrap_or(180);
                    RGBColor(r, g, b)
                } else {
                    RGBColor(31, 119, 180) // Default blue
                }
            } else {
                RGBColor(31, 119, 180) // Default blue
            };

            chart
                .draw_series(LineSeries::new(
                    spectrum_data
                        .iter()
                        .map(|(idx, val)| (*idx as f32, *val as f32)),
                    line_color.stroke_width(options.line_width as u32),
                ))
                .map_err(|e| anyhow::anyhow!("failed to draw line series: {e}"))?;
        }
    } // Backend is dropped here, pixel_buffer is now available

    // Convert to image and encode
    let img: RgbImage = image::ImageBuffer::from_vec(width, height, pixel_buffer)
        .ok_or_else(|| anyhow::anyhow!("plot image buffer had unexpected size"))?;

    let mut encoded_data = Vec::new();
    let mut encoder = image::codecs::jpeg::JpegEncoder::new_with_quality(&mut encoded_data, 85);
    encoder
        .encode(img.as_raw(), width, height, image::ExtendedColorType::Rgb8)
        .map_err(|e| anyhow::anyhow!("failed to JPEG encode plot: {e}"))?;

    Ok(encoded_data)
}

/// Render histogram plot to JPEG
///
/// Creates a 1D histogram (x = values, y = count/frequency) with optional fill
/// and support for overlaid multiple series.
pub fn render_histogram(
    data: crate::histogram_data::HistogramData,
    options: &crate::options::HistogramPlotOptions,
    _render_config: &mut RenderConfig,
) -> Result<PlotBytes> {
    use crate::histogram_data::{bin_values, BinnedHistogram, HistogramData, HistogramSeries};
    use crate::options::PlotOptions;
    use plotters::prelude::*;

    let base = options.base();
    let width = base.width;
    let height = base.height;
    let margin = base.margin;
    let x_label_area_size = base.x_label_area_size;
    let y_label_area_size = base.y_label_area_size;

    let x_min = *options.x_axis.range.start() as f64;
    let x_max = *options.x_axis.range.end() as f64;

    // Convert all data to binned series (list of (BinnedHistogram, gate_id))
    let series: Vec<(BinnedHistogram, u32)> = match data {
        HistogramData::RawValues(values) => {
            let binned = bin_values(
                &values,
                options.num_bins,
                *options.x_axis.range.start(),
                *options.x_axis.range.end(),
            );
            match binned {
                Some(b) => vec![(b, 0)],
                None => vec![],
            }
        }
        HistogramData::PreBinned { bin_edges, counts } => {
            let bin_centers: Vec<f64> = bin_edges
                .windows(2)
                .map(|w| (w[0] as f64 + w[1] as f64) / 2.0)
                .collect();
            let counts_f64: Vec<f64> = counts.iter().map(|&c| c as f64).collect();
            vec![(
                BinnedHistogram {
                    bin_centers,
                    counts: counts_f64,
                },
                0,
            )]
        }
        HistogramData::Overlaid(overlaid) => {
            let mut result = Vec::with_capacity(overlaid.len());
            for HistogramSeries { values, gate_id } in overlaid {
                if let Some(binned) = bin_values(
                    &values,
                    options.num_bins,
                    *options.x_axis.range.start(),
                    *options.x_axis.range.end(),
                ) {
                    result.push((binned, gate_id));
                }
            }
            result
        }
    };

    if series.is_empty() {
        // Empty plot - still render axes
        return render_empty_histogram(
            options, width, height, margin, x_label_area_size, y_label_area_size, x_min, x_max,
        );
    }

    // Optional: scale each series to its peak (max = 1.0)
    let series: Vec<(BinnedHistogram, u32)> = if options.scale_to_peak && series.len() > 1 {
        series
            .into_iter()
            .map(|(mut binned, gate_id)| {
                let max_count = binned.counts.iter().cloned().fold(0.0f64, f64::max);
                if max_count > 0.0 {
                    binned.counts.iter_mut().for_each(|c| *c /= max_count);
                }
                (binned, gate_id)
            })
            .collect()
    } else {
        series
    };

    // Compute y range
    let (y_min, y_max) = if options.baseline_separation > 0.0 && series.len() > 1 {
        // Stacked: each series gets baseline_separation offset
        let mut max_y = 0.0f64;
        let mut cumulative_offset = 0.0f64;
        for (binned, _) in &series {
            let peak = binned.counts.iter().cloned().fold(0.0f64, f64::max);
            max_y = max_y.max(cumulative_offset + peak);
            cumulative_offset += options.baseline_separation as f64;
        }
        (0.0, (max_y * 1.05).max(0.1))
    } else {
        let global_max = series
            .iter()
            .flat_map(|(b, _)| b.counts.iter())
            .cloned()
            .fold(0.0f64, f64::max);
        (0.0, (global_max * 1.05).max(0.1))
    };

    let mut pixel_buffer = vec![255; (width * height * 3) as usize];

    {
        let backend = BitMapBackend::with_buffer(&mut pixel_buffer, (width, height));
        let root = backend.into_drawing_area();
        root.fill(&WHITE)
            .map_err(|e| anyhow::anyhow!("failed to fill plot background: {e}"))?;

        let mut chart = ChartBuilder::on(&root)
            .margin(margin)
            .x_label_area_size(x_label_area_size)
            .y_label_area_size(y_label_area_size)
            .build_cartesian_2d(x_min..x_max, y_min..y_max)
            .map_err(|e| anyhow::anyhow!("failed to build histogram chart: {e}"))?;

        let mut mesh = chart.configure_mesh();
        mesh.x_max_light_lines(4).y_max_light_lines(4)
            .x_labels(10)
            .y_labels(10);

        if let Some(ref label) = options.x_axis.label {
            mesh.x_desc(label);
        } else {
            mesh.x_desc("Value");
        }
        mesh.y_desc("Count");

        mesh.draw()
            .map_err(|e| anyhow::anyhow!("failed to draw mesh: {e}"))?;

        let baseline_sep = options.baseline_separation as f64;
        let mut y_offset = 0.0f64;

        for (binned, gate_id) in &series {
            let (r, g, b) = options.gate_color(*gate_id);
            let color = RGBColor(r, g, b);
            let fill_color = RGBColor(r, g, b).mix(0.3);

            let points: Vec<(f64, f64)> = binned
                .bin_centers
                .iter()
                .zip(binned.counts.iter())
                .map(|(x, c)| (*x, y_offset + *c))
                .collect();

            if points.is_empty() {
                y_offset += baseline_sep;
                continue;
            }

            if options.histogram_filled {
                chart
                    .draw_series(AreaSeries::new(
                        points.iter().copied(),
                        y_offset,
                        fill_color,
                    ))
                    .map_err(|e| anyhow::anyhow!("failed to draw area series: {e}"))?;
                // Draw border line on top
                chart
                    .draw_series(LineSeries::new(
                        points.iter().copied(),
                        color.stroke_width(options.line_width as u32),
                    ))
                    .map_err(|e| anyhow::anyhow!("failed to draw histogram line: {e}"))?;
            } else {
                chart
                    .draw_series(LineSeries::new(
                        points.iter().copied(),
                        color.stroke_width(options.line_width as u32),
                    ))
                    .map_err(|e| anyhow::anyhow!("failed to draw histogram line: {e}"))?;
            }

            y_offset += baseline_sep;
        }
    }

    let img: RgbImage = image::ImageBuffer::from_vec(width, height, pixel_buffer)
        .ok_or_else(|| anyhow::anyhow!("plot image buffer had unexpected size"))?;

    let mut encoded_data = Vec::new();
    let mut encoder = image::codecs::jpeg::JpegEncoder::new_with_quality(&mut encoded_data, 85);
    encoder
        .encode(img.as_raw(), width, height, image::ExtendedColorType::Rgb8)
        .map_err(|e| anyhow::anyhow!("failed to JPEG encode plot: {e}"))?;

    Ok(encoded_data)
}

fn render_empty_histogram(
    options: &crate::options::HistogramPlotOptions,
    width: u32,
    height: u32,
    margin: u32,
    x_label_area_size: u32,
    y_label_area_size: u32,
    x_min: f64,
    x_max: f64,
) -> Result<PlotBytes> {
    use plotters::prelude::*;

    let mut pixel_buffer = vec![255; (width * height * 3) as usize];

    {
        let backend = BitMapBackend::with_buffer(&mut pixel_buffer, (width, height));
        let root = backend.into_drawing_area();
        root.fill(&WHITE)
            .map_err(|e| anyhow::anyhow!("failed to fill plot background: {e}"))?;

        let mut chart = ChartBuilder::on(&root)
            .margin(margin)
            .x_label_area_size(x_label_area_size)
            .y_label_area_size(y_label_area_size)
            .build_cartesian_2d(x_min..x_max, 0.0f64..1.0f64)
            .map_err(|e| anyhow::anyhow!("failed to build histogram chart: {e}"))?;

        let mut mesh = chart.configure_mesh();
        mesh.x_max_light_lines(4).y_max_light_lines(4);
        if let Some(ref label) = options.x_axis.label {
            mesh.x_desc(label);
        }
        mesh.y_desc("Count")
            .draw()
            .map_err(|e| anyhow::anyhow!("failed to draw mesh: {e}"))?;
    }

    let img: RgbImage = image::ImageBuffer::from_vec(width, height, pixel_buffer)
        .ok_or_else(|| anyhow::anyhow!("plot image buffer had unexpected size"))?;

    let mut encoded_data = Vec::new();
    let mut encoder = image::codecs::jpeg::JpegEncoder::new_with_quality(&mut encoded_data, 85);
    encoder
        .encode(img.as_raw(), width, height, image::ExtendedColorType::Rgb8)
        .map_err(|e| anyhow::anyhow!("failed to JPEG encode plot: {e}"))?;

    Ok(encoded_data)
}