//// Sept 25 2006
//T.J.Collins, MacBiophotonics, McMaster University, Hamilton, Canada.
//www.macbiophotonics.ca
//This example script shows how to identify different subcellular compartments for analysis 
//and how to calculate the intensity ratio between two compartments. e.g. nucleus to cytosol but could 
//plasmamembrane to cytoplasm etc.

//generate sample image
A1_0001()
rename(Ch1 = image)

//gnerate sample 'ch2' image
flip(yes, yes, image=Ch1)
rename(Ch2=image)


// find nucleus and cytosol 
nuclei_detection_B(Ch1)
rename(Nuclei_obj = Nuclei)
//Nuclei stencil is now Nucleis_obj.body

//detect whole cell objects
cytoplasm_detection_A(Ch1, nuclei=Nuclei_obj)
rename(Wholecells_obj = Wholecells)
//cytoplasm stencil is now Wholecells_obj.body


//detect plasmamembrane objects
CalcZone(objects=Wholecells_obj)
ZoneMask(-2, 2)
rename(Plasmamembrane_obj = objects)
//plasmamembrane stencil is now Plasmamembrane_obj.zonemask


//create cytoplasm only objects
//first create a 'not' nucleus image
not(image=Nuclei_obj.body.image)
//then create a image which is Not-nucleus AND whole cell
and(image, image = Wholecells_obj.body)
//calculate obects from this image
setAttr("body", image, objects = Wholecells_obj)
calcBorder(body)
rename(Cytoplasm_obj = objects)
//cytoplasm stencil is now Cytoplasm_obj.body


//create cytoplasm minus plasmamembrane objects
not(image=Plasmamembrane_obj.zonemask.image)
and(image, image = Cytoplasm_obj.body)
setAttr("body", image, objects = Wholecells_obj)
calcBorder(body)
rename(Cytoplasm_minus_PM_obj = objects)
//"cytoplasm minus plasmamembrane" stencil is now Cytoplasm_minus_PM_obj.body


//detect perinuclear objects
CalcZone(objects=Nuclei_obj)
ZoneMask(-4, 0)
rename(Perinuclear_obj = objects)
//Perinuclear stencil is now Perinuclear_obj.zonemask

//CALCUALTE INTENSITY OF A PARTICULAR CELL COMPARTMENT
//We can use the stencils to calculate mean intensity in a particular area of the image 
//using the calcIntensty() function:

//calcIntenisty(STENCIL NAME, image = IMAGE NAME, objects = OBJECTS, total = YES/NO)
//STENCIL NAME = the "short name" of the stencil to use, i.e. body or zonemask etc.
//IMAGE NAME = the name of the image to be anayled 
//OBJETCS= the name of the objedt list that holds your setncil, e.g Nuclei_obj or Wholecells_obj
//total = YES if you want to total intensities in the area; no if you want the mean intensity

//e.g.
calcIntensity(body, image=Ch2, Total=no, objects=Cytoplasm_obj)
rename(Ch2_Cytoplasm_obj = objects)
//you can find the mean intensity of the Ch2 cytoplasm objects here: Ch2_Cytoplasm_obj.intensity.mean


// CALCULATING THE RATIO OF TWO CELLS COMPARTMENTS
//calculate mean cytosolic intensity ratioed again mean nuclear intensity
//First create mean cytosolic intenisty image, the mean nuclear intensity image using 
//wholecells_obj.body stencil  - using this particular stencil is not important 
//so long as we use the same stencil for both compartments

//first calculate mean intensity for cytoplasm objects
calcIntensity(body, image= Ch1, Total = no, objects = Cytoplasm_obj)
//create new image of mean cytopalsm intensity
blank(Ch1.width, Ch1.height)
carryobjects(Wholecells_obj.body, objects.intensity, image=image)
rename(mCytoplasm = image)

//now calculate mean nuclear intensity 
calcIntensity(body, image= Ch1, Total = no, objects = Nuclei_obj)
//create new image of mean Nuclear intensity
blank(Ch1.width, Ch1.height)
carryobjects(Wholecells_obj.body, objects.intensity, image=image)
rename(mNuclei  = image)

//now create Nuc div cyto image
div(mNuclei, mCytoplasm, result_type="float", infinity=0)
rename(mNuc_div_mCyto = result)

//if you want individual intensities out..
calcIntensity(body, Total=no, objects=Cytoplasm_obj, Image=mNuc_div_mCyto)
rename(mNuc_div_mCyto_obj = objects)