Metis_xctg/Data/SurfacePlot/Volume.cs

using System;
using System.Drawing;
using System.Runtime.CompilerServices;

namespace PaintDotNet.Data.SurfacePlot
{
    class Volume
    {
        Image imp;

        private int widthX;
        private int heightY;
        private int depthZ;

        private int widthXm1;
        private int heightYm1;
        private int depthZm1;

        int[][] volumeLimits = new int[8][];
        int[][] cornerT = new int[8][];

        //double a = 0, b = 1;
        //private double min, max;
        //double oldMin, oldMax;
        //bool firstTime = true;

        private double zAspect = 1;
        private byte[][][] data3D;
        private int cutDist = -1000;

        private Lut lut = null;
        private int[] bufferPixels;
        private double[] zbufferPixels;
        private int bufferWidth;
        private int bufferHeight;
        private Transform tr;
        private int drawMode;
        private int threshold;
        //private int lutNr;
        private int dotsDeltaX = 1;  // subsampling factor in x direction (used by dots drawing)
        private int dotsDeltaY = 1;  // subsampling factor in y direction (used by dots drawing)
        private int dotsDeltaZ = 1;  // subsampling factor in z direction (used by dots drawing)
        private int renderDepth = 6;

        public Volume(Image imp)
        {
            this.imp = imp;
            //this.ip = imp.getProcessor();

            widthX = imp.Width;
            heightY = imp.Height;
            //这个应该是获取最亮的值
            depthZ = 255; //imp.getStackSize();

            volumeLimits[0][0] = 0; volumeLimits[0][1] = 0; volumeLimits[0][2] = 0;         // left  front down [0] 
            volumeLimits[1][0] = 0; volumeLimits[1][1] = heightY; volumeLimits[1][2] = depthZ;   // left  back  top  [1]
            volumeLimits[2][0] = widthX; volumeLimits[2][1] = 0; volumeLimits[2][2] = depthZ;   // right front top  [2]
            volumeLimits[3][0] = widthX; volumeLimits[3][1] = heightY; volumeLimits[3][2] = 0;      // right back  down [3]
            volumeLimits[4][0] = 0; volumeLimits[4][1] = 0; volumeLimits[4][2] = depthZ;   // left  front top  [4]
            volumeLimits[5][0] = 0; volumeLimits[5][1] = heightY; volumeLimits[5][2] = 0;       // left  back  down [5]
            volumeLimits[6][0] = widthX; volumeLimits[6][1] = 0; volumeLimits[6][2] = 0;        // right front down [6] 
            volumeLimits[7][0] = widthX; volumeLimits[7][1] = heightY; volumeLimits[7][2] = depthZ;   // right back  top  [7]


            widthXm1 = widthX - 1;
            heightYm1 = heightY - 1;
            depthZm1 = depthZ - 1;

            //data3D = new byte[depthZ][heightY][widthX]; //z y x
            data3D = new byte[depthZ][][];
            getMinMax();

            init();
        }

        public void setBuffers(int[] bufferPixels, double[] zbufferPixels, int bufferWidth, int bufferHeight)
        {
            this.bufferPixels = bufferPixels;
            this.zbufferPixels = zbufferPixels;
            this.bufferWidth = bufferWidth;
            this.bufferHeight = bufferHeight;
        }

        public void setTransform(Transform transform)
        {
            this.tr = transform;
        }

        private bool getMinMax()
        {
            /**
            ip.MinMaxLoc(out min, out max);
            Calibration cal = imp.getCalibration();
            if (cal != null)
            {
                if (cal.calibrated())
                {

                    min = cal.getCValue((int)min);
                    max = cal.getCValue((int)max);

                    double[] coef = cal.getCoefficients();
                    if (coef != null)
                    {
                        a = coef[0];
                        b = coef[1];
                    }
                }
            }

            if (firstTime)
            {
                if (zAspect == 1)
                    zAspect = (float)(cal.pixelDepth / cal.pixelWidth);
                if (zAspect == 0)
                    zAspect = 1;

                lut = new Lut();
                lut.readLut(imp);

                oldMin = min;
                oldMax = max;
                firstTime = false;
            }

            bool changed;
            if (oldMin == min && oldMax == max)
                changed = false;
            else
                changed = true;
            oldMin = min;
            oldMax = max;

            return changed;
            **/
            return false;
        }


        private void init()
        {
            /**
            ImageStack stack = imp.getStack();
            int bitDepth = imp.Depth();
            if (bitDepth == 8)
            {
                float scale = (float)(255f / (max - min));

                for (int z = 0; z < depthZ; z++)
                {
                    byte[] pixels = (byte[])stack.getPixels(z + 1);

                    int pos = 0;
                    for (int y = 0; y < heightY; y++)
                    {
                        for (int x = 0; x < widthX; x++)
                        {
                            int val = 0xff & pixels[pos++];
                            val = (int)((val - min) * scale);

                            if (val < 0f) val = 0;
                            if (val > 255) val = 255;

                            data3D[z][y][x] = (byte)(0xff & val);
                        }
                    }
                }
            }

            if (bitDepth == 16)
            {
                float scale = (float)(255f / (max - min));
                for (int z = 0; z < depthZ; z++)
                {
                    short[] pixels = (short[])stack.getPixels(z + 1);
                    int pos = 0;
                    for (int y = 0; y < heightY; y++)
                    {
                        for (int x = 0; x < widthX; x++)
                        {
                            int val = (int)((int)(0xFFFF & pixels[pos++]) * b + a - min);
                            if (val < 0f) val = 0;
                            val = (int)(val * scale);

                            if (val > 255) val = 255;
                            data3D[z][y][x] = (byte)(val);
                        }
                    }
                }
            }

            if (bitDepth == 32)
            {
                float scale = (float)(255f / (max - min));

                for (int z = 0; z < depthZ; z++)
                {
                    float[] pixels = (float[])stack.getPixels(z + 1);

                    int pos = 0;
                    for (int y = 0; y < heightY; y++)
                    {
                        for (int x = 0; x < widthX; x++)
                        {
                            float value = (float)(pixels[pos++] - min);
                            if (value < 0f) value = 0f;
                            int ivalue = (int)(value * scale);

                            if (ivalue > 255) ivalue = 255;
                            data3D[z][y][x] = (byte)(ivalue);
                        }
                    }
                }
            }
            **/
        }


        private void volumeDotsOrig()
        {
            int[] v = new int[3];
            byte[][] datay;
            byte[] datax;

            double step = dotsDeltaX / (widthX - 1);
            for (int z = 0; z < data3D.Length; z += dotsDeltaZ)
            {
                v[2] = z;
                datay = data3D[z];
                for (int y = 0; y < datay.Length; y += dotsDeltaY)
                {
                    v[1] = y;
                    datax = datay[y];

                    v[0] = 0;
                    tr.xyzPos(v);
                    double x1 = tr.X;
                    double y1 = tr.Y;
                    double z1 = tr.Z;

                    v[0] = datax.Length - 1;
                    tr.xyzPos(v);

                    double dx = (tr.X - x1) * step;
                    double dy = (tr.Y - y1) * step;
                    double dz = (tr.Z - z1) * step;

                    for (int x = 0; x < widthX; x += dotsDeltaX)
                    {
                        int val = 0xff & datax[x];

                        if (val >= threshold)
                        {

                            double z_ = z1;

                            if (z_ > cutDist)
                            {
                                int x_ = (int)(x1);
                                int y_ = (int)(y1);

                                if (x_ >= 0 && y_ >= 0 && x_ < bufferWidth && y_ < bufferHeight)
                                {
                                    int pos = y_ * bufferWidth + x_;

                                    if (z_ < zbufferPixels[pos])
                                    {
                                        zbufferPixels[pos] = z_;
                                        bufferPixels[pos] = val + 1;
                                    }
                                }
                            }
                        }
                        z1 += dz;
                        y1 += dy;
                        x1 += dx;
                    }
                }
            }
            for (int i = bufferPixels.Length - 1; i >= 0; i--)
            {
                int val = bufferPixels[i];
                if (val > 255)
                    val = 255;
                if (val > 0)
                    bufferPixels[i] = lut.colors[val];
            }
        }

        private void volumeDots()
        {
            int[] v = new int[3];

            byte[][] datay;
            byte[] datax;

            for (int z = 0; z < data3D.Length; z += dotsDeltaZ)
            {
                v[2] = z;
                datay = data3D[z];
                for (int y = 0; y < datay.Length; y += dotsDeltaY)
                {
                    v[1] = y;
                    datax = datay[y];

                    for (int x = 0; x < widthX; x += dotsDeltaX)
                    {
                        v[0] = x;

                        tr.xyzPos(v);
                        double x1 = tr.X;
                        double y1 = tr.Y;
                        double z1 = tr.Z;

                        int val = 0xff & datax[x];

                        if (val >= threshold)
                        {

                            if (z1 > cutDist)
                            {
                                int x_ = (int)(x1);
                                int y_ = (int)(y1);

                                if (x_ >= 0 && y_ >= 0 && x_ < bufferWidth && y_ < bufferHeight)
                                {
                                    int pos = y_ * bufferWidth + x_;

                                    if (z1 < zbufferPixels[pos])
                                    {
                                        zbufferPixels[pos] = z1;
                                        bufferPixels[pos] = lut.colors[val];
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }

        private void volumeSliceNearestNeighbor()
        {
            int xMin = bufferWidth - 1, xMax = 0, yMin = bufferHeight - 1, yMax = 0;

            for (int i = 0; i < 8; i++)
            {

                tr.xyzPos(volumeLimits[i]);

                int xi = (int)tr.X;
                if (xi < xMin) xMin = xi;
                if (xi > xMax) xMax = xi;
                int yi = (int)tr.Y;
                if (yi < yMin) yMin = yi;
                if (yi > yMax) yMax = yi;
            }

            xMin = (xMin < 0) ? 0 : xMin;
            xMax = (xMax >= bufferWidth) ? bufferWidth - 1 : xMax;
            yMin = (yMin < 0) ? 0 : yMin;
            yMax = (yMax >= bufferHeight) ? bufferHeight - 1 : yMax;

            int[] v = new int[3];

            v[2] = cutDist; // z
            for (int y = yMin; y < yMax; y++)
            {
                v[1] = y;
                int posOffset = y * bufferWidth;

                for (int x = xMin; x < xMax; x++)
                {
                    v[0] = x;
                    tr.invxyzPosf(v);

                    double x_ = tr.x;
                    double y_ = tr.y;
                    double z_ = tr.z;

                    if (x_ >= 0 && x_ < widthX && y_ >= 0 && y_ < heightY && z_ >= 0 && z_ < depthZ)
                    {

                        int val = 0xff & data3D[(int)z_][(int)y_][(int)x_];
                        int pos = posOffset + x;

                        int v_ = (cutDist * 20) + 128;
                        v_ = Math.Min(Math.Max(0, v_), 255);
                        v_ = (int)(0xFF000000 | (v_ << 8));
                        bufferPixels[pos] = lut.colors[val];
                        zbufferPixels[pos] = cutDist;
                    }
                }
            }
        }



        private void volumeSliceTrilinear()
        {
            int xMin = bufferWidth - 1, xMax = 0, yMin = bufferHeight - 1, yMax = 0;

            for (int i = 0; i < 8; i++)
            {

                tr.xyzPos(volumeLimits[i]);

                int xi = (int)tr.X;
                if (xi < xMin) xMin = xi;
                if (xi > xMax) xMax = xi;
                int yi = (int)tr.Y;
                if (yi < yMin) yMin = yi;
                if (yi > yMax) yMax = yi;
            }

            xMin = (xMin < 0) ? 0 : xMin;
            xMax = (xMax >= bufferWidth) ? bufferWidth - 1 : xMax;
            yMin = (yMin < 0) ? 0 : yMin;
            yMax = (yMax >= bufferHeight) ? bufferHeight - 1 : yMax;

            int[] v = new int[3];

            // draw slice at dist
            v[2] = cutDist; // z

            for (int y = yMin; y <= yMax; y++)
            {
                v[1] = y;
                int posOffset = y * bufferWidth;
                for (int x = xMin; x <= xMax; x++)
                {
                    v[0] = x;
                    tr.invxyzPosf(v);

                    double x_ = tr.x;
                    double y_ = tr.y;
                    double z_ = tr.z;

                    if ((x_ >= 0 && x_ < widthX) && (y_ >= 0 && y_ < heightY) && (z_ >= 0 && z_ < depthZ))
                    {
                        int pos = posOffset + x;

                        int val = trilinear(z_, y_, x_);

                        bufferPixels[pos] = lut.colors[val];
                        zbufferPixels[pos] = cutDist;
                    }
                }
            }
        }

        private bool isInside(int x, int y)
        {
            int[] p = new int[3];
            p[0] = x;
            p[1] = y;

            if (Misc.inside(p, cornerT[0], cornerT[1], cornerT[4])) return true;
            if (Misc.inside(p, cornerT[0], cornerT[1], cornerT[5])) return true;
            if (Misc.inside(p, cornerT[2], cornerT[3], cornerT[6])) return true;
            if (Misc.inside(p, cornerT[2], cornerT[3], cornerT[7])) return true;

            if (Misc.inside(p, cornerT[1], cornerT[3], cornerT[5])) return true;
            if (Misc.inside(p, cornerT[1], cornerT[3], cornerT[7])) return true;
            if (Misc.inside(p, cornerT[0], cornerT[2], cornerT[4])) return true;
            if (Misc.inside(p, cornerT[0], cornerT[2], cornerT[6])) return true;

            if (Misc.inside(p, cornerT[1], cornerT[2], cornerT[4])) return true;
            if (Misc.inside(p, cornerT[1], cornerT[2], cornerT[7])) return true;
            if (Misc.inside(p, cornerT[0], cornerT[3], cornerT[5])) return true;
            if (Misc.inside(p, cornerT[0], cornerT[3], cornerT[6])) return true;

            return false;
        }

        [MethodImplAttribute(MethodImplOptions.Synchronized)]
        public void volumeProjection_trilinear_front()
        {
            double width = widthX - 0.5f;
            double height = heightY - 0.5f;
            double depth = depthZ - 0.5f;

            transformVolumeLimits();

            int xMin = bufferWidth - 1, xMax = 0;
            int yMin = bufferHeight - 1, yMax = 0;
            int zMin = cutDist, zMax = cutDist;

            for (int i = 0; i < 8; i++)
            {
                tr.xyzPos(volumeLimits[i]);

                int xi = (int)tr.X;
                if (xi < xMin) xMin = xi;
                if (xi > xMax) xMax = xi;
                int yi = (int)tr.Y;
                if (yi < yMin) yMin = yi;
                if (yi > yMax) yMax = yi;
                int zi = (int)tr.Z;
                if (zi < zMin) zMin = zi;
                if (zi > zMax) zMax = zi;
            }

            xMin = (xMin < 0) ? 0 : xMin;
            xMax = (xMax >= bufferWidth) ? bufferWidth - 1 : xMax;
            yMin = (yMin < 0) ? 0 : yMin;
            yMax = (yMax >= bufferHeight) ? bufferHeight - 1 : yMax;

            int[] v1 = new int[3];
            int[] v2 = new int[3];

            if (cutDist > zMin)
                zMin = cutDist;

            int nd = zMax - zMin;

            if (nd > 0)
            {
                float scaleLum = (renderDepth > 1) ? (renderDepth * 255 / ((renderDepth - 1) * 255f)) : 2; ;
                v1[2] = zMin;
                v2[2] = zMax;

                double nd1 = 1 / nd;
                double d_z = (zMax - zMin) * nd1;

                for (int y = yMin; y <= yMax; y++)
                {
                    v1[1] = v2[1] = y;
                    int percent = 100 * (y - yMin) / (yMax - yMin);

                    for (int x = xMin; x <= xMax; x++)
                    {

                        if (isInside(x, y))
                        {
                            double z = zMin;

                            int pos = y * bufferWidth + x;
                            v1[0] = v2[0] = x;

                            int alpha = 255;
                            int V = 0;
                            int alpha2 = 255;
                            int V2 = 0;
                            tr.invxyzPosf(v1);
                            double xd1 = tr.x;
                            double yd1 = tr.y;
                            double zd1 = tr.z;

                            tr.invxyzPosf(v2);
                            double xd2 = tr.x;
                            double yd2 = tr.y;
                            double zd2 = tr.z;

                            double dx = (xd2 - xd1) * nd1;
                            double dy = (yd2 - yd1) * nd1;
                            double dz = (zd2 - zd1) * nd1;

                            int k = 0;

                            for (int n = nd; n >= 0; n--)
                            {
                                z += d_z;

                                if (xd1 >= 0 && xd1 < width &&
                                        yd1 >= 0 && yd1 < height &&
                                        zd1 >= 0 && zd1 < depth)
                                {
                                    int val1 = trilinear(zd1, yd1, xd1);

                                    int f = 255 - 10 * Math.Abs(val1 - threshold);
                                    if (f < 0)
                                        f = 0;
                                    if (f > 10)
                                    {
                                        alpha += f;
                                        V += f * 255;

                                        if (++k >= renderDepth)
                                            break;
                                    }
                                }

                                xd1 += dx;
                                yd1 += dy;
                                zd1 += dz;
                            }
                            int val = (int)(scaleLum * V / alpha);
                            if (val > 255)
                                val = 255;

                            int val2 = (int)(scaleLum * V2 / alpha2);
                            if (val2 > 255)
                                val2 = 255;
                            bufferPixels[pos] = lut.colors[val];
                            zbufferPixels[pos] = z;
                        }
                    }
                }
            }
        }

        private void volumeProjection_trilinear_back()
        {

            double width = widthX - 0.5f;
            double height = heightY - 0.5f;
            double depth = depthZ - 0.5f;

            transformVolumeLimits();

            int xMin = bufferWidth - 1, xMax = 0;
            int yMin = bufferHeight - 1, yMax = 0;
            int zMin = 10000, zMax = -10000;

            for (int i = 0; i < 8; i++)
            {
                tr.xyzPos(volumeLimits[i]);

                int xi = (int)tr.X;
                if (xi < xMin) xMin = xi;
                if (xi > xMax) xMax = xi;
                int yi = (int)tr.Y;
                if (yi < yMin) yMin = yi;
                if (yi > yMax) yMax = yi;
                int zi = (int)tr.Z;
                if (zi < zMin) zMin = zi;
                if (zi > zMax) zMax = zi;
            }

            xMin = (xMin < 0) ? 0 : xMin;
            xMax = (xMax >= bufferWidth) ? bufferWidth - 1 : xMax;
            yMin = (yMin < 0) ? 0 : yMin;
            yMax = (yMax >= bufferHeight) ? bufferHeight - 1 : yMax;

            xMin = (xMin < 0) ? 0 : xMin;
            xMax = (xMax > 511) ? 511 : xMax;
            yMin = (yMin < 0) ? 0 : yMin;
            yMax = (yMax > 511) ? 511 : yMax;

            int[] v1 = new int[3];
            int[] v2 = new int[3];

            if (cutDist < zMin)
                cutDist = zMin;

            int nd = zMax - cutDist;

            double dx, dy, dz;

            int[] vals = new int[zMax - cutDist];

            if (nd > 0)
            {
                v1[2] = cutDist;
                v2[2] = zMax;

                double nd1 = 1f / nd;
                for (int y = yMin; y <= yMax; y++)
                {
                    v1[1] = v2[1] = y;
                    int percent = 100 * (y - yMin) / (yMax - yMin);
                    for (int x = xMin; x <= xMax; x++)
                    {
                        if (isInside(x, y))
                        {
                            int pos = y * bufferWidth + x;

                            v1[0] = v2[0] = x;

                            int V = 0;
                            tr.invxyzPosf(v1);
                            double x1 = tr.x;
                            double y1 = tr.y;
                            double z1 = tr.z;

                            tr.invxyzPosf(v2);

                            dx = (tr.x - x1) * nd1;
                            dy = (tr.y - y1) * nd1;
                            dz = (tr.z - z1) * nd1;

                            int k = 0;
                            for (int n = nd; n >= 0; n--)
                            {
                                if (x1 >= 0 && x1 < width &&
                                        y1 >= 0 && y1 < height &&
                                        z1 >= 0 && z1 < depth)
                                {
                                    int val1 = trilinear(z1, y1, x1);

                                    if (val1 >= threshold)
                                    {
                                        vals[k++] = val1;

                                        if (k > renderDepth)
                                            break;
                                    }
                                }

                                x1 += dx;
                                y1 += dy;
                                z1 += dz;
                            }
                            for (int i = k - 1; i >= 0; i--)
                            {
                                int a = (vals[i] - threshold) + 50;

                                if (a > 255)
                                    a = 255;

                                V = (a * vals[i] + V * (255 - a)) / 255;
                            }

                            int val = V + 1;
                            if (val > 255) val = 255;
                            bufferPixels[pos] = lut.colors[val];
                            zbufferPixels[pos] = (int)tr.z;

                        }
                    }
                }
            }
        }


        private void transformVolumeLimits()
        {
            for (int i = 0; i < 8; i++)
            {
                tr.xyzPos(volumeLimits[i]);
                cornerT[i][0] = (int)tr.X;
                cornerT[i][1] = (int)tr.Y;
                cornerT[i][2] = (int)tr.Z;
            }
        }

        private int trilinear(double z, double y, double x)
        {
            int tx = (int)x;
            double dx = x - tx;
            int tx1 = (tx < widthXm1) ? tx + 1 : tx;
            int ty = (int)y;
            double dy = y - ty;
            int ty1 = (ty < heightYm1) ? ty + 1 : ty;
            int tz = (int)z;
            double dz = z - tz;
            int tz1 = (tz < depthZm1) ? tz + 1 : tz;

            int v000 = (int)(0xff & data3D[tz][ty][tx]);
            int v001 = (int)(0xff & data3D[tz1][ty][tx]);
            int v010 = (int)(0xff & data3D[tz][ty1][tx]);
            int v011 = (int)(0xff & data3D[tz1][ty1][tx]);
            int v100 = (int)(0xff & data3D[tz][ty][tx1]);
            int v101 = (int)(0xff & data3D[tz1][ty][tx1]);
            int v110 = (int)(0xff & data3D[tz][ty1][tx1]);
            int v111 = (int)(0xff & data3D[tz1][ty1][tx1]);

            return (int)(
                    (v100 - v000) * dx +
                    (v010 - v000) * dy +
                    (v001 - v000) * dz +
                    (v110 - v100 - v010 + v000) * dx * dy +
                    (v011 - v010 - v001 + v000) * dy * dz +
                    (v101 - v100 - v001 + v000) * dx * dz +
                    (v111 + v100 + v010 + v001 - v110 - v101 - v011 - v000) * dx * dy * dz + v000);
        }

        public int getWidth()
        {
            return widthX;
        }
        public int getHeight()
        {
            return heightY;
        }

        public double getZaspectRatio()
        {
            return zAspect;
        }


        public void setVolumeDrawMode(int volume_drawMode)
        {
            this.drawMode = volume_drawMode;
        }

        public void draw()
        {
            if (drawMode == JRenderer3D.VOLUME_DOTS)
                volumeDots();
            else if (drawMode == JRenderer3D.VOLUME_SLICE_NEAREST_NEIGHBOR)
                volumeSliceNearestNeighbor();
            else if (drawMode == JRenderer3D.VOLUME_SLICE_TRILINEAR)
                volumeSliceTrilinear();
            else if (drawMode == JRenderer3D.VOLUME_DOTS)
                volumeDots();
            else if (drawMode == JRenderer3D.VOLUME_PROJECTION_TRILINEAR_FRONT)
                volumeProjection_trilinear_front();
        }

        public void setVolumeThreshold(int volume_threshold)
        {
            this.threshold = volume_threshold;
        }

        public void setVolumeCutDist(int volume_cutDist)
        {
            this.cutDist = volume_cutDist;
        }

        public void setVolumeLut(int volume_lutNr)
        {
            if (lut != null)
                lut.setLut(volume_lutNr);
        }

        public void setVolumeDotsSubsampling(int volume_dotsDeltaX, int volume_dotsDeltaY, int volume_dotsDeltaZ)
        {
            this.dotsDeltaX = volume_dotsDeltaX;
            this.dotsDeltaY = volume_dotsDeltaY;
            this.dotsDeltaZ = volume_dotsDeltaZ;
        }
    }
}