From 796fba9165ed35d4896656f27e412a5024a2eafa Mon Sep 17 00:00:00 2001
From: Fredrik Orderud <fredrik.orderud@ge.com>
Date: Thu, 20 Feb 2020 09:56:01 +0100
Subject: [PATCH] Implement rudimentary CF simulation and display.

The important part here is the TestViewer GenerateBitmap function, which serves as a reference implementation of color-flow display.
---
 DummyLoader/Image3dSource.cpp |  67 +++++++++++++++++++++-
 DummyLoader/Image3dSource.hpp |   5 +-
 DummyLoader/Image3dStream.cpp |  18 ++++++
 TestViewer/MainWindow.xaml.cs | 101 +++++++++++++++++++++++++++++++---
 4 files changed, 181 insertions(+), 10 deletions(-)

diff --git a/DummyLoader/Image3dSource.cpp b/DummyLoader/Image3dSource.cpp
index d254c0d..133474f 100644
--- a/DummyLoader/Image3dSource.cpp
+++ b/DummyLoader/Image3dSource.cpp
@@ -37,6 +37,39 @@ Image3dSource::Image3dSource() {
         for (size_t i = 0; i < m_color_map_tissue.size(); ++i)
             m_color_map_tissue[i] = R8G8B8A8(static_cast<unsigned char>(i), static_cast<unsigned char>(i), static_cast<unsigned char>(i), 0xFF);
     }
+    {
+        // red to blue flow scale with green at high BW
+        assert(m_color_map_flow.size() == 256*256);
+        for (size_t bw = 0; bw < 256; ++bw) {
+            // increasing green for high bandwidth
+            uint8_t green = (bw >= 192) ? static_cast<unsigned char>(bw) : 0;
+
+            for (size_t freq = 0; freq < 256; ++freq) { // unsigned counter, so freq>127 corresponds to negative velocities
+                // increasing red for positive velocities
+                uint8_t red  = (freq < 128) ? 128+static_cast<unsigned char>(freq) : 0;
+                if (green)
+                    red = 0;
+
+                // increasing blue for negative velocities
+                uint8_t blue = (freq >= 128) ? 128+static_cast<unsigned char>(255-freq) : 0;
+                if (green)
+                    blue = 0;
+
+                m_color_map_flow[freq + bw*256] = R8G8B8A8(red, green, blue, 0xFF);
+            }
+        }
+    }
+    {
+        // flow arbitration scale
+        assert(m_flow_arb.size() == 256*256);
+        for (size_t bw = 0; bw < 256; ++bw) {
+            for (size_t freq = 0; freq < 256; ++freq) { // unsigned counter, so freq>127 corresponds to negative velocities
+                // show flow when |freq| >= 16
+                bool show_flow = std::abs(static_cast<int>(freq)) >= 16;
+                m_flow_arb[freq + bw*256] = show_flow ? 0xFF : 0x00;
+            }
+        }
+    }
     {
         // image geometry    X     Y       Z
         Cart3dGeom geom = { -0.1f, 0,     -0.075f,// origin
@@ -46,8 +79,9 @@ Image3dSource::Image3dSource() {
         m_bbox = geom;
     }
 
-    // a single tissue stream
+    // simulate tissue + color-flow data
     m_stream_types.push_back(IMAGE_TYPE_TISSUE);
+    m_stream_types.push_back(IMAGE_TYPE_BLOOD_VEL);
 }
 
 Image3dSource::~Image3dSource() {
@@ -73,8 +107,23 @@ HRESULT Image3dSource::GetStream(unsigned int index, Cart3dGeom out_geom, unsign
     CComPtr<IImage3dStream> stream_obj;
     {
         // on-demand stream creation
+        Cart3dGeom bbox = m_bbox;
+        if (m_stream_types[index] == IMAGE_TYPE_BLOOD_VEL) {
+            // shrink color-flow sector to make it more realistic
+            vec3f origin, dir1, dir2, dir3;
+            std::tie(origin,dir1,dir2,dir3) = FromCart3dGeom(bbox);
+
+            float SCALE_FACTOR = 0.8f;
+            origin += 0.5f*(1-SCALE_FACTOR)*(dir1 + dir2 + dir3);
+            dir1 *= SCALE_FACTOR;
+            dir2 *= SCALE_FACTOR;
+            dir3 *= SCALE_FACTOR;
+
+            bbox = ToCart3dGeom(origin, dir1, dir2, dir3);
+        }
+
         auto stream_cls = CreateLocalInstance<Image3dStream>();
-        stream_cls->Initialize(m_stream_types[index], m_bbox, out_geom, max_resolution);
+        stream_cls->Initialize(m_stream_types[index], bbox, out_geom, max_resolution);
         stream_obj = stream_cls; // cast class pointer to interface
     }
 
@@ -103,6 +152,20 @@ HRESULT Image3dSource::GetColorMap(ColorMapType type, /*out*/ImageFormat * forma
         memcpy(&color_map.GetAt(0), m_color_map_tissue.data(), sizeof(m_color_map_tissue));
         *map = color_map.Detach(); // transfer ownership
         return S_OK;
+    } else if (type == TYPE_FLOW_COLOR) {
+        *format = IMAGE_FORMAT_R8G8B8A8;
+        // copy to new buffer
+        CComSafeArray<uint8_t> color_map(4*static_cast<unsigned int>(m_color_map_flow.size()));
+        memcpy(&color_map.GetAt(0), m_color_map_flow.data(), sizeof(m_color_map_flow));
+        *map = color_map.Detach(); // transfer ownership
+        return S_OK;
+    } else if (type = TYPE_FLOW_ARB) {
+        *format = IMAGE_FORMAT_U8;
+        // copy to new buffer
+        CComSafeArray<uint8_t> color_map(static_cast<unsigned int>(m_flow_arb.size()));
+        memcpy(&color_map.GetAt(0), m_flow_arb.data(), sizeof(m_flow_arb));
+        *map = color_map.Detach(); // transfer ownership
+        return S_OK;
     }
 
     return E_NOTIMPL;
diff --git a/DummyLoader/Image3dSource.hpp b/DummyLoader/Image3dSource.hpp
index ba292e3..3f63717 100644
--- a/DummyLoader/Image3dSource.hpp
+++ b/DummyLoader/Image3dSource.hpp
@@ -40,7 +40,10 @@ class ATL_NO_VTABLE Image3dSource :
 private:
     ProbeInfo                m_probe;
     EcgSeries                m_ecg;
-    std::array<R8G8B8A8,256> m_color_map_tissue;
+
+    std::array<R8G8B8A8,256>     m_color_map_tissue;
+    std::array<R8G8B8A8,256*256> m_color_map_flow;
+    std::array<uint8_t,256*256>  m_flow_arb;
 
     Cart3dGeom               m_bbox = {};
     std::vector<ImageType>   m_stream_types;
diff --git a/DummyLoader/Image3dStream.cpp b/DummyLoader/Image3dStream.cpp
index 8256786..37381fb 100644
--- a/DummyLoader/Image3dStream.cpp
+++ b/DummyLoader/Image3dStream.cpp
@@ -55,6 +55,24 @@ void Image3dStream::Initialize (ImageType type, Cart3dGeom img_geom, Cart3dGeom
 
             m_frames.push_back(CreateImage3d(frameNumber*(duration/numFrames) + startTime, IMAGE_FORMAT_U8, dims, img_buf));
         }
+    } else if (type == IMAGE_TYPE_BLOOD_VEL) {
+        // velocity & bandwidth scale color-flow data
+        unsigned short dims[] = { 20, 15, 10 }; // matches length of dir1, dir2 & dir3, so that the image squares become quadratic
+        std::vector<byte> img_buf(2 * dims[0] * dims[1] * dims[2]);
+        for (size_t frameNumber = 0; frameNumber < numFrames; ++frameNumber) {
+            for (unsigned int z = 0; z < dims[2]; ++z) {
+                for (unsigned int y = 0; y < dims[1]; ++y) {
+                    for (unsigned int x = 0; x < dims[0]; ++x) {
+                        int8_t & freq = reinterpret_cast<int8_t&>(img_buf[0 + 2*(x + y*dims[0] + z*dims[0] * dims[1])]);
+                        byte & bw   = img_buf[1 + 2*(x + y*dims[0] + z*dims[0] * dims[1])];
+
+                        freq = static_cast<int8_t>(255*(0.5f - y*1.0f/dims[1])); // [+127, -128] along Y axis
+                        bw   = static_cast<uint8_t>(256*(x*1.0f/dims[0]));       // [0,255] along X axis
+                    }
+                }
+            }
+            m_frames.push_back(CreateImage3d(frameNumber*(duration/numFrames) + startTime, IMAGE_FORMAT_FREQ8POW8, dims, img_buf));
+        }
     } else {
         abort(); // should never be reached
     }
diff --git a/TestViewer/MainWindow.xaml.cs b/TestViewer/MainWindow.xaml.cs
index 9dba342..7b54a3e 100644
--- a/TestViewer/MainWindow.xaml.cs
+++ b/TestViewer/MainWindow.xaml.cs
@@ -44,8 +44,11 @@ public partial class MainWindow : Window
         IImage3dSource     m_source;
 
         IImage3dStream     m_streamXY;
+        IImage3dStream     m_streamXYcf;
         IImage3dStream     m_streamXZ;
+        IImage3dStream     m_streamXZcf;
         IImage3dStream     m_streamZY;
+        IImage3dStream     m_streamZYcf;
 
         public MainWindow()
         {
@@ -69,8 +72,11 @@ void ClearUI()
             ImageZY.Source = null;
 
             m_streamXY = null;
+            m_streamXYcf = null;
             m_streamXZ = null;
+            m_streamXZcf = null;
             m_streamZY = null;
+            m_streamZYcf = null;
 
             ECG.Data = null;
 
@@ -302,6 +308,8 @@ private void InitializeSlices()
             bboxXY.dir3_y = 0;
             bboxXY.dir3_z = 0;
             m_streamXY = m_source.GetStream(0, bboxXY, new ushort[] { HORIZONTAL_RES, VERTICAL_RES, 1 });
+            if (stream_count >= 2)
+                m_streamXYcf = m_source.GetStream(1, bboxXY, new ushort[] { HORIZONTAL_RES, VERTICAL_RES, 1 }); // assume 2nd stream is color-flow
 
             // get XZ plane (assumes 1st axis is "X" and 3rd is "Z")
             Cart3dGeom bboxXZ = bbox;
@@ -315,6 +323,8 @@ private void InitializeSlices()
             bboxXZ.dir3_y = 0;
             bboxXZ.dir3_z = 0;
             m_streamXZ = m_source.GetStream(0, bboxXZ, new ushort[] { HORIZONTAL_RES, VERTICAL_RES, 1 });
+            if (stream_count >= 2)
+                m_streamXZcf = m_source.GetStream(1, bboxXZ, new ushort[] { HORIZONTAL_RES, VERTICAL_RES, 1 }); // assume 2nd stream is color-flow
 
             // get ZY plane (assumes 2nd axis is "Y" and 3rd is "Z")
             Cart3dGeom bboxZY = bbox;
@@ -331,30 +341,68 @@ private void InitializeSlices()
             bboxZY.dir3_y = 0;
             bboxZY.dir3_z = 0;
             m_streamZY = m_source.GetStream(0, bboxZY, new ushort[] { HORIZONTAL_RES, VERTICAL_RES, 1 });
+            if (stream_count >= 2)
+                m_streamZYcf = m_source.GetStream(1, bboxZY, new ushort[] { HORIZONTAL_RES, VERTICAL_RES, 1 }); // assume 2nd stream is color-flow
         }
 
-        private void DrawSlices (uint frame)
+        private void DrawSlices(uint frame)
         {
             Debug.Assert(m_source != null);
 
             ImageFormat image_format;
-            byte[] color_map = m_source.GetColorMap(ColorMapType.TYPE_TISSUE_COLOR, out image_format);
+            byte[] tissue_map = m_source.GetColorMap(ColorMapType.TYPE_TISSUE_COLOR, out image_format);
             if (image_format != ImageFormat.IMAGE_FORMAT_R8G8B8A8)
                 throw new Exception("Unexpected color-map format");
 
-            // retrieve image slices
+            byte[] cf_map = m_source.GetColorMap(ColorMapType.TYPE_FLOW_COLOR, out image_format);
+            if (image_format != ImageFormat.IMAGE_FORMAT_R8G8B8A8)
+                throw new Exception("Unexpected color-map format");
+
+            byte[] arb_table = m_source.GetColorMap(ColorMapType.TYPE_FLOW_ARB, out image_format);
+            if (image_format != ImageFormat.IMAGE_FORMAT_U8)
+                throw new Exception("Unexpected color-map format");
+
+            uint cf_frame = 0;
+            if (m_streamXYcf != null) {
+                // find closest corresponding CF frame
+                double t_time = m_streamXY.GetFrame(frame).time;
+                double[] cf_times = m_streamXYcf.GetFrameTimes();
+
+                int closest_idx = 0;
+                for (int i = 1; i < cf_times.Length; ++i) {
+                    if (Math.Abs(cf_times[i] - t_time) < Math.Abs(cf_times[closest_idx] - t_time))
+                        closest_idx = i;
+                }
+
+                cf_frame = (uint)closest_idx;
+            }
 
             // get XY plane (assumes 1st axis is "X" and 2nd is "Y")
             Image3d imageXY = m_streamXY.GetFrame(frame);
-            ImageXY.Source = GenerateBitmap(imageXY, color_map);
-
+            if (m_streamXYcf != null) {
+                Image3d imageXYcf = m_streamXYcf.GetFrame(cf_frame);
+                ImageXY.Source = GenerateBitmap(imageXY, imageXYcf, tissue_map, cf_map, arb_table);
+            } else {
+                ImageXY.Source = GenerateBitmap(imageXY, tissue_map);
+            }
+        
             // get XZ plane (assumes 1st axis is "X" and 3rd is "Z")
             Image3d imageXZ = m_streamXZ.GetFrame(frame);
-            ImageXZ.Source = GenerateBitmap(imageXZ, color_map);
+            if (m_streamXZcf != null) {
+                Image3d imageXZcf = m_streamXZcf.GetFrame(cf_frame);
+                ImageXZ.Source = GenerateBitmap(imageXZ, imageXZcf, tissue_map, cf_map, arb_table);
+            } else {
+                ImageXZ.Source = GenerateBitmap(imageXZ, tissue_map);
+            }
 
             // get ZY plane (assumes 2nd axis is "Y" and 3rd is "Z")
             Image3d imageZY = m_streamZY.GetFrame(frame);
-            ImageZY.Source = GenerateBitmap(imageZY, color_map);
+            if (m_streamZYcf != null) {
+                Image3d imageZYcf = m_streamZYcf.GetFrame(cf_frame);
+                ImageZY.Source = GenerateBitmap(imageZY, imageZYcf, tissue_map, cf_map, arb_table);
+            } else {
+                ImageZY.Source = GenerateBitmap(imageZY, tissue_map);
+            }
 
             FrameTime.Text = "Frame time: " + imageXY.time;
         }
@@ -387,6 +435,45 @@ private WriteableBitmap GenerateBitmap(Image3d t_img, byte[] t_map)
             return bitmap;
         }
 
+        private WriteableBitmap GenerateBitmap(Image3d t_img, Image3d cf_img, byte[] t_map,  byte[] cf_map, byte[] arb_table)
+        {
+            Debug.Assert(t_img.dims.SequenceEqual(cf_img.dims));
+            Debug.Assert(t_img.format == ImageFormat.IMAGE_FORMAT_U8);
+            Debug.Assert(cf_img.format == ImageFormat.IMAGE_FORMAT_FREQ8POW8);
+
+            WriteableBitmap bitmap = new WriteableBitmap(t_img.dims[0], t_img.dims[1], 96.0, 96.0, PixelFormats.Rgb24, null);
+            bitmap.Lock();
+            unsafe {
+                for (int y = 0; y < bitmap.Height; ++y) {
+                    for (int x = 0; x < bitmap.Width; ++x) {
+                        byte t_val = t_img.data[x + y * t_img.stride0];
+                        ushort cf_val = BitConverter.ToUInt16(cf_img.data, 2*x + y*(int)cf_img.stride0);
+
+                        uint rgba = 0;
+                        if (arb_table[cf_val] > t_val) {
+                            // display color-flow overlay
+                            rgba = BitConverter.ToUInt32(cf_map, 4*cf_val);
+                        } else {
+                            // display tissue underlay
+                            rgba = BitConverter.ToUInt32(t_map, 4*t_val);
+                        }
+
+                        byte[] channels = BitConverter.GetBytes(rgba);
+
+                        // assign red, green & blue
+                        byte* pixel = (byte*)bitmap.BackBuffer + x * (bitmap.Format.BitsPerPixel / 8) + y * bitmap.BackBufferStride;
+                        pixel[0] = channels[0]; // red
+                        pixel[1] = channels[1]; // green
+                        pixel[2] = channels[2]; // blue
+                        // discard alpha channel
+                    }
+                }
+            }
+            bitmap.AddDirtyRect(new Int32Rect(0, 0, bitmap.PixelWidth, bitmap.PixelHeight));
+            bitmap.Unlock();
+            return bitmap;
+        }
+
         static void SwapVals(ref float v1, ref float v2)
         {
             float tmp = v1;