Stereo SEM image generation

[Adam Hacking]

Dear list,

Stereo SEM can be a powerful, accessible and high resolution tool for analysis of surface morphology.

Is anyone aware of any ImageJ plugins that can generate 3D models from stereo SEM images ? Is there any interest in development of such a plugin ?

Thanks

Adam Hacking


________________________________________________________________         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: [phone]

[emersonlucena]

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


----- Original Message -----
From: "Adam Hacking, Dr." <[hidden email]>
To: <[hidden email]>
Sent: Thursday, November 22, 2007 2:31 AM
Subject: Stereo SEM image generation

> Dear list,
>
> Stereo SEM can be a powerful, accessible and high resolution tool for
> analysis of surface morphology.
>
> Is anyone aware of any ImageJ plugins that can generate 3D models from
> stereo SEM images ? Is there any interest in development of such a plugin
> ?
>
> Thanks
>
> Adam Hacking
>
>
> ________________________________________________________________
> 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: [phone]
>
>
>
> --
> No virus found in this incoming message.
> Checked by AVG Free Edition.
> Version: 7.5.503 / Virus Database: 269.16.4/1146 - Release Date:
> 22/11/2007 18:55
>
>

[Orkun Ersoy]

Dear Adam,
Gary Chinga developed such a plugin (StereoJ)...
You can calculate height maps from stereo pairs.
You can find it in his website...
http://www.gcsca.net/IJ/StereoJ.html

Dr.Orkun ERSOY

Hacettepe University

Department of Geological Engineering

Beytepe Ankara Turkey

[phone] / 126


-----Original Message-----
From: ImageJ Interest Group [mailto:[hidden email]] On Behalf Of
Adam Hacking, Dr.
Sent: Thursday, November 22, 2007 6:32 AM
To: [hidden email]
Subject: Stereo SEM image generation

Dear list,

Stereo SEM can be a powerful, accessible and high resolution tool for
analysis of surface morphology.

Is anyone aware of any ImageJ plugins that can generate 3D models from
stereo SEM images ? Is there any interest in development of such a
plugin ?

Thanks

Adam Hacking


________________________________________________________________
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: [phone]

[Gary Chinga]

Dear all...

This is for information to all of you that have been interested in  
the StereoJ plugin. I developed a beta version of the StereoJ plugin  
last year. However, I was not satisfied with the performance of the  
plugin and decided not to release it yet. Due to serious time  
problems during the last months I have not had time to work with it  
until recently. Some bugs have also been fixed. A new version will be  
released, most probably in February/March next year, as I will be  
taking vacations in Chile  from December 3 to January 29.

Regards,

Gary.
http://www.gcsca.net




On Nov 24, 2007, at 12:46 PM, Orkun Ersoy wrote:

> Dear Adam,
> Gary Chinga developed such a plugin (StereoJ)...
> You can calculate height maps from stereo pairs.
> You can find it in his website...
> http://www.gcsca.net/IJ/StereoJ.html
>
> Dr.Orkun ERSOY
>
> Hacettepe University
>
> Department of Geological Engineering
>
> Beytepe Ankara Turkey
>
> [phone] / 126
>
>
> -----Original Message-----
> From: ImageJ Interest Group [mailto:[hidden email]] On Behalf Of
> Adam Hacking, Dr.
> Sent: Thursday, November 22, 2007 6:32 AM
> To: [hidden email]
> Subject: Stereo SEM image generation
>
> Dear list,
>
> Stereo SEM can be a powerful, accessible and high resolution tool for
> analysis of surface morphology.
>
> Is anyone aware of any ImageJ plugins that can generate 3D models from
> stereo SEM images ? Is there any interest in development of such a
> plugin ?
>
> Thanks
>
> Adam Hacking
>
>
> ________________________________________________________________
> 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: [phone]
>
>

[emersonlucena]

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
 [phone]

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: [phone]



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


  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

[emersonlucena]

Noel

This is email with Stereo plugin code.  The processing  is very fast due
matching process to be made with unidimensional transform and only one point
of frequency domain. So the idea is improve epipolar conditions of stereo
pair  we work here (decrease stage rotated angle below +-1 degrees and image
resolution above 640x480 pixels) to apply this plugin.
As stability of plugin, this is made using ideal pair stereo (one image
displacement few pixels to a copy of the same image). The elevation map must
be a plan. This technique is describe in  Hein, LRO, Quantitative
fractography by digital image processing: NIH Image macro tools for stereo
pair analysis and 3-D reconstruction, Journal of Microscopy, Vol. 204, pp
17-28, (2001).


Regards


Lucena

----- Original Message -----
From: "emersonlucena" <[hidden email]>
To: <[hidden email]>
Sent: Saturday, November 24, 2007 9:45 PM
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
> [phone]
>
> 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: [phone]
>
>
>
> ------------------------------------------------------------------------------
>
>
>  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
>
>
>
> --
> No virus found in this incoming message.
> Checked by AVG Free Edition.
> Version: 7.5.503 / Virus Database: 269.16.6/1150 - Release Date:
> 24/11/2007 17:58
>
>

[emersonlucena]

Ersoy, Noel


Hi. Everything ok?
Try this new version. I changed code at line 242 and it obtained better results on pair results from Adam ([hidden email]).

if(matching==trunca) displax=searchx;       // older version.
if(matching<=trunca) displax=searchx;       // new one.

Any questions, please contact me.


Lucena

[Jarek Grodek]

I apologise for my mistake, here is the text of the macro:

run("Set Scale...", "distance=21.4 known=1 pixel=1 unit=mm global");
run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00
show=Masks exclude include add");
rename("mask");
run("Erode");
run("Erode");
run("Erode");
run("Erode");
run("Erode");
run("Invert LUT");
run("Invert");
imageCalculator("Add", "ori.tif","mask");
//run("Image Calculator...", "image1=ori.tif operation=Add image2=mask");
selectWindow("mask");
close;
selectWindow("ori.tif");
run("Select None");


roiManager("Select", 0);print("1");run("Clear Results");
run("Set Measurements...", "area redirect=None decimal=3");
run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00
show=Nothing display exclude");
getStatistics(area);
        if (nResults==0) {
                print("sound kernel");
        } else {
                print("infested");
        }
        if (area>=0.050 && area<=0.051)
                print("egg");
        else if (area>=0.052 && area<=1.50)
                print("1st instar");
        else if (area>=1.51 && area<=3.50)
                print("2nd instar");
        else if (area>=3.51 && area<=5.50)
                print("3rd instar");
        else if (area>=5.51 && area<=10.50)
                print("4th instar");
        else if (area>=10.51 && area<=10.55)
                print("prepupa");
        else if (area>=10.56 && area<=10.59)
                print("pupa");


roiManager("Select", 1);print("2");run("Clear Results");
run("Set Measurements...", "area redirect=None decimal=3");
run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00
show=Nothing display exclude");
getStatistics(area);
        if (nResults==0) {
                print("sound kernel");
        } else {
                print("infested");
        }
i=nResults;
for (i=0; i<nResults; i++) {
getResult("Area",i)
        if (i++>=0.050 && i++<=0.051)
                print("egg");
        else if (i++>=0.052 && i++<=1.50)
                print("1st instar");
        else if (i++>=1.51 && i++<=3.50)
                print("2nd instar");
        else if (i++>=3.51 && i++<=5.50)
                print("3rd instar");
        else if (i++>=5.51 && i++<=10.50)
                print("4th instar");
        else if (i++>=10.51 && i++<=10.55)
                print("prepupa");
        else if (i++>=10.56 && i++<=10.59)
                print("pupa");

}



roiManager("Select", 2);print("3");run("Clear Results");
run("Set Measurements...", "area redirect=None decimal=3");
run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00
show=Nothing display exclude");
getStatistics(area);
        if (nResults==0) {
                print("sound kernel");
        } else {
                print("infested");
        }
        if (area>=0.050 && area<=0.051)
                print("egg");
        else if (area>=0.052 && area<=1.50)
                print("1st instar");
        else if (area>=1.51 && area<=3.50)
                print("2nd instar");
        else if (area>=3.51 && area<=5.50)
                print("3rd instar");
        else if (area>=5.51 && area<=10.50)
                print("4th instar");
        else if (area>=10.51 && area<=10.55)
                print("prepupa");
        else if (area>=10.56 && area<=10.59)
                print("pupa");

roiManager("Select", 3);print("4");run("Clear Results");
run("Set Measurements...", "area redirect=None decimal=3");
run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00
show=Nothing display exclude");
getStatistics(area);
        if (nResults==0) {
                print("sound kernel");
        } else {
                print("infested");
        }
        if (area>=0.050 && area<=0.051)
                print("egg");
        else if (area>=0.052 && area<=1.50)
                print("1st instar");
        else if (area>=1.51 && area<=3.50)
                print("2nd instar");
        else if (area>=3.51 && area<=5.50)
                print("3rd instar");
        else if (area>=5.51 && area<=10.50)
                print("4th instar");
        else if (area>=10.51 && area<=10.55)
                print("prepupa");
        else if (area>=10.56 && area<=10.59)
                print("pupa");

roiManager("Select", 4);print("5");run("Clear Results");
run("Set Measurements...", "area redirect=None decimal=3");
run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00
show=Nothing display exclude");
getStatistics(area);
        if (nResults==0) {
                print("sound kernel");
        } else {
                print("infested");
        }
getStatistics(area);
        if (area>=0.050 && area<=0.051)
                print("egg");
        else if (area>=0.052 && area<=1.50)
                print("1st instar");
        else if (area>=1.51 && area<=3.50)
                print("2nd instar");
        else if (area>=3.51 && area<=5.50)
                print("3rd instar");
        else if (area>=5.51 && area<=10.50)
                print("4th instar");
        else if (area>=10.51 && area<=10.55)
                print("prepupa");
        else if (area>=10.56 && area<=10.59)
                print("pupa");

run("Select None");



thanks i advance for help
Jarek Grodek

[Jarek Grodek]

once again i feel sorry, my firts message contained image for analysis, in
tiff format and NIH server rejected it.

i cite it : <<Dear IJ Users!

I would like to ask you for help with macro designing.
what i intend to do is to create macro wich will show number of objects
with internal damage and will print which object(grain) is
damaged(infested) in which extent.

below is presented what i created and artificial image is attached to this
e-mail.
hope there is someone who may help me to go one step forward?>>


if anyone could help contact me on [hidden email] and i will reply
with image for analysis.

thanks in advance
Jarek Grodek
On Pn Grudnia 3 2007, 9:50, Jarosław Grodek napisał(a):

> I apologise for my mistake, here is the text of the macro:
>
> run("Set Scale...", "distance=21.4 known=1 pixel=1 unit=mm global");
> run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00
> show=Masks exclude include add");
> rename("mask");
> run("Erode");
> run("Erode");
> run("Erode");
> run("Erode");
> run("Erode");
> run("Invert LUT");
> run("Invert");
> imageCalculator("Add", "ori.tif","mask");
> //run("Image Calculator...", "image1=ori.tif operation=Add image2=mask");
> selectWindow("mask");
> close;
> selectWindow("ori.tif");
> run("Select None");
>
>
> roiManager("Select", 0);print("1");run("Clear Results");
> run("Set Measurements...", "area redirect=None decimal=3");
> run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00
> show=Nothing display exclude");
> getStatistics(area);
> if (nResults==0) {
> print("sound kernel");
> } else {
> print("infested");
> }
> if (area>=0.050 && area<=0.051)
> print("egg");
> else if (area>=0.052 && area<=1.50)
> print("1st instar");
> else if (area>=1.51 && area<=3.50)
> print("2nd instar");
> else if (area>=3.51 && area<=5.50)
> print("3rd instar");
> else if (area>=5.51 && area<=10.50)
> print("4th instar");
> else if (area>=10.51 && area<=10.55)
> print("prepupa");
> else if (area>=10.56 && area<=10.59)
> print("pupa");
>
>
> roiManager("Select", 1);print("2");run("Clear Results");
> run("Set Measurements...", "area redirect=None decimal=3");
> run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00
> show=Nothing display exclude");
> getStatistics(area);
> if (nResults==0) {
> print("sound kernel");
> } else {
> print("infested");
> }
> i=nResults;
> for (i=0; i<nResults; i++) {
> getResult("Area",i)
> if (i++>=0.050 && i++<=0.051)
> print("egg");
> else if (i++>=0.052 && i++<=1.50)
> print("1st instar");
> else if (i++>=1.51 && i++<=3.50)
> print("2nd instar");
> else if (i++>=3.51 && i++<=5.50)
> print("3rd instar");
> else if (i++>=5.51 && i++<=10.50)
> print("4th instar");
> else if (i++>=10.51 && i++<=10.55)
> print("prepupa");
> else if (i++>=10.56 && i++<=10.59)
> print("pupa");
>
> }
>
>
>
> roiManager("Select", 2);print("3");run("Clear Results");
> run("Set Measurements...", "area redirect=None decimal=3");
> run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00
> show=Nothing display exclude");
> getStatistics(area);
> if (nResults==0) {
> print("sound kernel");
> } else {
> print("infested");
> }
> if (area>=0.050 && area<=0.051)
> print("egg");
> else if (area>=0.052 && area<=1.50)
> print("1st instar");
> else if (area>=1.51 && area<=3.50)
> print("2nd instar");
> else if (area>=3.51 && area<=5.50)
> print("3rd instar");
> else if (area>=5.51 && area<=10.50)
> print("4th instar");
> else if (area>=10.51 && area<=10.55)
> print("prepupa");
> else if (area>=10.56 && area<=10.59)
> print("pupa");
>
> roiManager("Select", 3);print("4");run("Clear Results");
> run("Set Measurements...", "area redirect=None decimal=3");
> run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00
> show=Nothing display exclude");
> getStatistics(area);
> if (nResults==0) {
> print("sound kernel");
> } else {
> print("infested");
> }
> if (area>=0.050 && area<=0.051)
> print("egg");
> else if (area>=0.052 && area<=1.50)
> print("1st instar");
> else if (area>=1.51 && area<=3.50)
> print("2nd instar");
> else if (area>=3.51 && area<=5.50)
> print("3rd instar");
> else if (area>=5.51 && area<=10.50)
> print("4th instar");
> else if (area>=10.51 && area<=10.55)
> print("prepupa");
> else if (area>=10.56 && area<=10.59)
> print("pupa");
>
> roiManager("Select", 4);print("5");run("Clear Results");
> run("Set Measurements...", "area redirect=None decimal=3");
> run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00
> show=Nothing display exclude");
> getStatistics(area);
> if (nResults==0) {
> print("sound kernel");
> } else {
> print("infested");
> }
> getStatistics(area);
> if (area>=0.050 && area<=0.051)
> print("egg");
> else if (area>=0.052 && area<=1.50)
> print("1st instar");
> else if (area>=1.51 && area<=3.50)
> print("2nd instar");
> else if (area>=3.51 && area<=5.50)
> print("3rd instar");
> else if (area>=5.51 && area<=10.50)
> print("4th instar");
> else if (area>=10.51 && area<=10.55)
> print("prepupa");
> else if (area>=10.56 && area<=10.59)
> print("pupa");
>
> run("Select None");
>
>
>
> thanks i advance for help
> Jarek Grodek
>
>

MSc Eng. Jaroslaw Grodek
Institute of Agrophysics
Polish Academy of Sciences in Lublin
Doswiadczalna 4
20-290 Lublin
Poland

[Z.J. Zhang]

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
 [phone]

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: [phone]



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


  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