/*

   Mip_.java
   
   12/12/01
   28/03/02
   24/05/02    Modifica per un problema di compilazione segnalato da Paola Bonetto
               (linee 309 e 317)
   28/05/02    Aggiunta del vero algoritmo MIP (denominazione dell'algoritmo
               precedente: TIP)

   Questo plugin genera una sequenza di proiezioni laterali (partendo da uno stack 
   di sezioni transassiali, implementando il diffuso algoritmo MIP (Maximal Intensity 
   Projection) o, in alternativa, una sua variazione da noi denominata TIP (Total
   Intensity Projection).

   Lo stack ottenuto puo' essere facilmente trasformato in file AVI (con il plugin
   AVI_Writer), ed eventualmente inserito in una presentazione PowerPoint (c).

   Le operazioni che consentono di ottenere questo sono:

   - verifica che l'immagine corrente sia uno stack;
   - apertura di una dialog box con la quale si richiede il numero di proiezioni;
   - creazione di uno stack vuoto per le proiezioni;

   - per ogni proiezione i-esima da costruire:

     - rotazione dell'intero stack di un angolo pari a i*(360/n.proiezioni)° secondo
       l'asse z;
     - proiezione del solo valore massimo di tutte le slice "frontali" (secondo
       l'asse x) su un'unica slice (MIP), oppure somma e media di tutte le slice 
       "frontali" (secondo l'asse x) su un'unica slice (TIP);
     - assegnazione della slice media come i-esima slice dello stack delle somme;
   
   - assegnazione allo stack delle somme di una tavola di lookup "fire";
   - visualizzione dello stack.

   BIBLIOGRAFIA
   
   - Bailer W, "Writing ImageJ Plugins - A Tutorial";
   - Sun Y, Parker DL, Performance Analysis of Maximum Intensity Projection algorithm for
     display of MRA images, IEEE Trans Med Imaging 1999 Dec; 18(12):1154-69;
   - Mroz L, Real-Time Volume Visualization on Low-End Hardware, PhD Thesis, 
     http://www.cg.tuwien.ac.at/~mroz/diss/html/diss.html;

*/

import ij.*;
import ij.gui.*;
import ij.io.*;
import ij.measure.*;
import ij.process.*;
import ij.plugin.*;
import ij.plugin.frame.*;
import ij.plugin.PlugIn;
import ij.plugin.filter.PlugInFilter;

import java.io.*;
import java.util.*;
import java.lang.String;

import java.awt.*;
import java.awt.image.*; 
import java.awt.event.*;
import java.awt.Color;
import javax.swing.*;

/*****************************************/
public class Mip_ implements PlugInFilter {
/*****************************************/

   static ImagePlus img;
   private static int nProj = 120;
   private double angCurr, angStep;
   private int sectCurr, nSect;

//---------------------------------------------
   public int setup(String arg, ImagePlus imp){
//---------------------------------------------
            
      this.img = imp;

      if (imp == null){
         IJ.showMessage("There are no images open");
         return DONE;
      }

 		if (arg.equals ("about")) {
         showAbout(); 
         return DONE;
      }

//    Si abilita l'elaborazione di stack a 8 o 16 bit ...
      return DOES_8G+DOES_16+STACK_REQUIRED+NO_CHANGES;
   }                                  

//-------------------------------------     
   public void run(ImageProcessor ip) { 
//-------------------------------------     
      
      long[] sumVector;
      ImageProcessor origSliceIp, modifSliceIp, projSliceIp;
      ImageStack origStack = img.getStack();
      int currProj;
      int i, j, k, x, y, z;
      double maxVol, minVol;
      double maxSlice, minSlice;
      double valToSubtract;
      double valDepthCorr = 1.0;
      boolean lDepthCorr = true;
      boolean lTIP = false;

//    Si calcola il minimo e il massimo valore di conteggio presente nell'intero volume ...

      IJ.showStatus("Calculating minimum in volume ...");

      maxVol = Double.MIN_VALUE;
      minVol = Double.MAX_VALUE;
      nSect = img.getStackSize();

      for (z = 1; z <= nSect; ++z) {

         origSliceIp = origStack.getProcessor(z);
         minSlice = origSliceIp.getMin();
         maxSlice = origSliceIp.getMax();

         if (minSlice < minVol)
            minVol = minSlice;

         if (maxSlice > maxVol)
            maxVol = maxSlice;
      }

      if (minVol < 0.0) {
         IJ.showMessage ("Mip", "Error: current stack contains negative values.");
         return;
     }

     if (minVol >= 30000)
        valToSubtract = 32768.0;
     else if (minVol >= 15000)
        valToSubtract = 16384.0;
      else if (minVol >= 7500)
        valToSubtract = 8192.0;
     else if (minVol >= 3250)
        valToSubtract = 4096.0;
     else
        valToSubtract = 0.0;

//    Si apre una dialog box per la richiesta del numero di proiezioni da costruire ...

      GenericDialog gd = new GenericDialog("Mip", IJ.getInstance());
      gd.addNumericField("Projections              ", nProj, 1);
      gd.addNumericField("Depth correction factor  ", valDepthCorr, 1);
      gd.addNumericField("Counts to subtract       ", valToSubtract, 1);
      gd.addCheckbox("Depth correction         ", true);
      gd.addCheckbox("Total Intensity Proj     ", false);
      gd.addMessage("(Contents: " + minVol + " - " + maxVol+ ")");
      gd.showDialog();
      if (gd.wasCanceled())
         return;
      nProj = (int)gd.getNextNumber();
      if (gd.invalidNumber()) {
         IJ.error("Number of projections invalid.");
         return;
      }
      valDepthCorr = (int)gd.getNextNumber();
      if (gd.invalidNumber()) {
         IJ.error("Depth correction factor invalid.");
         return;
      }
      valToSubtract = (int)gd.getNextNumber();
      if (gd.invalidNumber()) {
         IJ.error("Number of counts to subtract invalid.");
         return;
      }                                               
      lDepthCorr = gd.getNextBoolean();
      lTIP = gd.getNextBoolean();

      if (IJ.debugMode)
         IJ.showMessage("Mip", "lDepthCorr: " + lDepthCorr);

      if (nProj < 1) {
         IJ.error("At least one projection required.");
         return;
      }
      if ((valDepthCorr < 0.1 || valDepthCorr > 10.0) && lDepthCorr) {
         IJ.error("Depth correction factor must lie in the range 0.1 ... 1.0.");
         return;
      }
      if (valToSubtract < 0.0 || valToSubtract > maxVol) {
         IJ.error("Number of counts to subtract must lie in the range 0.0 ... " +
                  maxVol + ".");
         return;
      }
               
//    Verifica ulteriore che l'immagine sia ad 8 o 16 bit ...

      int imgType = img.getType();

      if (imgType != ImagePlus.GRAY8 && imgType != ImagePlus.GRAY16) {
         IJ.error("Image types allowed: only 8-bit and 16-bit grayscale.");
         return;
      }

//    Si passa dal numero di proiezioni all'angolo di rotazione per proiezione ...

      angStep = 360.0 / nProj;

      if (IJ.debugMode)
         IJ.showMessage("Mip - Debug", "Ready to create " + nProj + 
                        " projections (step angle " + angStep + ") " +
                        "from a stack of " + nSect + "sections. " +
                        "Type " + img.getType() + ".");

//    Si creano due ImagePlus, una per lo stack ruotato (modifImg) e una per lo stack
//    delle proiezioni (projImg) ...

      ImagePlus modifImg = new ImagePlus();
      ImagePlus projImg = new ImagePlus();

      if (modifImg == null || projImg == null) {
         IJ.showMessage("Mip", "Unable to create a new image.");
         return;
      }
    
      String origTitle = img.getTitle();

      if (imgType == ImagePlus.GRAY8) {
         modifImg = NewImage.createByteImage(origTitle + " - Rotation 0", 
                                             img.getWidth(), img.getHeight(),
                                             nSect, NewImage.FILL_BLACK);
         projImg = NewImage.createByteImage(origTitle + " - Projections", 
                                           img.getWidth(), nSect,
                                           nProj, NewImage.FILL_BLACK);
      }
      else if (imgType == ImagePlus.GRAY16) {
         modifImg = NewImage.createShortImage(origTitle + " - Rotation 0", 
                                              img.getWidth(), img.getHeight(),
                                              nSect, NewImage.FILL_BLACK);
         projImg = NewImage.createShortImage(origTitle + " - Projections", 
                                             img.getWidth(), nSect,
                                             nProj, NewImage.FILL_BLACK);
      }
      else {
         IJ.error("Internal error\n\nImage types allowed: " +
                  "only 8-bit and 16-bit grayscale.");
         return;   
      }

      if (modifImg == null || projImg == null) {
         IJ.showMessage("Mip", "Unable to create a new image.");
         return;
      }

//    Si costruisce lo stack delle proiezioni prima ruotando di un angolo dato lo stack
//    originale e copiandolo in modifImg (modifStack), poi creando il nuovo stack dalla
//    media dei pixel aventi stessa x e stessa z, sommati lungo l'asse y. L'ordine dei 
//    pixel secondo x deve essere invertito (verificato empiricamente per confrontto con
//    il plugin "VolumeJ / Volume Reconstructor") ...

      ImageStack modifStack = modifImg.getStack();
      ImageStack projStack = projImg.getStack();

      int modifSliceDimension = modifStack.getWidth() * modifStack.getHeight();
      int pixVal;
      double pixSum, pixNormVal, pixMax, pixAbsMax;
      double weightFact = 0.0;
      double volMaxSum = Double.MIN_VALUE;

      int index;
      int xDim = modifStack.getWidth();
      int yDim = modifStack.getHeight();

      IJ.showStatus("Calculating minimum counts per pixel in global volume ...");

      maxVol = Double.MIN_VALUE;
      minVol = Double.MAX_VALUE;

      for (z = 1; z <= nSect; ++z) {

         origSliceIp = origStack.getProcessor(z);
         minSlice = origSliceIp.getMin();
         maxSlice = origSliceIp.getMax();

         if (minSlice < minVol)
            minVol = minSlice;

         if (maxSlice > maxVol)
            maxVol = maxSlice;
      }

      if (minVol < 0.0) {
         IJ.showMessage ("Mip", "Error: current stack contains negative values.");
         return;
      }

      if (IJ.debugMode)
         IJ.showMessage("Mip", "Max counts: " + maxVol + "\nMin counts: " + minVol);

//    Se si e' selezionato il metodo TIP (Total Intensity Projection), si effettua una prima 
//    scansione completa del volume (ruotandolo a step e calcolando le somme per ciascun 
//    pixel delle future proiezioni) per cercare il massimo assoluto delle somme ...

      if (lTIP) {

         double localMax = 0.0;

         if (IJ.debugMode)
            IJ.showMessage("valToSubtract " + valToSubtract + " - valDepthCorr " + valDepthCorr);

         for (angCurr = 0, currProj = 1; angCurr < 360.0; angCurr += angStep, currProj++) {

            IJ.showStatus("Calculating max counts per pixel in projection # " + currProj + " ...");
            IJ.showProgress((double)angCurr / 360.0);
         
            projSliceIp = projStack.getProcessor(currProj);

            for (sectCurr = 1; sectCurr <= nSect; sectCurr++) {
  
               origSliceIp = origStack.getProcessor(sectCurr);

               modifSliceIp = modifStack.getProcessor(sectCurr);
               modifSliceIp.setInterpolate(true);
               modifSliceIp.copyBits(origSliceIp, 0, 0, Blitter.COPY);
               modifSliceIp.rotate(angCurr);
               modifSliceIp.add(-((int)valToSubtract));

//    Ciascuna sezione routata viene corretta sottraendo il minimo valutato sull'
//    intero volume ...

               for (x = 0; x < xDim; x++) {
               
                  pixSum = 0.0;
                  weightFact = 0.0;
            
                  for (y = 0; y < yDim; y++) {
                  
                     weightFact = ((double)y + 1.0) * valDepthCorr;
                     pixVal = modifSliceIp.getPixel(x, y);
                     if (pixVal < 0)
                        pixVal = 0;
                     
                     if (imgType == ImagePlus.GRAY16) {
                        if (lDepthCorr)
                           pixSum += ((double)pixVal / weightFact);           
                        else
                           pixSum += (double)pixVal;
                     }
                     else {
                        if (lDepthCorr)
                           pixSum += ((double)(0xff & pixVal) / weightFact);
                        else
                           pixSum += (double)(0xff & pixVal);
                     }
                  }
            
                  if (pixSum > volMaxSum)
                     volMaxSum = pixSum;
               }
            }
         }

         if (IJ.debugMode)
            IJ.showMessage("Mip", "Overall max counts: " + volMaxSum + ".");
      }

//    Se si e' selezionato il metodo TIP (Total Intensity Projection), una volta
//    trovato il massimo delle somme, si costruiscono effettivamente le proiezioni
//    impiegando il valore della somma normalizzato a 255 per le immagini di tipo
//    GRAY8, e normalizzato a 65535 per le immagini di tipo GRAY16 ...

      if (lTIP) {

         for (angCurr = 0, currProj = 1; angCurr < 360.0; angCurr += angStep, currProj++) {

            IJ.showStatus("Building projection # " + currProj + " ...");
            IJ.showProgress((double)angCurr / 360.0);
         
            projSliceIp = projStack.getProcessor(currProj);
   
            for (sectCurr = 1; sectCurr <= nSect; sectCurr++) {
         
               origSliceIp = origStack.getProcessor(sectCurr);
               modifSliceIp = modifStack.getProcessor(sectCurr);
               modifSliceIp.setInterpolate(true);
               modifSliceIp.copyBits(origSliceIp, 0, 0, Blitter.COPY);
               modifSliceIp.rotate(angCurr);
               modifSliceIp.add(-((int)valToSubtract));
           }
 
           for (z = 1; z <= nSect; z++) {

               modifSliceIp = modifStack.getProcessor(z);

               for (x = 0; x < xDim; x++) {
            
                  pixSum = 0.0;
                  weightFact = 0.0;

                  for (y = 0; y < yDim; y++) {

                     weightFact = ((double)y + 1.0) * valDepthCorr ;
                     pixVal = modifSliceIp.getPixel(x, y);
                     if (pixVal < 0)
                        pixVal = 0;                     
                  
                     if (imgType == ImagePlus.GRAY16) {
                        if (lDepthCorr) {
                           pixSum += ((double)pixVal / weightFact);      
                        }
                        else
                           pixSum += (double)pixVal;
                     }
                     else {
                        if (lDepthCorr)
                           pixSum += ((double)(0xff & pixVal) / weightFact);
                        else
                           pixSum += (double)(0xff & pixVal);
                     }
                  }
      
                  if (imgType == ImagePlus.GRAY16) {

                     pixNormVal = (pixSum / volMaxSum) * 32767.0;

                     if (IJ.debugMode) {
                        if (pixNormVal <= 32767.0)
                           projSliceIp.putPixel(xDim-x-1, z-1, (short)pixNormVal);
                        else {
                           IJ.showMessage("Mip", "Overflow. Pixel counts: " + pixNormVal +
                                          " - Angle: " + angCurr + " - z: " + z + " - x:" +
                                          x + " - y: " + y);
                           return;
                        }
                     }

                     else {
                        projSliceIp.putPixel(xDim-x-1, z-1, (short)pixNormVal);                                               
                     }
                  }

                  else {

                     pixNormVal = (pixSum / volMaxSum) * 255.0;

                     if (IJ.debugMode) {
                        if (pixNormVal <= 255.0)
                           projSliceIp.putPixel(xDim-x-1, z-1, (byte)pixNormVal);
                        else {
                           IJ.showMessage("Mip", "Overflow. Pixel counts: " + pixNormVal);
                           return;
                        }
                     }

                     else {
                        projSliceIp.putPixel(xDim-x-1, z-1, (byte)pixNormVal & 0xff);
                     }
                  }
               }
            }
         }
      }

//    Se non si e' selezionato il metodo TIP (Total Intensity Projection), viene implementato
//    l'algoritmo MIP (Maximal Intensity Projection) classico, che consiste nel selezionare,
//    per ogni pixel da proiettare, il solo valore massimo, normalizzandolo al massimo assoluto
//    di tutto il volume, prventivamente ricercato ...

      else {

//    Si ricerca il massimo assoluto di tutto il volume ...

         pixAbsMax = 0.0;
   
         for (z = 1; z <= nSect; z++) {

            origSliceIp = origStack.getProcessor(z);

            for (x = 0; x < xDim; x++)
               for (y = 0; y < yDim; y++)
                  if ((pixVal = origSliceIp.getPixel(x, y)) > pixAbsMax)
                     pixAbsMax = pixVal;        
         }

         if (IJ.debugMode)
            IJ.showMessage ("Debug", "Massimo assoluto: " + pixAbsMax);

         for (angCurr = 0, currProj = 1; angCurr < 360.0; angCurr += angStep, currProj++) {

            IJ.showStatus("Building projection # " + currProj + " ...");
            IJ.showProgress((double)angCurr / 360.0);
         
            projSliceIp = projStack.getProcessor(currProj);

            for (sectCurr = 1; sectCurr <= nSect; sectCurr++) {
  
               origSliceIp = origStack.getProcessor(sectCurr);
               modifSliceIp = modifStack.getProcessor(sectCurr);
               modifSliceIp.setInterpolate(true);
               modifSliceIp.copyBits(origSliceIp, 0, 0, Blitter.COPY);
               modifSliceIp.rotate(angCurr);
               modifSliceIp.add(-((int)valToSubtract));
           }
 
           for (z = 1; z <= nSect; z++) {

               modifSliceIp = modifStack.getProcessor(z);

               for (x = 0; x < xDim; x++) {
            
                  pixSum = 0.0;               
                  pixMax = 0.0;
                  weightFact = 0.0;

                  for (y = 0; y < yDim; y++) {

                     weightFact = ((double)y + 1.0) * valDepthCorr ;
                     pixVal = modifSliceIp.getPixel(x, y);
                     if (pixVal < 0)
                        pixVal = 0;                     
                  
                     if (imgType == ImagePlus.GRAY16) {
                        if (lDepthCorr) {
                           pixSum = ((double)pixVal / weightFact);      
                        }
                        else
                           pixSum = (double)pixVal;
                     }
                     else {
                        if (lDepthCorr)
                           pixSum = ((double)(0xff & pixVal) / weightFact);
                        else
                           pixSum = (double)(0xff & pixVal);
                     }
               
                     if (pixSum > pixMax)
                        pixMax = pixSum;
                  }

                  pixNormVal = (pixMax / pixAbsMax) * 32767.0;
                  projSliceIp.putPixel(xDim-x-1, z-1, (short)pixNormVal);
               }
            }
         }
      }

//    Si normalizza e si visualizza lo stack delle proiezioni dopo aver applicato
//    la lookup table "fire" ...

      projImg.setTitle(origTitle + " - " + nProj + " Projections");      

      projImg.updateAndDraw();
      projImg.show();

      IJ.runPlugIn("ij.plugin.LutLoader","fire");

      projSliceIp = projImg.getProcessor();
      projSliceIp.resetMinAndMax ();
      projImg.updateAndDraw();
      projImg.show();
   }

//-------------------     
   void showAbout() {
//-------------------     

		IJ.showMessage("About Mip ...",
			"'Mip' stands for 'Maximal Intensity Projection'. A set of projections\n" +
         "is generated starting from a stack of transaxial sections. The number\n" +
         "of the projections is settled by the user");
   }

}

// *** EOF ***