[Analysis] Channel Heights (or Levels)

[discoleo]

Channel Heights (or Levels)

Function to compute the height (levels) of linear channels.

Note:

• the function is written in C++ (Rcpp);
• should wait for R 4.3 & RTools 4.3;

R Code

### Channel Height
height.channel = function(m, id, val=-1, verbose = TRUE) {
    tmp = select.subgrid(m, id=id, pad.val = val - 1);
    colnames(tmp) = NULL;
    nrows0  = attr(tmp, "nr");
    nStart  = which(tmp[,2] == id);
    idStart = nrows0[1] + nStart - 1; # Original rows;
    isChannel = idStart %% 2 == 1;
    nStart = nStart[isChannel];
    if(verbose) print(nStart);
    tmp[tmp > 0] = 0;
    tmp = t(tmp);
    dim = dim(tmp);
    # nStart = (nStart - 1) * dim[1] + 1;
    tmp = heightChannel(nStart, tmp, dim[1], dim[2]);
    tmp = array(tmp, dim);
    tmp = t(tmp);
    tmp = tmp[- c(1, nrow(tmp)), - c(1, ncol(tmp)), drop=FALSE];
    attr(tmp, "nr") = nrows0;
    invisible(tmp);
}

height.channel.all = function(m, verbose = FALSE, val=-1) {
    m[m > 0] = 0;
    m = t(m);
    dim  = dim(m);
    padx = rep(val - 1, dim[2]);
    pady = rep(val - 1, dim[1] + 2);
    m = rbind(padx, m, padx);
    m = cbind(pady, m, pady);
    dim = dim + 2;
    nStart = seq(3, dim[2] - 1, by=2);
    m = heightChannel(nStart, m, dim[1], dim[2], verbose=verbose);
    dim = dim - 1;
    m = m[2:dim[1], 2:dim[2]];
    m = t(m);
    invisible(m);
}

# [very slow]
height.channel.debug = function(m, verbose = FALSE, val=-1) {
    ids = unique(m[,1]);
    ids = ids[ids > 0];
    if(verbose) print(ids);
    r = m;
    r[r > 0] = 0; # required for Ideal Pores;
    for(id in ids) {
        tmp = height.channel(m, id=id, val=val, verbose=verbose);
        nrs = attr(tmp, "nr");
        isVal = tmp > 0;
        r[seq(nrs[1], nrs[2]), ][isVal] = tmp[isVal];
    }
    invisible(r);
}

[discoleo]

C++ Code

// Ideal Pores: considered to NOT communicate with adjacent pores; NumericVector heightChannel(IntegerVector xStart, NumericVector m, const int nRows, const int nCols, const bool idealPores = false, const bool verbose = false) { // m = transposed Matrix, with Rows & Cols padded; // xStart = Pore entry {nr, nc} std::vector<int> xyPoreEntry; // Pore offset {nr, nc} for(int i = 0; i < xStart.length(); i++) { xyPoreEntry.push_back(1); xyPoreEntry.push_back((int) xStart[i] - 1); // Column offset } // Start std::vector<int> xyNewPores; // Pore entry {x, y} double idLevel = 1; if(verbose) Rcout << "Starting: Inflows = " << xyPoreEntry.size() << ";\\n"; while(true) { for(size_t nv = 0; nv < xyPoreEntry.size(); ) { int nr0 = xyPoreEntry[nv]; nv++; int nc0 = xyPoreEntry[nv]; nv++; int offsetCol = nc0*nRows; int valE = m[offsetCol + nr0]; if(valE != 0) continue; // Sub-Channel: int nStart = nr0; m[offsetCol + nStart] = idLevel; while(nStart > 0) { if(m[offsetCol + nStart - 1] == 0) { nStart --; m[offsetCol + nStart] = idLevel; } else break; } int nEnd = nr0; while(nEnd < nRows) { if(m[offsetCol + nEnd + 1] == 0) { nEnd ++; m[offsetCol + nEnd] = idLevel; } else break; } // search Pores: // Previous Column int offPrevW = (nc0 - 1)*nRows; int offPrevCh = (nc0 - 2)*nRows; int npos = nStart; while(npos <= nEnd) { int nposNew = offPrevW + npos; if(m[nposNew] == 0) { m[nposNew] = idLevel; // to check or not to check? if(m[offPrevCh + npos] == 0) { // if(verbose) Rcout << " Up = " << (nc0 - 2) << ","; xyNewPores.push_back(npos); xyNewPores.push_back(nc0 - 2); } } npos ++; } if( ! idealPores) { int nposNew = offPrevW + nStart - 1; if((m[nposNew] == 0) && (m[nposNew + 1] > 0)) { m[nposNew] = idLevel; xyNewPores.push_back(npos); xyNewPores.push_back(nc0 - 2); } nposNew = offPrevW + nEnd + 1; if((m[nposNew] == 0) && (m[nposNew - 1] > 0)) { m[nposNew] = idLevel; xyNewPores.push_back(npos); xyNewPores.push_back(nc0 - 2); } } // Next Column int offNextW = (nc0 + 1)*nRows; int offNextCh = (nc0 + 2)*nRows; npos = nStart; while(npos <= nEnd) { int nposNew = offNextW + npos; if(m[nposNew] == 0) { m[nposNew] = idLevel; if(m[offNextCh + npos] == 0) { xyNewPores.push_back(npos); xyNewPores.push_back(nc0 + 2); } } npos ++; } if( ! idealPores) { int nposNew = offNextW + nStart - 1; if((m[nposNew] == 0) && (m[nposNew + 1] > 0)) { m[nposNew] = idLevel; xyNewPores.push_back(npos); xyNewPores.push_back(nc0 + 2); } nposNew = offNextW + nEnd + 1; if((m[nposNew] == 0) && (m[nposNew - 1] > 0)) { m[nposNew] = idLevel; xyNewPores.push_back(npos); xyNewPores.push_back(nc0 + 2); } } } xyPoreEntry.clear(); if(xyNewPores.size() == 0) break; xyPoreEntry.insert(xyPoreEntry.begin(), xyNewPores.begin(), xyNewPores.end()); xyNewPores.clear(); idLevel ++; } if(verbose) Rcout << "\\n"; return m; }

[adrian0010]

solved in ec806855fe8f565396d18909e34df23bb8cbb73f