setBatchMode(true); /* * This macro code will take a multi-channel TIFF stack of a Tetrahymena cell stained with a robust basal body marker * (i.e., Centrin) and identify the location of all basal bodies within the image as well as each basal body's most * logical partner within a ciliary row. The code is split up into * */ initMacroReturn=initMacro();//Initializes the macro by getting the stage of the macro that needs to be run and the folder structure for what the macro should process. /* * Part 1-Find BBs */ if(initMacroReturn[0]=="1. Find BBs"){ initPart1Return=initPart1(initMacroReturn[1]);//Initializes the macro by getting the: 1) channel names 2) the channel number for the robust BB marker of interest and 3) whether the macro is to be run on a single image or a directory of images. This information is returned as an array. fileList=createFileList(initMacroReturn[1],initPart1Return[8]);//Creates a list of all the files (directory+file name) that the macro will run on. Returns the file list as an array. for(i=0; i0){ modifyBinaryImage(bbCentroidsTitle,"delete",deleteBBs,1,"square"); selectImage(mergeImage); close(); } /* * Add missed BBs */ run("Merge Channels...", "c1=&bbImage c2=&bbCentroidsTitle create keep"); mergeImage=getImageID(); setBatchMode("show"); addBBs=interactiveCursor(0.5,1,mergeImage,bbImage,1);//input for width and depth are in microns setBatchMode("hide"); if(addBBs.length>0){ modifyBinaryImage(bbCentroidsTitle,"add",addBBs,1,"square"); selectImage(mergeImage); close(); } /* * Rewrite the text file if the user updated the BBs */ if(deleteBBs.length>0 || addBBs.length>0){ roiManager("Reset"); selectImage(bbCentroidsTitle); makeBinaryImage(); bbs=convertSelectionsToRois(0,1); finalBbsCoords=convertRoisToCoords(1); File.delete(saveDir+saveNameNoExt+".txt"); writeTo=File.open(saveDir+saveNameNoExt+".txt");//Creates a .txt file in the parent directory with the same name as the image file. All image information will be written to this file. The file is assigned the varaible writeTo. writeToOpenFile(writeTo,dialogText,"None");//Writes the dialog information from line 5 to the writeTo file. writeToOpenFile is a custom function that adds text to an open .txt file. writeToOpenFile(writeTo,cellOutline,"Largest Feature"); writeToOpenFile(writeTo,cropBox,"Cropped Feature"); writeToOpenFile(writeTo,oaFeature,"OA"); writeToOpenFile(writeTo,finalBbsCoords,"Final BBs"); writeToOpenFile(writeTo,poleFeature,"Poles");//Writes the coordinates of the poles to the writeTo file. File.close(writeTo);//The write to file must be closed since only a single writeTo file can be open at any time } run("Close All");//Closes all the images that are currently open roiManager("Reset");//Resets the ROI manager so it is blank before the next file is run } } /* * Part 3-Make Connections */ if(initMacroReturn[0]=="3. Make Connections"){ initPart2Return=initPart2(initMacroReturn[1]);//Initializes Part 2 of the macro fileList=createFileList(initMacroReturn[1],initPart2Return[0]);//Creates a list of all the files (directory+file name) that the macro will run on. Returns the file list as an array. for(i=0; i0){//If the user inputed a string into the Channel 1 dialog, the increase images by 1. images++; } channel2=Dialog.getString;//Retrieve the name for channel 2 if(lengthOf(channel2)>0){//If the user inputed a string into the Channel 1 dialog, the increase images by 2. images++; } channel3=Dialog.getString;//Retrieve the name for channel 3 if(lengthOf(channel3)>0){//If the user inputed a string into the Channel 1 dialog, the increase images by 3. images++; } channel4=Dialog.getString;//Retrieve the name for channel 4 if(lengthOf(channel4)>0){//If the user inputed a string into the Channel 1 dialog, the increase images by 4. images++; } zStep=Dialog.getNumber; dialogArray=newArray(totalChannels,robustChannel,channel1,channel2,channel3,channel4,images,zStep);//The dialogArray will be filled with the information gathered from the dialog. return dialogArray; } /* * Creates a file list for the main part of the macro to process if it is in batch mode. */ function createFileList(fileStructure,file){//Arguments 1)"fileStructure" is the file structure that the analysis will work on, 2)"file" is the path to the image, folder or folder of folders, 3) "images" is the number channels for each image in the folder structure filesToProcess=newArray(); if(fileStructure=="Single File"){ dir=File.getParent(file)+File.separator; name=File.getName(file); path=dir+name; filesToProcess=Array.concat(filesToProcess,path); } if(fileStructure=="Single Folder"){ dir=file; temp=getFileList(dir); filesToProcess=newArray(); for(i=0; i0){ getSelectionCoordinates(x,y); z=getSliceNumber(); selectionArray=Array.concat(selectionArray,"x",x,"y",y,"z",z); run("Select None"); roiManager("Reset"); } } return selectionArray; } /* * Convert selections into list of coordinates. * Requires a binary image or a non binary image with a selection. */ function convertSelectionsToRois(slice,returnArray){ arrayCoords=newArray(); if(slice==1 || nSlices==1){ getThreshold(min,max); threshold=min+max; if(is("binary")==true || threshold>-2){ if(threshold==-2){ setThreshold(255,255); } run("Create Selection"); roiManager("Add"); if(selectionType==9){ roiManager("Split"); roiManager("Select", 0); roiManager("Delete"); } } run("Select None"); }else{ for(i=0; i-2){ if(threshold==-2){ setThreshold(255,255); } run("Create Selection"); getStatistics(area,mean,min,max,std,histo); if(max>0){ index=roiManager("Count"); roiManager("Add"); if(selectionType==9){ roiManager("Split"); roiManager("Select", index); roiManager("Delete"); run("Select None"); } } run("Select None"); } resetThreshold(); } resetThreshold(); } run("Select None"); resetThreshold(); if(returnArray==1){ for(i=0; i0){ roiManager("Select", deleteArray); roiManager("Delete"); run("Select None"); } roiManager("Select", 0); getStatistics(area,mean,min,max,std,histo); run("Convex Hull"); roiManager("Update"); roiManager("Select", 0); roiManager("Remove Slice Info"); largestCoords=convertRoisToCoords(1); run("Select None"); if(area0){ deleteArray=newArray(); for(i=0; i2000 || feret>60){ deleteArray=Array.concat(deleteArray,i); } } if(feretcirc==1){ if(feret>47 && circ<0.85){ if(deleteArray.length>0){ if(deleteArray[deleteArray.length-1]!=i){ deleteArray=Array.concat(deleteArray, i); } } else{ deleteArray=Array.concat(deleteArray, i); } } } run("Select None"); } if(deleteArray.length>0){ roiManager("Select", deleteArray); roiManager("Delete"); run("Select None"); } } if(roiManager("Count")>0){ array=convertRoisToCoords(1); } return array; } /* * A general method to extract features from images. * Gaussian is the size of the blurring radius in microns (3 micron works well to find Tet) * Laplacian is the size of the Laplacian radius in microns (1 */ function roughFeatureExtraction(gaussian, laplacian){ getVoxelSize(width,height,depth,unit); run("Z Project...", "projection=[Max Intensity]"); temp=getImageID; getStatistics(area,mean,min,max,std); run("Subtract...", "value=&mean"); run("Gaussian Blur...", "sigma=&gaussian scaled"); factor=(0.12578/width)*laplacian; run("FeatureJ Laplacian", "compute smoothing=&factor"); temp2=getImageID; run("Invert"); run("Select All"); shrink=-1*gaussian/2; run("Enlarge...", "enlarge=&shrink"); getStatistics(area,mean,min,max,std,histogram); setMinAndMax(min,max); run("Select None"); run("8-bit"); run("Select All"); run("Enlarge...", "enlarge=&shrink"); run("Gaussian Blur...", "sigma=&gaussian scaled"); setAutoThreshold("Triangle dark"); getThreshold(threshold,max); getStatistics(area,mean,min,max,std,histogram); setThreshold(threshold, 255); if(max>=threshold){ convertSelectionsToRois(1,0); features=convertRoisToCoords(1); } selectImage(temp); close(); selectImage(temp2); close(); return features; } /* * Determine if image has a selection. If image does not have a selection * it returns 1. If image has a selection it returns 0. */ function queryForImageSelection(){ getStatistics(area,mean,min,max,std,histo); getVoxelSize(width,height,depth,unit); wide=getWidth()*width; high=getHeight()*height; area2=wide*high; if(round(area)==round(area2)){ return 1; } else{ return 0; } } /* * Find all local maxima within a 3D stack */ function findLocalMaxima(selection,constant){ dup="Duplicate"; run("Duplicate...", "title=[&dup] duplicate"); dup=getImageID(); dup2="Duplicate2"; run("Duplicate...", "title=[&dup2] duplicate"); dup2=getImageID(); getVoxelSize(height,width,depth,unit); blurRadius=0.125/width; if(blurRadius<1){ blurRadius=1; } //run("Gaussian Blur...", "sigma=&blurRadius stack"); maxXyRadius=0.5/width; if(maxXyRadius<1){ maxXyRadius=1; } maxZRadius=1.2/depth; if(maxZRadius<1){ maxZRadius=1; } if(selection.length>0){ convertArrayToRois(selection,0); run("Enlarge...", "enlarge=-1 pixel"); run("Clear Outside","stack"); run("Select None"); } run("3D Fast Filters","filter=MaximumLocal radius_x_pix=&maxXyRadius radius_y_pix=&maxXyRadius radius_z_pix=&maxZRadius Nb_cpus=4"); temp=getImageID(); noSlices=nSlices; makeBinaryImage(); convertSelectionsToRois(0,0); /* * first round */ selectImage(dup); intensity=getRoiIntensity(); Array.getStatistics(intensity,min,max,globalmean,globalstd); globalThreshold=globalmean+globalstd; slices=getRoiSlice(); rollingAvgArray=rollingSliceAverage(5,noSlices); thresholdArray=getRoiAvgWithinRange(intensity,slices,rollingAvgArray,noSlices,constant); for(i=0; i0 && tempArray[ii]+i=array3[index] && array2[ii]<=array3[index+1]){ if(count1Reset==0){ counter1=ii; count1Reset=1; } tempArray[ii]=array1[ii]; } else if(array2[ii]>array3[index+1]){ counter2=ii; ii=objects; } } if(array3[index+1]==slices){ counter2=array2.length-1; } tempArray2=Array.slice(tempArray,counter1,counter2); Array.getStatistics(tempArray2,min,max,mean,std); cv=std/mean;//This needs to be refined for EM gain images thresholdArray[counter]=((constant-(cv*constant))*std)+mean; if(cv>1){ thresholdArray[counter]=std+mean; } counter++; } return thresholdArray; } /* * This loop will go through and delete all ROI that are below * their slice range threshold. * array1 is a list of object intensities * array2 is a list of slice thresholds */ function deleteRoisBelowThreshold(array1,array2,singleThreshold){ objects=roiManager("Count"); coordArray=newArray(); for(i=0; ithreshold && array1[i]>threshold2){ coordArray=Array.concat(coordArray,slice); coordArray=Array.concat(coordArray,x[0]); coordArray=Array.concat(coordArray,y[0]); } } roiManager("Reset"); for(i=0; i0){ run("Clear","slice"); } run("Select None"); run("3D Convex Hull"); convex=getImageID(); run("3D Distance Map", "threshold=1"); distanceMap=getImageID(); temp=newArray(); for(i=0; idistance){ temp=Array.concat(temp,i); } run("Select None"); } if(temp.length>0){ roiManager("Select", temp); roiManager("Delete"); } selectWindow("Log"); run("Close"); selectImage(dup); close(); selectImage(binary); close(); selectImage(convex); close(); selectImage(distanceMap); close(); } /* * Filter maxima based on skew and kurtosis. Also returns the 3D coordinate of the OA. * The BB image must be selected heading into the macro. */ function skewKurtFilterHomog(enlargeBox){ getVoxelSize(width,height,depth,unit); title=getTitle(); dup="Duplicate"; run("Duplicate...", "title=[&dup] duplicate"); dup=getImageID(); for(i=0; i0){ deleteArray=Array.concat(deleteArray, i); } } run("Select None"); } selectImage(oaObjects); close(); if(deleteArray.length>0){ roiManager("Select", deleteArray); roiManager("Delete"); roiManager("Deselect"); run("Select None"); } coords=convertRoisToCoords(1); convertArrayToRois(coords,1); returnCoords=cleanUpCoordArray(coords); returnArray=Array.concat(returnCoords,oaArray); return returnArray; } /* * This plots BB average location */ function plotAverageBBLocation(objectArray,neighborArray,cluster,plot){ getVoxelSize(width,height,depth,unit); if(plot==1){ dup="Duplicate"; run("Duplicate...", "title=[&dup] duplicate"); dup=getImageID(); run("Select All"); run("Clear","stack"); run("8-bit"); } averageArray=newArray(objectArray.length); for(i=0; irankedMaxDistanceArray.length-(avgSize+1)){ index=counter*6; poleArray[index]=array[i*3]; poleArray[index+1]=array[(i*3)+1]; poleArray[index+2]=array[(i*3)+2]; poleArray[index+3]=clusterCoordinatesArray[i*3]; poleArray[index+4]=clusterCoordinatesArray[(i*3)+1]; poleArray[index+5]=clusterCoordinatesArray[(i*3)+2]; counter++; } } Array.getStatistics(xArray,min,max,mean,std); pole1X=0; pole1Y=0; pole1Z=0; pole2X=0; pole2Y=0; pole2Z=0; counter1=0; counter2=0; counter=0; poleStartX=poleArray[0]*width; poleStartY=poleArray[1]*height; poleStartZ=poleArray[2]*depth; for(i=0; imean2){ poleArray[0]=pole1X; poleArray[1]=pole1Y; poleArray[2]=pole1Z; poleArray[3]=pole2X; poleArray[4]=pole2Y; poleArray[5]=pole2Z; } else{ poleArray[0]=pole2X; poleArray[1]=pole2Y; poleArray[2]=pole2Z; poleArray[3]=pole1X; poleArray[4]=pole1Y; poleArray[5]=pole1Z; } return poleArray; } function crossProduct(array){ vecA1=array[0]-(array[3]); vecA2=array[1]-(array[4]); vecA3=array[2]-(array[5]); vecB1=array[0]-(array[6]); vecB2=array[1]-(array[7]); vecB3=array[2]-(array[8]); vecC1=vecA2*vecB3-vecA3*vecB2; vecC2=vecA3*vecB1-vecA1*vecB3; vecC3=vecA1*vecB2-vecA2*vecB1; cp=newArray(vecC1,vecC2,vecC3); return cp; } function createPlaneEquation(array){ xConstant=array[0]; yConstant=array[1]; zConstant=array[2]; offset=(array[0]*(-1*array[3]))+(array[1]*(-1*array[4]))+(array[2]*(-1*array[5])); planeEquation=newArray(xConstant,yConstant,zConstant,offset); return planeEquation; } function distPointPlane(array){ topTerm=abs((array[3]*array[0])+(array[4]*array[1])+(array[5]*array[2])+(array[6])); bottomTerm=sqrt(pow(array[3],2)+pow(array[4],2)+pow(array[5],2)); distance=topTerm/bottomTerm; return distance; } function planePenalty(objects,neighbors,poles){ getVoxelSize(width,height,depth,unit); clusterSize=neighbors.length/objects.length; arrayReturn=newArray(); counter=0; for(i=0; idistance2){ arrayReturn[counter]=1; } else{ arrayReturn[counter]=10000; } counter++; } } return arrayReturn; } /* * Sort scores to find the lowest score. * Return the components of the lowest score along with the xyz coordinates for the lowest score */ function sortScores(cluster,score1,score2,score3,neighbors,points){ getVoxelSize(width,height,depth,unit); arrayReturn=newArray((score1.length/cluster)*6);//This array will return the XYZ and individual score components for every object along with a binary value indicating whether the object scoreArray=newArray(score1.length); for(i=0; ipoleDist){ tempArray[ii]=10000000000000; } } Array.getStatistics(tempArray,min,max,mean,std); for(ii=0; ii0){ for(ii=0; ii1){ Array.getStatistics(tempDistanceArray,min,max,mean,std); for(ii=0; ii="1.37r"){ setOption("DisablePopupMenu", true); } getPixelSize(unit, pixWidth, pixHeight); shift=1; ctrl=2; rightButton=4; alt=8; leftButton=16; insideROI = 32; x2=-1; y2=-1; z2=-1; flags2=-1; stop=false; bbArray=newArray(); while (!stop) { selectImage(currentImage); resetMinAndMax(); getCursorLoc(x, y, z, flags); overlay=false; if (x!=x2 || y!=y2 || z!=z2 || flags!=flags2) { s = " "; if (flags&leftButton!=0){ hyperStack=Stack.isHyperstack; if(hyperStack==true){ Stack.getPosition(channel, slice, frame); }else{ slice=z; } selectImage(binaryImage); array=centerBox(x,y,slice,boxWidth,boxDepth,zSearch); centerx=parseFloat(array[0]); centery=parseFloat(array[1]); centerz=parseFloat(array[2]); upperXFinal=parseFloat(array[3]); upperYFinal=parseFloat(array[4]); boxWidthPix=parseFloat(array[5]); boxDepthPix=parseFloat(array[6]); overlay=true; bbArray=Array.concat(bbArray,centerx,centery,centerz); wait(250); } if (flags&rightButton!=0){ stop=true; } if (flags&shift!=0){ waitForUser("shit"); } } x2=x; y2=y; z2=z; flags2=flags; if(bbArray.length>1 && overlay==true){ selectImage(currentImage); if(hyperStack==true){ Stack.setSlice(centerz); }else{ setSlice(centerz); } setColor("green"); Overlay.drawRect(upperXFinal,upperYFinal,boxWidthPix,boxWidthPix); if(hyperStack==true){ Overlay.setPosition(1,centerz,1); }else{ Overlay.setPosition(centerz); } Overlay.setPosition(centerz); Overlay.show(); wait(100); } wait(100); } if (getVersion>="1.37r"){ setOption("DisablePopupMenu", false); } return bbArray; } /* * Center the box on the highest intensity voxel */ function centerBox(x,y,slice,boxWidth,boxDepth,zSearch){ title2=getImageID(); getVoxelSize(width,height,depth,unit); totalSlices=nSlices; boxWidthPix=round(boxWidth/width); boxDepthPix=round(boxDepth/depth); if(boxWidthPix%2==0){ boxWidthPix=boxWidthPix+1; } if(boxDepthPix%2==0){ boxDepthPix=boxDepthPix+1; } setSlice(slice); makeRectangle(x-((boxWidthPix-1)/2),y-((boxWidthPix-1)/2),boxWidthPix,boxWidthPix); getSelectionBounds(upperX, upperY, void, void); if(zSearch==1){ startSlice=slice-((boxDepthPix-1)/2); startRemainder=0; stopRemainder=0; frontTruncated=false; if(startSlice<=0){ startRemainder=abs(startSlice-1); startSlice=1; frontTruncated=true; } stopSlice=slice+((boxDepthPix-1)/2); backTruncated=false; if(stopSlice>totalSlices){ stopRemainder=abs(stopSlice-totalSlices); stopSlice=totalSlices; backTruncated=true; } rangeString=toString(startSlice)+"-"+toString(stopSlice); run("Duplicate...", "title=[Temp] duplicate range=&rangeString"); bits=bitDepth; temp=getImageID(); radiusX=getWidth; radiusY=getHeight; radiusZ=nSlices; if(frontTruncated==true){ if(bits==8){ newImage("TempAdd", "8-bit black", radiusX, radiusY, startRemainder); } if(bits==16){ newImage("TempAdd", "16-bit black", radiusX, radiusY, startRemainder); } run("Concatenate", " title=[Concatenated Stacks] image_1=TempAdd image_2=Temp"); rename("Temp"); temp=getImageID(); } if(backTruncated==true){ if(bits==8){ newImage("TempAdd", "8-bit black", radiusX, radiusY, stopRemainder); } if(bits==16){ newImage("TempAdd", "16-bit black", radiusX, radiusY, stopRemainder); } run("Concatenate", " title=[Concatenated Stacks] image_1=Temp image_2=TempAdd"); rename("Temp"); temp=getImageID(); } run("3D Fast Filters","filter=MaximumLocal radius_x_pix=&radiusX radius_y_pix=&radiusX radius_z_pix=&radiusZ Nb_cpus=4"); temp2=getImageID(); Stack.getStatistics(voxels,mean,min,max,std); coordsArray=newArray(); midX=((radiusX-1)/2); midY=((radiusY-1)/2); midZ=((radiusZ-1)/2)+1; for(i=0; i3){ distanceArray=newArray(); for(i=0; imidZ){ centerZ=slice+(coordsArray[(i*3)+2]-midZ); } } } }else{ centerX=upperX+coordsArray[0]; centerY=upperY+coordsArray[1]; if(coordsArray[2]==midZ){ centerZ=slice; } if(coordsArray[2]midZ){ centerZ=slice+(coordsArray[2]-midZ); } } selectImage(temp); close(); selectImage(temp2); close(); } if(zSearch==0){ run("Duplicate...", "title=[Temp]"); getRawStatistics(nPixels, mean, min, max, std, histogram); temp=getImageID(); radiusX=getWidth; radiusY=getHeight; intensityArray=newArray(); coordsArray=newArray(); midX=((radiusX-1)/2); midY=((radiusY-1)/2); for(i=0; i3){ distanceArray=newArray(); for(i=0; i