// @Auther : M.Hattori &丂H.Kunitomo 2008.9.23
// Dicom夋憸偺Multi_Frequency_Process(MFP)張棟
// 掇臣陛台倮傪巊偭偨懡廳夝憸搙暘夝偱偺MFP張棟


rename("Org_Image");
org=getImageID();

 w=getWidth() ;
 h=getHeight() ;

 f1=getNumber("G1:",1.0);
 f2=getNumber("G2:",0.8);
 f3=getNumber("G3:",0.6);
 f4=getNumber("G4:",0.4);
 f5=getNumber("G5:",0.2);

 
 xm=getNumber("Xm:",160);
 ym=getNumber("Ym:",1);

 run("Duplicate...","title=G3_1_Image") ;
 run("Gaussian Blur...", "sigma=1");
 imageCalculator("Subtract create", "Org_Image","G3_1_Image");
//run("Image Calculator...", "image1=Org_Image operation=Subtract image2=G3_1_Image creat");
run("Multiply...","value="+f1) ;
rename("Sus_ori_Ga1");

 selectImage("G3_1_Image");
 run("Duplicate...","title=G3_2_Image") ;
 Wt=round(w/2);
 He=round(h/2);
 run("Size...", "width="+Wt+" height="+He+" constrain");
 run("Gaussian Blur...", "sigma=1");
 run("Duplicate...","title=G3_3_Image") ;
 selectImage("G3_2_Image");
 run("Size...", "width="+w+" height="+h+" interpolate");
 imageCalculator("Subtract create", "G3_1_Image","G3_2_Image");
//run("Image Calculator...", "image1=G3_1_Image operation=Subtract image2=G3_2_Image creat");
 run("Multiply...","value="+f2) ;
 rename("Sus_Ga1_Ga2");

 selectImage("G3_3_Image");
 Wt=round(w/4);
 He=round(h/4);
 run("Size...", "width="+Wt+" height="+He+" constrain");
 run("Gaussian Blur...", "sigma=1");
 run("Duplicate...","title=G3_4_Image") ;
 selectImage("G3_3_Image");
 run("Size...", "width="+w+" height="+h+" interpolate");
 imageCalculator("Subtract create", "G3_2_Image","G3_3_Image");
//run("Image Calculator...", "image1=G3_2_Image operation=Subtract image2=G3_3_Image creat");
 run("Multiply...","value="+f3) ;
 rename("Sus_Ga2_Ga3");

 selectImage("G3_4_Image");
 Wt=round(w/8);
 He=round(h/8);
 run("Size...", "width="+Wt+" height="+He+" constrain");
 run("Gaussian Blur...", "sigma=1");
 run("Duplicate...","title=G3_5_Image") ;
 selectImage("G3_4_Image");
 run("Size...", "width="+w+" height="+h+" interpolate");
 imageCalculator("Subtract create", "G3_3_Image","G3_4_Image");
//run("Image Calculator...", "image1=G3_3_Image operation=Subtract image2=G3_4_Image creat");
 run("Multiply...","value="+f4) ;
 rename("Sus_Ga3_Ga4");
 
 selectImage("G3_5_Image");
 Wt=round(w/16);
 He=round(h/16);
 run("Size...", "width="+Wt+" height="+He+" constrain");
 run("Gaussian Blur...", "sigma=1");
 selectImage("G3_5_Image");
 run("Size...", "width="+w+" height="+h+" interpolate");
 imageCalculator("Subtract create", "G3_4_Image","G3_5_Image");
//run("Image Calculator...", "image1=G3_4_Image operation=Subtract image2=G3_5_Image creat");
 run("Multiply...","value="+f5) ;
 rename("Sus_Ga4_Ga5");



 imageCalculator("Add create", "Sus_ori_Ga1","Sus_Ga1_Ga2");
//run("Image Calculator...", "image1=Sus_ori_Ga1 operation=Add image2=[Sus_Ga1_Ga2] create");

 imageCalculator("Add create", "Sus_Ga2_Ga3","Result of Sus_ori_Ga1");
//run("Image Calculator...", "image1=Sus_Ga2_Ga3 operation=Add image2=[Result of Sus_Ga1_Ga2] create");

 imageCalculator("Add create", "Sus_Ga3_Ga4","Result of Sus_Ga2_Ga3");
//run("Image Calculator...", "image1=Sus_Ga3_Ga4 operation=Add image2=[Result of Sus_Ga2_Ga3] create");

 imageCalculator("Add create", "Sus_Ga4_Ga5","Result of Sus_Ga3_Ga4");
//run("Image Calculator...", "image1=Sus_Ga4_Ga5 operation=Add image2=[Result of Sus_Ga3_Ga4] create");

 rename("Sus_Image");

 for (y=0;y<h;y++){
     for(x=0;x<w;x++){
         sm=getPixel(x,y) ;
         if (sm==0) { pv=0; }
               else { pv=ym*abs(sm)*((exp(xm/sm)-1) / (exp(xm/sm)+1)) ; }
         if (isNaN(pv)) { pv=sm ; }
         setPixel(x,y,pv);
     }
  } 


 imageCalculator("Add create", "Org_Image","Sus_Image");
//run("Image Calculator...", "image1=Org_Image operation=Add image2=[Sus_Image] create");


 rename("MFP_Xm,Ym"+xm+","+ym);
