Making 2D cross-sectional slices of a 3D array

2019-05-25 20:24发布

I'm trying to make realtime mock-ultrasound images from volumetric CT data. The trick is that the user controls the position of the probe, which defines the plane they are seeing.

What I've done so far is read the pixel data from all the dicom images into a single 3D array of pixels, and now what I need to do is reslice that 3D array at different angles. Sorry if the following description comes off a little sloppy, but imagine a 3D rectangular box (say 100 pixels wide and deep [x,z], and 500 long [y]) and a 2D "viewing plane" (say 50 x 50 pixels). Say the starting position of the viewing plane (origin defined as the middle point at the close edge of the plane - [0,25]) is with the origin at [50,250,0] (dead center of the top surface, looking down), oriented left to right and piercing the rectangle straight down. Thus, the viewing plane has three parameters that can be changed - the location of the origin, the rotation around the vertical (the line running from the origin to the corresponding point on the opposite edge of the plane), and the "tilt" (rotation of the plane around the line where it intersects with the box). So the user can change those three parameters, and the output is an image built from the pixels "touched" by the viewing plane.

Again, I apologize if the description is sloppy, but I'm a med student without a strong mathematics background. Any help would be greatly appreciated.

2条回答
虎瘦雄心在
2楼-- · 2019-05-25 20:42

Sounds like an interesting problem and I started thinking about it but soon ran into some issues. It is not as simple or straightforward as you may first think! As a start I simplified it to a case of taking a 1D slice through a 2D array. It soon became clear that for some slices it was not obvious for all pixels which ones would form part of the slice. I have produced a pdf to show what I mean. This is the link to the pdf document Issues 2D. I or others will need more thought before coming up with a possible solution. Sorry I cannot be of more help at the moment.

查看更多
成全新的幸福
3楼-- · 2019-05-25 20:49

I would write the 2D equation for a line, solve for each value of x, and round the resulting y variable to the nearest integer– Edje09 yesterday

Sticking to the 2D case for the moment the method you suggest has two main issues

  1. If the line is steeper than gradient of 1 then some pixels can be missed.
  2. rounding can choose a pixel above the one you would want to choose.

This pdf shows the issues and a possible solution for the 2D case which can then be built on for the 3D case.

EDIT After further thought I may have produced an written pdf outline solution for the 3D case that can be turned into an algorithm and hence into code. This is as far as I have got I have done no checking and cannot guarentee its correctness but hopefully will take you a stage further.

EDIT CODE ADDED The following Javascript code seems to do what you require. It is quite slow so you need to wait after clicking SET. Also the 'pane' does not clear between views so you cannot tell anything is happening until the 'pane' is refilled. I have only tested by using 2 images to represent 100 pixels in the z direction. The first code line in the function getPixels deals with this limitation, remove for a full set of images in the z direction. The test I have carried out are fairly superficial but seem to pass OK. Better with a full set of images.

I have imagined the 3D array as a series of D images image(0) at the back running the z direction to image(D-1) at the front. Each image having width W in the x direction and height H in the y direction. Thanks for the challenge I enjoyed it.

Links to a zipped folder of images used is at end of code.

<!DOCTYPE HTML>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<!--
Copyright (c)  2013   John King
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-->
<title>3D Slicer</title>
<style type="text/css">
    div, canvas, img {
        position: absolute;
    }

    img {
        top:0px;
        left:0px;
        visibility:hidden;
    }
    input {
        text-align: right;
    }
    .seen {
        visibility: visible;
    }   
    #canvas3D {
        left:10px;
        top:10px;
        visibility:hidden;
    }
    #canvas2D {
        left:10px;
        top:50px;
        border:1px solid black;
    }
    #frame {
        left:650px;
        top:10px;
        border:1px solid black;
        background-color: #DDDDDD;
        width:600px;
        height:600px;
    }
    #framehead {
        left:0px;
        top:0px;
        height:25px;
        width:100%;
        border-bottom: 1px solid black;
        background-color: #999999;
    }
    #userdata {
        top:10px;
        left:10px;
    }
    #originins {
        top:10px;
        left:10px;
        width:260px;
    }
    #origintext {
        top:200px;
        left:10px;
        width:260px;
    }
    #origininput {
        top:225px;
        left:10px;
        width:260px;        
    }
    #originlimits {
        top:250px;
        left:10px;
        width:260px;
    }
    #thetaimg {
        top:10px;
        left:225px;
    }
    #thetatext {
        top:200px;
        left:225px;
        width:260px;
    }
    #thetainput {
        top:225px;
        left:225px;
        width:260px;
    }
    #thetalimits {
        top:250px;
        left:225px;
        width:260px;
    }
    #psiimg {
        top:10px;
        left:440px;
    }
    #psitext {
        top:200px;
        left:440px;
        width:260px;
    }
    #psiinput {
        top:220px;
        left:440px;
        width:260px;
    }
    #psilimits {
        top:250px;
        left:440px;
        width:260px;
    }
    #setButton {
        top:310px;
        left:10px;
        width:260px;
    }
    #axes {
        top:350px;
        left:10px;
    }
</style>
<script type="text/javascript">

    //add a trim function to string if not present - strips white space from start and end of string
    if(typeof String.prototype.trim !== 'function') {
        String.prototype.trim = function() {
            return this.replace(/^\s+|\s+$/g, ''); 
        }
    }

    // abbreviation function for getElementById
    function $(id) {
      return document.getElementById(id);
    }

    //parameters for 3D array of pixels set in code
    var W=100; //width of array in x direction, must be even
    var D=100; //depth of array in z direction, must be even
    var H=500; //height of array in y direction

    //parameters for the rectangular plane PQRS that will select the pixels for a 2D array by slicing through the 3D array
    //PQRS moves in such a way that PQ remains parallel to xz plane and PS remains parallel to yz plane
    //these parameters set in code
    var L=50; //length of rectangle PQ
    var B=50; //breadth of rectangle PS

    //Initialisation of parameters that can be changed by the user.
    var O=new Point(W/2,0,D/2); //O is middle of PQ
    var theta=0; //angle PQ is rotated after plane is rotated about a vertical axis through O, must be between -PI/2 and PI/2
    var psi=0; //angle PS is rotated after plane is rotated about PQ as an axis, must be between -PI/2 and PI/2

    //variable for canvases
    var c3D, c2D;


    /*getPixel gets an individual pixel from the 3D array of pixels formed by a stack of D (for depth) 2D images
     * numbered from 0 to D-1, with 0 being the image at the back.
     * Each image having width W and height H pixels.
     * 0<= x <W, 0<= y <H,  0<= z <D
     * each image is on the canvas canvas3D
     * 
     * for this test img0.jpg will be used for img0.jpg to img49.jpg  and img50.jpg will be used for img50 to img99
     */
    function getPixel(x,y,z) {
        // line below only required because just two images img0.jpg and img50.jpg are used for testing
        z=Math.floor(z/50)*50;          
        //Remove above line if full series of images used in z direction
        this.ctx.drawImage($("i"+z),0,0);
        var imdata=this.ctx.getImageData(0,0,this.width,this.height);
        var col=4*(y*this.width+x);
        var pix=new Pixel();
        pix.red=imdata.data[col++];
        pix.green=imdata.data[col++];
        pix.blue=imdata.data[col++];
        pix.alpha=imdata.data[col];
        return pix;
    }

    //Pixel Object
    function Pixel() {
        this.red;
        this.green;
        this.blue;
        this.alpha;
    }

    //Point Object
    function Point(x,y,z) {
        this.x=x;
        this.y=y;
        this.z=z;
    }

    function Point2D(a,d) {
        this.a=a;
        this.d=d;
    }

    function setValues() {
        c2D.ctx.clearRect(0,0,c2D.width,c2D.height);
        var Oobj=Ochecked($("Oin").value);
        if(!Oobj.OK) {
            $("Oin").style.backgroundColor="#F1B7B7";
            return
        }
        $("Oin").style.backgroundColor="#FFFFFF";
        O=Oobj.point;
        var th=parseInt($("thetain").value.trim());
        if(isNaN(th)) {
            $("thetain").style.backgroundColor="#F1B7B7";
            return
        }
        if(th<=-90 || th>90) {
            $("thetain").style.backgroundColor="#F1B7B7";
            return
        }
        $("thetain").style.backgroundColor="#FFFFFF";
        theta=th*Math.PI/180;
        var si=parseInt($("psiin").value.trim());
        if(isNaN(si)) {
            $("psiin").style.backgroundColor="#F1B7B7";
            return
        }
        if(si<=-90 || si>90) {
            $("psiin").style.backgroundColor="#F1B7B7";
            return
        }
        $("psiin").style.backgroundColor="#FFFFFF";
        psi=si*Math.PI/180;
        printPane();
    }

    function Ochecked(Ovalue) {
        Ovalue=Ovalue.trim();
        var V=Ovalue.split(",");
        if(V.length!=3) {return {OK:false}};
        var x=parseInt(V[0].trim());
        var y=parseInt(V[1].trim());
        var z=parseInt(V[2].trim());
        if(isNaN(x) || isNaN(y) || isNaN(z))  {return {OK:false}};
        if(x<0 || x>=W) {return {OK:false}};
        if(y<0 || y>=H) {return {OK:false}};
        if(z<0 || z>=D) {return {OK:false}};
        p=new Point(x,y,z);
        return {OK:true,point:p};
    }

    function printPane(){
        var p = new Point(O.x-Math.round((L/2)*Math.cos(theta)),O.y,O.z - Math.round((L/2)*Math.sin(theta)));
        var q = new Point(O.x+Math.round((L/2)*Math.cos(theta)),O.y,O.z + Math.round((L/2)*Math.sin(theta)));
        var s = new Point(p.x,p.y+Math.round((B)*Math.cos(psi)),p.z + Math.round((B)*Math.sin(psi)));
        var n = new Point2D(q.x-p.x,q.z-p.z);       
        var PQincVec=getIncVec(n.a,n.d);
        n = new Point2D(s.y-p.y,s.z-p.z);       
        var PSincVec=getIncVec(n.a,n.d);
        var pixel,col;
        var PSpoint =new Point(p.x,p.y,p.z); // points along PS initialised to start at P
        var PQpoint; //variable for points along line parallel to PQ
        var imdata=c2D.ctx.getImageData(0,0,c2D.width,c2D.height);          
        for(var ps=0;ps<PSincVec.length;ps++) {
            //increment along line PS
            PSpoint.y+=PSincVec[ps].a;
            PSpoint.z+=PSincVec[ps].d;
            PQpoint =new Point(PSpoint.x,PSpoint.y,PSpoint.z); // points along line parallel to PQ initialised to current point on PS
            for(var pq=0;pq<PQincVec.length;pq++) {
                //increment along line PQ
                PQpoint.x+=PQincVec[pq].a;
                PQpoint.z+=PQincVec[pq].d;
                //check that PQpoint is inside 3D array
                if(0<=PQpoint.x && PQpoint.x<W && 0<=PQpoint.y && PQpoint.y<H && 0<=PQpoint.z && PQpoint.z<D) {
                    pixel=c3D.getPixel(PQpoint.x,PQpoint.y,PQpoint.z);
                    //write pixel from point along line parallel to PQ onto plane
                    col=4*(ps*c2D.width+pq);
                    imdata.data[col++]=pixel.red;
                    imdata.data[col++]=pixel.green;
                    imdata.data[col++]=pixel.blue;
                    imdata.data[col]=pixel.alpha;
                }               
            }
        }
        c2D.ctx.putImageData(imdata,0,0);
    }

    function getIncVec(a,d) {
        var r,t;
        if(a>Math.abs(d)) {
            var incVec=getIncs(a,Math.abs(d));
        }
        else {
            var incVec=getIncs(Math.abs(d),a);
            for(var i=0;i<incVec.length;i++) {
                r=incVec[i];
                t=r.a;
                r.a=r.d;
                r.d=t;
            }
        }
        if(d<0) {
            for(var i=0;i<incVec.length;i++) {
                incVec[i].d*=-1;
            }
        }
        return incVec;
    }

    function getIncs(a,d) {
        var p=new Point2D(0,0);
        var vec=[];
        vec.push(p);
        for(var i=0;i<a;i++) {
            p=new Point2D(1,Math.floor((i+1)*d/a) - Math.floor(i*d/a));
            vec.push(p);
        }
        return vec;
    }

    function main() {
        //set limits and values for user input.
        $("Oin").value=O.x+","+O.y+","+O.z;
        $("thetain").value=theta;
        $("psiin").value=psi;
        $("originlimits").innerHTML="0&lt;= x &lt;"+W+"<br>0&lt;= y &lt;"+H+"<br>0&lt;= z &lt;"+D;
        //set canvas3D so that pixels are readable
        c3D=$("canvas3D");
        c3D.width=W;
        c3D.height=H;
        c3D.ctx=c3D.getContext('2d');
        c3D.getPixel=getPixel;

        //set canvas2D so that pixels are settable
        c2D=$("canvas2D");
        c2D.width=L;
        c2D.height=B;
        c2D.ctx=c2D.getContext('2d');
        c2D.initialise=initialise;

        $("hide").style.width=L+"px";
        $("hide").style.height=B+"px";
    }
</script>
</head>
<body onload="main()">
    <!-- list of images for 3D array -->
    <img id="i0" src="images/img0.jpg">
    <img id="i50" src="images/img50.jpg">
    <!-- end of list of images for 3D array -->

    <canvas id="canvas3D"></canvas>
    <div id="frame">
        <div id="framehead">&nbsp;&nbsp;&nbsp;View of Slicing Pane</div>
        <canvas id="canvas2D"></canvas>
    </div>
    <div id="userdata">
        <div id="originins">Enter in form x,y,z </br> eg 40,27,83</div>
        <div id="origintext">Position for Origin O</div>
        <div id="origininput"><input id="Oin"></div>
        <div id="originlimits">limits</div>
        <img class="seen" id="thetaimg" src="images/theta.png">
        <div id="thetatext">Theta in degrees</div>
        <div id="thetainput"><input id="thetain"></div>
        <div id="thetalimits">-90 &lt; theta &lt;=90</div>
        <img class="seen" id="psiimg" src="images/psi.jpg">
        <div id="psitext">Psi in degrees</div>
        <div id="psiinput"><input id="psiin"></div>
        <div id="psilimits">-90 &lt; psi &lt;=90</div>
        <div id="setButton"><input type="button" value="SET" onclick="setValues()"></div>
        <img class="seen" id="axes" src="images/axes.jpg">
    </div>
<div id="msg"></div>    
</body>
</html>

images used in code

查看更多
登录 后发表回答