Login  Register

Re: Stereo SEM image generation

Posted by Z.J. Zhang on Jan 16, 2008; 7:58pm
URL: http://imagej.273.s1.nabble.com/Stereo-SEM-image-generation-tp3697564p3697566.html

Hi Lucena:

When I tried to compile the plugin, I receive the following error message. I am a biologist and don't know anything about programming. Could you help?


-----This is the error message I received-------
Note: sun.tools.javac.Main has been deprecated.
C:\ImageJ\plugins\SEM-pair-1j.java:97: '}' expected.
showtable(cal, ip3, imp,width,height,scalex,scaley);//***********************************************linha 132 }}
                                                    ^
C:\ImageJ\plugins\SEM-pair-1j.java:99: Statement expected.
public int[] s2ints(String s) {
^
C:\ImageJ\plugins\SEM-pair-1j.java:195: '}' expected.
}//reconstruction end
 ^
C:\ImageJ\plugins\SEM-pair-1j.java:197: Statement expected.
public void showtable (Calibration cal,ImageProcessor ip3, ImagePlus imp, int width, int height, float scalex, float scaley){
^
C:\ImageJ\plugins\SEM-pair-1j.java:5: Public class Reconstruction_3Dversao5 must be defined in a file called "Reconstruction_3Dversao5.java".
public class Reconstruction_3Dversao5 implements  PlugIn {
             ^
5 errors, 1 warning
---------end of error message----------------

Thank you,

Zhaojie Zhang

-----Original Message-----
From: emersonlucena [mailto:[hidden email]]
Sent: Saturday, November 24, 2007 4:45 PM
To: [hidden email]
Subject: Re: Stereo SEM image generation

Adam

Copy the code int to plugin folder  of ImageJ and compile.
When opening the Reconstruction 3D_versao5  plugin
1-Choose Reconstruction
2-Angle: enter angle value (degrees).
   Unit: scale of stereo pair.
   Searching-x: displacement used during matching process
   Stop Searching: use different values between 0 and 1. There is no  rule in in this item.
3-Enter matrix dimension will be used to obtain elevation map and ok.
    You will have a elevation map of surface. Use Surface Plot of ImageJ to obtain 3D image or Plot Profile to analyse.

Lucena


Reconstruction 3D_versao5   Code
*****************************************************************************************************************************************
import ij.*;
import ij.plugin.*;
import ij.plugin.filter.*;//PlugInFilter
import ij.process.*;
import ij.gui.*;
import java.awt.*;
import java.util.*;
import ij.measure.*;

public class Reconstruction_3Dversao5 implements  PlugIn {

private ImagePlus imp;
private float[] C;
private float[] S;
private float[] fht;
private static final int TABLE=0,RECONSTRUCTION=1;
private static String[] showStrings = {"Table","Reconstruction"};
private static int option;
ResultsTable rt = new ResultsTable();


int displax=0;
int disx=0;
double matching=0;
boolean yesTable=false;

public void run(String arg) {

IJ.showMessage("About Reconstruction_3D",
     "This plug-in make elevation maps from stereo pair\n" +
     "The  left projection must be open first\n");

ImagePlus imp= WindowManager.getCurrentImage();
int[] wList = WindowManager.getIDList();

GenericDialog gd1=new GenericDialog ("Stereo_Reconstruction");
gd1.addChoice("Options:", showStrings, showStrings[option]);
gd1.showDialog();
if (gd1.wasCanceled()) return ;
option = gd1.getNextChoiceIndex();
Calibration cal = imp.getCalibration();
float scalex=(float)cal.pixelWidth;
float scaley=(float)cal.pixelHeight;
 

if (option==1){
ImageProcessor ip = imp.getProcessor();
ImagePlus img1 = WindowManager.getImage(wList[0]);
ImagePlus img2 = WindowManager.getImage(wList[1]);
ImageProcessor ip3;
ImageProcessor ip1 = img1.getProcessor();

int width = ip1.getWidth();
int height = ip1.getHeight();
ImageStack stack3 = img1.getStack();
ImageStack stack4 = new ImageStack(width, height);
String label;
label = stack3.getSliceLabel(1);
ip3 = ip1.convertToFloat();
stack4.addSlice(label, ip3);
img1.setStack(null, stack4);
ImageProcessor ip4;
ImageProcessor ip2 = img2.getProcessor();

int width2 = ip2.getWidth();
int height2 = ip2.getHeight();
ImageStack stack1 = img2.getStack();
ImageStack stack2 = new ImageStack(width2, height2);
label = stack1.getSliceLabel(1);
ip4 = ip2.convertToFloat();
stack2.addSlice(label, ip4);
img2.setStack(null, stack2);

GenericDialog gd=new GenericDialog ("Parameters of Stereo_Reconstruction");
  gd.addNumericField("Angle:",0.0,3);
  gd.addStringField("Unit:","µm");
  gd.addNumericField("Searching-x:",(30),3);
  gd.addNumericField("Stop Searching:",0.07,3);
  gd.showDialog();
    if (gd.wasCanceled())return ;
    short theta=(short) gd.getNextNumber();
    String unit=gd.getNextString();
    int matchx=(int)(gd.getNextNumber());
    float trunca= (float)gd.getNextNumber();
    int maxN;
    float[] pixels1 =(float[])ip3.getPixels();
    float[]pixels2 =(float[])ip4.getPixels();
 
String sizes = "2,4,8,16,32,64";
int[] boxSizes;
String s = IJ.getString("Dimensions used:", sizes);
  if (s.equals(""))return;
  boxSizes = s2ints(s);
  if (boxSizes==null || boxSizes.length<1)return;
float dimensions=boxSizes.length;
for (int i=0; i<dimensions; i++){
  maxN=boxSizes[i];
  String title3=("Working with kernel of "+maxN+"x"+maxN);
  IJ.showStatus(title3);
  reconstruction( boxSizes,i,dimensions , cal, imp,ip2, unit,theta,maxN, pixels1, pixels2, width, height,matchx,trunca,scalex,scaley);
}
GenericDialog gd4 = new GenericDialog("More Operations?");
gd4.addCheckbox("Table?", false);
gd4.showDialog();
if (gd4.wasCanceled()) return;
yesTable=gd4.getNextBoolean();
if(yesTable)option=0;
}

if (option==0){
ImageProcessor ip3;
if (imp==null){
  IJ.error("One grayscale image required");
  return;
}
if(yesTable){
  int[] wList1 = WindowManager.getIDList();
  ImagePlus img3 = WindowManager.getImage(wList1[2]);
  ip3 = img3.getProcessor();
}
else{
  ImagePlus imp3 = WindowManager.getCurrentImage();
  ip3 = imp3.getProcessor();
}
int width = ip3.getWidth();
int height = ip3.getHeight();
showtable(cal, ip3, imp,width,height,scalex,scaley);//***********************************************linha 132
}}
 
public int[] s2ints(String s) {
  StringTokenizer st = new StringTokenizer(s, ", \t");
  int nInts = st.countTokens();
  int[] ints = new int[nInts];
  for(int i=0; i<nInts; i++) {
    try {ints[i] = Integer.parseInt(st.nextToken());}
    catch (NumberFormatException e) {IJ.write(""+e); return null;}
    }
    return ints;
}

public void reconstruction(int[] boxSizes,int i,float dimensions,Calibration cal,ImagePlus imp,ImageProcessor ip2,String unit,short theta,int maxN, float[] pixels1,float[] pixels2, int width, int height, int matchx, float trunca, float scalex, float scaley){
 
  int slice=0;
  int width1 = 0;
  int height1 = 0;
  double f[]=new double [(maxN*maxN)];
  double g[]=new double [(maxN*maxN)];
  float H[]=new float [(width-matchx-(maxN/2))*(height-maxN)];
  float xx[]=new float [(width-maxN)];
  float yy[]=new float [(height-maxN)];

  ImageWindow win = imp.getWindow();
  win.running = true;
  long startTime = System.currentTimeMillis();
  IJ.showStatus("Building the elevation map...");

  for (int y=(maxN/2)  ; y <height-(maxN/2); y++){      
     IJ.showProgress((double)y/(height-(maxN/2)));
     if(!win.running)break;
     for (int x=(maxN/2)  ; x <width-matchx-(maxN/2); x++){
     
       for (int row=y-(maxN/2); row <y+(maxN/2); row++){
        for (int col=x-(maxN/2); col <x+(maxN/2); col++){
          f[(col-x+(maxN/2)+((row-y+(maxN/2))*maxN))]=pixels1[(col+(row*width))];
        }
       }
     double [] FFHT=fft(ip2,f,maxN);
     for (int searchx=x  ; searchx <x+(matchx); searchx++){
      for (int rowg=y-(maxN/2); rowg <y+(maxN/2); rowg++){
        for (int colg=searchx-(maxN/2); colg <searchx+(maxN/2); colg++){
         g[((colg-searchx+(maxN/2))+((rowg-y+(maxN/2))*maxN))]=pixels2[(colg+(rowg*width))];
        }        
       }
      double [] GFHT=fft(ip2,g,maxN);
      matching=(double)Math.abs((GFHT[0]/FFHT[0])-1);

//***************************
//slice++;
// rt.incrementCounter();
// rt.addValue("matching", matching);
// rt.addValue("(GFHT[0]/FFHT[0])",(GFHT[0]/FFHT[0]));
// if (slice==1)IJ.setColumnHeadings(rt.getColumnHeadings());
// IJ.write(rt.getRowAsString(rt.getCounter()-1));
//****************************


//Check trunca values
      if(matching==trunca) displax=searchx;//******************************************************************* line 242
     }//searchx end
   H[(x-(maxN/2))+(y-(maxN/2))*(width-matchx-(maxN/2))]=displacement(imp, displax, theta,x,y,scalex);//,display-(matchy/2)
  }
        //********************************************************************line 246
}

String title=("Depth Image with Kernel of "+maxN+"x"+maxN);
new ImagePlus(title, new FloatProcessor((width-matchx-(maxN/2)),(height-maxN),H, null)).show();
ImageProcessor ipp1;
ImageProcessor ipp2;
ImageProcessor ipp3;

int[] wListp = WindowManager.getIDList();
if (wListp.length>3){
 ImagePlus impp1 = WindowManager.getImage(wListp[wListp.length-1]);
 ImagePlus impp2 = WindowManager.getImage(wListp[wListp.length-2]);
 ipp1=impp1.getProcessor();
 ipp2=impp2.getProcessor();

int width2 = ipp2.getWidth();
int height2 = ipp2.getHeight();
width1 = ipp1.getWidth();
height1 = ipp1.getHeight();
float[] H1 =(float[])ipp1.getPixels();
float[] H2 =(float[])ipp2.getPixels();

for(int x = 0; x <width1; x++){
 for(int y = 0; y <height1 ; y++){
  if (i==(dimensions-1) ){
   H1[x+y*width1]=H1[x+y*width1]+H2[x+((boxSizes[i]/2)-(boxSizes[i-1]/2))+(y+((boxSizes[i]/2)-(boxSizes[i-1]/2)))*width2];
   H1[x+y*width1]=H1[x+y*width1]/dimensions;
  }
else
   H1[x+y*width1]=H1[x+y*width1]+H2[x+((boxSizes[i]/2)-(boxSizes[i-1]/2))+(y+((boxSizes[i]/2)-(boxSizes[i-1]/2)))*width2];
   H1[x+y*width1]=H1[x+y*width1];
}}// x and y end

String title1=("Result of add "+(boxSizes[i])+"x"+(boxSizes[i])+" with "+(boxSizes[i-1])+"x"+(boxSizes[i-1]));
new ImagePlus(title1, new FloatProcessor(width1,height1,H1, null)).show();
}

int[] wListp1 = WindowManager.getIDList();
if (i==(dimensions-1) ){
 long time = System.currentTimeMillis()-startTime;
 double seconds = time/1000.0;
 IJ.write("Reconstruction Time:"+IJ.d2s(seconds)+"seconds");
}
}//reconstruction end

public void showtable (Calibration cal,ImageProcessor ip3, ImagePlus imp, int width, int height, float scalex, float scaley){

 int slice=0;
 float xxx[]=new float [(width)];
 float yyy[]=new float [(height)];
 float[] E =(float[])ip3.getPixels();
 for(int x = 0; x <width; x++){
  for(int y = 0; y <height ; y++){
   xxx[x]=x*scalex;
   yyy[y]=y*scaley;
 slice++;
 rt.incrementCounter();
 rt.addValue("x", xxx[x]);
 rt.addValue("y", yyy[y]);
 rt.addValue("Elevation",E[x+y*(width)] );
 if (slice==1)IJ.setColumnHeadings(rt.getColumnHeadings());
 IJ.write(rt.getRowAsString(rt.getCounter()-1));
}}
}//show table end

public  float  displacement( ImagePlus imp,float displax, short theta, int x,int y, float scalex) {

   float dis=0;
   float PI = (float)( Math.PI);
   float delta= (float)(((theta*PI)/180));
   float max1=(float)(scalex*(displax-x));
   dis=(float) (max1/(2*Math.sin(delta/2.0)));
   return dis;
  }

/**
This file contains a Java language implementation of the
Fast Hartley Transform algorithm which is covered under
United States Patent Number 4,646,256.

This code may therefore be freely used and distributed only
under the following conditions:

 1)  This header is included in every copy of the code; and
 2)  The code is used for noncommercial research purposes only.

 Firms using this code for commercial purposes may be infringing a United
 States patent and should contact the

 Office of Technology Licensing
 Stanford University
 857 Serra Street, 2nd Floor
 Stanford, CA   94305-6225
 (415) 723 0651

This implementation is based on Pascal
code contibuted by Arlo Reeves.
*/

/** Fast Hartley Transform. */

 public double[] fft(ImageProcessor ip,double[] pixels,int maxN) {

  makeSinCosTables(maxN);
  double[] fht = pixels;
  rc2DFHT(fht, maxN);
  return fht;
 }

 void makeSinCosTables(int maxN) {

  int n = maxN/4;
  C = new float[n];
  S = new float[n];
  double theta = 0.0;
  double dTheta = 2.0 * Math.PI/maxN;
  for (int i=0; i<n; i++) {
   C[i] = (float)Math.cos(theta);
   S[i] = (float)Math.sin(theta);
   theta += dTheta;
  }
 }

 /** Row-column Fast Hartley Transform */
 void rc2DFHT(double[] x,  int maxN) {
  //IJ.write("FFT: rc2DFHT (row-column Fast Hartley Transform)");
  for (int row=0; row<maxN; row++)
   dfht3(x, row*maxN, maxN);
  transposeR(x, maxN);
  for (int row=0; row<1; row++)//maxN******************
   dfht3(x, row*maxN,  maxN);
  //transposeR(x, maxN);*************************

  int mRow, mCol;
  double A,B,C,D,E;
  for (int row=0; row<1; row++) { //*********maxN/2 Now calculate actual Hartley transform
   for (int col=0; col<1; col++) {//************maxN/2
    mRow = (maxN - row) % maxN;
    mCol = (maxN - col)  % maxN;
    A = x[row * maxN + col]; //  see Bracewell, 'Fast 2D Hartley Transf.' IEEE Procs. 9/86
    B = x[mRow * maxN + col];
    C = x[row * maxN + mCol];
    D = x[mRow * maxN + mCol];
    E = ((A + D) - (B + C)) / 2;
    x[row * maxN + col] = A - E;
    x[mRow * maxN + col] = B + E;
    x[row * maxN + mCol] = C + E;
    x[mRow * maxN + mCol] = D - E;
   }
  }
 }

 /* An optimized real FHT */
 void dfht3 (double[] x, int base,  int maxN) {
  int i, stage, gpNum, gpIndex, gpSize, numGps, Nlog2;
  int bfNum, numBfs;
  int Ad0, Ad1, Ad2, Ad3, Ad4, CSAd;
  double rt1, rt2, rt3, rt4;

  Nlog2 = log2(maxN);
  BitRevRArr(x, base, Nlog2, maxN); //bitReverse the input array
  gpSize = 2;     //first & second stages - do radix 4 butterflies once thru
  numGps = maxN / 4;
  for (gpNum=0; gpNum<numGps; gpNum++)  {
   Ad1 = gpNum * 4;
   Ad2 = Ad1 + 1;
   Ad3 = Ad1 + gpSize;
   Ad4 = Ad2 + gpSize;
   rt1 = x[base+Ad1] + x[base+Ad2];   // a + b
   rt2 = x[base+Ad1] - x[base+Ad2];   // a - b
   rt3 = x[base+Ad3] + x[base+Ad4];   // c + d
   rt4 = x[base+Ad3] - x[base+Ad4];   // c - d
   x[base+Ad1] = rt1 + rt3;      // a + b + (c + d)
   x[base+Ad2] = rt2 + rt4;      // a - b + (c - d)
   x[base+Ad3] = rt1 - rt3;      // a + b - (c + d)
   x[base+Ad4] = rt2 - rt4;      // a - b - (c - d)
   }

  if (Nlog2 > 2) {
    // third + stages computed here
   gpSize = 4;
   numBfs = 2;
   numGps = numGps / 2;
   //IJ.write("FFT: dfht3 "+Nlog2+" "+numGps+" "+numBfs);
   for (stage=2; stage<Nlog2; stage++) {
    for (gpNum=0; gpNum<numGps; gpNum++) {
     Ad0 = gpNum * gpSize * 2;
     Ad1 = Ad0;     // 1st butterfly is different from others - no mults needed
     Ad2 = Ad1 + gpSize;
     Ad3 = Ad1 + gpSize / 2;
     Ad4 = Ad3 + gpSize;
     rt1 = x[base+Ad1];
     x[base+Ad1] = x[base+Ad1] + x[base+Ad2];
     x[base+Ad2] = rt1 - x[base+Ad2];
     rt1 = x[base+Ad3];
     x[base+Ad3] = x[base+Ad3] + x[base+Ad4];
     x[base+Ad4] = rt1 - x[base+Ad4];
     for (bfNum=1; bfNum<numBfs; bfNum++) {
     // subsequent BF's dealt with together
      Ad1 = bfNum + Ad0;
      Ad2 = Ad1 + gpSize;
      Ad3 = gpSize - bfNum + Ad0;
      Ad4 = Ad3 + gpSize;

      CSAd = bfNum * numGps;
      rt1 = x[base+Ad2] * C[CSAd] + x[base+Ad4] * S[CSAd];
      rt2 = x[base+Ad4] * C[CSAd] - x[base+Ad2] * S[CSAd];

      x[base+Ad2] = x[base+Ad1] - rt1;
      x[base+Ad1] = x[base+Ad1] + rt1;
      x[base+Ad4] = x[base+Ad3] + rt2;
      x[base+Ad3] = x[base+Ad3] - rt2;

     } /* end bfNum loop */
    } /* end gpNum loop */
    gpSize *= 2;
    numBfs *= 2;
    numGps = numGps / 2;
   } /* end for all stages */
  } /* end if Nlog2 > 2 */
 }

 void transposeR (double[] x, int maxN) {
  int   r, c;
  double  rTemp;

  for (r=0; r<maxN; r++)  {
   for (c=r; c<maxN; c++) {
    if (r != c)  {
     rTemp = x[r*maxN + c];
     x[r*maxN + c] = x[c*maxN + r];
     x[c*maxN + r] = rTemp;
    }
   }
  }
 }

 int log2 (int x) {
  int count = 15;
  while (!btst(x, count))
   count--;
  return count;
 }


 private boolean btst (int  x, int bit) {
  //int mask = 1;
  return ((x & (1<<bit)) != 0);
 }

 void BitRevRArr (double[] x, int base, int bitlen, int maxN) {
  int    l;
  double[] tempArr = new double[maxN];
  for (int i=0; i<maxN; i++)  {
   l = BitRevX (i, bitlen);  //i=1, l=32767, bitlen=15
   tempArr[i] = x[base+l];
  }
  for (int i=0; i<maxN; i++)
   x[base+i] = tempArr[i];
 }

 private int BitRevX (int  x, int bitlen) {
  int  temp = 0;
  for (int i=0; i<=bitlen; i++)
   if ((x & (1<<i)) !=0)
    temp  |= (1<<(bitlen-i-1));
  return temp & 0x0000ffff;
 }

 private int bset (int x, int bit) {
  x |= (1<<bit);
  return x;
 }
//End of routine
}

*****************************************************************************************************************


  ----- Original Message -----
  From: Adam Hacking, Dr.
  To: [hidden email]
  Sent: Saturday, November 24, 2007 2:24 PM


  Lucena,

  I would be interested in trying it. Can you provide some more information ? Will the model produce a 3D image (similar to an AFM image) suitable for analysis ?

  Adam

  emersonlucena <[hidden email]> wrote: Adam

  I development a plugin for 3D models from stereo SEM ideal images (one image
  displacement 10 pixels to a copy of the same image). The aim of this
  research is determined only influence of imaging conditions system on
  reconstruction quality. Anyway, you can try use it if your real pair stereo
  have better epipolar conditions (decrease stage rotated angle below +-1
  degrees and image resolution above 640x480 pixels) than we work here. Are
  you interested?
  And please, sorry my english. I need improved it too.

  Cheers

  Lucena


  ________________________________________________________________ Dr. S Adam Hacking, Post Doctoral Fellow
        JTN Wong Laboratories for Mineralized Tissue Research,
        Center for Bone and Periodontal Research,
        McGill University   740 Dr. Penfield Ave. Rm. 2300A
        Montreal, QC, Canada H3A 1A4
        Ph.:  514-398-5112
        Fax: 514-398-4020



------------------------------------------------------------------------------


  No virus found in this incoming message.
  Checked by AVG Free Edition.
  Version: 7.5.503 / Virus Database: 269.16.5/1149 - Release Date: 24/11/2007 10:06