requires("1.54p");

macro "Batch_yH2AX_C3_Faster_withBatchReports_v4" {

    // =========================================================
    // PARAMETERS
    // =========================================================
    rollingBallSpots   = 8;
    gaussianSigmaSpots = 1.0;
    prominenceSpots    = 20;

    minSpotArea        = 1;
    maxSpotArea        = 150;

    localBoxRadius     = 10;
    borderMargin       = 1;

    saveOverlay        = true;
    fontSizeLabel      = 18;

    // =========================================================
    // SELECT FOLDERS
    // =========================================================
    c3Dir = getDirectory("Cartella immagini C3");
    if (c3Dir == "") exit("Nessuna cartella C3 selezionata.");

    roiDir = getDirectory("Cartella ROI nuclei");
    if (roiDir == "") exit("Nessuna cartella ROI selezionata.");

    outDir = getDirectory("Cartella output");
    if (outDir == "") exit("Nessuna cartella output selezionata.");

    if (!endsWith(c3Dir, "/")) c3Dir = c3Dir + "/";
    if (!endsWith(roiDir, "/")) roiDir = roiDir + "/";
    if (!endsWith(outDir, "/")) outDir = outDir + "/";

    fileList = getFileList(c3Dir);

    // velocizza molto la macro
    setBatchMode(true);

    // =========================================================
    // GLOBAL BATCH REPORTS
    // =========================================================
    batchPerNucleusCsv = outDir + "Batch_C3_yH2AX_per_nucleus.csv";
    batchSummaryCsv    = outDir + "Batch_C3_yH2AX_summary.csv";

    batchPerNucleusHeader =
        "Sample_ID," +
        "Nucleus_ID," +
        "Area," +
        "Mean," +
        "StdDev," +
        "IntDen," +
        "Mean_over_area," +
        "Spot_candidates_per_nucleus," +
        "Spot_count_per_nucleus," +
        "Nucleus_centroid_X," +
        "Nucleus_centroid_Y\n";
    File.saveString(batchPerNucleusHeader, batchPerNucleusCsv);

    batchSummaryHeader =
        "Sample_ID," +
        "Nuclei_count," +
        "Mean_nucleus_area," +
        "Mean_yH2AX_mean," +
        "Mean_yH2AX_sd," +
        "Mean_yH2AX_IntDen," +
        "Total_spot_candidates," +
        "Total_spots," +
        "Mean_spots_per_nucleus," +
        "SD_spots_per_nucleus," +
        "Find_Maxima_prominence," +
        "Min_spot_area," +
        "Max_spot_area," +
        "Border_margin\n";
    File.saveString(batchSummaryHeader, batchSummaryCsv);

    // =========================================================
    // LOOP FILES
    // =========================================================
    for (f = 0; f < fileList.length; f++) {

        name = fileList[f];
        path = c3Dir + name;

        if (File.isDirectory(path))
            continue;

        if (!(endsWith(name, "_C3.tif") || endsWith(name, "_C3.tiff") ||
              endsWith(name, "_C3.TIF") || endsWith(name, "_C3.TIFF")))
            continue;

        open(path);
        originalTitle = getTitle();

        run("Select None");
        roiManager("Reset");

        // -----------------------------------------------------
        // base name
        // -----------------------------------------------------
        baseName = originalTitle;
        if (endsWith(baseName, ".tif") || endsWith(baseName, ".TIF"))
            baseName = substring(baseName, 0, lengthOf(baseName) - 4);
        else if (endsWith(baseName, ".tiff") || endsWith(baseName, ".TIFF"))
            baseName = substring(baseName, 0, lengthOf(baseName) - 5);

        if (endsWith(baseName, "_C3"))
            sampleBase = substring(baseName, 0, lengthOf(baseName) - 3);
        else
            sampleBase = baseName;

        roiZip = roiDir + sampleBase + "_C1_Nuclei_ROIs.zip";

        if (!File.exists(roiZip)) {
            print("ROI non trovate: " + sampleBase);
            close();
            continue;
        }

        perNucleusCsv = outDir + sampleBase + "_C3_yH2AX_per_nucleus.csv";
        summaryCsv    = outDir + sampleBase + "_C3_yH2AX_summary.csv";
        overlayTif    = outDir + sampleBase + "_C3_yH2AX_overlay.tif";

        // -----------------------------------------------------
        // prepare image
        // -----------------------------------------------------
        Stack.getDimensions(w, h, c, z, t);

        if (z > 1) {
            run("Z Project...", "projection=[Max Intensity]");
            rename("C3_MAX");
            selectWindow(originalTitle);
            close();
        } else {
            rename("C3_MAX");
        }

        selectWindow("C3_MAX");
        run("Select None");
        run("Duplicate...", "title=C3_OVERLAY");

        // -----------------------------------------------------
        // load nuclei ROIs once
        // -----------------------------------------------------
        roiManager("Reset");
        roiManager("Open", roiZip);
        nucleiCount = roiManager("count");

        if (nucleiCount == 0) {
            close("*");
            roiManager("Reset");
            continue;
        }

        // -----------------------------------------------------
        // arrays
        // -----------------------------------------------------
        nucArea      = newArray(nucleiCount);
        nucMean      = newArray(nucleiCount);
        nucSD        = newArray(nucleiCount);
        nucIntDen    = newArray(nucleiCount);
        nucX         = newArray(nucleiCount);
        nucY         = newArray(nucleiCount);
        spotCount    = newArray(nucleiCount);
        spotCand     = newArray(nucleiCount);
        normMean     = newArray(nucleiCount);

        totalSpots   = 0;
        totalCand    = 0;
        totalMean    = 0;
        totalSD      = 0;
        totalIntDen  = 0;
        totalNucArea = 0;

        spotXs = "";
        spotYs = "";

        // =====================================================
        // LOOP NUCLEI
        // =====================================================
        for (i = 0; i < nucleiCount; i++) {

            // measurements on original image
            selectWindow("C3_MAX");
            roiManager("select", i);

            run("Set Measurements...", "area mean standard integrated centroid decimal=3");
            if (isOpen("Results")) {
                selectWindow("Results");
                run("Clear Results");
            }
            run("Measure");

            nucArea[i]   = getResult("Area", 0);
            nucMean[i]   = getResult("Mean", 0);
            nucSD[i]     = getResult("StdDev", 0);
            nucIntDen[i] = getResult("IntDen", 0);
            nucX[i]      = getResult("X", 0);
            nucY[i]      = getResult("Y", 0);

            if (nucArea[i] > 0)
                normMean[i] = nucMean[i] / nucArea[i];
            else
                normMean[i] = 0;

            totalMean += nucMean[i];
            totalSD += nucSD[i];
            totalIntDen += nucIntDen[i];
            totalNucArea += nucArea[i];

            // create processed nucleus-only image
            selectWindow("C3_MAX");
            run("Select None");
            run("Duplicate...", "title=TMP_NUC");

            selectWindow("TMP_NUC");
            roiManager("select", i);
            run("Clear Outside");

            run("8-bit");
            run("Subtract Background...", "rolling=" + rollingBallSpots);
            if (gaussianSigmaSpots > 0)
                run("Gaussian Blur...", "sigma=" + gaussianSigmaSpots);
            run("Enhance Contrast...", "saturated=0.35");

            // build inner ROI directly from the nucleus ROI
            selectWindow("C3_MAX");
            roiManager("select", i);
            run("Enlarge...", "enlarge=-" + borderMargin);

            hasInnerRoi = 1;
            if (selectionType() == -1)
                hasInnerRoi = 0;

            candidates = 0;
            accepted   = 0;

            if (hasInnerRoi) {

                // apply inner ROI directly on TMP_NUC and find maxima there
                selectWindow("TMP_NUC");
                run("Clear Outside");
                run("Find Maxima...", "prominence=" + prominenceSpots + " output=[Point Selection]");

                if (selectionType() != -1) {
                    getSelectionCoordinates(xpts, ypts);

                    for (p = 0; p < xpts.length; p++) {

                        xg = xpts[p];
                        yg = ypts[p];

                        // explicit control: point must stay in original nucleus ROI
                        selectWindow("C3_MAX");
                        roiManager("select", i);
                        if (selectionContains(xg, yg) == 0)
                            continue;

                        // local crop
                        selectWindow("TMP_NUC");

                        x0 = xg - localBoxRadius;
                        y0 = yg - localBoxRadius;
                        bw = 2 * localBoxRadius + 1;
                        bh = 2 * localBoxRadius + 1;

                        if (x0 < 0) x0 = 0;
                        if (y0 < 0) y0 = 0;
                        if (x0 + bw > w) bw = w - x0;
                        if (y0 + bh > h) bh = h - y0;

                        makeRectangle(x0, y0, bw, bh);
                        run("Duplicate...", "title=TMP_SPOT");

                        selectWindow("TMP_SPOT");
                        run("Select None");
                        run("8-bit");

                        setAutoThreshold("Default");
                        setOption("BlackBackground", false);
                        run("Convert to Mask");

                        localX = xg - x0;
                        localY = yg - y0;

                        keepArea = -1;

                        // wand directly on the particle containing the maximum
                        if (localX >= 0 && localX < bw && localY >= 0 && localY < bh) {
                            doWand(localX, localY);

                            if (selectionType() != -1) {
                                getSelectionBounds(rx, ry, rw, rh);

                                touchesBorder = 0;
                                if (rx <= 0 || ry <= 0 || (rx + rw) >= bw || (ry + rh) >= bh)
                                    touchesBorder = 1;

                                if (touchesBorder == 0) {
                                    getStatistics(areaLoc);
                                    keepArea = areaLoc;
                                }
                            }
                        }

                        candidates++;

                        if (keepArea >= minSpotArea && keepArea <= maxSpotArea) {
                            accepted++;

                            if (spotXs == "") {
                                spotXs = "" + xg;
                                spotYs = "" + yg;
                            } else {
                                spotXs = spotXs + "," + xg;
                                spotYs = spotYs + "," + yg;
                            }
                        }

                        if (isOpen("TMP_SPOT")) {
                            selectWindow("TMP_SPOT");
                            close();
                        }
                    }
                }
            }

            spotCand[i]  = candidates;
            spotCount[i] = accepted;
            totalCand += candidates;
            totalSpots += accepted;

            if (isOpen("TMP_NUC")) {
                selectWindow("TMP_NUC");
                close();
            }
        }

        // =====================================================
        // PER NUCLEUS CSV
        // =====================================================
        if (isOpen("Results")) {
            selectWindow("Results");
            run("Clear Results");
        }

        for (i = 0; i < nucleiCount; i++) {
            setResult("Image", i, sampleBase);
            setResult("Nucleus_ID", i, i + 1);
            setResult("Area", i, nucArea[i]);
            setResult("Mean", i, nucMean[i]);
            setResult("StdDev", i, nucSD[i]);
            setResult("IntDen", i, nucIntDen[i]);
            setResult("Mean_over_area", i, normMean[i]);
            setResult("Spot_candidates_per_nucleus", i, spotCand[i]);
            setResult("Spot_count_per_nucleus", i, spotCount[i]);
            setResult("Nucleus_centroid_X", i, nucX[i]);
            setResult("Nucleus_centroid_Y", i, nucY[i]);
        }

        updateResults();
        saveAs("Results", perNucleusCsv);

        for (i = 0; i < nucleiCount; i++) {
            batchLineNucleus =
                sampleBase + "," +
                (i + 1) + "," +
                d2s(nucArea[i], 6) + "," +
                d2s(nucMean[i], 6) + "," +
                d2s(nucSD[i], 6) + "," +
                d2s(nucIntDen[i], 6) + "," +
                d2s(normMean[i], 6) + "," +
                spotCand[i] + "," +
                spotCount[i] + "," +
                d2s(nucX[i], 3) + "," +
                d2s(nucY[i], 3) + "\n";
            File.append(batchLineNucleus, batchPerNucleusCsv);
        }

        // =====================================================
        // SUMMARY
        // =====================================================
        Array.getStatistics(spotCount, minS, maxS, meanS, sdS);

        meanNucArea = totalNucArea / nucleiCount;
        meanNucMean = totalMean / nucleiCount;
        meanNucSD   = totalSD / nucleiCount;
        meanIntDen  = totalIntDen / nucleiCount;

        if (isOpen("Results")) {
            selectWindow("Results");
            run("Clear Results");
        }

        setResult("Image", 0, sampleBase);
        setResult("Nuclei_count", 0, nucleiCount);
        setResult("Mean_nucleus_area", 0, meanNucArea);
        setResult("Mean_yH2AX_mean", 0, meanNucMean);
        setResult("Mean_yH2AX_sd", 0, meanNucSD);
        setResult("Mean_yH2AX_IntDen", 0, meanIntDen);
        setResult("Total_spot_candidates", 0, totalCand);
        setResult("Total_spots", 0, totalSpots);
        setResult("Mean_spots_per_nucleus", 0, meanS);
        setResult("SD_spots_per_nucleus", 0, sdS);
        setResult("Find_Maxima_prominence", 0, prominenceSpots);
        setResult("Min_spot_area", 0, minSpotArea);
        setResult("Max_spot_area", 0, maxSpotArea);
        setResult("Border_margin", 0, borderMargin);
        updateResults();
        saveAs("Results", summaryCsv);

        batchLineSummary =
            sampleBase + "," +
            nucleiCount + "," +
            d2s(meanNucArea, 6) + "," +
            d2s(meanNucMean, 6) + "," +
            d2s(meanNucSD, 6) + "," +
            d2s(meanIntDen, 6) + "," +
            totalCand + "," +
            totalSpots + "," +
            d2s(meanS, 6) + "," +
            d2s(sdS, 6) + "," +
            prominenceSpots + "," +
            minSpotArea + "," +
            maxSpotArea + "," +
            borderMargin + "\n";
        File.append(batchLineSummary, batchSummaryCsv);

        // =====================================================
        // OVERLAY
        // =====================================================
        if (saveOverlay) {

            selectWindow("C3_OVERLAY");
            run("RGB Color");

            roiManager("Deselect");
            roiManager("Show All with labels");
            run("Flatten");
            rename("OVERLAY");

            if (spotXs != "") {
                xs = split(spotXs, ",");
                ys = split(spotYs, ",");

                setForegroundColor(255, 0, 0);
                for (p = 0; p < xs.length; p++) {
                    xg = parseFloat(xs[p]);
                    yg = parseFloat(ys[p]);
                    makeOval(xg-2, yg-2, 5, 5);
                    run("Draw");
                }
            }

            setForegroundColor(255, 255, 0);
            setFont("SansSerif", fontSizeLabel, "antialiased");
            for (i = 0; i < nucleiCount; i++) {
                drawString("" + spotCount[i], nucX[i] + 5, nucY[i] - 5);
            }

            saveAs("Tiff", overlayTif);
            close();
        }

        close("*");
        roiManager("Reset");
        run("Collect Garbage");
    }

    setBatchMode(false);

    print("=== DONE ===");
    print("Batch per nucleus salvato in: " + batchPerNucleusCsv);
    print("Batch summary salvato in: " + batchSummaryCsv);
}