add code to read MRI data

diff --git a/src/edu/berkeley/qfat/Main.java b/src/edu/berkeley/qfat/Main.java
index 56eeb83..f404144 100644 (file)
--- a/src/edu/berkeley/qfat/Main.java
+++ b/src/edu/berkeley/qfat/Main.java
@@ -318,6 +318,88 @@ public class Main extends InteractiveMeshViewer {
         //fixupTile();
     }
 
+    public void marchingCubes2() {
+        try {
+
+            final float[][][] samples = new float[256][256][124];
+            double total = 0;
+            int count = 0;
+            for(int i=1; i<=124; i++) {
+                String ix = ""+i;
+                while(ix.length() < 3) ix = "0"+ix;
+                FileInputStream fis = new FileInputStream("/Users/megacz/Desktop/projects/mri brain images/spgr/I."+ix);
+                DataInputStream dis = new DataInputStream(new BufferedInputStream(fis));
+                dis.skip(7904);
+                for(int x=0; x<256; x++)
+                    for(int y=0; y<256; y++) {
+                        short s = dis.readShort();
+                        samples[x][y][i-1] = (float)s;
+                        total += samples[x][y][i-1];
+                        count++;
+                    }
+            }
+            System.out.println("done reading samples; average sample is " + (total/count));
+
+            /*
+            PrintWriter pw = new PrintWriter(new FileOutputStream("/tmp/out"));
+            pw.println("new short[][][] {");
+            for(int x=0; x<256; x++) {
+                if (x>0) pw.println("  ,");
+                pw.println("  {");
+                for(int y=0; y<256; y++) {
+                    if (y>0) pw.println("    ,");
+                    pw.print("    {");
+                    for(int z=0; z<124; z++) {
+                        if (z>0) pw.print(", ");
+                        if ((z % 20) == 0) {
+                            pw.println();
+                            pw.print("      ");
+                        }
+                        pw.print(shorts[x][y][z]);
+                    }
+                    pw.println();
+                    pw.println("    }");
+                }
+                pw.println("  }");
+            }
+            pw.println("};");
+            pw.flush();
+            pw.close();
+            */
+
+            Mesh mesh = new Mesh(false);
+            mesh.coalesce = true;
+            MarchingCubes.march(new VoxelData() {
+                    float radius = 1.0f;
+                    public float getMaxX()        { return  2.4f; }
+                    public float getMinX()        { return    0f; }
+                    public int   getNumSamplesX() { return   256/10; }
+                    public float getMaxY()        { return  2.4f; }
+                    public float getMinY()        { return    0f; }
+                    public int   getNumSamplesY() { return   256/10; }
+                    public float getMaxZ()        { return 1.86f; }
+                    public float getMinZ()        { return    0f; }
+                    public int   getNumSamplesZ() { return   124/10; }
+                    public float getSample(Point p) {
+                        int x = (int)Math.floor(p.x / ((double)getMaxX() / (double)256));
+                        int y = (int)Math.floor(p.y / ((double)getMaxY() / (double)256));
+                        int z = (int)Math.floor(p.z / ((double)getMaxZ() / (double)124));
+                        if ( (x<0) || (x>=samples.length) ) return 0;
+                        if ( (y<0) || (y>=samples[x].length) ) return 0;
+                        if ( (z<0) || (z>=samples[x][y].length) ) return 0;
+                        return samples[x][y][z];
+                    }
+                },
+                200,
+                mesh);
+            setTile(mesh);
+            //fixupTile();
+
+        } catch (Exception e) {
+            e.printStackTrace();
+        }
+    }
+
     public void hexBrick(boolean offset, boolean rotated) {
         setTile(new Mesh(false));
                 float width  = (float)0.8;
@@ -948,6 +1030,9 @@ public class Main extends InteractiveMeshViewer {
             tileMenu.add(new MyMenuItem("Marching Cubes") { public void hit() {
                 marchingCubes();
             }});
+            tileMenu.add(new MyMenuItem("Brain") { public void hit() {
+                marchingCubes2();
+            }});
 
             // Finally, add all the menus to the menu bar.
             add(tileMenu);