//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// "Standard Feature Output" version 0.14 January 2011                                                                                                                                                                                  
// Developed by Jarkko Ylanko at McMaster Biophotonics Facility 2009-2010
// Additional MBF procedures developed by Tony J Collins and Sean Cianflone at McMaster Biophotonics Facility 2008-2009         
// Report questions, bugs and concerns to info@macbiophotonics.ca
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 1. START
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 1.1 - USER INPUTS
input(IN_mode, "Segmentation adjustment", "MBF_FeatureExtraction mode selection", "s", "Select analysis mode. Optimize object detection in 'Segmentation adjustment', analyze image data in 'Feature extraction', or review classification results in 'Image gallery' (only available after image data analyzed in MBF_Classify). Possible values: Segmentation adjustment, Feature extraction, Image Display. Default: Segmentation adjustment.", 0.5)
//input(IN_adjust, 1, "MODE 1. Adjust segmentation parameters", "y", "Browse images and adjust object segmentation parameters. Refer to Nuclei, Cytoplasm and Spot Detection options below. It is strongly recommended that multiple fields of view and/or wells are considered when confirming segmentation accuracy. Default: ON.",1)
//input(IN_fx, 0, "MODE 2. Feature extraction", "y", "Extract selected features from image data. Available features and file-saving options are selected below. Default: OFF.", 2)
//input(IN_gallery, 0, "MODE 3. Image gallery", "y", "Display object classification results on select images. Currently compatible with classification data obtained by MBF_Classify MatLab procedures only. Default: OFF.", 3)
input(IN_ImageONOFF, 0, "Show object reference images in Player Window", "y", "View object borders from primary segmentation procedures in Player window. Default: off.", 4) // if one would like to see the object_detection procedure progress in the Player (note: this will display object borders as determined by nuclei_- and cytoplasm_detection prior to any object filtering coded below)
input(IN_fov, 0, "Select field of view for analysis", "i", "Field of view to analyze (1...total field count); if 0, all fields will be analyzed. Default: All Fields (0).", 5) // indicate field of interest
//input(IN_fieldselect, "1", "Field selection mode", "i", "Field range to analyze. 1 = all fields, 2 = every 2nd field, etc. Default: All Fields.") // indicate field of interest
input(IN_Stencil, "All", "Select object masks for analysis", "s", "Object stencils to use for analysis below. Possible values: All, Nucleus and Cytoplasm, Nucleus, Cytoplasm, Cell. Default: All", 6)
//input(IN_ChDefault, 1, "Automatic channel assignment: Channel assignment", "i", "Automatic channel assignments, assuming nuclear reference in Camera 3. Default: on.") // default assumes DRAQ5 for nuclear reference
input(IN_Ch1, 0, "Reference channel for nuclei detection: Nuclei Detection", "i", "Indicate channel corresponding to object reference image according to slice order of multi-channel image file. When set to zero, the reference will be set to the channel presumed to contain far red signal in Opera Camera 3, if present, or else the first channel in the series. By convention, the reference channel will be reported as 'IM_1' (image) or 'Ch1' (channel) in the output data, and additional '#' channels as 'IM_#' or 'Ch#'. Default: 0.", 7) 
input(IN_RGBCh1, "Not used", "RGB reference channel: Nuclei Detection", "s", "Indicate RGB channel corresponding to object reference image. By convention, the reference channel will be reported as 'Image 1' or 'Channel 1' in the output below. Possible values: Red, Green, Blue, Not used. Default: Not used.",8) 
input(IN_cyto, 1, "Cytoplasm detection: Cytoplasm detection", "y", "Segmentation of cytoplasm. Detection parameter adjustments below. Default: on.", 9) // toggle for cyto detection
input(IN_CytoRef, 1, "Cytoplasm reference image: Cytoplasm detection", "i", "Indicate channel corresponding to cytoplasm reference image, relative to nuclear reference image specified by convention as 'IM_1' or 'Ch1' above. Value '#' selects the #th image in the multi-channel series as cytoplasm reference. Default: 1, use nuclear reference image.", 9.1)
input(IN_spots, 0, "Small spot detection: Spot Detection", "y", "Segmentation of small (<30 pixel) and large (>30 pixel) foci of intensity. Default: off.",14) // toggle spot detection on/off
input(IN_LSpots, 0, "Large spot detection: Spot Detection", "y", "BETA VERSION. Segmentation of large (>30 pixel) subcellular aggregates. Default: off.",10.1) // toggle spot detection on/off
input(IN_LSpotThresh, 3, "Large spot intensity threshold adjustment: Spot Detection", "i", "BETA VERSION. Adjustment for intensity threshold applied to segmenting large spots; 0..oo, high values select for brighter spots. Default: 3.",11) // intensity threshold for scoring large spots
input(IN_LSpotRound, 0.7, "Large spot roundness threshold adjustment: Spot Detection", "i", "BETA VERSION. Adjustment for roundness threshold applied to segmenting large spots; values approaching 0 tolerate irregular spots. Default: 0.7",12) // roundness threshold for large spots
input(IN_RefSEG, 0, "Include segmentation on reference image: Spot Detection", "y", "Attempt segmentation of sub-cellular structures using reference image. Default: off.",10) // toggle 1/2 for IM_1
input(IN_intensity, 0, "[INT] General intensity measurements: Features", "y", "Calculate intensity features on reference and signal channel(s). Default: off.", 15) // toggle intensity calculation suite on/off
input(IN_morph, 0, "[MOR] General morphology measurements: Features", "y", "Calculate area, perimeter, and compactness features on reference image. Default: off.",16) // toggle morphology calculation suite on/off
input(IN_Zone, 0, "Perinuclear zone segmentation: Other Detection", "y", "Define zone radiating out from nuclear border for intensity measurements. Only applies if intensity measurements are being calculated. Default: off.", 14.1) // define perinuclear space
input(IN_ZoneCount, 3, "Perinuclear zone width: Other Detection", "i", "Pixel width of perinuclear zone. Default: 3 pixels.",14.2) // tune size of perinuclear zone for intensity measurements
input(IN_Haralick, 0, "[TXT] Haralick texture features: Features", "y", "Haralick texture features.", 19) // toggle texture analysis on/off
input(IN_HarStepNum, 1, " ..Number of steps in Haralick calculation: Features", "i", "Number of steps to use in Haralick calculations. Default: 1.",20) // adjust number of Haralick steps
input(IN_TAS, 0, "[TXT] Threshold adjacency statistics: Features", "y", "Threshold adjacency statistics.", 21) // toggle texture analysis on/off
input(IN_TRangeMode, 1, " ..TAS range mode: Features", "i", "Threshold adjacency statistics calculation range mode. 1, two range points; 2, four range points; 3, nine range points. Default: 1", 22) // TAS calculation restriction
input(IN_TASmode, 1, " ..TAS output mode: Features", "i", "Threshold adjacency statistics calculation mode. 1, output all TAS steps; 2, every second step; 3, every third. Default: 1.", 23) // TAS feature output restriction
input(IN_RadMom, 0, "[TXT] Radial moments: Features", "y", "Calculates the kth radial moment (given by the integral of (r^k)*I(r) where I(r) is the intensity of the object based on object radius) of detected objects highlighting either the center or border as defined by the user.", 24) // toggle texture analysis on/off
input(IN_RefTXT, 0, "[TXT] Include texture analysis of reference image: Features", "y", "Measure texture features from reference image. Default: off.", 25) // toggle 2/2 for IM_1
input(IN_Pearsons, 0, "[CLC] Pearsons coefficient: Features", "y", "Calculates Pearson's coefficient, see Reference: Press, W.H., Flannery, B.P., Teukolsky, S.A. & Vetterling, W.T. (1988) Numerical Recipes in C. Cambridge University Press, Cambridge.", 26) // toggle colocalisation analysis on/off
input(IN_Manders, 0, "[CLC] Manders coefficient: Features", "y", "Calculates Manders colocalisation values. See Reference:Manders, EEM, Verbeek, FJ., and Aten, JA. (1993). Measurement of co-localisation of objects in dual-colour confocal images. Journal of Microscopy 169, 375-382.", 27) // toggle colocalisation analysis on/off
input(IN_ICQ, 0, "[CLC] ICQ coefficient: Features", "y", "Calculates ICQ colocalisation values. See Reference: Li, Qi, Lau, Anthony, Morris, Terence J., Guo, Lin, Fordyce, Christopher B., and Stanley, Elise F. (2004). A Syntaxin 1, G{alpha}o, and N-Type Calcium Channel Complex at a Presynaptic Nerve Terminal: Analysis by Quantitative Immunocolocalization. Journal of Neuroscience 24, 4070-4081.", 28) // toggle colocalisation analysis on/off
input(OP_data, 0, "Write data to file: Output", "y", "Write object-level data as .txt file.", 29) // specify whether or not to write data files, default = no, in order to avoid accidental redundancy in OP files
input(OP_destination, "default", "Target directory: Output", "s", "Destination directory for data output. By default, a new directory named 'output' will be created in the same location as the image data selected for analysis.", 30) // default given by dir from MBF_getinfo3()
input(OP_treat, 0, "Save data by treatment: Output", "y", "Save object-level data in separate files according to treatment (defined below). Default: no separation, all data in single output file.", 31) // save data in multiple files according to multiple treatments
input(OP_treatID, "NA", "Filename treatment: Output", "s", "Treatment name to be used in output data filenames, if applicable. Possible values: Treatment_Sum, Treatment01, Treatment02, Treatment03, NA. Default: NA.", 32)
input(OP_imagecount, "5", "Number of object images to output (Mode 3): Output", "i", "Select the maximum number of object images per classification per field to be output. Used only if in Image Display mode. If set to 0, all object images will be generated. Default: 5.", 33)
input(OP_FieldResult, "Not used", "Output field image (Mode 3): Output", "s", "Optional: output classification results as stencils on field image. Used only in Image Display mode. Possible values: Not used, All fields, Informative fields. A Field is considered informative only if multiple classification events are present within its object population. Default: off.", 34)
set(IN_cleanup=1) // toggle on/off (1/0) to remove clutter from Data Explorer window
set(IN_NAN=0) // toggle on/off (1/0) to remove NAN from output data
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
if(IN_mode=="Segmentation adjustment")
set(IN_adjust=1, IN_fx=0, IN_gallery=0)
end()
if(IN_mode=="Feature extraction")
set(IN_adjust=0, IN_fx=1, IN_gallery=0)
end()
if(IN_mode=="Image Display")
set(IN_adjust=0, IN_fx=0, IN_gallery=1)
end()
set(ModeSelect=IN_adjust+IN_fx+IN_gallery)
// 1.2 - OBLIGATE CODE
singlewell(compact=yes)
// 1.2.1 - EXPERIMENT INDEX (from MBF_GetInfo3 in MBF_ScriptBits.proc by Tony Collins)
set(platemap_tab="null")
if(length(wellIndex)==7)
        set(wellIndexLng="0"&wellIndex)
        else()
        set(wellIndexLng=wellIndex)
end()
// wellindex
set(row=substr(wellIndexLng, 1,3 )+0)
set(column=substr(wellIndexLng, 4,3 )+0)
eval(set(tmp=sourcedata.sublayout[0]))
if(errorcode>0)
        if(defined("sourcedata.field"))
                set(sourcedata.sublayout=sourcedata.field)
                else()
                set(sourcedata.sublayout=1)
        end()
        // defined(sourcedata.field)
end()
// errorcode
if(defined("sourcedata.field"))
        TableFilter("Field==1", table=sourcedata)
        else()
        TableFilter("sublayout==1", table=sourcedata)
        set(Table.Field=Table.sublayout)
        delete(table.sublayout)
end()
// defined(sourcedata.field)
if(defined("sourcedata.Channel"))
        set(nCh=sourcedata.Channel.max)
        else()
        set(nCh=sourcedata.sourceimage.length)
end()
// defined(sourcedata.channel)
if(defined("sourcedata.field"))
        set(fov=sourcedata.Field.max)
        else()
        set(fov=sourcedata.sublayout.max)
end()
// defined(sourcedata.field)
set(nImages=sourcedata.@rowcount/fov)
set (directoryStop=at(wellIndex, sourcedata.sourcefilename[0], -1))
set (dir=substr (sourcedata.sourcefilename[0], 0, directoryStop-1))
set(lastBackSlash=at("/",sourcedata.sourcefilename[0],-2))
set(exptFolder=substr(sourcedata.sourcefilename[0], 0, lastBackSlash))
set (firstSlash=at("/", sourcedata.sourcefilename[0], 1))
set (secondSlash=at("/", sourcedata.sourcefilename[0], 2))
if(defined("sourcedata.barcode"))
        set(barcode=sourcedata.Barcode[0])
        else()
        set(barcode="NoBarcode")
        set(sourcedata.barcode=barcode)
end()
// defined(sourcedata.barcode)
set(skip=0)
eval(readfile(filename=dir&"\\platemap.txt"))   
if(warningcode!=2050)
        readtable(filename=dir&"\\platemap.txt", fileformat="ascii", columns="firstrow", separators="\t")
        tablefilter("wellindex=="&wellindex, table=table)
        if(table.@rowcount>0)
                foreach(0..table.@columncount-1, "c")
                output(table.[c][0], table.@columns[c])
                end()
                // c
        Insert_At(table, sourcedata.sourcefilename[0], "Path", 1)
        Insert_At(table, sourcedata.Barcode[0], "Barcode", 1)
        set(platemap_tab=table) 
        else() // ie. if no "platemap.txt" is defined
        warning("Platemap is missing from image directory or "&wellindex&" is not accounted for!", confirm=no)
        set(skip=1)
        if(defined("table"))
                delete(table)
        end()
        create("table")
        insert_at(table, "temp", "temp")
        Insert_At(table, "Treatment03", "Treatment03", 1)
        Insert_At(table, "Treatment02", "Treatment02", 1)
        Insert_At(table, "Treatment01", "Treatment01", 1)
        Insert_At(table, "Treatment_Sum", "Treatment_Sum", 1)
        Insert_At(table, Column, "Column", 1)
        Insert_At(table, Row, "Row", 1)
        Insert_At(table, "Dye", "Dye", 1)
        Insert_At(table, "Cells", "Cells", 1)
        Insert_At(table, sourcedata.sourcefilename[0], "Path", 1)
        Insert_At(table, sourcedata.Barcode[0], "Barcode", 1)
        Insert_At(table, wellIndex, "WellIndex", 1)
        delete(table.temp, table.sourceimage, table.sourcefilename, table.frameindex, table.field)
        set(platemap_tab=table) 
        end()
        // table.@columncount>0
end()
//warningcode!2050
if(skip!=1)
delete(barcode,directoryStop,errorcode,errors,errors_xml,firstSlash,lastBackSlash,secondSlash,table,warningcode,wellIndexLng)
set(extensionpt=at(".", sourcedata.sourcefilename[0],-1))
set(exe=substr(sourcedata.sourcefilename[0],extensionpt+1))
if(exe!="flex")
        if(IN_RGBCh1=="Not used")
                warning("Object reference channel for RGB image has not been assigned! Confirm channel assignment and try again.", confirm=yes)
                stop()
        else()
                if(IN_RGBCh1=="Blue")
                        push(_RGB, "", "red", "green")
                end()
                // blue
                if(IN_RGBCh1=="Red")
                        push(_RGB, "", "blue", "green")
                end()
                // red
                if(IN_RGBCh1=="Green")
                        push(_RGB, "", "blue", "red")
                end()
                // green
        end()
        // not used
        set(image=sourcedata.sourceimage[0])
        RGBsplit(image)
        set(IM_1=_[&IN_RGBCh1])
        set(rc2=2)
        set(nCh=1)
        foreach(1..2, "rc")
                set(sigIM=_RGB[rc])
                set(IMtest=_[&sigIM])
                if(IMtest.max>0)
                        set(_["IM_"&rc2]=_[&sigIM])
                        set(rc2=rc2+1)
                        set(nCh=nCh+1)
                end()
                // image.max>0
        end()
        // "rc"
        if(in_cleanup==1)
        delete(_RGB, rc, rc2,IMtest,sigIM)
        end()
        // cleanup
end()
// image is RGB ie. non-flex
delete(exe, extensionpt)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 1.3 - CREATE OUTPUT DIRECTORIES
if(!defined("dir"))
        set(dir="\\\mercury\operaimages\\"&OP_destination)
        else()
        set(date=__date__)
        set(datebreak1=at(".", date,1))
        set(datebreak2=at(".", date,2))
        set(dateyear=substr(date,datebreak2+1,4))
        set(datemonth=substr(date,datebreak1+1,datebreak2-datebreak1-1))
        set(dateday=substr(date,0,datebreak1))
        set(date=dateyear&"."&datemonth&"."&dateday)
        set(dir2=dir&"\\output\\"&date&"_output")
delete(datebreak1, datebreak2, dateyear, datemonth, dateday)
end()
// if defined dir
if(OP_data==1 AND IN_mode=="Feature extraction")
        MkDir(dir2)
end()
//if OP_data == 1
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 1.4a - IMAGE DISPLAY
if(IN_gallery==1)
eval(readfile(filename=dir&"\DirTemp.txt"))
//eval(readfile(filename=dirIM&"\ClassificationCoordinates.txt"))
if(warningcode!=2050)
readfile(filename=dir&"\DirTemp.txt")
set(content=substr(content,0,at("/",content,-1)))
set(dirIM=content&"/ImageGallery_"&date)
readtable(filename=dirIM&"\CoordinatesTemp.txt")
set(filter=WellIndex*1)
tablefilter("wellindex==filter", table=table)
set(CoordinateTable=table)
delete(table)
else()
getfiles()
set(dirIM=substr(files[0],0,at("/",files[0],-1)))
write(dirIM&"/", dir&"\DirTemp.txt", "ascii")
set(dirIM=dirIM&"/ImageGallery_"&date)
MkDir(dirIM)
foreach(0..files.length-1, "fl")
set(_["filename"&fl+1]=Files[fl])
set(IN_filename=Files[fl])
set(filetypebreak=at(".", IN_filename,-1))
set(filetype=substr(IN_filename,filetypebreak,4))
set(fileshortname=substr(IN_filename,0,filetypebreak))
if(filetype!=".dat")
        readtable(filename=IN_filename)
        //writetable(table, dirIM&"\CoordinatesTemp.txt")
        else()
        eval(readfile(filename=fileshortname&".txt"))
                if(errorcode==0)
                readtable(filename=fileshortname&".txt")
                end()
                // 0
        ReadTable(filename=IN_filename, separators="\t")
        push(_columnheaders, "wellindex", "barcode", "dir", "platemap", "fov", "Xcoord", "Ycoord", "classification")
        foreach(2..table.rowcount-1, "t")
        //foreach(0..table.rowcount-1, "t")
                set(coord=table.A[t])
                foreach(0.._columnheaders.length-1, "s")
                        set(header=_columnheaders[s])
                        set(space=(8-s)*-1)
                        set(headerspace=at(" ",coord,space))
                                if(s==7)
                                        push(_["index_"&s],substr(coord, headerspace))
                                        //set(headervalue=1*substr(coord, headerspace))
                                        else()
                                        set(spacediff=at(" ", coord, space+1)-headerspace)
                                        push(_["index_"&s],substr(coord, headerspace, spacediff))
                                        //set(headervalue=1*substr(coord, headerspace, spacediff))
                                end()
                                // s==7
                end()
                // s
        end()
        // foreach t
        create("table")
        foreach(0.._columnheaders.length-1,"s")
                Insert_At(table, _["index_"&s],_columnheaders[s],table.length)
                delete(_["index_"&s])
        end()
        // s
        delete(_columnheaders)
        eval(readfile(filename=dirIM&"\CoordinatesTemp.txt"))
        if(errorcode==2050)
                //writetable(table, fileshortname&".txt")
                //writetable(table, dir2&"\ClassificationCoordinates.txt")
                //writetable(table, dirIM&"\CoordinatesTemp.txt")
                writetable(table,dirIM&"\CoordinatesTemp.txt")
                else()
                write(table,dirIM&"\CoordinatesTemp.txt","ascii", append=yes)
                end()
                // 2050
end()
// if filetype!=.DAT
set(_["table"&(fl+1)]=table)
if(fl>0)
        appendtable(tableX, _["table"&(fl+1)])
        set(tableX=table)
        else()
        set(tableX=table)
end()
//n>1
set(filter=WellIndex*1)
tablefilter("wellindex==filter", table=tableX)
set(CoordinateTable=table)
delete(table)
end()
// "fl"
end()
// error2050
end()
// IN_gallery

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 1.4 - INPUT FIELD OF VIEW
if(!defined("fov"))
        set(fov=1)
end()
// if !defined fov
if(IN_fov>fov)
        set(IN_fov=fov)
end()
// IN_fov>fov
if(IN_fov>0) // analyze one field
        set(StartField = IN_fov)
        set(EndField = IN_fov)
        else()
        set(StartField=1)
        set(EndField=fov)
end()
        // if IN_fov==0
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 1.5 - DETERMINE NUCLEAR REFERENCE ASSIGNMENTS
if(IN_RGBCh1=="Not used")
push(ChannelSeries1, "", "0", "1", "2", "3")
push(ChannelSeries2, "", "0", "-1", "1", "2")
push(ChannelSeries3, "", "0", "-2", "-1", "1")
push(ChannelSeries4, "", "0", "-3", "-2", "-1")
if(IN_Ch1==0)
        if(nCh>1)
                foreach(1..nCh, "j") // "j" used for same nCh loop below
                        set(cam=substr(SourceData.name[j-1], 8,1 )+0)
                        set(cam3=no)
                        if(cam==3)
                                set(nuc_ch=(j-1)) // if reference channel=ch_1, then (j-1)=0; if 2, then (j-1)=1, etc.
                                set(ChannelSeries=_["ChannelSeries"&j])
                                delete("ChannelSeries1", "ChannelSeries2", "ChannelSeries3", "ChannelSeries4")
                                set(cam3=yes)
                        end()
                        // if cam==3
                        if(cam3==no)
                                set(IN_Ch1=1)
                        end()
                        //cam3=no
                        delete(cam3)
                end()
                //foreach 1..nCh
        end()
        //if nCh>1
        else()
        set(ChannelSeries=_["ChannelSeries"&IN_Ch1])
        set(nuc_ch=IN_Ch1-1)
end()
// if IN_Ch1==1
end()
// if IN_RGBCh1=="Not used"
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 1.6 MISCELLANEOUS INPUTS
if(IN_ImageONOFF==1)
        set(_imageview="yes")
        else()
        set(_imageview="no")
end()
// if imageonoff =1
push(_signal, "", "Ch1", "Ch2", "Ch3", "Ch4") // by convention, Ch1 = IM_1, ie nuclear reference
if(IN_cyto==0)
set(IN_Stencil="Nucleus")
end()
//if(IN_cyto==0)
if(IN_Stencil=="All")
push(_stencil, "", "Cell", "Nucleus", "Cytoplasm")
end()
// IN_Stencil
        if(IN_Stencil=="Nucleus and Cytoplasm")
        push(_stencil, "", "Nucleus", "Cytoplasm")
        end()
        // IN_Stencil
                if(IN_Stencil=="Nucleus")
                push(_stencil, "", "Nucleus")
                end()
                // IN_Stencil
                        if(IN_Stencil=="Cytoplasm")
                        push(_stencil, "", "Cytoplasm")
                        end()
                                // IN_Stencil
                                if(IN_Stencil=="Cell")
                                push(_stencil, "", "Cell")
                                end()
                                // IN_Stencil
set(fieldcount=0)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 2. ANALYSIS
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 2.1 - SIGNAL CHANNEL ASSIGNMENTS
foreach(StartField..EndField, "i")
if(IN_RGBCh1=="Not Used")
        if(IN_Ch1==0)
                if(nCh>1)
                        foreach(1..nCh, "j") // "j" used for same nCh loop below
                        set(_["IM_"&j]=SourceData.SourceImage[(i-1)*nCh+nuc_ch+ChannelSeries[j]]) // specify as signal of interest in channel "j"
                        end()
                        // for each nCh
                else()
                        set(IM_1=SourceData.SourceImage[i-1])
                end()
                // if nCh>1
                else()
                foreach(1..nCh, "j")
                        set(_["IM_"&j]=SourceData.SourceImage[(i-1)*nCh+nuc_ch+ChannelSeries[j]]) // specify as signal of interest in channel "j"
                end()
                // j
        end()
        // if IN_Ch1==0
end()
// if IN_RGBCh1=="Not used"
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 2.2 - IMAGE QUALITY CHECKPOINT - OPTIONAL
if(IN_gallery==1)
        set(IN_adjust=1)
end()
// IN_gallery
if(IN_adjust==1)
nuclei_detection_select(IM_1, _imageview)
if(IN_cyto==0)
        cytoplasm_detection_a(IM_1, no)
        else()
        if(IN_CytoRef!=1 AND IN_CytoRef<(nCh+1))
        set(CytoRef=_["IM_"&IN_CytoRef])
        else()
        set(CytoRef=IM_1)
        end()
        // IN_CytoRef==1
        cytoplasm_detection_a(CytoRef, _imageview)
        delete(CytoRef)
end()
//IN_cyto
objectfilter(area>0, objects=wholecells) 
removeborderobjects(2, objects=objects)
if(objects.count>0)
CalcMassCentre(objects=objects)
SetAttr("aaFieldofView", filledvec(objects.@count, i))
SetAttr("aaYcoord", objects.MassCentre.y.target)
SetAttr("aaXcoord", objects.MassCentre.x.target)
set(fieldcount=fieldcount+1)
calcattr(centers_area, nuclei.area[objectindex], objects=objects)
calcattr(centers_contrast, nuclei.contrast[objectindex], objects=objects)
calcattr(centers_intensity, nuclei.intensity[objectindex], objects=objects)
set(_["obj_"&i]=objects)
if(defined("all_cells"))
        addobjects(objects, DeleteGeometry=yes, objects=all_cells)
        set(all_cells=objects)
        else()
        set(all_cells=objects)
end()
//if defined all_count
        set(IM_blue=(255*IM_1/IM_1.max))
                if(nCh>1)
                set(IM_green=(255*IM_2/IM_2.max))
                        if(nCh>2)
                        set(IM_red=(255*IM_3/IM_3.max))
                        RgbJoin(red=IM_red, blue=IM_blue, green=IM_green)
                        set(_["RGB_"&i]=image)
                        delete("IM_red")
                        else()
                        RgbJoin(blue=IM_blue, green=IM_green)
                        set(_["RGB_"&i]=image)
                        end()
                        // (nCh>2)
                        delete("IM_green")
                else()
                RgbJoin(blue=IM_1)
                set(_["RGB_"&i]=image)
                end()
                // (nCh>1)
        delete("IM_blue")
end()
//object.count>0
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 2.2b - IMAGE DISPLAY CONT.
if(objects.count>0)
if(IN_gallery==1)
set(dir3=dir2&"\\ImageGallery")
//      if(OP_data==1)
//      MkDir(dir3)
//      end()
//      //OP_data==1
tablefilter("fov*1==i", table=CoordinateTable)
if(table.@rowcount!=0)
set(table.Xcoord=table.Xcoord*1)
set(table.Ycoord=table.Ycoord*1)
set(CoordinateTableX=table)
set(Categories=CoordinateTableX.classification*1)
set(CategoriesCount=Categories.max)
        if(categories.max>3)
        set(CategoriesCount=3)
        end()
        // categories.max>3
set(objects=_["obj_"&i])
foreach(0..objects.count-1, "oc")
tablefilter(Xcoord==objects.aaXcoord[oc] AND Ycoord==objects.aaYcoord[oc], table=CoordinateTableX)
if(table.@rowcount!=1)
        if(table.@rowcount>1)
                warning("Error in classified data set: multiple objects at same X and Y coordinate. Inspect input data set.", confirm=yes)
                stop()
        end()
        // multiple objects at a given X and Y
        if(table.@rowcount==0)
                push(_vecClass, 0)
//              warning("Error in classified data set: at least one object is missing an X and Y coordinate. Confirm segmentation parameters and try again.", confirm=yes)
//              stop()
        end()
        // no objects with the same X and Y
else()
push(_vecClass, table.classification[0])
end()
// catch for multiple or no objects at input X and Y
end()
// oc
calcattr(Category, _vecClass[objectindex], objects=_["obj_"&i])
delete(_vecClass)
calcerosion(-1, cytoplasm_border, objects=objects)
set(objects_all=objects)
// FIELD OUTPUT
if(OP_FieldResult!="Not used")
set(OP_Field=0)
if(OP_FieldResult=="Informative fields" AND categories.max>categories.min)
set(OP_Field=1)
else()
if(OP_FieldResult=="All fields")
set(OP_Field=1)
end()
// all fields
end()
// informative fields
if(OP_Field==1)
MkDir(dirIM&"\FieldCollection")
        //set(IM_category=IM_2*255/(IM_2.mean)+3*(IM_2.stddev))
        //set(IM_category=IM_2*255/IM_2.mean)
        //set(IM_category=_["RGB_"&i])  
        //set(IM_category=IM_2)
split(5000,image=IM_2)
delete(bright_image)
        gamma(image=image)
        set(IM_category=image)
        push(_colour, "", "rgb(255,255,0)", "rgb(255,0,0)", "rgb(238,130,238)")
        foreach(1..Categories.max, "c")
                if(c<4)
                set(colour=_colour[c])
                objectfilter(Category==c, objects=objects_all)
                if(objects.count>0)
                        set(_["objects_category_"&c]=objects)
                        //carrypixels(objects.Cytoplasm_border_eroded,200, image=IM_category)   
                        //carrypixels(objects.Cytoplasm_border_eroded, colour, image=image)
                        carrypixels(objects.Cytoplasm_border, IM_category.max, image=IM_category)       
                        carrypixels(objects.Cytoplasm_border, colour, image=image)
                set(IM_category=image)
                end()
                // objects.count>0
                end()
                // c<4
        end()
        // c 
        ImageView(IM_category, "Signal image showing classification results for Field "&i, image=IM_category)
//set(IM_temporary=IM_category)
        eval(readfile(dirIM&"\FieldCollection\\"&wellindex&"_F"&i))
        if(warningcode==2050)
                if(OP_data==1)
                        writeimage(dirIM&"\FieldCollection\\"&wellindex&"_F"&i&".jpg", image=IM_category, imageformat="jpg")
                end()
                // OP_data
        end()
        // warning code
        delete(_colour, colour)
        end()
// OP_field==1
end()
// OP_field result != not used

// OBJECT-LEVEL IMAGE OUTPUT
set(imagecount=OP_imagecount)
//set(IM_category=IM_2*255/IM_2.max)
split(5000,image=IM_2)
delete(bright_image)
        gamma(image=image)
set(IM_category=image)
foreach(1..Categories.max, "c")
        if(c<4)
        objectfilter(Category==c, objects=objects_all)
        if(objects.count>0)
        set(_["objects_category_"&c]=objects)
        set(objectstmp=_["objects_category_"&c])
        if(OP_data==1)
                set(dir4=dirIM&"\\class_"&c)
                MkDir(dir4)
        end()
        // OP_data==1
        if(objectstmp.count>imagecount)
                set(ObjectStart=floor((objectstmp.count-1)/2)-floor(imagecount/2))
                set(ObjectEnd=floor((objectstmp.count-1)/2)+floor(imagecount/2))
                else()
                set(ObjectStart=0)
                set(ObjectEnd=objectstmp.count-1)
        end()
        // ObjectLoop definition
                foreach(ObjectStart..ObjectEnd, "d")
                        set(X=objectstmp.aaXcoord[d])
                        set(Y=objectstmp.aaYcoord[d])
                        objectfilter(aaXcoord==X AND aaYcoord==Y, objects=objectstmp)
                                carrypixels(objects.Cytoplasm_border, IM_category.max, image=IM_category)
                                carrypixels(objects.Cytoplasm_border, "white", image=image)
                                set(FieldImage=image)
                                crop(X-50, Y-50,X+50, Y+50, image=image)
                                set(CropImage=image)
                                convelems(image, size=1)
                                if(OP_data==1)
                                        writeimage(dir4&"\\"&wellindex&"_F"&i&"_x"&X&"_y"&Y&".jpg", image=result, imageformat="jpg")
                                end()
                                // OP_data
                end()
                // d
        end()
        //objects.count>0
        end()
        // c<4
end()
// c
delete(image,IM_category,Categories,CategoriesCount,CoordinateTableX,objects,objects_all,objects_category_1,objects_category_2,objectstmp)
if(c==3)
delete(objects_category_3)
end()
//c==3
delete(c, objectstart, objectend)
end()
// present well is in the coordinate table
end()
// IN_gallery
end()
// objects.count>0 (line407)

else() // do all the work below
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 2.3 - OBJECT SEGMENTATION AND INITIAL MEASUREMENTS
// 2.3.1 - NUCLEI DETECTION
nuclei_detection_select(IM_1, _imageview)
calcattr(roundness, objects=nuclei)
renameattr(compact=roundness)
calcarea(border, objects=objects)
CalcWidthLength()
set(Nuclei=objects)
// 2.3.2 - CYTOPLASM DETECTION
if(IN_cyto==0)
        cytoplasm_detection_a(IM_1, no)
        else()
        if(IN_CytoRef!=1 AND IN_CytoRef<(nCh+1))
        set(CytoRef=_["IM_"&IN_CytoRef])
        else()
        set(CytoRef=IM_1)
        end()
        // IN_CytoRef==1
        cytoplasm_detection_a(CytoRef, _imageview)
        delete(CytoRef)
end()
//IN_cyto
objectfilter(area>0, objects=wholecells)
deleteattr(area)
RemoveBorderObjects(2, objects=objects) 
if(objects.count>0)
set(fieldcount=fieldcount+1)
CalcMassCentre(objects=objects)
SetAttr("aaFieldofView", filledvec(objects.@count, i))
SetAttr("aaYcoord", objects.MassCentre.y.target)
SetAttr("aaXcoord", objects.MassCentre.x.target)
deleteattr("NumberOfNuclei", objects=objects)
set(wholecells=objects)
//if(IN_RefSEG==1)
//      set(ChStart=1)
//      set(ChEnd=nCh)
//      else()
//      set(ChStart=2)
//      set(ChEnd=nCh)
//end()
//// IN_RefSEG==1
if(IN_RefSEG!=1)
        if(nCh>1)
        set(ChStart=2)
        set(ChEnd=nCh)
set(ChEnd=2) // REMOVE THIS -- 20110323
        else()
        set(ChStart=1)
        set(ChEnd=nCh)
        end()
        // nCh>1
else()
set(ChStart=1)
set(ChEnd=nCh)
end()
// IN_RefSEG!=1

// 2.3.3 - SPOT DETECTION
if(IN_spots==1 OR IN_Lspots==1)
set(wholecellstemp=wholecells)
set(stencilcount=1)
        if(IN_Stencil=="Nucleus and Cytoplasm" OR IN_Stencil=="All")
        Stencil2Objects(wholecells.centers)
        calcarea()
        set(WholeCells_1=objects)
        Stencil2Objects(wholecells.cytoplasm)
        calcarea()
        set(WholeCells_2=objects)
        set(stencilcount=2)
        end()
        // IN_Stencil = "Nucleus and Cytoplasm"
                if(IN_Stencil=="Nucleus")
                Stencil2Objects(wholecells.centers)
                calcarea()
                set(WholeCells_1=objects)
                end()
                // IN_Stencil = "Nucleus"
                        if(IN_Stencil=="Cytoplasm")
                        Stencil2Objects(wholecells.cytoplasm)
                        calcarea()
                        set(WholeCells_1=objects)
                        end()
                        // IN_Stencil = "Cytoplasm"
                                if(IN_Stencil=="Cell")
                                Stencil2Objects(wholecells.body)
                                calcarea()
                                set(WholeCells_1=objects)
                                end()
                                // IN_Stencil = "Cell"
foreach(ChStart..ChEnd, "j")
set(signal=_signal[j])
set(FieldImage=_["IM_"&j])
foreach(1..stencilcount, "e")
//if(stencilcount==2)
        //set(stencilx=_stencil[e+1])
        //else()
        set(stencilx=_stencil[e])
//end()
////stencilcount
set(objects_temp=_["WholeCells_"&e])
calcintensity(image=FieldImage, objects=objects_temp)
set(objects_temp=objects)
// 2.3.3.2 - LARGE SPOTS // adapted from MBFACA001 made for Sandra Nelson in March 2010
if(IN_Lspots==1)
foreach(0..objects_temp.count-1, "c") // detect "large spots" one object at a time; c = cell counter
        set(body=objects_temp[c].body) // 
        Mask2Stencil(body) // convert mask of object #c to stencil..
        Stencil2Objects(Stencil) // .. and stencil to objectlist
        calcarea() // calculating area for the cell mask
        calcintensity(image=FieldImage, Total=no) // calculating intensity for the cell mask
        //set(object_c=objects) // rename objectlist for later
        set(image=FieldImage*body.image) // define image to include signal found within current cell mask
        set(thld=IN_LSpotThresh) // adjust threshold as needed via input
        cluster(30, thld*objects_temp.intensity.mean, image) // identify bright regions above given threshold, ie large spots
        calcattr("wholecell_index", &c) // assign index for corresponding cell mask
        calcintensity(image=image, total=no) // spot intensity
        calcattr(roundness) // spot roundness -- may be used for filter below if only the roundest objects are of interest
        calcarea(border) // spot perimeter
        renameattr(perimeter=border_area)
        objectfilter(area<objects_temp[c].area) // area filter
        objectfilter(roundness>IN_LSpotRound) // roundness filter
        if(objects.count>0) // collect large spot attributes into vectors
                push(_LScount, objects.count)
                push(_LSarea, objects.area.mean)
                push(_LSintensity, objects.intensity.mean)
                push(_LSperimeter, objects.perimeter.mean)
                push(_LSroundness, objects.roundness.mean)
        //      push(_LSareaofcell, objects.area.sum/wholecells[c].area)
                        if(defined("LargeSpotsTemp"))
                                addobjects(objects, CheckOverlap=no, DeleteGeometry=no, objects=LargeSpotsTemp) // collect large spot attributes into objectlist
                                set(LargeSpotsTemp=objects)
                                else()
                                set(LargeSpotsTemp=objects)
                        end()
                        // if defined LargeSpots
                else() // ie no large spots scored
                push(_LScount, 0)
                push(_LSarea, 0)
                push(_LSintensity, 0)
                push(_LSperimeter, 0)
                push(_LSroundness, 0)
        //      push(_LSareaofcell, 0)
        end()
        //if attrtable.rowcount>0
end()
// foreach c
set(_["wholecellstemp."&signal&"_MOR_"&stencilx&"_Spots_Large_Count"]=_LScount)
set(_["wholecellstemp."&signal&"_MOR_"&stencilx&"_Spots_Large_AvgArea"]=_LSarea)
set(_["wholecellstemp."&signal&"_MOR_"&stencilx&"_Spots_Large_AvgPerimeter"]=_LSperimeter)
set(_["wholecellstemp."&signal&"_MOR_"&stencilx&"_Spots_Large_AvgRoundness"]=_LSroundness)
set(_["wholecellstemp."&signal&"_INT_"&stencilx&"_Spots_Large_AvgIntensity"]=_LSintensity)
if(defined("LargeSpotsTemp"))
        set(_["LargeSpots_"&stencilx]=LargeSpotsTemp)
end()
//LargeSpotsTemp
end()
// if IN_Lspots == 1

// 2.3.3.1 - SMALL SPOTS
if(IN_spots==1)
spot_detection_a(FieldImage, "body", no, ShowOutputParameters=no, WholeCells=objects_temp)
if(defined("LargeSpotsTemp"))
Inverse(image=LargeSpotsTemp.body.mask.image) // adjust small spot count from Spot_Detection procedure above to remove overlap with large spots
set(image=image-254)
set(SpotCorrection=WholeCells.SpotCenters.mask.image*image)
Mask(image=SpotCorrection)
Mask2Stencil()
Stencil2Objects()
CalcMassCentre()
set(SpotCorrection=objects.MassCentre.mask.image)
calcintensity(body, SpotCorrection, objects=WholeCells, Total=yes)
RenameAttr(SmallSpotCount=intensity)
RenameAttr(SmallSpotAverageIntensity=SpotAverageCenterIntensity)
DeleteAttr(SpotsArea)   // adjust NumberOfSpots to include Large Spots but not Small Spots that coincide with the Large Spot masks
//calcattr(NumberOfSpots, LargeSpots.count+SmallSpotCount) // adjust NumberOfSpots to include Large Spots but not Small Spots that coincide with the Large Spot masks
set(WholeCellsSpotCorrected=objects) // WholeCell objectlist including spot data (but not all spot stencils)
else()
RenameAttr(SmallSpotCount=NumberOfSpots, objects=wholecells)
RenameAttr(SmallSpotAverageIntensity=SpotAverageCenterIntensity)
end()
// defined LargeSpots
set(SmallSpots=spots)
delete(integratedspotsignalperarea)
set(_["SmallSpots_"&stencilx]=SmallSpots)
end()
// IN_spots==1


// 2.3.3.3 - ADD SPOT DATA TO WHOLECELL OBJECT LISTS

// 2.3.3.3.3 - RENAMING OF ATTRIBUTES
//calcattr(spots_countestimate, objects.Spots_small_count[objectindex]+objects.Spots_large_count[objectindex])
if(IN_Lspots==1)
        if(IN_spots==1)
        set(_["wholecellstemp."&signal&"_MOR_"&stencilx&"_Spots_Small_Count"]=objects.SmallSpotCount)   
        set(_["wholecellstemp."&signal&"_INT_"&stencilx&"_Spots_Small_AvgIntensity"]=objects.SmallSpotAverageIntensity)
        foreach(0..objects.count-1, "c")        
                push(_totalspot, objects.SmallSpotCount[c]+_LScount[c])
        end()
        // c
        delete(c)
        set(_["wholecellstemp."&signal&"_MOR_"&stencilx&"_Spots_Total_Count"]=_totalspot)
        end()
        // IN_Spots==1
        else()
        set(_["wholecellstemp."&signal&"_MOR_"&stencilx&"_Spots_Small_Count"]=wholecells.NumberOfSpots) 
        set(_["wholecellstemp."&signal&"_INT_"&stencilx&"_Spots_Small_AvgIntensity"]=wholecells.SpotAverageCenterIntensity)
end()
// IN_Lspots==1
delete(_LSarea, _LScount, _LSintensity, _LSperimeter, _LSroundness, _totalspot)

// PER FIELD OUTPUTS //
if(IN_ImageONOFF==1)
foreach(ChStart..ChEnd, "j2")
        if(defined("SmallSpots"))
                set(spotstemp1=_["SmallSpots_"&stencilx])
                convert(2, image=spotstemp1.border.mask.image) // output image of small spots into Acapella Player window
                set(IM_SpotMask=image)
                gamma(image=_["IM_"&j2])
                set(IM_Jadjusted=image)
                paint(IM_SpotMask, image=IM_Jadjusted)
                set(_["IM_view_small"&j2]=image)
                imageview(_["IM_view_small"&j2], "SmallSpots_"&stencilx&"_Ch"&j2, delay@s=1)
                delete(spotstemp1)
        end()
        // defined smallspots
        if(defined("LargeSpotsTemp"))
                set(spotstemp2=_["LargeSpots_"&stencilx])
                convert(2, image=spotstemp2.border.mask.image) // output image of large spots into Acapella Player window
                paint(image, image=IM_Jadjusted)
                set(_["IM_view_large"&j2]=image)
                imageview(_["IM_view_large"&j2], "LargeSpots_"&stencilx&"_nCh"&j2, delay@s=1)
                delete(_["IM_view_large"&j2], image, spotstemp2)
        end()
        // defined LargeSpots
                delete(image, _["IM_view_small"&j2], IM_SpotMask, IM_Jadjusted, SpotCorrection)
end()
// j
delete(j2)
end()
// if(IN_ImageONOFF==1)
if(defined("LargeSpotsTemp"))
delete(LargeSpotsTemp)
end()
// defined LargeSpotsTemp
if(defined("SmallSpots"))
delete(SmallSpots)
end()
// defined SmallSpots
end()
// foreach stencil "e"
delete(FieldImage, e)
end()
// foreach ch
delete(j, stencilcount)
set(WholeCells=wholecellstemp)
// 2.3.3.3.4 - EXPLORER CLEANUP
if(IN_cleanup==1)
Delete(image, body, objects, objects_temp, wholecellstemp, WholeCells_1, WholeCells_2, WholeCells_3, Stencil, SpotCorrection, mask, SpotCandidates, Spots, NumberOfSpotCandidates, NumberOfSpots,SpotsPerArea,SpotsPerObject,Spots_large_counter,Spots_small_counter,spots_large_area,spots_small_area, integratedspotssignalperarea, integratedspotsignalperarea_backgroundsubtracted, integratedspotsignalpercellularsignal, integratedspotsignalpercellularsignal_backgroundsubtracted, distance_vector, thld, ItemCount)
if(IN_Spots==1)
delete(WholeCellsSpotCorrected)
end()
//IN_Spots==1
end()
// cleanup

end()
// IN_spots
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 2.4 - INTENSITY MEASUREMENTS
// 2.4.1 - INTENSITY BY COMPARTMENT
if(IN_intensity==1)
foreach(1..nCh, "j") // foreach channel
        set(signal=_signal[j])
        set(image=_["IM_"&j])
        if(IN_Stencil == "All" OR IN_Stencil == "Cell")
        calcintensity(body, image, CalcStdDev=yes, objects=WholeCells)
        set(_["wholecells."&signal&"_INT_Cell_intensity"]=objects.intensity)
        set(_["wholecells."&signal&"_INT_Cell_intensity_stddev"]=objects.intensity_stddev)
        deleteattr(intensity_stddev, objects=objects)
        end()
        // IN_Stencil = All or Cell
                if(IN_Stencil == "All" OR IN_Stencil == "Nucleus and Cytoplasm" OR IN_Stencil == "Nucleus")
                calcintensity(centers, image, CalcStdDev=yes, objects=WholeCells)
                set(_["wholecells."&signal&"_INT_Nucleus_intensity"]=objects.centers_intensity)
                set(_["wholecells."&signal&"_INT_Nucleus_intensity_stddev"]=objects.centers_intensity_stddev)
                end()
                // IN_Stencil = All, Nucleus & Cytoplasm, or Nucleus
                        if(IN_Stencil == "All" OR IN_Stencil == "Nucleus and Cytoplasm" OR IN_Stencil == "Cytoplasm")
                        calcintensity(Cytoplasm, image, CalcStdDev=yes, objects=WholeCells)
                        set(_["wholecells."&signal&"_INT_Cytoplasm_intensity"]=objects.cytoplasm_intensity)
                        set(_["wholecells."&signal&"_INT_Cytoplasm_intensity_stddev"]=objects.cytoplasm_intensity_stddev)
                        end()
                        // IN_Stencil = All, Nucleus & Cytoplasm, or Cytoplasm
        set(objects=wholecells)
        foreach(1..9, "u")
                set(quant=u/10)
                if(IN_Stencil == "All" OR IN_Stencil == "Cell")
                CalcStat("quantile", quant, body, Image=image)
                set(_["wholecells."&signal&"_INT_Cell_quantile"&(u*10)]=_["objects.quantile"])
                end()
                        if(IN_Stencil == "All" OR IN_Stencil == "Nucleus and Cytoplasm" OR IN_Stencil == "Nucleus")
                        CalcStat("quantile", quant, centers, Image=image)
                        set(_["wholecells."&signal&"_INT_Nucleus_quantile"&(u*10)]=_["objects.quantile"])
                        end()
                                if(IN_Stencil == "All" OR IN_Stencil == "Nucleus and Cytoplasm" OR IN_Stencil == "Cytoplasm")
                                CalcStat("quantile", quant, cytoplasm, Image=image)
                                set(_["wholecells."&signal&"_INT_Cytoplasm_quantile"&(u*10)]=_["objects.quantile"])
                                end()
        end()
        //foreach u
        if(IN_Zone==1) // create perinuclear zone
                Stencil2Objects(wholecells.centers)
                calczone(IN_ZoneCount, objects=objects, ZoneType="standard")
                zonemask(filt_low=-1, filt_high=-3)
                renameattr(Nucleus_outerzone=zonemask)
                set(WholeCells_CentersZone=objects)
                calcintensity(nucleus_outerzone, image, CalcStdDev=yes, objects=WholeCells_CentersZone)
                set(_["wholecells."&signal&"_INT_PerinuclearZone_intensity"]=objects.nucleus_outerzone_intensity)
                set(_["wholecells."&signal&"_INT_PerinuclearZone_intensity_stddev"]=objects.nucleus_outerzone_intensity_stddev)
                if(IN_cleanup==1)
                        delete(WholeCells_CentersZone)
                end()
                // cleanup
        end()
        // if IN_Zone=1
end()
// intensities of each channel
end()
// if(IN_intensity==1)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 2.5 - MORPHOLOGY MEASUREMENTS
if(IN_morph==1)
set(objects=wholecells)
// 2.5.1 - COMPACTNESS 
if(IN_Stencil == "All" OR IN_Stencil == "Cell")
calcattr(roundness, objects=objects)
renameattr(Ch1_MOR_Cell_compact=roundness)
end()
// IN_Stencil == "All" OR IN_Stencil == "Cell")
if(IN_Stencil == "All" OR IN_Stencil == "Nucleus and Cytoplasm" OR IN_Stencil == "Nucleus")
calcattr(Ch1_MOR_Nucleus_compact, nuclei.compact[objectindex])
end()
// IN_Stencil = All, or Nucleus & Cytoplasm, or Nucleus
// 2.5.2 - AREA AND PERIMETER
if(IN_Stencil == "All" OR IN_Stencil == "Cell")
calcarea(border)
renameattr(Ch1_MOR_Cell_area=area, Ch1_MOR_Cell_perimeter=border_area)
calcattr(Ch1_MOR_Cell_perim2area, objects.Ch1_MOR_Cell_perimeter[objectindex]/objects.Ch1_MOR_Cell_area[objectindex])
end()
// IN_Stencil == All, Cell
        if(IN_Stencil == "All" OR IN_Stencil == "Nucleus and Cytoplasm" OR IN_Stencil == "Nucleus")
        calcattr(Ch1_MOR_Nucleus_perimeter, nuclei.border_area[objectindex])
        calcarea(centers)
        renameattr(Ch1_MOR_Nucleus_area=centers_area)
        calcattr(Ch1_MOR_Nucleus_perim2area, objects.Ch1_MOR_Nucleus_perimeter[objectindex]/objects.Ch1_MOR_Nucleus_area[objectindex])
        end()
        // IN_Stencil = All, Nuc & Cyto, Nucleus
                if(IN_Stencil == "All" OR IN_Stencil == "Nucleus and Cytoplasm" OR IN_Stencil == "Cytoplasm")
                calcarea(cytoplasm)
                renameattr(Ch1_MOR_cytoplasm_area=cytoplasm_area)
                calcarea(cytoplasm_border)
                renameattr(Ch1_MOR_cytoplasm_perimeter=cytoplasm_border_area)
                calcattr(Ch1_MOR_cytoplasm_perim2area, objects.Ch1_MOR_cytoplasm_perimeter[objectindex]/objects.Ch1_MOR_cytoplasm_area[objectindex])
                end()
                // IN_Stencil = all, nuc & cyto, cytoplasm
// 2.5.3 - WIDTH AND LENGTH
if(IN_Stencil == "All" OR IN_Stencil == "Cell")
CalcWidthLength()
renameattr(Ch1_MOR_Cell_fulllength=full_length, Ch1_MOR_Cell_halfwidth=half_width)
end()
// if(IN_Stencil == "All" OR IN_Stencil == "Cell")
if(IN_Stencil == "All" OR IN_Stencil == "Nucleus and Cytoplasm" OR IN_Stencil == "Nucleus")
calcattr(Ch1_MOR_Nucleus_fulllength, nuclei.full_length[objectindex])
calcattr(Ch1_MOR_Nucleus_halfwidth, nuclei.half_width[objectindex])
end()
// if(IN_Stencil == "All" OR IN_Stencil == "Nucleus and Cytoplasm" OR IN_Stencil == "Nucleus")
set(wholecells=objects)
end()
// if(IN_morph==1)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 2.6 - TEXTURE FEATURE MEASUREMENTS
if(IN_RefTXT==1)
        set(ChStart=1)
        set(ChEnd=nCh)
        else()
        set(ChStart=2)
        set(ChEnd=nCh)
end()
// IN_RefTXT==1
// 2.6.1 - STENCIL ASSIGNMENTS
if(IN_Haralick==1 OR IN_TAS==1 OR IN_RadMom==1)
if(IN_Stencil=="All")
Stencil2Objects(wholecells.body)
calcarea()
set(WholeCells_1=objects)
Stencil2Objects(wholecells.centers)
calcarea()
set(WholeCells_2=objects)
Stencil2Objects(wholecells.cytoplasm)
calcarea()
set(WholeCells_3=objects)
end()
// IN_Stencil = all
        if(IN_Stencil=="Nucleus and Cytoplasm")
        Stencil2Objects(wholecells.centers)
        calcarea()
        set(WholeCells_1=objects)
        Stencil2Objects(wholecells.cytoplasm)
        calcarea()
        set(WholeCells_2=objects)
        end()
        // IN_Stencil = "Nucleus and Cytoplasm"
                if(IN_Stencil=="Nucleus")
                Stencil2Objects(wholecells.centers)
                calcarea()
                set(WholeCells_1=objects)
                end()
                // IN_Stencil = "Nucleus"
                        if(IN_Stencil=="Cytoplasm")
                        Stencil2Objects(wholecells.cytoplasm)
                        calcarea()
                        set(WholeCells_1=objects)
                        end()
                        // IN_Stencil = "Cytoplasm"
                                if(IN_Stencil=="Cell")
                                Stencil2Objects(wholecells.body)
                                calcarea()
                                set(WholeCells_1=objects)
                                end()
                                // IN_Stencil = "Cell"
foreach(ChStart..ChEnd, "j") // foreach channel
set(signal=_signal[j])
set(image00=_["IM_"&j])
foreach(1.._stencil.size-1, "e")
set(stencilx=_stencil[e])
set(objects00=_["WholeCells_"&e])
set(signal=_signal[j])
// 2.6.2 - MOMENTS
if(IN_RadMom==1)
foreach(2..5, "v")
        CalcMoment(v, Image=image00, objects=objects00)
        set(_["wholecells."&signal&"_TXT_"&stencilx&"_moment"&v]=_["objects.moment"&v])
end()
//foreach v
// Radial Moments -- Created by Sean Cianflone, DWA lab, June 11, 2008; adapted from MBF_texture_suite.proc lines #495-553 by Jarkko Ylanko, MBF, December 9, 2010
foreach(1..2, "i3")
        if(i3==1)
                set(type="center")
                else()
                set(type="border")
        end()
        // if i3
set(original=image00)
set(ROI_obj=objects00)
set(k=2)
calcstat("sum", stencil=ROI_obj.body, image=original)
rename(ROI_obj_sum=objects)
calczone(objects=ROI_obj, Stencil=body, ZoneType="standard")
        zoneimage(objects=objects)
        set(radial=zoneimage)   
CalcStat("max", Image=radial, objects=ROI_obj)
        set(ROI_obj=objects)
blank(original.width, original.height, 1)
convelems(image, "float", nonnumeric=0)
carryobjects(ROI_obj.body, ROI_obj.max, image=result)
rename(maxzone=image)
set(stdradial=radial/maxzone)
set(stdborderrad=((-1)*stdradial)+1)
set(direction="C")
if(type=="border")
        set(direction="B")
        set(radimage=stdborderrad)
                        else()
                        if(type=="center")
                                set(radimage=stdradial)
                        end()
end()
// type definition
set(baseimage=radimage)
foreach(1..k, "i2")
        set(_["mult"&0]=1)
        set(_["mult"&i2]=baseimage*_["mult"&(i2-1)])
end()
//i2 - set
set(image=_["mult"&k])
foreach(1..k, "i2")
        delete(_["mult"&i2])
end()
//i2 - delete
delete(_["mult"&0], i2)
set(expedimage=image)
calcintensity(Image=expedimage*original, Total=yes, objects=ROI_obj)
blank(original.width, original.height, 0)
convelems(image, "float", nonnumeric=0)
carryobjects(ROI_obj.body, objects.intensity, image=result)
set(radmoment=image)
setAttr("Radial"&direction&k, objects.intensity/ROI_obj_sum.sum, objects=ROI_obj)
rename(radialmom_obj=objects)
DeleteAttr(max, objects=radialmom_obj)
delete(appliedshift,baseimage,direction,expedimage,k,maxzone,non_numerics_count,original,radial,radialmom_obj,radimage,radmoment,result,ROI_obj,ROI_obj_sum,stdborderrad,stdradial,v,ZoneImage)
        if(i3==1)
                set(_["WholeCells."&signal&"_TXT_"&stencilx&"_RadialC2"]=_["objects.RadialC2"])
                else()
                set(_["WholeCells."&signal&"_TXT_"&stencilx&"_RadialB2"]=_["objects.RadialB2"])
        end()
        // if i3==1
end()
// i3
delete(v,i3, type)
end()
// if IN_RadMom = 1
// 2.6.3 - HARALICK FEATURES
if(IN_Haralick==1)
foreach(1..IN_HarStepNum, "r")
// Haralick Features - Created by Tony Collins, MBF, 2008; adapted from MBF_texture_suite.proc lines #1-146 by Jarkko Ylanko, MBF, December 13, 2010
        set(objects_in=objects00)
        set(in_image=image00)
        set(im_ref=image00)
        set(step=r)
        set(orientation="east") 
                if(orientation=="north")
                        set(xstep=0)
                        set(ystep=-step)
                end()
                // north
                        if(orientation=="east")
                                set(xstep=step)
                                set(ystep=0)
                        end()
                        // east
                                if(orientation=="south")
                                        set(xstep=0)
                                        set(ystep=step)
                                end()
                                // south
                                        if(orientation=="west")
                                                set(xstep=-step)
                                                set(ystep=0)
                                        end()
                                        // west
        if(im_ref.elemtype==5) // if unsigned char 
                set(signalHARALICK=im_ref)
                set(signalHARALICK.factor=1)
                else()
                convelems(im_ref, size=1, shiftnegative=yes, clipoutofrange=no, neginfinity=0, nonnumeric=0, allowfactorchange=1 )
                set(signalHARALICK=result)
                set(signalHARALICK.factor=1)
        end()
        // elemtype==5
        set(maxshift=iif(abs(xstep)>abs(ystep), abs(xstep), abs(ystep)))
        if(maxshift>1)
                blank(im_ref.width,im_ref.height,1)
                rectmask((abs(xstep)-xstep)/2, (abs(ystep)-ystep)/2, im_ref.width-(abs(xstep)+xstep)/2, im_ref.height-(abs(ystep)+ystep)/2,image=image)
                and(image=objects_in.body.image,mask=mask.image)
                Stencil2Objects(image)
                set(objects_corrected=objects_in)
                else()
                set(objects_corrected=objects_in) 
        end()
        // maxshift>1
        set(kernelwidth=iif(abs(xstep)>abs(ystep), 2*abs(xstep)+1, 2*abs(ystep)+1))
        blank(kernelwidth,kernelwidth)
        set(image[(kernelwidth-1)/2+xstep,(kernelwidth-1)/2+ystep]=1)
        set(image.type="mask")
        set(kernel1=image.vector)
        convolution(image=signalHARALICK,convolutionkernel=kernel1)
        set(IM_neighbour=image) 
        mul(signalHARALICK, 256,result_type="unsigned,short")
        plus(result,IM_neighbour,result_type="unsigned,short")
        and(image=result,mask=objects_corrected.body.mask.image)
        set(result=image)
        set(result.type="stencil")
        stencil2objects(result) 
        calcarea()
        objectfilter(area>0) 
        calcstat("max", attrname="I_pixel",stencil=body, image=signalHARALICK) 
        calcstat("max", attrname="I_neighbour",stencil=body, image=IM_neighbour) 
        set(IM_IntensityPairs=objects.body.image) // hold on - from line 82 of RC57, aka line 85 from RC117     
        deleteattr(area,body,border,index)  
        set(OL_IntensityPairs=objects) 
        Foreach(0..objects_corrected.count-1,"i2") 
                objectfilter(objectindex==i2,objects=objects_corrected)
                set(ObjectArea=objects.area.sum)
                and(image= IM_IntensityPairs,mask=objects.body.mask.image)
                setattr(single_object,image.vector,objects=OL_IntensityPairs)
                calcarea(single_object)
                deleteattr(single_object)
                objectfilter(single_object_area) 
                calcattr(f,1.0*single_object_area/ObjectArea)
                calcattr(asm,f*f)
                calcattr(contrast, f*(i_pixel-I_neighbour)*(i_pixel-I_neighbour))
                calcattr(IDM, f/(1+(i_pixel-I_neighbour)*(i_pixel-I_neighbour)))
                calcattr(entropy,-f*ln(f))
                push(asm_vec, 1.0*objects.asm.sum)
                push(contrast_vec, 1.0*objects.contrast.sum)
                push(IDM_vec, 1.0*objects.IDM.sum)
                push(entropy_vec, 1.0*objects.entropy.sum)
                push(sum_vec, 1.0*objects.f.sum)
        end()
        set(Haralick_table=tbl("Haralick_ASM_"& step=asm_vec, "Haralick_contrast_" & step=contrast_vec, "Haralick_IDM_"&step=IDM_vec, "Haralick_entropy_"&step=entropy_vec, "Haralick_sum_"&step=sum_vec))
        set(_["WholeCells."&signal&"_TXT_"&stencilx&"_ASM_"&r]=_["haralick_table.Haralick_ASM_"&r])
        set(_["WholeCells."&signal&"_TXT_"&stencilx&"_contrast_"&r]=_["haralick_table.Haralick_contrast_"&r])
        set(_["WholeCells."&signal&"_TXT_"&stencilx&"_entropy_"&r]=_["haralick_table.Haralick_entropy_"&r])
        set(_["WholeCells."&signal&"_TXT_"&stencilx&"_IDM_"&r]=_["haralick_table.Haralick_IDM_"&r])
        delete("haralick_table",appliedshift,asm_vec,contrast_vec,entropy_vec,i2,IDM_vec,IM_IntensityPairs,IM_neighbour,im_ref,in_image,kernel1,kernelwidth,maxshift,non_numerics_count,ObjectArea,objects_corrected,objects_in,OL_Intensity,Pairs,orientation,result,step,sum_vec,xstep,ystep, signalHARALICK,OL_IntensityPairs)
end()
// foreach haralick step value r
delete(r)
end()
// if IN_Haralick=1
// 2.6.4 - TEXTURE ADJACENCY STATISTICS
if(IN_TAS==1)
        if(IN_TRangeMode==1)
        push(_TRange, "", "0.2", "0.4") // percentage of starting point pixels with neighbours of 0/1 = 1
        end()
        // TRangeMode 1
        if(IN_TRangeMode==2)
        push(_TRange, "", "0.2", "0.4", "0.6", "0.8") // percentage of starting point pixels with neighbours of 0/1 = 1
        end()
        // TRangeMode 2
        if(IN_TRangeMode==3)
        push(_TRange, "", "0.1", "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9") // percentage of starting point pixels with neighbours of 0/1 = 1
        end()
        // TRangeMode 3
                if(IN_TASmode==1)
                push(_TASmode,"",1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27) // keep each TAS statistic
                end()
                //mode 1
                if(IN_TASmode==2)
                push(_TASmode,"",2,4,6,8,10,12,14,16,18,20,22,24,26) // keep every other statistic
                end()
                //mode 2
                if(IN_TASmode==3)
                push(_TASmode,"",3,6,9,12,15,18,21,24,27) // keep every third stat
                end()
                //mode 3
foreach(1.._TRange.size-1, "s")
        set(TRange=_TRange[s])
        set(TNum=100*TRange)
// Texture Adjacency Statistic - Created by Tony Collins, MBF, 2008; adapted from MBF_texture_suite.proc lines #239-355 by Jarkko Ylanko, MBF, December 13, 2010
set(in_objects=objects00)
set(in_image=image00)
set(ratio_method="mean")
set(removeNAN="false")
set(range=TRange)
calcStat(ratio_method, Stencil=body, Image=in_image, objects=in_objects)
blank(in_image.width, in_image.height)
convelems(image, sign = "signed")
if(ratio_method=="mean")
        carryobjects(in_objects.body, objects.mean, image=result)
        else()
        carryobjects(in_objects.body, objects.median, image=result)
end()
// ratio_method=mean
set(img_avg=image)
div(in_image, image, result_type="float")
mask("inf", image=result)
mul(not(mask), result, result_type= "float")
set(ratio_image=result)
create("table")
foreach(0..2, "i3")
        if(i3==0) 
                set(range_min=range)
                set(range_max=range)
        end()
        if(i3==1)
                set(range_min=range)
                set(range_max=500)
        end()
        if(i3==2)
                set(range_min=0)
                set(range_max=500)      
        end()
        mask(1-range_min, image=ratio_image)
        set(mask1=mask)
        mask(1+range_max, image=ratio_image)
        set(mask2=mask)
        minus(mask1, mask2)
        rename(mask_all=result)
        blank(3,3,9)
        rename(kernel=image)
        convolution(image=mask_all, ConvolutionKernel=kernel, Faster=yes, ConvolutionKernelFactor=9)
        rename(adj_img=image)
        calcstat("sum", AttrName="tas_mTotal"&i3+1, Stencil=body, Image=mask_all, objects=objects)
        if(removeNAN)
                objectfilter("tas_mTotal"&i3+1&">0", objects=objects)
        end()
        // remove NAN
        mul(adj_img, mask_all)
        set(_["adj_mask"&i3]=result)
        foreach(1..9, "p")
                mask(p, image=_["adj_mask"&i3])
                set(mask3=mask)
                mask(p+1, image=_["adj_mask"&i3])
                rename(mask4=mask)
                minus(mask3, mask4, neg_method="zero")
                set(_["m"& p]=result)
                if((i3*9)+p<10) 
                        set(pad="0")
                        else()
                        set(pad="")
                end()
                calcstat("sum", AttrName="tas"& pad & (i3*9)+p , Stencil=body, Image=_["m"& p], objects=objects)
                div(_["objects.tas"&pad & (i3*9)+p], _["objects.tas_mTotal"&i3+1], result_type="float")
                set(_["objects.tas"&pad & (i3*9)+p]=result)
        end()
//p
end()
// i3
foreach(1..27, "t")
        if(t<10) 
                set(pad="0")
                else()
                set(pad="")
        end()
        // t<19
        Insert_At(table, _["objects.TAS"& pad & t], "TAS"&pad & t & "_" & (range*100), table.columncount)
end()
// t
foreach(1..3, "m")
        Insert_At(table, _["objects.tas_mTotal"&m], "tas_mTotal"&m , table.columncount)
end()   
// m
set(TAS_table=table)
        foreach(1.._TASmode.size-1, "m2")
                set(TASmode=_TASmode[m2])
                if(TASmode<10)  
                        set(_["WholeCells."&signal&"_TXT_"&stencilx&"_TAS_0"&TASmode&"_"&TNum]=_["TAS_table.TAS0"&TASmode&"_"&TNum])
                        else()
                        set(_["WholeCells."&signal&"_TXT_"&stencilx&"_TAS_"&TASmode&"_"&TNum]=_["TAS_table.TAS"&TASmode&"_"&TNum])
                end()
                //if m2<10
        end()
        // foreach TAS#, ie m2
        delete("TAS_table")
end()
// foreach TAS range value s
delete(s, m2, TNum, TRange, _TRange, _TASmode, TASmode, _TRangeMode)
end()
// if(IN_TAS==1)
end()
//foreach 1..3 e
end()
//foreach 1..nCh j
delete(adj_img,adj_mask0,adj_mask1,adj_mask2,appliedshift,i3,in_image,in_objects,kernel,m,m1,m2,m3,m4,m5,m6,m7,m8,m9,mask1,mask2,mask3,mask4,mask_all,non_numerics_count,p,pad,range,range_max,range_min,ratio_image,ratio_method,removeNAN,result,t,table,img_avg,image00,objects00)
end()
// if IN_textfeat=1

// 2.6.5 - CLEANUP
if(IN_cleanup==1)
delete(e,stencilx)
delete(WholeCells_1,WholeCells_2,WholeCells_3)
end()
// cleanup
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 2.7 - SIGNAL COLOCALIZATION MEASUREMENTS
if(IN_Pearsons==1 OR IN_Manders==1 OR IN_ICQ==1)
/////////////////////////////////////
// 2.7.1 - TWO SIGNAL ANALYSIS  
// Rollingball background correction: Created by Tony Collins, MBF, 25 Sept 2006; adapted from MBF_BG_Correction.proc by Jarkko Ylanko, MBF, 14 December 2010
foreach(1..nCh, "cc0")
        set(im_cur=_["IM_"&cc0])
        set(resizefactor=4)
        set(ballradius=15)
        resize(1/resizefactor, image=im_cur, yfactor=1/resizefactor)
        eros(edge=ballradius)
        resize(resizefactor, image=image, yfactor=resizefactor)
        set(backval = image.mean)
        set(xslip = (im_cur.width - image.width)/2)
        set(yslip = (im_cur.height - image.height)/2)
        redimension(im_cur.width, im_cur.height, xslip, yslip, backval, image=image, BackgroundMethod="mirror")
        set(BGimage = image)
        minus(im_cur, BGimage, neg_method="zero")
        set(BGCimage=result)
        set(_["IM_"&cc0&"bgc"]=BGCimage)
end()
// cc0
delete(cc0, im_cur, resizefactor, ballradius, backval, xslip, yslip, BGimage, BGCimage, result)

foreach(1..nCh, "cc1")
foreach(1..nCh, "cc2")
if(cc1!=cc2 AND cc2>cc1)
set(Ch1=_["IM_"&cc1&"bgc"])
set(Ch2=_["IM_"&cc2&"bgc"])
if(IN_Pearsons==1)
// Pearsons Coefficient: created by Tony Collins, MBF, 25 Sept 2006; adapted from MBF_ColocalisationCoefficients b03.proc lines #1-128 by Jarkko Ylanko, MBF, 14 December 2010
        set(ROI_obj=WholeCells)
        CalcIntensity(Image=Ch1, objects=ROI_obj, CalcStdDev=yes)
        rename(ROI_obj = objects)
        blank(Ch1.width, Ch1.height)
        convelems(image, sign = "signed")
        carryobjects(ROI_obj.body, ROI_obj.intensity, image=result)
        rename(Ch1m = image)
        minus(Ch1, Ch1m, result_type="long, signed")
        rename(Ch1subCh1m = result)
        CalcIntensity(Image=Ch2, objects=ROI_obj, CalcStdDev=yes)
        rename(Ch2_obj = objects)
        blank(Ch1.width, Ch1.height)
        convelems(image, sign = "signed")
        carryobjects(ROI_obj.body, Ch2_obj.intensity, image=result)
        rename(Ch2m = image)
        minus(Ch2, Ch2m, result_type = "long, signed")
        rename(Ch2subCh2m = result)
        mul(Ch2subCh2m, Ch1subCh1m, result_type = "long, signed")
        rename(PDM = result)
        blank(Ch1.width, Ch1.height)
        convelems(image, "floating", sign="unsigned")
        minus(PDM, 0, neg_method = "zero", result_type="float")
        set(PDM_pos_image = result)
        calcIntensity(Image = result, objects = ROI_obj, Total=yes)
        carryObjects(ROI_obj.body, objects.intensity, image = result)
        set(positive_PDM=image)
        mul(PDM, -1, result_type = "long")
        set(PDM_neg_image = result)
        minus(result, 0, neg_method = "zero", result_type="float")
        calcIntensity(Image = result, objects = ROI_obj, Total=yes)
        carryObjects(ROI_obj.body, objects.intensity, image = result)
        set(negative_PDM=image)
        minus(positive_PDM, negative_PDM, result_type = "float, signed")
        set(Pnumerator=result)
        mul(Ch1subCh1m, Ch1subCh1m, result_type = "float, signed")
        rename(Ch1subCh1mSqrd = result)
        calcIntensity(body, Ch1subCh1mSqrd,objects = ROI_obj, total=yes)
        rename( Ch1subCh1mSqrd_obj =objects)
        mul(Ch2subCh2m, Ch2subCh2m, result_type = "float, signed")
        rename(Ch2subCh2mSqrd = result)
        calcIntensity(body, Ch2subCh2mSqrd,objects = ROI_obj, total=yes)
        rename( Ch2subCh2mSqrd_obj =objects)
        blank(Ch1.width, Ch1.height)
        convelems(image, "floating", sign="signed")
        carryObjects(ROI_obj.body, Ch1subCh1mSqrd_obj.intensity, image = result)
        rename(sumCh1subCh1mSqrd = image)
        blank(Ch1.width, Ch1.height)
        convelems(image, "floating", sign= "signed")
        carryObjects(ROI_obj.body, Ch2subCh2mSqrd_obj.intensity, image = result)
        rename(sumCh2subCh2mSqrd = image)
        mul(sumCh2subCh2mSqrd, sumCh1subCh1mSqrd, result_type = "float, signed")
        CompAny("sqrt(result)") 
        rename(Pdenominator = result)
        div(Pnumerator, Pdenominator, result_type="float")
        rename(Pearsons =result)
        calcIntensity(body,Pearsons, objects = ROI_obj)
        renameattr(PearsonsR=intensity, objects=objects)
        delete(appliedshift,Ch1m,Ch1subCh1m,Ch1subCh1mSqrd,Ch1subCh1mSqrd_obj,Ch2m,Ch2subCh2m,Ch2subCh2mSqrd,Ch2subCh2mSqrd_obj,Ch2_obj,negative_PDM,non_numerics_count,Pdenominator,PDM,PDM_neg_image,PDM_pos_image,Pearsons,Pnumerator,positive_PDM,ROI_obj,sumCh1subCh1mSqrd,sumCh2subCh2mSqrd)
        set(_["WholeCells.Ch"&cc1&"_CLC_Ch"&cc2&"_PearsonsR"]=objects.PearsonsR)
        end()
        // IN_Pearsons==1
if(IN_Manders==1)
        // Mander's Coefficients: created by Tony Collins, MBF, 25 Sept 2006; adapted from MBF_ColocalisationCoefficients b03.proc lines #206-325 by Jarkko Ylanko, MBF, 14 December 2010
        set(ROI_obj=WholeCells)
        calcIntensity(body,Ch1, objects = ROI_obj, total=yes)
        rename(Ch1sum_obj = objects)
        calcIntensity(body,Ch2, objects = ROI_obj, total=yes)
        rename(Ch2sum_obj = objects)
        blank(Ch1.width, Ch1.height)
        convelems(image, "floating", sign="signed")
        carryObjects(ROI_obj.body, Ch1sum_obj.intensity, image = result)
        rename(Ch1sum=image)
        blank(Ch1.width, Ch1.height)
        convelems(image, "floating", sign="signed")
        carryObjects(ROI_obj.body, Ch2sum_obj.intensity, image = result)
        rename(Ch2sum=image)
        set(maskthresh=1)
        mask(image=Ch1, threshold=maskThresh, userealvalues=yes)
        rename(Ch1mask=mask)
        mask(image=Ch2, threshold=maskThresh, userealvalues=yes)
        rename(Ch2mask=mask)
        mul(Ch1, Ch2mask)
        rename(Ch1coloc = result)
        mul(Ch2, Ch1mask)
        rename(Ch2coloc = result)
        calcIntensity(body,Ch1coloc, objects = ROI_obj, total=yes)
        rename(Ch1coloc_obj = objects)
        calcIntensity(body,Ch2coloc, objects = ROI_obj, total=yes)
        rename(Ch2coloc_obj = objects)
        blank(Ch1.width, Ch1.height)
        convelems(image, "floating", sign="signed")
        carryObjects(ROI_obj.body,Ch1coloc_obj.intensity, image = result)
        rename(Ch1colocsum=image)
        blank(Ch1.width, Ch1.height)
        convelems(image, "floating", sign="signed")
        carryObjects(ROI_obj.body, Ch2coloc_obj.intensity, image = result)
        rename(Ch2colocsum=image)
        div(Ch1colocsum, Ch1sum, result_type="float")
        rename(M1=result)
        calcIntensity(body, M1, objects = ROI_obj)
        rename(M1_obj=objects)
        div(Ch2colocsum, Ch2sum, result_type="float")
        rename(M2=result)
        calcIntensity(body,M2, objects = ROI_obj)
        rename(M2_obj=objects)
        mul(Ch1, Ch2)
        rename(Ch1Ch2=result)
        calcIntensity(body, Ch1Ch2, objects=ROI_obj, Total=yes)
        rename(Ch1Ch2_obj=objects)
        blank(Ch1.width, Ch1.height)
        convelems(image, "floating", sign="signed")
        carryObjects(ROI_obj.body,Ch1Ch2_obj.intensity, image = result)
        rename(Rnumerator=image)
        mul(Ch1,Ch1)
        rename(Ch1sqrd=result)
        calcIntensity(body, Ch1sqrd, objects=ROI_obj, Total=yes)
        rename(Ch1sqrd_obj=objects)
        blank(Ch1.width, Ch1.height)
        convelems(image, "floating", sign="signed")
        carryObjects(ROI_obj.body,Ch1sqrd_obj.intensity, image = result)
        rename(Ch1sqrdsum=image)
        mul(Ch2,Ch2)
        rename(Ch2sqrd=result)
        calcIntensity(body, Ch2sqrd, objects=ROI_obj, Total=yes)
        rename(Ch2sqrd_obj=objects)
        blank(Ch1.width, Ch1.height)
        convelems(image, "floating", sign="signed")
        carryObjects(ROI_obj.body,Ch2sqrd_obj.intensity, image = result)
        rename(Ch2sqrdsum=image)
        mul(Ch2sqrdsum,Ch1sqrdsum, result_type = "float, signed")
        CompAny("sqrt(result)") 
        rename(Rdenominator = result)
        div(Rnumerator, Rdenominator,result_type="float")
        rename(MandersR =result)
        calcIntensity(body,MandersR, objects = ROI_obj)
        rename(MandersR_obj = objects)
        setAttr("M1", M1_obj.intensity, objects=ROI_obj)
        setAttr("M2", M2_obj.intensity, objects=objects)
        setAttr("MandersR", MandersR_obj.intensity, objects=objects)
        delete(appliedshift,Ch1Ch2,Ch1Ch2_obj,Ch1coloc,Ch1colocsum,Ch1coloc_obj,Ch1mask,Ch1sqrd,Ch1sqrdsum,Ch1sqrd_obj,Ch1sum,Ch1sum_obj,Ch2coloc,Ch2colocsum,Ch2coloc_obj,Ch2mask,Ch2sqrd,Ch2sqrdsum,Ch2sqrd_obj,Ch2sum,Ch2sum_obj,M1,M1_obj,M2,M2_obj,MandersR,MandersR_obj,maskthresh,non_numerics_count,Rdenominator,Rnumerator,ROI_obj)
        set(_["WholeCells.Ch"&cc1&"_CLC_Ch"&cc2&"_MandersR"]=objects.MandersR)
        set(_["WholeCells.Ch"&cc1&"_CLC_Ch"&cc2&"_M1"]=objects.M1)
        set(_["WholeCells.Ch"&cc1&"_CLC_Ch"&cc2&"_M2"]=objects.M2)
end()
// IN_Manders==1
if(IN_ICQ==1)
//131-203
// Quantitative Immunocolocalization (ICQ): created by Tony Collins, MBF, 25 Sept 2006; adapted from MBF_ColocalisationCoefficients b03.proc lines #131-203 by Jarkko Ylanko, MBF, 15 December 2010
        set(ROI_obj=WholeCells)
        CalcIntensity(Image=Ch1, objects=ROI_obj, CalcStdDev=no)
        rename(ROI_obj = objects)
        blank(Ch1.width, Ch1.height)
        convert(4,image=image)
        convelems(image, sign = "signed")
        carryobjects(ROI_obj.body, ROI_obj.intensity, image=result)
        rename(Ch1m = image)
        minus(Ch1, Ch1m, result_type="long, signed")
        rename(Ch1subCh1m = result)
        CalcIntensity(Image=Ch2, objects=ROI_obj, CalcStdDev=no)
        rename(Ch2_obj = objects)
        blank(Ch1.width, Ch1.height)
        convelems(image, "floating", sign = "signed")
        carryobjects(ROI_obj.body, Ch2_obj.intensity, image=result)
        rename(Ch2m = image)
        minus(Ch2, Ch2m, result_type = "long, signed")
        rename(Ch2subCh2m = result)
        mul(Ch2subCh2m, Ch1subCh1m, result_type = "long, signed")
        rename(PDM = result)
        minus( PDM, 0, result_type="automatic", neg_method="zero")
        rename(positive_images = result)
        mask(0.5, image=positive_images, userealvalues=yes)
        rename(positivePixels= mask)
        calcIntensity(body, positivePixels, objects=ROI_obj, Total=yes)
        rename(positive_obj = objects)
        blank(Ch1.width, Ch1.height)
        convelems(image, sign = "signed")
        carryObjects(ROI_obj.body, positive_obj.intensity, image = result)
        rename(ICQnumerator = image)
        blank(Ch1.width, Ch1.height)
        convelems(image, sign = "signed")
        calcAttr("area", objects=ROI_obj)
        carryObjects(ROI_obj.body, objects.area, image = result)
        rename(ICQdenominator = image)
        div(ICQnumerator, ICQdenominator, result_type="float")
        rename(ICQ = result)
        calcIntensity(body, ICQ, objects = ROI_obj)
        renameattr(ICQ=intensity, objects=objects)
        set(_["WholeCells.Ch"&cc1&"_CLC_Ch"&cc2&"_ICQ"]=objects.ICQ)
        delete(appliedshift,Ch2m,Ch2subCh2m,Ch2_obj,foreach_index,ICQ,ICQdenominator,ICQnumerator,non_numerics_count,PDM,positivePixels,positive_images,positive_obj,ROI_obj,Ch1m,Ch1subCh1m)
end()
// IN_ICQ==1
delete(Ch1,Ch2)
end()
// cc1!=cc2 AND cc2>cc1
end()
// cc2
end()
// cc1
foreach(1..nCh, "cc0")
delete(_["IM_"&cc0&"bgc"])
end()
// cc0
delete(cc0,cc1,cc2)
/////////////////////////////////////
end()
// if IN_Pearsons==1 OR IN_Manders==1 OR IN_ICQ==1
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 2.8 - DATA ACCUMULATION
if(defined("all_cells"))
        addobjects(wholecells, CheckOverlap=no, DeleteGeometry=yes, objects=all_cells)
        set(All_Cells=objects)
        else()
        set(All_Cells=wholecells)
end()
// if defined all_cells
end()
//if objects.count>0
end()
// IN_adjust
end()
// start..endfield
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 3. DATA OUTPUT
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 3.1 - OUTPUT IF NO ANALYSIS PERFORMED
if(defined("all_count"))
        output(all_count.count, "Preliminary object count")
end()
// 3.2 - OUTPUT IF ANALYSIS PERFORMED
if(defined("all_cells"))
CreateScalarAttrTable(all_cells)
set(table=attrtable)
insert_At(table, filledvec(table.@rowcount, __date__), "DateOfAnalysis", 3)
insert_At(table, table.aaFieldofView, "FieldofView", 0)
insert_At(table, table.aaXcoord, "Xcoord", 1)
insert_At(table, table.aaYcoord, "Ycoord", 2)
delete(table.aaFieldofView,table.aaXcoord,table.aaYcoord)
set(attrtable=table)
// 3.2.1 - INCLUSION OF PLATE INFORMATION
if(defined("platemap_tab"))
        set(table_first=platemap_tab)
        set(table_second=attrTable)
        if(table_first.@columncount==0)
                set(table=table_second)
                return()
        end()
        // empty table
        if(table_first.@rowcount<table_second.@rowcount)
                create("table")
                set(table_first_new=table)
                foreach(0..table_first.@columncount-1, "tf")
                        set(newVecName = table_first.columns[tf])
                        set(newVec  =_["table_first."&newVecName])
                        set(newFilledVec = filledvec(table_second.@rowcount, newvec[0]))
                        Insert_At(table_first_new, newfilledvec, newVecName, table_first_new.@columncount)
                        set(table_first_new=table)
                end()
                // tf
                set(table_first=table_first_new)
        end()
        // row count check
        foreach(0..table_second.@columncount-1, "ts")
                set(newVecName = table_second.columns[ts])
                set(newVec  =_["table_second."&newVecName])
                Insert_At(table_first, newvec, newVecName, table_first.@columncount)
                set(table_first=table)
        end()
        // ts
        delete(tf, ts, table_first, table_second, newVecName, newVec, newFilledVec, table_first_new)
        else()
        set(wellindex=1*wellindex)
        set(table.aaWellIndex=filledvec(attrTable.rowcount, wellindex))
end()
// if defined platemap_tab
// 3.2.2 - OBJECT-LEVEL DATA OUTPUT
if(IN_gallery==1)
set(OP_data=0)
end()
// pre-empt output of non-image data
if(OP_data==1)
        if(OP_treat==1)
                if(OP_treatID!="NA")
                        set(_OPtreatID=_["table."&OP_treatID])
                        set(_OPtreatID=_OPtreatID[0])
                        else()
                        set(_OPtreatID="NA")
                end()
                // if OP_treatID defined
                set(fileOut=dir2&"\\_"&_OPtreatID&"_output.txt")
                else()
                set(fileOut=dir2&"\\_output.txt")
        end()
        // save by OP_treat = 1
        set(table=table)
        set(filename=fileOut)
        set(h1="")
        set(keys=table.columns)
        eval(readfile(filename))        
        if(warningcode==2050)
                foreach(keys, "k") 
                        typeinfo(table[k])
                        if(typeinfo.elemclass == "memblock")
                                delete(table[k])
                                else()
                                set(h1=h1 & unquote|text(k,"name") & "\t")
                        end()
                        // if memblock
                end()
                // foreach k
                        set(h1 = h1 & "\n")
        end()
        // warningcode=2050
        write(table, "<h1", "ascii", append=yes)
        write(h1,filename,"ascii", append=yes)
end()
//if OP_data = yes
if(IN_cleanup==1)
delete(OP_treat, OP_treatID, _OPtreatID)
end()
// cleanup
// 3.2.3 - WELL-LEVEL DATA OUTPUT
if(IN_fx==1)
set(table=table)
set(stat="mean")
foreach(1..table.@columncount-1, "s")
if(s==15)
        if(IN_fov==0)
                output(fov, "Number of fields selected")
                output(fov-fieldcount, "Number of fields skipped")
                else()
                output(1, "Number of fields selected")
        end()
        // IN_fov
        output(table.@rowcount,"Object count")
        end()
        //s=15
        if(s<11 OR s>13)
        set(elem=_["table." & table.@columns[s] & ".elemtype"])
        if(elem==13)    
                set(value = table.[s][0])
                output(value,table.@columns[s])
                else()
                output(_["table."& table.@columns[s]&"."& stat],table.@columns[s])      
        end()
        // elem=13
        end()
        // ignore X, Y, fov
end()
// s
delete(stat, value, elem, s)    
end()
// IN_fx==1
set(All_Data=table)     
end()
// if defined all_cells
// 3.2.4 - FIELD INFORMATION

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 4. EXPLORER CLEANUP
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
if(IN_cleanup==1)
if(IN_gallery==1)
delete(coordinatetable, cropimage, d, dir3, dir4, dirIM, appliedshift, errorcode, errors, errors_xml, fieldimage, filter, imagecount, non_numerics_count, oc, OP_Field, result, warningcode, X, Y)
end()
//in_gallery==1
        //set(_["Field_"&i&"_Nuclei"]=Nuclei)
        //set(_["Field_"&i&"_WholeCells"]=WholeCells)
        //set(_["Field_"&i&"_Ch1"]=IM_1)
        //if(defined("IM_2"))
        //      set(_["Field_"&i&"_Ch2"]=IM_2)
                delete(IM_2)
        //end()
        //// if defined IM_2
        //if(defined("IM_3"))
        //      set(_["Field_"&i&"_Ch3"]=IM_3)
                delete(IM_3)
        //end()
        //// if defined IM_3
        delete(Nuclei,WholeCells,IM_1)
        delete(AttrTable,platemap_tab,blue,green,red,image,objects,_signal,_stencil)
        delete(IN_Ch1,IN_ChDefault,IN_colocal,IN_fov,IN_ImageONOFF,IN_RGB,IN_adjust,IN_spots,IN_textfeat,_imageview,signal,TotalCellArea,StartField,EndField,quant,OP_data,OP_destination,foreach_index,fieldcount,exptFolder,distance_vector,dir2,dir,u,nCh,nImages,i,j,IN_Haralick,IN_HarStepNum,IN_ICQ,IN_intensity,IN_LSpotThresh,IN_Manders,IN_morph,IN_Pearsons,IN_RadMom,IN_TAS,IN_TASmode,IN_Zone,IN_ZoneCount,c)
        delete(cam, content, IN_cyto, IN_LSpotRound, IN_LSpots, IN_NAN, IN_RefSEG, IN_RefTXT, IN_Stencil, IN_TRangeMode, nuc_ch)
        delete(IN_cleanup, fov, IN_cytoRef, IN_fx, IN_gallery, IN_RGBCh1, OP_FieldResult, OP_imagecount, row,column, table, date, ChEnd, ChStart)
        delete(ChannelSeries)
end()
// cleanup
end()
// skip!=1